## Abstract

We present a theory for the coupled flow of ice, subglacial water and subglacial sediment, which is designed to represent the processes which occur at the bed of an ice sheet. The ice is assumed to flow as a Newtonian viscous fluid, the water can flow between the till and the ice as a thin film, which may thicken to form streams or cavities, and the till is assumed to be transported, either through shearing by the ice, squeezing by pressure gradients in the till, or by fluvial sediment transport processes in streams or cavities. In previous studies, it was shown that the dependence of ice sliding velocity on effective pressure provided a mechanism for the generation of bedforms resembling ribbed moraine, while the dependence of fluvial sediment transport on water film depth provides a mechanism for the generation of bedforms resembling mega-scale glacial lineations. Here, we combine these two processes in a single model, and show that, depending largely on the granulometry of the till, instability can occur in a range of types which range from ribbed moraine through three-dimensional drumlins to mega-scale glacial lineations.

## 1. Introduction

Drumlins are small, rounded hills which occur in swarms. On the ground they give rise to a typically rolling landscape, which can, for example, be found in large parts of Ireland (where they were first distinctively observed). Where they become inundated, they are particularly remarkable, as the individual hills then stand out, as shown in figure 1, which is a view from Westport harbour in Co. Mayo of the drumlins of Clew Bay.

The aspect of pattern which is evident on the ground in many drumlin fields is particularly notable in satellite photographs, or even better in digital elevation models such as that in figure 2, which shows that the drumlins of north central Ireland blanket the entire landscape.

The origin of these landforms has been a subject of discussion for well over a hundred years. Early investigators such as Bryce [1] correctly perceived of drumlins as being formed by the transport of sediments. Bryce was able to determine, by examination of the kinds of rocks in the till which formed drumlins in Northern Ireland, that they had formed through the action of a current which had swept down from the northwest. Those being diluvial times, Bryce assumed this was a current of water. It was not until Agassiz's glacial theory that opinion changed, and Kinahan & Close [2], describing the drumlins in the vicinity of Clew Bay, associated them squarely with the flow of the great ice sheets of the ice age. Over the following century, the debate on origin largely centred on the alternative concepts of erosional or depositional origin (e.g. [3]), but this period is remarkable by its absence of any quantitative theory. It was with the work by Smalley & Unwin [4], Menzies [5] and Boulton [6] that the role of deforming till was explicitly considered; but still there was no descriptive theory. Only with the landmark paper of Hindmarsh [7] was the interaction between ice flow and deformable sediment formulated mathematically, with the result that instability could be predicted to occur with wavelengths of several hundred metres, corresponding to those which are observed.

Hindmarsh's numerical study was followed by an analytic study by Fowler [8], which yielded similar conclusions. This promising beginning has since run into theoretical obstacles, as described in the studies by Schoof [9,10]. There are three issues that arise. The first is a practical one: as the bedforms grow, they soon reach an amplitude where the ice separates from the till bed, forming lee-side cavities. Cavitation is not in itself a problem, but it raises severe difficulties in continuing the numerical evolution of the bedforms. Schoof circumvented this by seeking travelling wave solutions, and suggested that the resulting bedforms would be of insignificant elevation compared to the tens of metres which can be obtained.

A second difficulty is that the instability described by Hindmarsh [7] and Fowler [8] is fundamentally two-dimensional in nature. The difficulty here is that drumlins are three-dimensional objects, although it is not at all necessary that the initial instability be three-dimensional. Despite this, the search for three-dimensional waveforms has been elusive. Schoof [9] thought that the two-dimensional nature of the instability represented a serious obstacle to the theory, a view taken on trust by Pelletier [11]. Consideration of this issue represents one of the purposes of the present paper.

A third consideration concerns the stratigraphy of drumlins. While many consist of an essentially homogeneous mass of till, there are certain fields of drumlins which provide evidence of stratification, and at first sight this seems incompatible with a formation mechanism which involves deformation of the till; the point has been emphasized by Hart [12, p. 194] and Hiemstra *et al*. [13], for example. A little thought shows that this is not a real objection in net-erosive environments, and indeed most of the principal stratigraphic constituents of drumlins can be plausibly explained [14].

