## Abstract

The instability theory of drumlin formation has been very successful in predicting the existence of ribbed moraine, as well as its amplitude and wavelength. However, the theory as it stands has not yet been shown to have the capability of predicting the existence of three-dimensional bedforms—drumlins—or their more extreme cousins, mega-scale glacial lineations. We extend the instability theory to include a dynamic description of the local subglacial drainage system, and in particular, we show that a uniform water-film flow between ice and deformable subglacial till is unstable, and that as a consequence, lineations will form. Predictions of the transverse wavelengths are consistent with observations.

## 1. Introduction

Mega-scale glacial lineations (MSGLs; Clark 1993) are very long topographic ridge-groove corrugations, up to tens of metres in height, hundreds of metres in width and up to 100 km in length (figure 1); much longer than drumlins, which are usually less than a kilometre in length. They were discovered on the now-exposed beds of former ice sheets, and their presence was inferred to record the location of ice streams (Clark 1993). This hypothesis has recently been validated by King *et al.* (2009), who have found MSGLs under the Rutford ice stream in Antarctica. While it is generally agreed that MSGLs are formed through the action of ice (although see Shaw *et al.* (2008) and the reply by Ó Cofaigh *et al.* (2010)), it is as yet unclear how they do so. The question is of practical, as well as scientific interest, since it has a bearing on the drag that the bed offers the ice, and hence on how rapidly the ice streams discharge. In a warming climate, these are important ingredients of the process whereby ice sheets melt or collapse and the sea level rises.

There have been a number of suggestions as to how such lineations could be formed. Clark *et al.* (2003) suggest that they are formed through the ‘combing’ action of basal ice as it ploughs through sediment. Schoof & Clarke (2008), following Shaw & Freschauf (1973), suggest that they may be formed through a transverse secondary flow in basal ice, which is caused by a non-Newtonian rheology. Many have argued that ribbed moraine, drumlins, megaflutes and MSGLs are part of a continuous spectrum of bedforms (Sugden & John 1976; Aario 1977; Boulton 1987; Rose 1987; Clark 1993), with MSGLs as the upper-end member in the scale. In this view, the bed-forming mechanism provides the observed shapes, depending on the local ice-flow conditions. Most obviously, one might suppose that the transition from ribs to drumlins to lineations occurs as the ice flow increases.

The instability theory of drumlin formation (Hindmarsh 1998; Fowler 2000, 2009; Schoof 2007*a*; Dunlop *et al*. 2008) is the first theory that is able to successfully predict qualitative characteristics of drumlinoid topography, such as amplitude and wavelength, on the basis of a physically based mathematical model. As such, its adherents might expect that the theory would provide the building blocks for all the main types of bedforms that are observed. However, this appears not to be so. In particular, Schoof (2007*a*) pointed out three difficulties with the theory. The first is that almost as soon as the bedforms begin to grow, cavities form on their lee. In a later paper, Schoof (2007*b*) finds that there are finite-amplitude travelling wave solutions, whose elevation scales with the deforming till thickness. Since this is most probably of the order of centimetres to metres, it apparently mitigates against observations of drumlins of tens of metres elevation (although, in fact, the average elevation of drumlins is actually about 7 m; Clark *et al.* 2009). This difficulty was potentially resolved by Fowler (2009), who showed that finite-amplitude bedforms of appropriate amplitude could indeed be obtained, even when cavities were present. However, Fowler’s formulation of the problem is different from Schoof’s, and relies on a plausible but heuristic assumption that cavities are infilled by sediments, in the manner of crag and tail landforms.

The second difficulty raised by Schoof is the observation that some drumlins consist of stratified sands and gravels (e.g. Menzies & Brand 2007). While such an observation is consistent with the instability theory, provided there is net erosion of the underlying sediments, it seems difficult to reconcile it to the notion that the bed-forming instability leads to travelling waves. However, in Fowler’s (2009) formulation, he finds that, in fact, the bedforms that emerge at finite amplitude are actually stationary, and this would remove the immediate difficulty concerning stratification.

