## Abstract

It is proposed that the formation of the subglacial bedforms known as drumlins occurs through an instability associated with the flow of ice over a wet deformable till. We pose a mathematical model that describes this instability, and we solve a simplified version of the model numerically in order to establish the form of finite-amplitude two-dimensional waveforms. A feature of the solutions is that cavities frequently form downstream of the bedforms; we allow the model to cater for this possibility and we provide an efficient numerical method to solve the resulting free boundary problem.

## 1. Introduction

Drumlins are small ovoid hills that form ubiquitously under ice sheets and glaciers. They litter the landscape of North America, Great Britain, Northern Europe and other formerly glaciated areas. It is only recently, with the advent of high-quality digital elevation models (DEMs), that the extent of their coverage has become apparent, and such terrains cover an estimated 70, 50, 40 and 15 per cent of the areas of Canada, Ireland, Scandinavia and Great Britain, respectively (Clark *et al*. 2009). Figure 1 shows a view of the drumlins of Clew Bay near Westport in Ireland, whereas figure 2 shows a DEM of the drumlins of County Clare, Ireland.

Although drumlins have been described for well over 100 years (Kinahan & Close 1872), the cause of their formation has escaped quantitative explanation for much of that time. In a landmark paper, Hindmarsh (1998) proposed a mathematical model to explain such bedforms. In his theory, ice flows over a layer of wet, deformable till, and the resultant shear flow in the ice and the till causes an instability to occur at the interface, much in the way that air flow over water causes water waves or (a better analogy) water (or air) flow over sand causes fluvial (or aeolian) dunes to be formed. Fowler (2000) addressed essentially the same theory, but solved the linear stability problem associated with perturbations of the uniform flat state analytically. He found, as did Hindmarsh, that instability occurred commonly for realistic kinds of assumptions concerning the rheology of till. More recently, the theory has been extended by Schoof (2007*a*,*b*), who in the latter paper addresses the difficult issue of cavitation.

This promising beginning to a comprehensive theory of subglacial bedforms has run into two obstacles. Although the theory appears to predict waveforms of the correct wavelength, there are a number of reasons to doubt its veracity; we discuss these in more detail below. There are three apparent problems in the theory: only two-dimensional instabilities are predicted; the theory seems unable to explain evidence of internal stratification; and inappropriate amplitudes are apparently predicted because cavities are formed at low elevation, and these may limit further drumlin growth (Schoof 2007*b*).

The other obstacle arises because of the use by both Hindmarsh and Fowler (and also Schoof) of a Boulton–Hindmarsh rheology for till, that is to say, one with a power law relating strain rate to stress τ, in the form
1.1
where *p*_{e} is the effective pressure, equal to the difference between overburden pressure and pore water pressure, and *A*, *a* and *b* are positive material constants. The apparent observations of Geoffrey Boulton at Breidamerkurjökull in 1977, reported by Boulton & Hindmarsh (1987), which motivated the choice of this rheology, were seemingly discredited by later laboratory experiments of Kamb (1991). Working with samples of till from West Antarctica, Kamb showed that till behaved approximately plastically, with deformation occurring essentially without limit at a particular yield stress. The disagreement between these two points of view has led to controversy (Tulaczyk *et al*. 2000; Iverson & Iverson 2001; Fowler 2003) and has, to some extent, obscured the development of the drumlin instability theory. This is unfortunate because the theory does not at all rely on details of the till rheology, but assumes only that there is basal ice motion (sliding), that till deforms, is thus transported and that the rate of ice slip and till transport depends on both basal stress τ and on the effective pressure at the ice–till interface *N*. This point was made clear by Schoof (2007*a*). In the present paper, we do not presume a till rheology, but just the two assumptions mentioned earlier.