The first difficulty alluded to above, and in effect also the third, was addressed by Fowler [15], who devised a method to compute the evolution of drumlins after cavitation occurred. He incidentally found that the resulting waveforms became stationary, thus allowing the erosional mechanism of Stokes *et al*. [14] to provide an explanation for layered stratigraphy. The present paper addresses the second issue discussed above, that of the formation of three-dimensional bedforms.

It is a long-standing concept in theories of subglacial bedform genesis (e.g. [16]) that drumlins form part of a continuum traversed as some critical parameter varies, between the extremes of ribbed moraine [17] and mega-scale glacial lineations (MSGLs) [18]. This attractive idea is predicated on the notion that one single mechanism can encompass the whole variety of bedforms, but it is by no means obvious that this is the case. Indeed, apart from the Hindmarsh–Fowler instability theory, a number of alternative explanations have been put forward to explain ribbed moraine, drumlins and MSGLs (e.g. [19–26]). An early review of drumlin formation theories is that of Menzies [27], and a more recent one is that of Clark [28].

Schoof's [9] comments drove Dunlop *et al*. [29] to the defensive position of advocating the instability theory only in connection with ribbed moraine formation, with the implicit recognition that ribbed moraine and drumlins were somehow different, and this informal distinction has gained a certain currency. It is not, however, one that we espouse, and indeed the principal purpose of the present paper is to show that a suitable model describing combined ice, sediment and water flow can indeed plausibly explain the whole spectrum of bedforms, from ribbed moraine through drumlins to MSGLs. Before we commence this, however, we need to describe and distinguish the various components of the instability theory as it has been developed thus far.

Hindmarsh [7], Fowler [8] and Schoof [9] all developed a theory based on the Stokes flow of ice near the bed (i.e. the ice surface is far away, so that the ice thickness is effectively infinite), with the ice assumed to be Newtonian viscous, and an underlying motion of water-saturated deformable till, described by an Exner equation for the ice/till interface, in which the till flux *q* is assumed to depend on basal shear stress *τ* and effective pressure *N*. The effective pressure is defined by
*p*_{i} is ice pressure, which can be assumed known, and *p*_{w} is the basal water pressure at the ice/till interface, which must be determined. Hindmarsh, Fowler and Schoof all assume this water pressure is in hydrostatic equilibrium, presumably with a subglacial drainage system at a (locally) constant effective pressure.

Fowler [30] extends the theory to ice of finite depth and also treats three-dimensional perturbations, something also done by Schoof [9]. An evolution of the theory is offered by Fowler [31], who finally considers the water flow as an active contributor to the model. That edition of the theory concerns MSGLs, which are found via an instability analogous to the Smith–Bretherton theory of sub-aerial rivulet formation [32]. Fowler describes water flow either as a groundwater flow through porous till, or as a separated Poiseuille-like flow. When the water layer is absent, till deforms as before, but when the water layer is present, sediment transport occurs as interfacial bedload, and it is this latter that drives the rill-like instability.

In order to formulate a meaningful instability problem, Fowler adopted the stance of Creyts & Schoof [33] that a possible subglacial drainage mechanism is that of a thin film trickling at the base of a glacier. The concept circumvents the Walder [34] instability problem by allowing clast sizes larger than the film thickness. In retrospect, this seems the obvious basic flow, which takes account of the roughness of the bed. In this paper, we adopt this viewpoint from the outset, thus allowing a sediment layer and a water layer to exist everywhere under the ice. The Fowler [31] MSGL model is thus modified by replacing through-till porous flow by a distributed film flow. The hydraulic potential in this water layer is to be determined, and in particular, the stream system becomes part of the problem to be solved; and, in fact, there is no theoretical distinction to be made between stream and cavity; they are both simply regions of separated flow.

The structure of the remainder of the paper is as follows. In §2, we re-expound the Hindmarsh instability theory in its present form. Our presentation is a compact summary of the successive constituent sub-models of ice flow, water flow, sediment transport and ice closure, the last of these being a new development consequent to our assumption of a Creyts–Schoof water film. These sub-models are written dimensionlessly, dimensionless parameters are estimated, and an asymptotically simplified model is then presented, combining the four sub-models.

Section 3 provides a uniform steady solution, and its linear stability is analysed; it is shown that the growth rate is determined by the solution of a quadratic equation, whose roots are plotted in terms of horizontal and transverse wavenumbers as functions of a critical parameter *Π*. Section 4 forms the conclusion of the study, in which also the meaning and interpretation of the results are discussed.