The third and final difficulty is that the instability is fundamentally two dimensional. It forms transverse ridges that we interpret as the precursors of ribbed moraine, but purely longitudinal rolls that we might expect to be the precursors of MSGLs are weakly stable (Fowler 2010). Fowler raised the possibility that the two-dimensional bedforms might be subject to a secondary transverse instability, but three-dimensional numerical calculations (to be reported elsewhere) suggest that this is not the case. In this paper, we therefore focus on this issue, and provide a generalization of the instability model that is able to predict the formation of longitudinal rolls, which we suppose represent the onset of glacial lineations.

## 2. Mathematical model

The key to the formation of MSGLs is the basal hydraulic system. Water is already an essential component of the instability model, insofar as it is assumed that ice slides over wet till, where the water pressure is sufficiently high that the till can deform. In previously presented versions of the theory (e.g. Schoof 2007*a*; Fowler 2009), it was assumed that the basal hydraulic system passed through the bed-forming environment without interacting with it, other than providing a local channel effective pressure. This view is similar to that portrayed schematically by Lliboutry’s (1968) cartoon of a linked cavity system.

It is, however, obvious that such a view is overly simplistic, and a more sophisticated description is desirable for several reasons. The development of the ribbed moraine causes a barrier to water flow, which is also driven by the ice-surface slope in the direction of ice flow. Thus, the ice-till surface must have an effect on the water system, and this suggests that effective pressure must respond dynamically to the coupled ice-till flow.

Furthermore, in the absence of the bed-forming instability, a uniform film of water flowing downstream can be expected to be unstable (Smith & Bretherton 1972; Walder 1982) owing to the interaction of the water flow with the underlying sediment, provided the flow is sufficient to evacuate sediment. That this is likely is because there is now evidence that drainage under ice sheets sometimes occurs through the discharge of subglacial lakes (Wingham *et al.* 2006; Fricker *et al*. 2007), with the flow being sufficient to transport sediment. It is plausible to suppose that this same rilling instability will occur below an ice sheet, and, just as it does subaerially, cause the underlying topography to evolve into a number of streams separated by ridges or lineations.

The difference between the theory provided here and that of earlier work lies in the conceptual picture that we have of the subglacial stream system. When the streams are small, then the description of the local stream flow in terms of a locally defined effective pressure makes sense (Röthlisberger 1972), and we regain the earlier theory (Schoof 2007*a*; Fowler 2009). A different viewpoint, which we might call ‘wide stream theory’, does not distinguish the scale of the stream flow from that of the ice, and thus the streams become indistinguishable from cavities, with their depth being controlled by the criterion that the effective pressure at the ice–water interface is zero. The only difference is that cavities are formed through the deficit pressures generated by ice flow, whereas streams are caused by the excess pore pressures generated by permeable flow through the till. As we shall see, from a mathematical point of view, small stream theory is a particular case of wide stream theory. Our view of the subglacial stream system bears resemblance to the experimental results of Catania & Paola (2001), and our theoretical description of subglacial water flow is similar to that of Le Brocq *et al.* (2009), although they do not include a dynamic dependence of the water pressure on the ice flow.

The geometry of the flow we consider is indicated in figure 2. We take axes (*x*,*y*,*z*), with *z* vertically upwards and *x* pointing in the direction of ice flow. The ice surface is denoted by *z*=*z*_{i}. We define *N* to be the effective stress at the base of the ice, which we take to be at *z*=*s*. If *N*>0, then ice is assumed to rest on the underlying wet till, while if *N* reaches zero, then a water film or stream or cavity of depth *h* is formed. In general, the bed will thus consist of two parts: attached, where *N*>0 and *h*=0; and separated, where *N*=0 and *h*>0. We now wish to describe the water flow in each regime.