Drumlins represent one variant of the subglacial bedforms induced (we suppose) by basal sliding. They are three-dimensional, with their long axes aligned with the direction of ice flow. Two other forms are noteworthy. The first is ribbed moraine, which takes the form of approximate two-dimensional parallel ridges aligned transverse to the direction of ice flow (Dunlop & Clark 2006; Lindén *et al*. 2008). The second is that of mega-scale glacial lineations (MSGLs) that are extremely long sinuous ridges aligned parallel to the ice flow, which are thought to have underlain former ice streams (Clark 1993). The different bedforms are plausible different members of a continuum (Sugden & John 1976), in which the aspect ratio of the bedform may be caused by the effect of a controlling parameter, the most likely candidate for which is the ice velocity. This leads us to a hypothetical picture of the way in which these bedforms develop sequentially as the ice velocity increases.

At very low velocities, a level ice–till interface may be stable. As the ice velocity increases, instability sets in first as transverse two-dimensional rolls, much as thermal convection in a fluid initially emerges as two-dimensional rolls. We surmise that the resultant finite-amplitude waveforms, which represent ribbed moraine, are themselves unstable to a transverse secondary instability, just as in thermal convection. We suppose then that the resulting waveforms become first drumlins, and then as the ice velocity increases further, the drumlins become elongated, and finally at ice stream speeds, they form MSGLs.

This sequence of events is entirely hypothetical, but plausible, and the papers of Hindmarsh (1998) and Fowler (2000) are evidence for the initial steps of the hypothesis. In his paper, Schoof (2007*a*) is less optimistic. He considers the two-dimensionality of the instability unsupportive of the instability hypothesis. Arguably, this may be misplaced pessimism, since while the linear instability may be suggestive of finite-amplitude bedform shape, it does not necessarily constrain it precisely. Schoof’s other criticisms are based on the observation in his numerical solutions that cavities typically begin to form in the lee of bedforms when they have a small amplitude that is comparable with the thickness of the deforming till layer. To extend his results, Schoof (2007*b*) considered the ice flow problem when cavities are present. He found a family of travelling wave solutions, whose amplitudes scale with the deforming layer thickness, which is quite small. However, the dimensionless amplitude in his fig. 4.1 is about 17, which indicates that these travelling waves can indeed obtain reasonable amplitudes. The stability, and hence realizability, of these travelling wave solutions remains unresolved. Schoof’s final criticism that drumlins can consist of stratified sands and gravels is one that we will not address specifically in this paper.

Our principal purpose here is to address the issue of including cavitation in the model, but with the aim of finding fully time-dependent solutions.

## 2. Mathematical model

The mathematical model we propose describes the two-dimensional flow of ice, considered as a constant viscosity fluid, over a layer of wet, deformable till. The horizontal and vertical axes are (*x*,*z*), and the objective is to calculate the elevation of the ice–till surface *z*=*s*(*x*,*t*). The geometrical framework is indicated in figure 3. The model was first proposed by Hindmarsh (1998) and analysed by him and Fowler (2000). Schoof (2007*a*) also presented the model and, in particular, rendered it dimensionless, and we follow his presentation, with some minor modifications that are described in appendix A.

A major simplification to this model is obtained by assuming that the aspect ratio ν of the bed is small. If we put ν=0, the basal boundary conditions (A 40) can be applied on *z*=0 and are
2.1
all applied at *z*=0. In these (dimensionless) equations, τ_{nn} is the deviatoric normal stress, Ψ is a stream function for the perturbation of the ice velocity from the basic shear flow, τ is the basal shear stress, is the (leading-order) sliding velocity, *N* is the basal effective pressure, *Π* is the reduced pressure, *q* is the sediment flux, *a* is the deforming till thickness and *V* is the mean till velocity. See appendix A for further details. It should be noted that the sliding velocity is a function of time only and is to be determined as described below.

We can use the Fourier transform of a function
2.2
to solve the ice flow equations. The result is that
2.3
where and are functions of *k*. Applying the approximate boundary conditions (2.1), we find
2.4
where *v*=Ψ_{z} and ℱ also denotes the Fourier transform.

We can invert these results using the fact that the Fourier transform of the Hilbert transform
2.5
satisfies
2.6
We then find that *s* is found by solving the system
2.7