## 2. The instability theory

The Hindmarsh instability theory is built around a description of ice, water and sediment flow on a local horizontal scale of kilometres. The geometry of the flow is shown in figure 3. The coordinates are *x*, *y*, *z*, with *z* vertical, and *x* in the downstream ice flow direction. The ice surface is at *z*=*z*_{i}, the ice base and water surface at *z*=*s*, and the water base and sediment surface is at *z*=*b*. The water depth is thus

### (a) Ice flow

The Stokes flow equations for the flow of ice, treated as a Newtonian viscous fluid, are given by (2.5) of Fowler [30], and are
*P* is reduced pressure (ice pressure minus cryostatic pressure), **u** is ice velocity, *ρ*_{i} is ice density, *g* is the acceleration due to gravity and *η*_{i} is the assumed constant ice viscosity. Defining the horizontal shear stress vector as
*τ*_{1}=*τ*_{13} and *τ*_{2}=*τ*_{23} are the two horizontal components of shear stress), the surface boundary conditions are the stress balance and kinematic conditions
*t*, *x* and *y* denote partial derivatives, *a* is the accumulation rate, and **u**=(*u*,*v*,*w*). The basal boundary conditions are then the sliding law and the kinematic boundary condition (ignoring the relatively small melting term), thus
**u**_{b}=(*u*,*v*), evaluated at *z*=*s*, and called the sliding velocity; *N* is the effective pressure introduced in (1.1). To be specific, we will assume a generalized Weertman sliding law of the form
*R* is a roughness coefficient [36].

The Stokes equations require two boundary conditions at top and bottom, so that the kinematic condition in (2.4) determines *z*_{i}, and we see from (2.5) that if *s*, *s*_{t} and *N* are known then the ice flow is fully determined. Before proceeding to water and sediment flow, we complete an approximate solution for the ice flow, based on the assumption of small basal ice slope, |**∇***s*|≪1, which is an accurate assumption.

#### (i) Non-dimensionalization

Following Fowler [30], we define the bedform elevation scale *d*_{D}, the till deformation scale *d*_{T} and the bedform length scale *l*_{D}, by
*u*_{0} and *τ*_{b} are velocity and stress scales (defined below), *ϕ* is sediment porosity, and
*ρ*_{s} being the sediment particle density and *ρ*_{w} the water density. This is equivalent to Fowler [30] if his effective pressure scale *N*_{c}=*τ*_{b}.

The aim of the instability theory is to study perturbations about a uniform state, which we take to be the shear profile
*d*_{i} is the ice depth; obviously, this varies over the scale of the ice sheet but can be taken as a constant on the much smaller bedform length scale *l*_{D}. The values of *u*_{0} and the downstream ice surface slope *S* (or *τ*_{b}) are determined from the ice sheet model at a much larger length scale and can be taken to be given in the present context.

We now scale the ice flow equations by writing
*ν* and corrugation parameter *σ* are defined, respectively, by
**u*** and **x*** are omitted in the sequel. The dimensionless mean basal velocity *R** is defined by

For small *ν*, the bed conditions at *z*=*νs* can be linearized to apply at *z*=0, and because *S*≪1, the ice surface conditions can be linearized to apply at *z*=1. The perturbed ice flow thus satisfies the biharmonic equation in the strip 0<*z*<1 and can be solved by Fourier transform using the prescribed stress conditions at the ice surface, together with the sliding law and the normal stress at the bed. Specifically, if we define
*Φ*=−*Θ*), then the kinematic conditions at the bed and the ice surface take the dimensionless form
*w*=*K*), where
*w* and *Ξ* can be written as linear combinations of *H*, *F* and *Φ*; additionally, **k**=(*k*_{1},*k*_{2}), and *k*_{1} and *k*_{2} are downstream and transverse wavenumbers, respectively, we obtain
^{−2} is small, we can put *Ξ*=0, and simplification of (2.21) yields an expression for *s* of the form
*J* and *L* are, respectively,

### (b) Water flow

In dimensional terms, the hydraulic potential can be shown to be
*ψ*; this is possible, but as we shall see, unlikely.

Slow water flow through a thin film or stream of depth *h* is given by a local Poiseuille flow, and takes the form
*η*_{w} is water viscosity and *Γ* is the local water source due to internal melting in the basal ice, and is given approximately by
*G* is geothermal heat flux, *q*_{T} is heat flux into the ice and *L* is latent heat. The dependence of *Γ* on ice flow may have a dynamical effect on a large scale [37], but in fact *Γ* is negligible on the small length scale of concern here.