We define a number of different pressures: *p*_{w} is the water pressure at the ice base *z*=*s*; *p*_{T} is the pressure in the till; and *p*_{i} is the ice pressure, and in the usual way for ice flow, we subtract the cryostatic component by writing
2.1
where *p*_{a} is atmospheric pressure, *ρ*_{i} is the density of ice, *g* is the acceleration due to gravity and we call *Π* the reduced pressure. Fowler (2010) showed that the relaxation time of the pore pressure in the till was rapid, and so we suppose that the pore pressure *p* in the till is hydrostatic,
2.2
where *ρ*_{w} is the density of water. If ** σ** is the stress tensor in the ice, then the effective stress at the ice base is
2.3
where

**is the unit normal (upwards, say) at the interface. The effective pressure in the till is denoted by**

*n**p*

_{e}, and can be determined to be (using the fact that

*p*

_{T}=−

*σ*

_{nn}at

*z*=

*s*) 2.4 where 2.5

*ρ*

_{s}is the sediment density and

*ϕ*is the till porosity. Strictly, equation (2.4) assumes that

*ϕ*is constant, which is probably a reasonable approximation; relation (2.4) applies at both attached and separated parts of the base.

The hydraulic head at the ice base is *p*_{w}+*ρ*_{w}*gs*, and it is convenient to add a constant, thus we define the hydraulic head to be
2.6
where *d*_{i} is the local ice depth (at *x*=0). This is not the same as *z*_{i} because we allow for a small ice-surface slope to drive both ice and water flow. The deviatoric stress tensor is defined by
2.7
where ** δ** is the unit tensor, and it follows from this that
2.8
where
2.9
This determines

*N*in terms of ice flow and hydraulic gradient.

### (a) Groundwater flow

In general, water is formed at the base through basal melting at typical rates of millimetres per year, and this water provides a flux that must be evacuated, essentially in the direction of ice flow. If *h*=0, then we suppose water flows through the till and underlying sediments by Darcian flow.

The porosity should be a function of the local effective pressure in the till, thus *ϕ*=*ϕ*(*p*_{e}), while the local water flux is given by −(*k*/*η*_{w})**∇***ψ*, where *k* is the permeability and *η*_{w} is the viscosity of water. Here and below, **∇** is the horizontal gradient operator (∂/∂*x*,∂/∂*y*). If we assume a vertical hydrostatic equilibrium for the pore water (Fowler 2010), then *ψ* does not vary vertically. Then, the water flux integrated over a column of till is
2.10
Conservation of water mass then takes the form
2.11
if we ignore the local melt production term. The justification for this is that the water flux arises through basal melting over ice-sheet length scales *l*_{i}∼10^{6} m, whereas we are concerned with the bedform length scales *l*∼10^{2}–10^{3} m, so that the relative size of the melting term in equation (2.11) to the flux divergence term will be ∼(*l*/*l*_{i})≪1.

To find a simple form of this equation, suppose that the variations of *s*, *ϕ* and *k* are small. Then, *ϕ*≈*ϕ*(*N*)+*ϕ*′Δ*ρ*_{T}*g*(*s*−*z*), and if the variation in *s* is a good deal less than the depth of the permeable sediments, then equation (2.11) can be written approximately as
2.12
If, however, the depth of the aquifer was constant, for example being the depth of deformable till, then the term in *s*_{t} would be absent. More generally, we might include this term, but with a smaller coefficient than in equation (2.12).

### (b) Stream/cavity flow

When the ice separates from the till, water can flow through the resulting channel, whether this be a stream or a cavity. Assuming laminar flow, the water flux is (Batchelor 1967) 2.13 and mass conservation then takes the form 2.14 In fact, there will still be water flow through the till, and the total mass-flow equation should include this, but as we shall see, the flux through a stream is very much larger than that in the till, so that one can reasonably ignore it.

### (c) Sediment flux