The determination of follows from the requirement that the mean of *v* is zero. Thus, is determined by the requirement that
2.8

The model in equation (2.7) is equivalent to that presented by Schoof (2007*a*), with some minor modifications. In particular, we identify a drumlin elevation scale *d*_{D} that is distinct from the till thickness scale *d*_{T}, and their ratio *d*_{T}/*d*_{D} defines the (small) parameter α. Schoof (2007*a*) defines the sediment flux
2.9
as a function of τ and *N*; equivalently, we would have *a* and *V* as functions of τ and *N*. Simple assumptions are to take, following equation (A 36),
2.10
where *c*,*c*′≤1. There is some redundancy here, since in practice, it is the product *c**c*′ that is important. Note that since *A*_{τ}>0 and *A*_{N}<0, as well as *U*_{τ}>0 and *U*_{N}<0, then also *Q*_{τ}>0 and *Q*_{N}<0, as assumed in equation (A 22).

If we adopt equation (2.10), then the entrainment rate is defined by equation (2.7)_{4}, and the model is equivalent to that of Schoof (2007*a*). Alternatively, in detachment-limiting conditions, we might suppose that entrainment of sediment by erosion is rate limiting, in which case, a suitable generalization of equation (2.10) might be
2.11
where the transport-limited result in equation (2.10) is regained if ϵ→0.

It should be emphasized that the derivation of the expressions in equation (2.7) relies on the smallness of ν and *not* on any linearization procedure; thus the model in equation (2.7) is a fully nonlinear model for the evolution of drumlins.

### (a) Instability

Linear stability of the uniform state has been studied variously by Hindmarsh (1998), Fowler (2000) and Schoof (2007*a*). The summary below follows Schoof’s discussion. We take *Q*=*Q*(τ,*N*)=*A*(τ,*N*)*V* (τ,*N*). The basic uniform state is *s*=0, *N*=1, τ=θ, and *q*=*q*_{0}=*c**c*′(θ/μ−1), as determined by equation (2.10), since we can take
2.12
by choice of the velocity scale *u*_{0}. For greater flexibility, we retain the symbol , which will allow us to pinpoint the destabilizing term in the model. Linearizing equation (2.7) about this state and taking the Fourier transform as before, we find, after some algebra, that solutions are proportional to e^{σt}, where
2.13
the wave speed is
2.14
in which
2.15
and the growth rate is
2.16
It follows from this that the uniform flow is unstable if *R*>0 (presuming ).

There is a simple characterization for this instability (Schoof 2007*a*). Since , then is a function of *N* only and *R*≈*Q*′(*N*). Thus instability occurs if *Q*′>0, i.e. the sediment flux increases with effective pressure (although *Q*_{N}<0).

Because α is small, the maximum value of *r* occurs when *k* is large, and the relevant scaling for the linear stability analysis when α≪1 is obtained by rescaling time and space in equation (A 38) to (A 40) via
2.17
The dimensional time scale for growth of the fastest growing wavelength is
2.18
and the dimensional wavelength of this fastest growth rate is
2.19
To estimate these directly, we use values η=6 bar yr, *u*_{0}=100 m yr^{−1}, *d*_{T}=5 m, *N*_{c}=0.4×10^{5} Pa, *R*=1 and ; then,
2.20
The suggestion implicit in our scaling is that, over a longer time scale, the drumlins will coarsen and evolve to a longer length scale (formally; in fact , where *l* is the drumlin length scale defined in equation (A 37)).

## 3. An extended model

We now consider the problem of solving the nonlinear set of equations (2.7). This problem was studied by Schoof (2007*a*). So long as *Q*′(*N*)>0, the instability of the uniform state will continue to grow. There are two physical constituents of the model, which may act to prevent this. The first (which we call capping) is that, for sufficiently large *N* (), the till becomes immobile, and *Q* must decrease to zero. Thus, *Q*′<0 for sufficiently large *N*, and this may stabilize the growth of *s*.

The second is that when *N* reaches zero, cavitation will occur. In effect, this can prevent *s* decreasing indefinitely, because equation (2.7)_{3} still applies and (roughly) *N*∼1+*s*. While the effect of capping is easily effected by making *Q* decrease at large *N*, the effect of cavitation is less easily dealt with. This is discussed in the following section.

### (a) Cavitation

When the model described earlier is solved numerically, it is found that the interface evolves unstably, but that, fairly early in the computation, the interfacial effective pressure reaches zero (Schoof 2007*a*). Cavities form when this happens, and the ice separates from the underlying sediments. When cavities are present, we identify *s* with the cavity roof, and then its determination is made differently to the part where the ice is attached, where *N*>0. To be precise, the boundary conditions at *z*=0 are
3.1
but when a cavity is formed, these are replaced by the conditions
3.2

In terms of the model (2.7), the single change is to replace the Exner equation by the extra condition *N*=0; assuming that *f*(*u*,0)≡0, the shear stress is automatically equal to zero in this case.

The question arises as to how this can be solved numerically. The mixed boundary conditions provide a challenge for any method, because the cavity endpoints are not known, and there are potential stress singularities at these points. This problem has been studied in depth by Schoof (2007*b*). When the ice separates from the bed, the prescription of the bed *z*=*b* must be determined separately. In Schoof’s formulation, he supposes that no sediment transport occurs in the cavity (since there is no shear stress), and thus he poses the Exner equation,
3.3
This equation is supplemented by requiring that *b* is continuous at the downstream ends of the cavities, but *b* is not necessarily continuous at the upstream ends. In general, we can expect the sediment flux *q* to be positive upstream of a cavity, and the resulting jump in *q* causes a ‘shock’ to form in the till bed, and this propagates forwards at a finite rate. The shock will correspond physically to a slip face. Schoof then seeks travelling wave solutions of the problem using complex variable methods and finds a one-parameter family of such solutions in terms of their period.

Schoof’s generalization of the model equation (2.7) (in which we take *q*=*Q*(*N*)), can be succinctly written in the form
3.4
together with the alternative contact conditions
3.5
where *C* denotes the cavitated bed and *C*′ denotes the uncavitated bed; this assumes *Q*(0)=0, as is assumed here and also by Schoof (2007*b*). The ice base *s* is assumed to be continuous, but neither *b* nor *N* is continuous at both cavity endpoints.

Here, we propose a different way of formulating the problem with cavitation. Our strategy is to find a kind of weak formulation, such that the differing boundary conditions can be replaced by a single set of conditions. The resulting (still discontinuous) formulation will then be solved by approximating the discontinuous conditions by a suitable continuous set of boundary conditions.

Coalescence of the zero shear stress condition τ=0 with the sliding law *u*=*U*(τ,*N*) is easily done by suitable selection of the sliding law. For example, the generalized Weertman law (with coefficients such that *f*(1,1)=θ)
3.6
(suitably modified when τ<μ*N*) allows τ=0 automatically when *N*=0.

We now come to the gist of the matter. We allow for cavitation by extending the definition of *A* in equation (2.10) to the case *N*=0 by specifying
3.7
A typical form of *A*(*N*) is shown in figure 4. By extending the definition in this way, we retain the Exner equation and enable the determination of the cavity roof by allowing the ‘flux’ *q* to vary in such a way that *N*=0.

This way of extending the model can be written succinctly in the form (taking *c*′=1)
3.8
together with the alternative contact conditions
3.9
This can be compared directly with equations (3.4) and (3.5).

The fundamental mathematical distinction between Schoof’s formulation and that here is that, in Schoof’s formulation, discontinuities in *N* cause discontinuities in the sediment flux, and this leads to a discontinuity in *b*. In the present formulation, we will find that the sediment depth (and thus also the flux) is continuous, even though *N* is discontinuous: there is no tendency to form shocks.

Schoof’s formulation also has a simple and appealing physical basis, which lies in the assumption that sediment in cavities does not deform. Indeed, the situation appears comparable to the lee faces of fluvial and aeolian dunes, where the slip faces are manifestations of sediment shocks. However, a distinction in the subglacial case is that there is a lateral effective pressure gradient in the till between uncavitated parts of the bed and cavitated parts. This suggests the possibility that till may be squeezed into cavities from the sides, and indeed this is evidenced by the existence of crag and tail features. If such infilling can occur, then it suggests a formulation as given here, where the effective sediment flux (which arises through this infilling) is just such as to keep the till interface at the cavity roof. However, it should be pointed out that the present model formulation does not include cavity infilling in an explicit way.

### (b) Numerical solutions

The presence of the Hilbert transforms in equation (3.8) suggests the use of a spectral method, but this method also has difficulties. We used a spectral method to solve the equations, in the form
3.10
in which we take ϵ≪1, and where *A*(*N*) is given by the graph
3.11
assuming the generalized sliding law (3.6). We take and *c*=1 for convenience. Then the choice *c*′=μ/(θ−μ) allows the steady value *A*=1 when *N*=1. This is a unimodal function if *b*<1, for which *A*′(1)>0 if *b*θ>μ, which thus gives the instability criterion for this particular choice of sliding law and sediment flux.

However, we cannot use equation (3.11) directly, and we must approximate *A*. In our first simulations, we thus chose
3.12
with small values of δ corresponding to the approach to equation (3.11), see figure 5.

The solutions behave in the following way when started with initial data close to *N*=1 and when *A*′(1)>0. There is an initial phase where the solution for *s* grows in amplitude, propagating forwards as it grows. During this initial phase, *N* also grows in amplitude, oscillating spatially about *N*=1 on the unstable branch of the *A*(*N*) curve (figure 5). When the maximum of *N* reaches the maximum of *A*, there is a rapid transition to a new regime in which *N* resembles a square wave, oscillating back and forth between values, such that *A* is constant. In this regime, the solution reaches a steady, finite-amplitude state, consisting of stationary drumlins.

Figures 6 and 7 show the final steady states that were achieved in one particular run. We interpret the result in figure 7 as follows. The wiggles in the figure are a consequence of the Gibbs phenomenon, in which a finite Fourier truncation aims to approximate a piecewise continuous function. We thus infer that an accurate solution would portray the steady *N* profile as piecewise constant. Second, the minimum value of *N* represents the attempt of the effective pressure to reach zero; it does not get there because of the approximation in equation (3.12). We thus interpret the lower value of *N* as representing cavitation.

In order to test these interpretations, we modify the output for *N* to suppress the spurious Gibbs wiggles by taking a moving average of the form
3.13
equivalent to filtering the Fourier transform by multiplication by the windowing filter (i.e. we put ). In addition, we replace the approximation to *A* in equation (3.12) by
3.14
where
3.15
and the limit we seek is obtained when δ→0. This has the effect of putting the lower limit of *N* near zero, evidently attractive for cosmetic reasons (figure 8).

Figures 9 and 10 show the results of a computation with these alterations in place. We see that cavities form ubiquitously downstream of obstacles and that the drumlins themselves remain fairly regular in shape. We have tried other methods of smoothing or improving the results; we mention two. In one method, we filtered with *r*=1/(1+λ*k*^{2}), with the idea of suppressing high-frequency oscillations in *N*. This improves the numerics, which can now run with a time step of 10^{−5}, and to some extent smoothes *N*, but at the price of introducing spurious wiggles in *s*. The other method changes the equation for *a* to
3.16
where λ is small. In this case, we find that, as we would expect, the solution for *N* is indeed smoothed, but that now the waves travel at non-zero speed, to the right.

## 4. Discussion

We have discussed the onset of cavitation in the instability theory of drumlin formation and have developed a numerical technique to follow the solutions for the evolving bedforms past this onset. In our formulation of the model, it is then necessary to specify what happens at the bed of the cavity. Implicitly, we assume that the till beneath has no strength at *N*=0, so that it can freely flow into the cavity. In effect, there is no cavity, or, at least, the cavity is infilled with weak till. Since this is essentially what a crag and tail is, it seems a reasonable assumption to make, at least in some cases. We must contrast our model to that of Schoof (2007*b*). In his model, he assumes that the sediment is immobile in cavities, so that the cavities are water filled. It is difficult to compare results of the two versions of the theory, because Schoof restricts his attention to travelling wave solutions, in which consequently the sediment flux is discontinuous, whereas we do not find a tendency to form shocks (in *a*).

A further discrepancy between the models lies in the fact that Schoof does not consider till locking at high *N* and thus he effectively takes *A*=*N* (his eqn (2.6)). It can be seen, in our solutions, that the decreasing part of the graph of *A*(*N*) is essential in order that we can have *A*= constant in a stable steady state (the stability being associated with values of *A* in which the slope of the graph of *A*(*N*) is decreasing).

Non-trivial steady states are not possible in the Schoof formulation (3.4) and (3.5) because, in that case, the flux *q* must be continuous (otherwise a shock would propagate), and as *q*= constant in a steady state, we would infer that *Q*=1 is constant, and thus also *N*, if *Q* is monotonic. There is a possibility for non-trivial steady states either for the Schoof formulation with till locking or with the present formulation without till locking, but in both cases, we might expect such possible solutions, even if they exist, to be unstable, because of the positive-sloping portion of the *Q*(*N*) curve at one of the two values of *N* for which *Q*=1. Certainly, we have seen no numerical evidence for such steady states in our computations.

An interesting consequence of till locking is that, as the effective pressure increases at the summit, the shear stress at the bed will become concentrated there, and thus (since the far-field stress owing to the ice depth and slope does not change) the stress lower down the flanks will decrease. If the decrease is to a value below the yield stress, then the lower parts of the drumlin may become immobile, and in this way, one may obtain a situation in which the drumlin becomes capped by a layer of immobile till, with only the low-lying area between drumlins having low stress and also low effective pressure. The sediment in these areas is then marshy. Stream flow will be concentrated at the base and can erode the sediments and, as it does so, we might imagine that the capped drumlins can gently subside, as the sediment at their base is squeezed out. There is no evidence to support this idea, but it is consistent with what is known of till deformation and erosion.

In an extreme version of this scenario, the sediment removal is entirely by stream flow, and the till is virtually undeformable. This has been suggested by C. Schoof (2008, personal communication) as a possible mechanism for the formation of stratified drumlins in Washington State. The drumlins can retain previous fluvially stratified structure as they erode because they become stationary (the till is at high effective pressure), and the drumlins emerge as the landscape is sculpted by the water flow.

In an even more extreme version, Shaw *et al.* (Shaw 1983, 1994; Shaw & Kvill 1984; Shaw & Sharpe 1987; Fisher & Shaw 1992) have suggested that drumlins are sculpted in massive subglacial floods. An attraction of this radical theory lies in the more recent discoveries of subglacial lakes in Antarctica and of the fact that they can flow from one to another in subglacial floods. It has even been suggested that massive floods of the type Shaw seeks might have occurred below the former Laurentide ice sheet (Evatt *et al*. 2006). However, it is by no means simple to generate such large floods. The two basic problems with the Shaw mechanism are that there is no predictive way to get floods of the right magnitude in the right place, and, worse, there is no theory that predicts what putative bedforms might actually be produced. In actual fact, it may well be that basal water flow plays a key role in certain aspects of the present theory, but it seems unlikely to involve water in the biblical quantities imagined by Shaw. For further comments on Shaw’s theory, see Eyles (2006), who also reviews more generally the effects of water in the subglacial environment.

## 5. Conclusions

In this paper, we have extended the solution of Hindmarsh’s instability model of drumlin formation to allow for the presence of cavities. We do this by firstly representing till flux as *q*=*a**V* , where *a* is the deforming till depth and *V* is the mean till velocity, which we simply assume is some constant fraction of the basal ice velocity. We suppose that the depth of the deforming till is determined by the position where the yield stress is attained, and we suppose that the sliding velocity (and also till flux) satisfies a generalized Weertman law of the form given in equation (3.6).

The consequence of these fairly general assumptions is that till flux will be a hump-like function of effective pressure *N*, and the uniform state in which the bed is flat is then unstable on the part of this relation where *q* increases with *N*. Our numerical investigations suggest that the resultant instabilities cause drumlins to grow to a finite amplitude and that they become stationary. Estimates of height and length of drumlins are roughly consistent with observations, but a detailed comparison will await 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. Particular thanks to Chris Clark for the prize of a Barbour jacket. Thanks to the referees, Christian Schoof and Anders Schomacker, whose comments have resulted in a much improved text.

## Appendix A

#### (a) Effective pressure

In this appendix, we describe more fully the instability model for drumlin formation. The model was first described by Hindmarsh (1998). The version presented here stems from the work of Schoof (2002), and more specifically follows Schoof (2007*a*), with some minor modifications.

If the pore water pressure at the interface *z*=*s* is and the overburden normal stress there is *P*^{s}, then, assuming hydrostatic and lithostatic balance,
A 1
where ρ_{w} and ρ_{s} are the densities of water and sediment, respectively, and ϕ is the sediment porosity. Within the till, the effective pressure *p*_{e} is defined as
A 2
and thus
A 3
where
A 4
and we define *N* to be the effective normal stress at the ice–till interface, i.e.
A 5

The interfacial normal stress *P*^{s} is related to the stress in the ice by
A 6
where σ_{nn} is the normal stress in the ice, τ_{nn} is the deviatoric normal stress in the ice and is the ice pressure at the bed. We define a reduced pressure *Π* in the ice by
A 7
where *p*_{a} is the atmospheric pressure, and we define the effective pressure in the drainage system as
A 8
where *p*_{c} is the water pressure in the local drainage system, which we presume known.

It then follows that A 9 where A 10

There are two important consequences of equations (A 3) and (A.9). The effective pressure at the ice–till interface increases with its elevation. This causes decreased till deformation at higher elevation. Additionally, the effective pressure also increases with till depth; in effect, only a finite layer of till will deform.

#### (b) Ice flow

The equations describing slow two-dimensional flow of ice, considered to be an incompressible fluid of viscosity η, are
A 11
Here, *Π* is the reduced pressure defined in equation (A 7) and *z*_{i}′=∂*z*_{i}/∂*x*. The velocity components in the (*x*,*z*) directions are (*u*,*w*). We take the ice surface slope *z*_{i}′ to be constant because the relevant horizontal length scale for drumlin formation is much less than that over which ice sheets vary.

#### (c) Boundary conditions

Appropriate boundary conditions in the far field are those of Schoof (2007*a*),
A 12
τ_{b} is the ‘basal shear stress’, defined by
A 13
The validity of the conditions (A.12) depends on the value of the dimensionless parameter
A 14
in which *l* is the drumlin horizontal length scale; equation (A 12) is appropriate if σ≪1.

At the bed *z*=*s*, one of the boundary conditions is equation (A 9); the normal deviatoric stress is defined by
A 15
The shear stress at the ice–till interface (note that this is not the same as τ_{b}) is
A 16
The sliding law at the bed is then taken to be
A 17
where *U*(τ,*N*) is the sliding velocity, which is discussed further below. We will also use the sliding law in the more common form
A 18

The final condition at the base is the kinematic condition
A 19
wherein we ignore basal melting. We now consider the determination of *s*.

#### (d) Sliding and deformation of till

The evolution of the bed *s* is determined by the Exner equation
A 20
in which *q* is the till flux, taken to depend, like the sliding velocity *U*, on bed stress and effective pressure,
A 21
Since τ and *N* are already given by equations (A 16) and (A 9), this completes the specification of the ice flow problem, once *U* and *Q* are given, we assume that the partial derivatives satisfy
A 22
Equivalently to the assumptions on *U*, we suppose
A 23

We suppose that till is ‘plastic’ in the sense that it has a yield stress τ_{c}; for shear stresses below this, no deformation will occur. We take
A 24
where μ is an *O*(1) coefficient of friction. We define a depth scale
A 25
For an elevation change of *d*_{D}, the effective pressure at the ice–till interface will change by *N*_{c}, which we suppose is comparable to τ_{c}/μ for flow over wet, deformable sediment. Thus, this depth scale gives a drumlin height sufficient to cause the till at the summit to cease deforming. We have in mind that the ‘locking’ of the till at the summit will provide a self-limiting mechanism to limit growth of drumlins and thus that *d*_{D} may be a representative drumlin height. Choosing *N*_{c}=0.4×10^{5} Pa gives a value *d*_{D}≈50 m.

We define a second depth scale
A 26
then, we can see that the effective pressure increases by *N*_{c}∼τ_{c}/μ in a depth *d*_{T} below the ice–till interface. This suggests that deforming till will be limited to a finite depth of *O*(*d*_{T}) below the ice–till interface, and in particular, it suggests that in order of magnitude, *Q*∼*d*_{T}*U*.

Actually, it is appropriate to include the depth of deforming till explicitly. This is because there can be no till flux if there is no deforming till. We propose a prescription based on similar considerations that arise when considering sediment transport in rivers (Howard 1994), where, in principle, one must consider conservation of deformable sediment as a separate equation (Tucker & Slingerland 1994; Fowler *et al.* 2007). We retain the Exner equation (A 20), but now we prescribe sediment flux as
A 27
where
A 28
is the mean till velocity and *a* is the deforming till thickness; in uniform conditions, we suppose that
A 29
and *Q*=*A**V* still satisfies equation (A 22). Conservation of deforming sediment requires that, in addition to equation (A 20), *a* should satisfy the conservation law
A 30
where represents the rate of entrainment of sediment from the undeforming (B horizon) sediment. This equation may be useful in situations in which sediment supply is limited (detachment-limited kinetics, in fluvial parlance (Howard 1994)).

#### (e) Non-dimensionalization

We define a length scale *l*, yet to be chosen, which is to be a representative length scale for drumlins. We define a stream function ψ via
A 31
We take the basic shear flow without bed perturbations to be
A 32
and when the aspect ratio *d*_{D}/*l* of the emerging bedforms is small (as is the case in practice), then the perturbation to the stream function is correspondingly small. This corresponds to the situation that occurs in sliding theory (Fowler 1981).

We now scale the model by choosing A 33 the dimensionless velocity is introduced here because the development of the topography causes a change in the mean sliding velocity. We scale the depth of the till by writing A 34 Thus, the dimensionless effective pressure in the till is A 35 and the yield criterion (A 24) implies that deformation will occur if the dimensionless till depth ζ satisfies A 36

The value of *u*_{0} is determined by the magnitude of the sliding velocity, and the horizontal length scale is defined by balancing the stress and strain rates, thus
A 37

If we choose *u*_{0}=100 m yr^{−1}, η=2×10^{13} Pa s, *N*_{c}=0.4×10^{5} Pa, *d*_{D}=50 m, for example, then *l*=274 m. Other typical values, with ρ_{s}=2.5×10^{3} kg m^{−3} and ϕ=0.4, are *d*_{T}=4.6 m, and the time scale is 30 yr.

With this choice of scaling, the dimensionless model for the ice flow is
A 38
with far-field matching condition
A 39
The basal conditions take the form
A 40
and these are all applied at *z*=ν*s*.

The dimensionless parameters σ, θ, ν and α are defined by
A 41
With *l*=300 m, *z*_{i}=1500 m, thus τ_{b}=0.15×10^{5} Pa with an assumed ice surface slope of 10^{−3}, *d*_{D}=50 m, *N*_{c}=0.4 bar and *d*_{T}=5 m, typical values are
A 42

The far-field matching condition (A 39) is appropriate if σ≪1, and so we adopt the formal limit σ→0. This simply removes the gravitational term from the momentum equation in (A 38).

## Footnotes

- Received November 26, 2008.
- Accepted May 21, 2009.

- © 2009 The Royal Society