We scale these equations using values in (2.11), and in addition
*h*_{0} is a typical size of the water film depth, prescribed below; this leads us to the dimensionless forms,
*q*_{w} of the basal water flux, given by
*l*_{i} is a relevant horizontal length scale for the large-scale ice flow. Equation (2.31) thus defines a mean film thickness
*t*_{D} was defined in (2.11).

### (c) Sediment flow

Next, we write a suitable Exner equation (describing conservation of the mobile sediment) for the evolution of the sediment surface *z*=*b*. This amalgamates the till flux terms described in the original instability theory, e.g. Fowler [30], and the sediment transport terms appropriate for fluvial flow described by Fowler [31]. We take it in the form
*h*_{A} is the thickness of the deforming till, which is limited by the yield stress *μN*, thus [15]
*μ* is the coefficient of friction; here *η*_{s}, corresponding to a yield stress material with the normal flow rule, and as always, this does not imply the till is viscous, but the application of more probable rheologies (e.g. [38]) is an uncertain distraction to the present purpose.

The terms on the right represent net suspended sediment erosion (*E*), which we take to be zero, indicating a balance between suspended sediment erosion and deposition. The term inside square brackets represents bedload sediment transport *Q* and is generally non-zero and increasing downstream, as the glacier strips its bed away. The second equation defines the effective stress *τ*_{e} transmitted by the water flow to the bed, consisting of the actual viscous stress in the water plus a term representing the propensity of sediment grains to roll downhill: *D*_{s} is the sediment grain size. *b* in this equation was written as *s* in eqn (2.19) of Fowler [31], on the basis that *h* is small; (2.34) is the correct form, however; cf. Fowler *et al*. [39], eqn (2.3). The bedload *Q* will generally be a nonlinear function of this effective stress, but the slow conditions underneath a glacier do not necessarily suggest the use of a sub-aerial correlation. However, Fowler [31] found that the streaming instability relied on the existence of a yield stress, something that is essentially also necessary in the sub-aerial case. Our more general results here belie this conclusion.

We non-dimensionalize the equations as previously, with the extra choices
*q*_{b} is a typical value of bedload sediment flux. It follows that the dimensionless sediment model is

### (d) Ice closure

Counting equations, we see that (since *b*=*s*−*δh*) (2.30) and (2.38) provide two equations for the three unknowns *N*, *s* and *h*, bearing in mind that from (2.30), (2.21), (2.16) and (2.15), *ψ*=*ψ*(*N*,*s*,*s*_{t}). An extra equation is necessary to close the system. In the classical stability theory [7,8], this is taken to be the condition of constant hydraulic potential, *ψ*=−*N*_{c}, but in the present context, that is not enough, because we have explicitly introduced the water layer thickness *h*. Fowler [15] accommodated *h*>0 (the formation of cavities) with the extra statement:

In the present case, we need an extra equation, because we allow non-zero film thickness even when *N*>0. This extra condition is the film closure condition, and can be proposed by analogy with the corresponding condition in Röthlisberger's [40] drainage theory. The necessity for an extra condition occurs because the presence of a partially clast-supported film implies that at the clast separation scale (say metres), a uniform normal (effective) stress in the ice above the bed is partitioned between supporting clasts and the intervening film where it is zero. Thus [33] the ice will sink into the film at a rate (for a Glen law model) *Kl*_{c}*N*^{n}, where *K* would be proportional to the Glen flow law constant, and *l*_{c} is clast spacing. This closure is offset by the production of melt, which thus leads to the closure model
*l*_{c} is an increasing function of *h*, which we may suppose tends to infinity at a critical value *h*_{c}, which properly defines a stream (or a cavity), where *h*>*h*_{c}. It may be perverse to assume a nonlinear rheology here when ice is considered Newtonian elsewhere, and thus it may be preferable instead to pose the closure condition

To non-dimensionalize this, we define a typical clast spacing when *h*∼*h*_{0} to be *l**_{c}, and then we write
*h* is dimensionless and *N**∼*O*(1), so that (2.46) becomes in dimensionless terms
*N**(*h*) is largely one of convenience. We should have *N**≈0 for *h*≫1, *N**∼*O*(1) for *h*∼1 and *h*=0. In fact, it is best to have *h*→0, since in reality the small neglected melting term in (2.30)_{2} precludes the film completely disappearing as long as *ω*>0. A simple choice which satisfies these constraints is