The sediment flux owing to the deformation of till has been described by Fowler (2010). It has two elements. The motion of the ice at the interface exerts a drag on the till, and this will cause a mechanical deformation and consequent sediment flux. Quite how this depends on the sliding velocity *u*_{b} will depend on the choice of till rheology. As a saturated granular material, till has a plastic yield stress, but this statement is not sufficient to determine the flow behaviour (Fowler 2003). Plausible flow rules (e.g. Iverson & Iverson 2001; Jop *et al.* 2006) suggest an effective viscous behaviour, and a resulting Couette-type flow would have the mean till velocity being proportional to *u*_{b}. In addition, the squeezing of till via horizontal gradients of the effective pressure *p*_{e} leads to an additional, Poiseuille-type flux. Combining these yields an expression for the sediment flux ** q** of the form
2.15
where
2.16
is the depth of deformable till;

^{1}here, [

*x*]

_{+}means , 2.17 and

*τ*is the basal shear stress.

^{2}We suppose that

*c*<1, and that

*b*is given, following Fowler (2010), by 2.18 where the various constituent quantities are defined below.

### (d) Bedload transport

When the ice separates from the till, the ice can no longer drag the till, and we suppose that the sediment can move instead owing to overlying water flow. The effective shear stress *τ*_{e} exerted by a laminar flow of depth *h* and hydraulic gradient −**∇***ψ* is
2.19
where the second term on the right represents the downslope force on the particles, taken to have diameter *D*_{s}. This effective stress generates a sediment transport flux
2.20
one popular choice of transport law is the Meyer-Peter & Müller (1948) law
2.21
where
2.22
and *K*=8; *τ**_{c}=0.05 is the critical Shields stress.

### (e) Stream- and sediment-flow models

Summarizing, we have two equations for stream flow and sediment flow. The second of these is the Exner equation. For attached flow, *N*>0 and *h*=0, we have
2.23
where *ϕ*′=*ϕ*′(*p*_{e}) is the derivative of *ϕ* with respect to effective pressure and *r*_{a} is a number representing the aquifer depth: *r*_{a}=1 for a full aquifer in 0<*z*<*s*, while *r*_{a}=0 for an aquifer of fixed depth; generally we may suppose that 0≤*r*_{a}≤1.

For separated flow, where *N*=0 and *h*>0, we have
2.24
(note that the till–water interface is at *z*=*s*−*h*). These equations determine *s* and *N*, and *s* and *h*, respectively; *ψ* is defined by equation (2.8), and the stream- and sediment-flow models are coupled to the ice flow through the normal stress term in equation (2.8).

### (f) Ice-flow model

The ice-flow model has been described frequently before (Hindmarsh 1998; Fowler 2000; Schoof 2002, 2007*a*; Fowler 2009, 2010), and is only summarized here. We follow Fowler’s (2010) exposition, which assumes flow of Newtonian viscous ice over the bed, as shown in figure 2. The ice base is at *z*=*s* and the (gently sloping) ice surface is at *z*=*z*_{i}. When written in terms of the reduced pressure *Π*, the Stokes equations take the form
2.25
where ** u** is the ice velocity, and

*η*

_{i}is the ice viscosity. Kinematic boundary conditions are applied at the top and bottom (ignoring both surface accumulation/ablation and basal melting), which describe the time evolution of each surface in terms of the normal velocity of the ice. These conditions specify that ice particles in the surface remain there. In addition, we prescribe boundary conditions of continuous stress across the top surface, and no normal velocity at the base, together with a sliding law there of the form 2.26 where

**is the basal traction.**

*τ*### (g) Non-dimensionalization

The equations are non-dimensionalized in the manner devised by Schoof (2007*a*), as modified by Fowler (2009). We define an effective pressure scale *N*_{c} (which will be chosen later). Three relevant length scales are then
2.27
which is a measure of bedform elevation;
2.28
which is a measure of deforming till depth; and
2.29
which is a suitable horizontal length scale for bedforms; here, *u*_{0} is the horizontal velocity scale, which is determined from the sliding law. Since MSGLs have typical widths of the order of hundreds of metres, it is reasonable to expect that *l* (which has a similar order of magnitude) is appropriate. It should be emphasized that this choice of scale does not affect the following analysis.

In a basic uniform state where *s*=0, *Π*=0 and the ice surface has a downward slope *S* in the *x*-direction, there is a basic shear flow for the ice given by
2.30
where
2.31
is the basal shear stress of the underlying ice flow; ** i** is the unit vector in the

*x*-direction

^{3}. The equations are scaled by choosing 2.32 the film-thickness scale is chosen to balance a presumed upstream water flux

*Q*

_{0}, itself derived from a typical melt rate

*m*over a typical catchment length

*l*

_{i}. If we take

*m*∼1 mm yr

^{−1}and

*l*

_{i}∼1000 km, then

*Q*

_{0}∼10

^{3}m

^{2}yr

^{−1}. The choice of

*h*

_{0}is then 2.33

A number of parameters are introduced in non-dimensionalizing. In particular, we define
2.34
These represent the aspect ratio of the lower surface (*ν*) and the ‘corrugation’ of the lower surface (*σ*). The value of the parameter *θ* depends on the choice of the effective pressure scale *N*_{c}. Previously (Fowler 2009), this was chosen as the prescribed channel effective pressure scale. In the current context, there may be no prescribed *N*_{c}, and it is then natural to choose *N*_{c}=*τ*_{b}, and thus *θ*=1, which we now do. The dimensionless ice depth *h*_{i} is defined in terms of a surface perturbation *H* by writing
2.35
which points out the role of the mean slope of the ice surface.

Central to the solution of the ice-flow problem is the assumption that the aspect ratio is small. This allows us to linearize the geometry to a half space, as well as linearizing the geometric nonlinearities in the stress and velocity components. With this simplification, we find that the ice velocity is uniform at leading order, , and thus the leading order form of the dimensionless sliding law is 2.36 and is defined by the requirement that the average stress be equal to one, i.e. 2.37 where the overbar denotes a space average.

The dimensionless normal compressive stress at the base is written as
2.38
and then the dimensionless hydraulic gradient is
2.39
The surface perturbation *H* satisfies the approximate equation derived from the kinematic condition
2.40
which equates the ice-surface motion to the vertical ice velocity (=*Ξ*), and the quantities *Θ* and *Ξ* are determined in terms of the quantities *H*, *K* and *F*, where the latter two are defined by
2.41
The determination of *Θ* and *Ξ* is described by Fowler (2010). The Fourier transforms of *Θ* and *Ξ* are given by linear combinations of those of *H*, *K* and *F*.

The additional dimensionless parameters arising from the ice-flow solution are 2.42 Typical values are given below using the estimates of dimensional quantities in table 1.

### (h) Stream and sediment flow

The equations for stream and sediment flow, equations (2.23) and (2.24), take the following dimensionless form, taking into account the small aspect ratio. For attached flow where *N*>0,
2.43
where
2.44

For separated flow, where *N*=0 and *h*>0, we have
2.45

The parameters *Γ*, *r*′, *β*, *ε*, *κ*, *δ*, *γ* and *τ** are defined by
2.46

In table 1 we give values of the dimensional constants that appear in the model, and in table 2, we derive the corresponding scales and typical values of the dimensionless parameters. In comparison with the scales of Fowler (2009), the primary difference is the choice of stress scale of 10^{4} Pa, compared with 4×10^{4} Pa. More importantly, we retain the ice viscosity value of 2×10^{13} Pa s used by Fowler (2009) and Schoof (2007*a*). Most easily, this is explained by using Paterson’s (1994) recommended values for Glen’s flow law of *n*=3 and *A*=6.8×10^{−24} Pa^{−3} s^{−1}, together with a stress of 6×10^{4} Pa. However, there is nothing sacrosanct about these values. One might suppose that for a stress of 10^{4} Pa, one should choose a higher viscosity of 0.7×10^{15} Pa s, but it is generally thought that at lower stress, the Glen exponent is reduced (or *A* increases), which would result in no such increase. In addition, it is known that basal ice flow over small slope topography causes an enhancement of stress (Fowler 1981) by a factor inverse to the slope magnitude. For both these reasons, we suppose that the lower viscosity may be the more appropriate one.

### (i) A reduced model