### (e) Parameter estimation

The values of the constants which are assumed in the model are given in table 1. Some of these are typical values, and some are intelligent guesses, so that the constant values of the scales in tables 2 and 3 are also representative.

The value of the bedload transport *q*_{b} is estimated as follows. Values of *q*_{b}/*q*_{w} for a turbulent Alpine proglacial stream are of order 10^{−4} [41]. We might suppose a similar or lower fraction for the (mostly) slower flow under an ice sheet. The assumption

The simplest way to estimate *D*_{s} is to determine it by the critical Shields stress for sediment transport [35, p. 273]. The Shields stress is given by
*τ**>*τ**_{c}, where *τ**_{c} has a typical value of 0.06, but it depends on the particle Reynolds number
*D*_{s} based on *τ**=*τ**_{c}≈0.06, we find *D*_{s}≈43 μm, but then *Re*_{p}∼0.14, and at such low values, *τ**_{c}≈0.1/*Re*_{p}, which indicates sediment transport if *τ*_{w} is just 0.04 Pa.

The simple conclusion is that sediment transport need not occur in the millimetric water film, but rather in the cavities which are formed behind growing ribbed moraine. This is essentially the picture portrayed by Fowler [15] and Chapwanya *et al*. [42] and is consistent with an absence of streams. In this case, a suitable value of *γ* is zero, and the specification of *D*_{s} is irrelevant. On the other hand, we know that streams do occur, and indeed that floods occur between subglacial lakes (e.g. [43,44]), so some mechanism to form streams must exist. This can be due either to the Walder [34] thermal mechanism, or due to the pre-emptive formation of cavities and their evolution to a linked cavity system. However, there is also another mechanism which does not rely on fluvial transport. A local thickening of the water film causes a reduction of *N*, and thus (in a developed rib system, [15]) an increase in till deformation depth and thus till flux. This also appears to be a destabilizing mechanism. Fowler's [31] stream formation mechanism relied on bedload transport, but in that paper the sliding law dependence was explicitly excluded. In either event, we will choose a value of grain size commensurate with transport by a developed stream of depth *h*_{s}=20 cm. This gives a water stress of approximately 1.3 Pa, and the Shields criterion suggests *D*_{s}=80 μm.

The choice of cooling rate *q*_{T} is based on a conductive heat flux *k*Δ*T*/*d*_{i}, enhanced by horizontal advection; the value in table 1 is enhanced by 1.5. The viscosity *η*_{i} would have a normal value of 2×10^{13} Pa s at a shear stress of 10^{5} Pa, but two orders of magnitude higher at the lower stress of 10^{4} Pa; however, this is offset by the enhancement of shear stress induced by flow over the induced bedforms, so we choose an intermediate value. The ‘viscosity’ *η*_{s} of till is an effective viscosity, if we suppose a normal flow rule. Our estimate is based on the definition

The values in table 2 are computed from expressions given in the text. The choice of typical clast spacing *Π*≈1. In reality, it should depend on the granulometry of the till, being small for fine-grained sediments, and larger for coarser till. However, it should also depend on the small-scale roughness of the bed, since if streams form, it would seem that *l**_{c} should be comparable to stream width.

### (f) A reduced model

The model consists of equations (2.15), (2.18), (2.16), (2.21), (2.30), (2.37), (2.38), (2.42) and (2.48), which determine *b*, *h*, *s*, *ψ*, *N* together with the subsidiary variables *A*, *Φ*, *Ξ*, *w*, *F*, *H*. For simplicity, we define
*ω*, we obtain a reduced model in the form
*β* and *γ* to facilitate the discussion below.

This reduced model can be sensibly compared with that of Fowler [15]. It is in fact the same if *σ*=0 (the apparently inconsequential limit of deep ice) whence we can take *ψ*=−1 as constant, and neglect *β*, *γ* and *δ*. Consulting (2.59), it seems this is a perfectly good solution. It allows water to trample roughshod over hills, the Creyts–Schoof water film allowing easy passageway.