Suppose that *h*=0 everywhere, so that the ice remains adjacent to the till. In a steady state, we must have the water flux ** Q**=

*Q*

_{0}

**, and when this is written in dimensionless terms, it is 2.47 where 2.48 Using the estimates in tables 1 and 2, we find**

*i**Ω*∼1.3×10

^{5}. This large value implies that it is not possible for permeable drainage through till to evacuate subglacial meltwater. Consequently, we suppose that some kind of stream system must always exist.

It is a matter of conjecture what form such a stream system should take. One simple possibility is that the high pressure generated by forcing the water through the till will simply evacuate a sufficient number of fine particles to raise the permeability to the level where it can transmit the water. This would occur through piping in the till. However, it is much more probable that such ‘piping’ would occur at the weak ice–till interface, thus forming a film of water some millimetres deep. This causes a conceptual difficulty, if we think of till as a fine-grained continuum. For one thing, Walder (1982) showed that such a film would be unstable (at least for ice flow over hard beds). This is not in itself a problem, insofar as we are endeavouring to find a comparable instability when the ice flows over a deformable bed, but the typical generalized Weertman sliding law of the dimensionless form
2.49
implies that the shear stress becomes zero when the effective pressure reaches zero, and this raises a difficulty, insofar as we assume that the ice delivers a finite basal stress to the bed. We can resolve this conundrum if we relinquish the idea that the granularity of the till is finer than the water-film depth. If the effect of the high water pressures is to wash out particles of silt size, then we may conceptually suppose that the residual clast fraction of the till is able to provide a resistance to the ice. In this view, the till provides a resistance even when the effective pressure is zero, and the stress only decreases to zero when the water-film thickness reaches a certain value, perhaps of the order of centimetres, depending on the granulometry of the till. At present, we assume that this is the case. This is basically the view of Creyts & Schoof (2009), who have explicitly studied models of water-film flow below ice supported by larger clasts. An alternative point of view would allow the effective permeability of the till to increase dramatically as *N* approaches zero.

In this case, a uniform state consists of a water film of constant thickness above a flat bed, and we can consider the basic model to consist of a separated flow, given by equations (2.45). Now we adopt the limits *ε*→0, *δ*→0, so that the reduced model for stream and sediment flow is (with *N*=0)
2.50

## 3. Stability analysis

We now wish to consider the stability of a uniform state in which there is a constant-thickness water film. By choice of the scale for *h*, the dimensionless water flux upstream is −*h*^{3}**∇***ψ*=** i**, and thus the basic uniform water flow over a flat till bed is given by
3.1
The quantities

*Θ*in equation (2.50) and

*Ξ*in equation (2.40) are determined by the perturbed ice flow:

*Θ*is (minus) the perturbed normal stress at the ice–till interface, while

*Ξ*is the normal velocity of the ice surface. Both quantities are zero in the basic state of the uniform flow. Since it is known (Schoof 2007

*a*) that the ribbing instability is explicitly due to the dependence of the till flux on

*N*, and more particularly through the increasing dependence of the basal shear stress on

*N*in the sliding law, we choose to suppress this dependence in order to highlight the rilling instability present in equation (2.50). In keeping with our discussion in the preceding section, we may do this in the present instance by taking a sliding law of the dimensionless Weertman form 3.2 which explicitly suppresses the dependence on

*N*. It then follows from equation (2.37) that , and from equation (2.41), that

*F*=0. We assume that

*λ*≪1 in equation (2.40), so that

*Ξ*=0. We then find from Fowler (2010) that the Fourier transform is given by 3.3 where the Fourier transform of

*g*(

*x*,

*y*) is defined by 3.4

*k*=|

**|=|(**

*k**k*

_{1},

*k*

_{2})|, and 3.5 the coefficient functions in equation (3.5) are defined by Fowler (2010). Figure 3 shows a graph of

*M*(

*ξ*),

*ξ*=(

*k*/

*σ*).

*M*is positive, with ,

*M*∼(4/

*ξ*) as

*ξ*→0.