It can also be compared to Fowler's [31] analysis of stream and MSGL formation, which takes *δ*=*A*=0, and replaces (2.51) with the assumption *N*=0, which would follow from the assumption that *Π*≫1. It is reasonable to anticipate that the model without both sets of simplifications will admit both forms of the instability.

## 3. Linear stability

The steady uniform solution of (2.59) is given by
*h* and *b* are not fixed by the reduced equations. In practice, therefore, it is best to fix the uniform state by means of two specific upstream boundary conditions, for example for two of *h*, *Ψ* and *b*.

Linear stability analysis proceeds along the same lines as in Fowler [15,31]. We write
*A*, *A*′=*A*′(*N*) and *f*_{N}=∂*f*/∂*N* are evaluated at the basic state (3.1), *κ* was defined in (2.41), and
*Q* is the dimensionless bedload flux, evaluated at the base state *τ*_{e}=*σ*. Defining the Fourier transform *Σ*:
*Σ*_{+}.

Simplification of (3.9) follows from neglecting the small coefficients *β*, *κ* (thus also *T*) and *δ*; then the quadratic takes the simpler form:

Because also *ε*≪1, there is a fast root (which is *Σ*_{−}) of *O*(1/*ε*), which is easily found to be stable (it represents the rapid relaxation of the film closure equation to equilibrium), and a slow root, which for *ε*→0, is given by
*σΠ*=0 and *U*=1 and *V* =0), under the identification
*D* in the 2010 paper is absent as it arises from the neglected terms in *T* in (3.10).)

A question of some interest is how the maximal growth rate depends on *k*_{1} and *k*_{2}. Remember, this was one of the issues concerning the instability theory: there seemed to be no genuine three-dimensional instabilities. Figure 4 shows contour maps of growth rate in conditions of instability, when *Π*. Specifically, transverse ridges (ribbed moraine) are preferred for small values of *Π*, whereas longitudinal ridges (lineations) are preferred for large *Π*; and for intermediate values, the maximal growth rate occurs for genuinely three-dimensional drumlinoid shapes.

A comment on these results is necessary. Firstly, while the approximate result (3.14) is indeed valid as *ε*→0, it is actually not accurate for our estimated value *ε*=2.7×10^{−2}, apparently because the maximal growth rates occur for *εk*^{2} is not small). With the selected parameter values, accuracy ensues for *Π* make the approximate result more robust. One noteworthy comment is that the MSGL instability criterion of Fowler [31] required *C*>3*B* in order that *E*>0. In the present combined theory, this criterion is not necessary. Indeed, approximation of the values of *Σ*_{+} taking *β*=*γ*=*δ*=*κ*=0 has minimal effect on the results.

Plotting of these results in terms of length and width (defined as half-periods) is of potential use in comparison with observed data, although it must be borne in mind that linear stability results do not necessarily bear close relation to observed finite amplitude results. An example of such a plot is in figure 5, where the regions of maximal growth lie close to observed ranges of actual bedforms [46–48], though the mega-ribs of Greenwood & Kleman [49] and the similar-sized traction ribs of Sergienko & Hindmarsh [50] are somewhat larger. In particular, the suggestion of our theory that ribs, drumlins and lineations all form through a single process appears to be consistent with observed length–width distributions of these various bedforms, which appear to form overlapping clusters in parameter space [48].

Of possible interest, or possible concern, is the retention of the ribbing mode as *Π* increases, although it becomes weaker. By comparison with figure 4, this may be associated with the fact that *Re* *Σ*_{+} decreases slowly to zero as *k* increases, and does not become negative at a large wavenumber. We suspect this may be associated with oversimplified physics, which does not provide a stabilizing mechanism at a large wavenumber, and we expect this may be associated with an over-simple closure relation, which lacks any spatial derivatives.

However, we should also point out the danger of trying to over-interpret figure 5. It is merely a linear stability diagram and can only give hints of finite amplitude observations. Ribbed moraine has a fairly definite preferred width to length ratio of about three to four [48]. Infinitely wide ribs do not occur, and this may be due to the necessity to provide a passageway for the downstream flow of water.

## 4. Conclusion

In this paper, we have provided a theoretical model for the coupled evolution of ice, water and sediment at the base of an ice sheet. The theory is a development of earlier more rudimentary versions, and the principal theoretical development presented here lies in the consideration of the water in a Creyts–Schoof water film. This concept securely allows a water film at the base without the associated unrealities of flotation and vanishing stress.

The model essentially consists of four partial differential equations, representing evolution of basal ice surface (2.24), conservation of water mass (2.30), conservation of sediment (2.38) and ice closure (2.48). The first three of these are securely based on physical principles, but the provenance of the closure equation is less clear, despite its familiarity in the case of Röthlisberger channels, since it seems that the ice surface kinematic equation (2.24) ought to represent the same physical process. The argument that both equations are necessary is based on scale: the closure equation represents ice motion at the microscale of the water film, while the kinematic equation represents the macroscopic ice motion associated with flow over the bed. The distinction is similar to the scale separation used in determining the sliding law [51].

The earliest versions of the instability theory assumed a passive, hydrostatic water system. In this latest edition of the model, we have included a fully integrated water system in the form of a Creyts–Schoof film, with the included ability to describe both streams and cavities. Analysis of the instability of the uniform state of the system shows that the preferred mode of spatial growth (and thus what might correspond to the actual observed pattern) ranges from ribbed moraine pattern through drumlin pattern to glacial lineation pattern as a single critical dimensionless parameter *Π* increases from 0.7 to 10. In fact, consultation of (3.10) suggests (since *δ* is small) that it is the product *σΠ* which is important.

The fact that this parameter is a product suggests a complicated combination of stability-determining circumstances. In particular, the fact that *σ* is the ratio of the drumlin length scale *l*_{D} to the ice depth scale *d*_{i} suggests that, in contradiction to earlier results (e.g. [15]), ice depth is an important quantity. Actually, this appears not to be the case, and it is better to view *σΠ* as a single controlling parameter. From (2.7), (2.12), (2.49) and (2.10), it is defined as
*Γ* is defined in (2.28), and we see that it is in fact independent of ice depth. A simple rescaling in (2.10) and (2.11) using *l*_{D} instead of *d*_{i} would lead to the parameter defined above appearing in the closure equation.

The principal dependence of *σΠ* is on the clast spacing *l**_{c}, and if we had to make a prediction based on this definition, it would be that fine-grained sediments (large *l**_{c}) promote ribs or drumlins. The dependence on ice velocity is relatively weak; we might take *Γ*∼*τ*_{b}*u*_{0} and *σΠ* is rather small.

The dependence of the results of this analysis on the clast spacing is not only intriguing but also frustrating. The closure equation is the least satisfactory part of the model, and a sensible definition of *l**_{c} necessarily hard to tie down. Earlier we chose *l**_{c}=1.2 m, ostensibly so that *Π*∼1. A more intelligent choice would be such that the dimensional equilibrium value *N*=*τ*_{b}/*Π* corresponds to observed values, commonly thought to be 0.5 bar =5×10^{4} Pa beneath Whillans Ice Stream [52]. Using the value of *τ*_{b} in table 2, this suggests *Π*=0.18, corresponding to a choice *l**_{c}=0.21 *m*. Note, however, that we require *τ*_{b}>*μN* is the till is to deform, corresponding to *Π*>*μ*. The issue is clouded in practice by the observation that much of the resistance in Antarctic ice streams is provided by lateral drag, as opposed to the basal drag as assumed here [53].

This discussion is not very encouraging, because it suggests (figure 4) that high values of

So, we think that the role of the clast spacing, and the provision of the closure equation, will require an increasing sophistication in future development of this theory; in particular, while the present form may be reasonable for ice flow over a macroscopically flat bed, we suspect that a more elaborate version will be necessary in the presence of streams or bed topography. One possible development may involve the prescription of *l*_{c} as a non-local function of film thickness *h*, to allow for the situation where full stream separation occurs. A related point may be the slow decay of growth rate in figure 4. Using the approximation (3.14), the growth rate is
*Re* *Σ*∼*O*(1/*k*) as *Re* *Σ*∼−*k*^{2}, and it is possible that more sophisticated closure models will allow this. In particular, it is likely that this slow decay is a contributor to our present problems in providing numerical solutions of the model.

## Funding statement

A.C.F. acknowledges the support of the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie) funded by the Science Foundation Ireland PI grant 12/1A/1683.

## Acknowledgements

For continuing and enduring fruitful discussions, our thanks to Chris Clark, Chris Stokes, Matteo Spagnolo, Paul Dunlop and Richard Hindmarsh.

- Received March 7, 2014.
- Accepted August 8, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.