To study the stability of the basic state, we write
3.6
and consider *Ψ*, *η* and *s* to be small. Linearization of the model then leads to the equations
3.7
where
3.8
The equations (3.7) are obtained after some algebra, which is omitted here, but is very similar to that involved in studying the Smith–Bretherton overland flow model (Fowler *et al.* 2007). From equations (2.50), (3.3) and (2.41), we have
3.9
where we assume an exponential time dependence, . Substituting this in equation (3.7)_{2}, we find that
3.10
where *r* is the growth rate and *c* is the wave speed, and these are given by
3.11
where
3.12

We now wish to examine equation (3.11) for conditions of instability. We note firstly that if *E*<0, then the steady state is stable. This is, for example, true for the case where the critical Shields stress is zero, *τ**_{c}=0. Secondly, *D* increases and *E* decreases with *k*_{1}, so that transverse (rib-like) wave components are stabilizing. Indeed, purely transverse waves are stable. This is the antithesis of the result from the instability theory when the drainage is considered to be passive. Since longitudinal rolls are the most unstable, we now put *k*_{1}=0; then we have *k*=*k*_{2}, and we find that the wave speed *c*=0, and the growth rate is
3.13
where
3.14

Now let us take the Meyer-Peter and Müller choice in equation (2.45),
3.15
Inserting this in equation (3.14), we find
3.16
and thus *E**>0 and instability is possible if *τ**<*τ*_{0}<2*τ**. More precisely, we find that *E**>(*q*_{0}/*τ*_{0}) if
3.17
In this case, *r*_{||} is positive for 0<*k*<*k*_{⊥}, where the critical wave number is given by
3.18
this defines *k*_{⊥} implicitly (and uniquely), since *M*=*M*(*k*/*σ*), and *k*^{3}*M* is a monotonically increasing function. At *k*=*k*_{⊥}, the growth rate is infinite, and for *k*>*k*_{⊥}, it is negative.

If, on the other hand,
3.19
then *r* is positive for *k*>*k*_{⊥}. In both cases, the most rapid (infinite) growth rate occurs at *k*=*k*_{⊥}, and this then provides the natural expected wavenumber for the patterns we should expect to see.

## 4. Discussion

The equations (2.50) closely resemble those of the Smith–Bretherton model of overland flow (Smith & Bretherton 1972; Fowler *et al.* 2007), which is hardly surprising, since the physics of the two situations are very similar; the only difference in the present case (apart from the supposed laminarity of the water flow, as the Reynolds number is about 30 for a delivered flux of 10^{3} m^{2} yr^{−1}) being the effect of the overlying ice on the till elevation through the presence of the term *Θ* in equation (2.50).

The mechanism of instability closely resembles that of the Smith–Bretherton model, and is due to the positive feedback between water flow, sediment transport and erosion. Water-flow depth will be larger in an incipient stream, and this allows faster flow and thus higher erosion and sediment transport. Consequently, the stream deepens further, and this positive feedback causes the instability.

Mathematically, the instability arises from the relation for *η*_{x} in equation (3.7)_{1} that leads to the negative diffusion term −∇^{2}*ψ* in equation (3.7)_{2}. In the Smith–Bretherton model, the corresponding term is −∇^{2}*s*, and thus the instability when it occurs has a growth rate , and is unbounded at large wavenumbers. This ill-posedness in the rill model is regularized by the inclusion of the term corresponding to that in *δ* in equation (2.45), and this causes the maximally growing wavelength to be small but finite. In the present case, the ice acts as a dampener, causing large-wavenumber growth to be suppressed, and the infinite rapid growth is shifted to a finite wavenumber.

Inclusion of the term in *δ* in equation (2.45) modifies the stability analysis straightforwardly. Since *δ*≪1, we find that the term only becomes significant at small longitudinal wave numbers, *k*_{1}∼*δ*. In this case, we still have *k*≈*k*_{2}, and the expression for the growth rate is modified to
4.1
where
4.2
Figure 4 shows the variation of *r* with *k*_{2} for *k*_{1}=0 and *k*_{1}>0. We can see from equation (4.2) that the presence of a non-zero longitudinal wavenumber causes a stabilization of the system, and the removal of the singularity in equation (3.13). Although this analysis is only linear, it is natural to suppose that the preferred lateral wavelength will be *k*_{⊥}, while we expect that the typical longitudinal length scale of the lineations will be such that *L*∼1, and thus
4.3

We can now define the expected length scales of the lineations by 2*πl*/*k*, where *k* is the relevant wavenumber and *l* is the length scale. This gives us expected widths of
4.4
and typical longitudinal length scales of
4.5
Although it requires nonlinear calculation, the expected elevation is of *O*(*d*_{D}).

To estimate these, we suppose that in the unstable range *τ*_{0}≈(3/2)*τ**, so that
4.6
and from figure 3, we select the low *σ* limit *M*=2. Using the other parameter estimates in table 2, we then have the preferred width
4.7
while the preferred length scale is
4.8
assuming a nominal value of *L*=0.1. The elevation scale is *d*_{D}=12.3 m.

These values appear consistent with observation. It remains to assess the likelihood of instability occurring, that is, whether *τ*_{0}>*τ**. From table 2, we estimate *τ**=1.13, whereas *τ*_{0}=*σ*^{2/3}*γ*≈0.3. Uncertainties in parameter estimation mean that the fact that *τ*_{0} is comparable to *τ** means that instability will plausibly occur. The loosest choice of parameters in table 1 are perhaps the grain size *D*_{s}=100*μ* and the water flux *Q*_{0}=10^{3} m^{2} yr^{−1}. A reduction of grain size to 25*μ* increases *γ* and thus *τ*_{0} beyond the instability threshold. In fact, since sufficient water flow for sediment transport is the essential condition for instability, this is virtually guaranteed, since, as discussed above, the formation of a water film at the till surface essentially requires the removal of the finest size particles. In addition, the possible propensity for discharge to be episodic rather than uniform suggests that the effective value of *τ** may be lower than indicated here, where the value is that appropriate for a constant discharge.

## 5. Conclusions

We have generalized the instability theory of drumlin formation to allow for a dynamic hydraulic regime, rather than the passive, ‘small-stream’ version that has been studied thus far. This generalization may also be seen as a generalization of the Smith–Bretherton model of subaerial overland flow to include the effects of an overlying viscous substrate. We find that the rilling instability still occurs, but the effect of the overlying ice is to dampen the growth at very short wavelengths, and the resulting maximal growth rates occur at wavelengths that are comparable to those observed in MSGLs. For the particular choice of parameters used here, the width is 394 m. In our presentation, we explicitly exclude the mechanism that is responsible for the growth of ribbed moraine, that is, the dependence of the sliding velocity on the effective pressure. In this case, we find the most unstable modes are longitudinal rolls (lineations), but that downstream variations become important at very long length scales, as much as 50 km.

We anticipate that when the ribbing instability is put back in the model, and both attached and separated regions are treated consistently, that genuine three-dimensional instabilities will occur, and also that their finite-amplitude expressions will reveal the whole spectrum of behaviour, from ribbed moraine to drumlins to MSGLs, and that the stream/cavity system will be an integral part of the resulting flow. This anticipation awaits future work.

## Acknowledgements

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

## Footnotes

↵1 Since effective pressure increases with depth, the yield stress, which we take to be

*μp*_{e}, where*μ*is the coefficient of friction, is only exceeded for a finite depth.↵2 We distinguish the basal shear stress

*τ*, which is the shear stress at the ice–till interface, from what we might call the ice-sheet basal shear stress, defined below in equation (2.34), which is the shear stress delivered by the ice sheet to the basal region. The two would be equal if the bed was flat.↵3 It should be noted that although this does not exactly satisfy the top-surface boundary conditions, it is an accurate approximation for small surface slopes

*S*.

- Received January 7, 2010.
- Accepted March 22, 2010.

- © 2010 The Royal Society