Subglacial swamps.

The existence of both water and sediment at the bed of ice streams is well documented, but there is a lack of fundamental understanding about the mechanisms of ice, water and sediment interaction. We pose a model to describe subglacial water flow below ice sheets, in the presence of a deformable sediment layer. Water flows in a rough-bedded film; the ice is supported by larger clasts, but there is a millimetric water layer submerging the smaller particles. Partial differential equations describing the water film are derived from a description of the dynamics of ice, water and mobile sediment. We assume that sediment transport is possible, either as fluvial bedload, but more significantly by ice-driven shearing and by internal squeezing. This provides an instability mechanism for rivulet formation; in the model, downstream sediment transport is compensated by lateral squeezing of till towards the incipient streams. We show that the model predicts the formation of shallow, swamp-like streams, with a typical depth of the order of centimetres. The swamps are stable features, typically with a width of the order of tens to hundreds of metres.

Introduction

Glaciers and ice sheets commonly have basal temperatures that reach melting point, on account of the geothermal heat flux from the Earth, and the insulating effect of the ice cover; consequently, most glaciers and ice sheets have a basal water system that is formed from the resulting meltwater. In glaciers, this can be enhanced by percolation of surface meltwater to the bed, whereas in the interior of ice sheets the basal meltwater is formed entirely by subglacial melting.

It is widely acknowledged that the resulting subglacial hydrologic system exerts a fundamental control on ice sheet dynamics; the volume of evidence for an active subglacial water system in Antarctica and under ice masses is constantly growing (e.g. [1,2]). The large variation in ice velocities across Antarctica, and particularly the presence of ice streams, is commonly considered to be a consequence of this subglacial hydrologic system (e.g. [3,4]). Furthermore, it is important to know about water movement at the bed because of its critical role in phenomena such as potentially hazardous subglacial flooding (e.g. [5]), as well as in glacial erosion and landform evolution (e.g. [6,7]). An improved understanding of the drainage system will allow us to learn more about critical feedbacks between the water, ice and bed.

Subglacial drainage systems commonly discussed in the literature, such as linked cavities [8,9], thin patchy films [10] and Röthlisberger channels [11] are framed in the context of the classic hard glacial bed: cavities/thin patchy films form behind bedrock protrusions, while the channels are mostly thought of as rock-floored and incised into the ice, and the issue of sediment transport through them is rarely discussed (cf. Fowler [12, p. 699 ff].). However, evidence suggests that the observed rapid basal velocities in parts of Antarctica are facilitated by the presence of a layer of till at the base of the sheet [13-16], and therefore these hard-bed drainage systems may not be relevant to subglacial water flow beneath ice streams. A soft bed has rather become the paradigm assumption.

For probable permeabilities and till thickness, Darcian porous flow through the sediment is likely to be insufficient to evacuate all the meltwater present at the bed [17,18], and as such the till is water saturated. As a granular material, the till is expected to deform with Coulomb-plastic behaviour [16,19-22], and its behaviour depends on the effective pressure, which is determined by the drainage system. In an effort to understand the possible nature of subglacial water flow over a layer of saturated till, Walder & Fowler [18] developed the concept of canals, which are in essence inverted Röthlisberger channels, incised in the underlying sediment. In their work, they showed that these canals would form a distributed system owing to the fact that effective pressure decreased with increasing water flux. However, their model can be criticized on a number of counts. For example, they assumed a pseudo-viscous flow law for till, an assumption that has led to debate [23-25], and on which an effective conclusion has yet to be reached. More pertinently, Walder and Fowler's model was conceptually similar to that of Röthlisberger, insofar as it assumed a local channel surrounded by till and ice, and the resulting model was consequently essentially algebraic in nature. It also made assumptions concerning the shape of channels. One was that the canal depth was controlled by the necessity to enable bedload transport to remove sediment. The other was that whether a canal or channel formed depended on the relative softness of ice and till, and this depended on the effective pressure.

The aim of this paper is to develop a more sophisticated model of subglacial water flow interacting with the sediment in the presence of a soft bed. In particular, we abandon the restrictive view of channels or canals as pre-formed, and we aim to describe their evolution from an initial spatially extensive rough water film. As an eventual aim, we hope to characterize sub-ice-stream water flow sufficiently to allow for a useful parametrization of its effect on the ice flow. This requires a suitable averaging over the presumably small-scale dimensions of hydraulic stream flow to describe the larger scale ice stream flow. This would be a significant improvement on recent work looking at ice-stream formation over a Coulomb-plastic till, which simply uses a basic diffusion equation to describe the evolution of a water layer over saturated till, assuming the water diffuses from high-to-low effective thickness [26-28]. Furthermore, other water flow models that currently exist on the larger scale generally consist of a simple distributed system for the water, where the water flows down gradients in hydraulic potential (e.g. [29,30]). While such models are able to demonstrate where we expect the water to flow, and the sensitivities of the flow paths to small changes in the system, they do not begin to include any of the intricacies associated with the subglacial hydrologic system. There have been various efforts made at forming more detailed hydrology models in the literature, some considering simple flow-band models [31,32], and others with two-dimensional models including the details of interaction between a distributed drainage system and a channellized drainage system (e.g. [33-35]). For these latter models, the channellized drainage is on the grid scale. In this paper, however, we look at modelling on the small scale, with a view to upscale later to model the water flow on the ice-sheet scale.

In §2, we present the coupled model for water and sediment transport that forms the basis of our study. We non-dimensionalize the system in §3, before simplifying and solving it in §3a. A discussion of the results and their context follows in §4, and the conclusions of the study are in §5.

Set-up and governing equations

We conceive of water at the base of an ice stream in the following way. Water is formed by the melting of ice at grain boundaries, whence it seeps down to the bed in a steady drizzle of some 5 mm yr−1. The ice is underlain by sediments, or till, which thus become saturated with water. Typically, the sediments are relatively impermeable, and this has two effects: the water cannot drain through the sediments and must consequently form some kind of hydraulic system between the ice and the bed. The second consequence is that this causes the water to be at high pressure, so that the till is deformable. The simplest kind of water system is a uniformly distributed thin film at the interface between the ice and the till. With a clean separation, as imagined by Weertman [36], the resulting flow is always unstable [37], but we can suppose such a film could exist if it is thinner than supporting clast sizes; we call this a Creyts–Schoof film [38]. The idea is that a thin water film whose thickness is measured in millimetres will generally not drown the coarser particles of the till, so that despite the partially lubricating effect of the water, the ice still grips the coarse till and can transmit a stress to it. Evidently, the degree to which this happens depends on the thickness of the water film. The film is subject to Walder's [37] instability mechanism, and we can expect some kind of streams to form. Our model will show this.

We consider a cartesian set of axes (x,y,z) in which x points in the ice flow direction, y is across flow and z is vertically upwards. As shown in figure 1, the till surface is denoted by z=b, the lower ice surface is denoted by z=s and the upper ice surface is denoted by z=zi, all of these being functions of x, y and time t. The depth of the water film is thus

Figure 1.

Geometry of the flow under consideration. A three-dimensional schematic and a cross section with notation labelled.

The hydraulic head driving water flow at the ice/till interface is where pw is water pressure, pa is atmospheric pressure (assumed constant), ρw and ρi are the densities of water and ice, g is the acceleration due to gravity, and di is the local ice depth, a constant and distinct from zi which varies in the downstream ice flow direction; the constant term in parenthesis is added for convenience. The ice pressure pi is approximately cryostatic, and it is common in analysing ice flow dynamics (as also in numerous other fluid flows) to define a reduced pressure Π which is the ice pressure minus its cryostatic component. At the base, this takes the form [6], eqn (2.1), and thus (cf. [6], eqn (2.8)) where τnn is the deviatoric normal stress of the ice, and is the effective normal stress at the ice/sediment interface. Equation (2.4) reflects the well-known fact that basal water flow is driven primarily by the ice surface slope, and the basal slope only contributes 9% to the flow direction (because Δρwi≈0.09ρi). In the following subsections, we describe briefly the constituent parts of the model involving ice, water and sediment motion, together with the important concept of a closure relation.

Ice flow

Fowler [39,6] derived expressions for the solution of the problem of ice flow over locally varying small amplitude topography. If x denotes the direction of flow, then we split the ice surface elevation into the sum of a local average value and the perturbation from this; specifically, we write where is the (small) average ice upper surface slope, Z∼O(1) is the dimensionless local average depth and H is a dimensionless measure of the local ice surface perturbation due to flow over any short wavelength basal topography. Typically, we suppose Z varies on length scales of hundreds of kilometres, while H varies on scales of hundreds of metres. It follows that the hydraulic potential is where Φ is defined as and the (average) basal stress is (note that Φ=−Θ in the notation of Fowler [6]). The hydraulic potential is defined up to the addition of an arbitrary function of time, and we will use this later.

We are able to solve the ice flow problem over relatively small length scales, typical of those which will concern us here. The method is as follows: the ice flow equations are those of slow flow, and if we assume a constant viscosity, they reduce to the equations of Stokes flow, which are linear. In addition, the ice upper and lower surfaces are small perturbations of level surfaces, so that the boundary conditions on them can be linearized about the unperturbed flat surfaces. The consequence of this is that we obtain a linear problem with inhomogeneous boundary conditions on upper and lower surfaces, which can therefore be explicitly solved by taking a Fourier transform in the x-direction, and solving the resulting linear differential equations analytically in the vertical.

The boundary conditions are those of no stress at the upper surface, while at the base the sliding law provides a relation between basal downstream velocity, shear stress and effective stress N. If we suppose that the normal stress at the base, and more specifically Φ, in (2.10), is given then the Fourier transforms of all other quantities at the upper and lower surfaces can be written ultimately in terms of the Fourier transforms of Φ and N.

Specifically, the upper and lower surface kinematic conditions take the form where in both conditions, we can ignore small source terms due to surface accumulation and basal melting, respectively. Here, u is the ice velocity in the x-direction (essentially constant), and v,w are the vertical velocities at surface and base, respectively.

It turns out that the ice surface perturbation H relaxes rapidly to a quasi-equilibrium, so that H can be determined, and w depends only on Φ and N and primarily the former. Specifically, we have, approximately, where where is the Fourier transform, and k=|k|. Note that the quantities, particularly x and k, are dimensional in these definitions (which have been translated from the dimensionless equivalents in [6]), but M is dimensionless and is shown in fig. 3 of [6]. Particularly, for wavelengths less than about three times the ice depth (, thus ), M≈2 is constant. Note that (2.13) ignores a similar integral term due to sliding and thus N, which is in fact negligible for . If we presume the short wavelength limit in which M≈2, then a direct calculation shows that A rather better approximation which takes account of the fact that is finite (and specifically ) as k→0 is to choose since then This means that the kernel J(x) in (2.13) decays exponentially rather than algebraically at large wavenumber, so that w is insensitive to short wavelength variations of Φ (i.e. they are damped).

Water flow

The base of the ice rests on a cobbled bed of subglacial till, typically having a range of particle sizes varying from coarse clasts to fine sediments. Between the coarser cobbles, meltwater from the ice, produced by the combined effects of viscous heating in the basal ice and geothermal heat flux from below, flows in a thin film which can drown the finer particles, but not the coarsest clasts (unless it becomes quite deep); this is the Creyts–Schoof water film [38].

Mass conservation of the water flow in the film takes the form where, on a regional scale, G is the geothermal heat flux, τb is the basal ice stress, u is the basal ice velocity, qT is the sensible heat flux into the overlying ice and L is the latent heat. In our model, we take Γ to be constant because we are concerned with small length scales; Γ will of course depend on ice depth and velocity, with its variation occurring over larger length scales. It is important to realize that the frictional heat source τbu arises from an integration of the ice viscous dissipation term over the basal sliding region, and τb and u here refer to conditions near the base of the far-field ice sheet flow, but far from the actual interface itself [40]; they are the quantities used in the ice sheet sliding law.

We assume a local Poiseuille flow in the water film, which implies where ηw is the viscosity of water; thus h satisfies

Film closure

Equation (2.22) must be supplemented with a further closure equation, which is analogous to the closure equation in Röthlisberger's [11] channel flow model. The reason for this is that the effective normal stress in (2.6) varies rapidly along the Creyts–Schoof film. Clearly, the actual normal stress is equal to the water pressure on those parts of the film which are unsupported, with the excess stress being provided by the supporting clasts. As a consequence, a Röthlisberger-type closure condition acts on the film due to the differential stress along the ice surface. A suitable prescription for this closure equation is of the form where l* is a measure of the distance between supporting clasts. It will depend on film thickness h because as the film thickens, the distance between the supporting clasts increases, so that l* is an increasing function of h. When the water film becomes deep enough for separation of the bed and ice to occur, the water flows in streams or cavities, and it seems that l* should then become equal to the stream width, or more generally be a non-local functional of h. Here we simply define l* as where l0 is a length scale that represents the typical clast spacing in the absence of water, and N*(h) (denoted Λ(h) in [41]) is some decreasing function of h, with N*(0)=1.

Sediment transport

Sediment transport is most usually discussed in the context of fluvial or aeolian transport, where it is due to the action of the water or wind stress on the sand grains [12]. The resulting evolution of the sediment surface is then described by the Exner equation, which expresses conservation of sediment by relating sediment surface elevation to the particle flux. The same situation can occur subglacially, but now the sediment can be transported both by shear in the water film and by shear due to the overlying ice; in addition, there is a squeezing effect caused by horizontal differences in the ice normal stress (and thus also the effective stress).

The evolution of the sediment surface z=b in our model is thus given by the Exner equation where ϕ is the bed porosity, and the volumetric sediment flux has two constituents, a bedload volume flux qb and a creep flux qc; E represents net erosion of the bed as suspended sediment and can also be written as a divergence (of minus the suspended sediment flux). We suppose here i is the unit vector in the ice flow direction x and A is the depth of the deformable till layer. The first term represents the shear-driven deformation by the ice, and the second represents till squeezing by effective pressure gradient. The depth A is determined by the horizon below which the effective pressure (which increases with depth) is sufficiently large that the yield stress is not exceeded. It is given by Fowler [6] as where ρs being the sediment density and μ is the coefficient of interparticle friction. Note that τ in (2.27) is the actual shear stress at the interface and not the same as the ‘far-field’ basal stress τb used in the sliding law. The two are related [42], if we suppose that the actual sliding law at the interface is by the fact that where the overline denotes a local spatial average. In particular, if we choose then, since u is not perturbed at leading order at the interface,

The bedload transport is defined in terms of an effective stress in the water layer, which combines the stress exerted by the water flow, with the effect of the bed slope. The effective stress is [6] where Ds is the grain size of the sediment. We then write where B is generally a function of τe. Typical forms (e.g. the Meyer-Peter/Müller Law) involve measurements on turbulent flows, together with absence of transport below some critical yield stress, thus B∝[τe−τc]1/2+, but such sub-aerial relations may be inappropriate, since for thin films and relatively slow flows, we can expect the flow to be laminar and not well described by sub-aerial formulae. Note, however, the specific requirement that B=0 when τw=0, i.e. h=0.

Normally, one assumes the bedload layer to be a very thin layer, one or two grains thick, at the base of a deeper water flow. More precisely, the bedload is a constituent of the fluid flow, so that the water flux is really considered to comprise both water and mobile sediment. In particular, we must have qb

The model for coupled ice, water and sediment evolution thus consists of the seven equations (2.1), (2.9), (2.12), (2.13), (2.22), (2.23) and (2.25) for h, s, b, ψ, Φ, w and N.

Boundary conditions

To pose suitable boundary conditions, we consider a rectangular domain in which the upstream end x=0 might be a divide, or a cold-temperate transition line; in either case, it is reasonable to put In this case, there can be no hydraulic sediment transport, and so we must choose B=0. If the till is not deformable, then A=0 and no further condition is necessary. If it is, it seems reasonable that the till flux is just , and then we have Alternatively, we might prescribe an inlet water mass flow together with (2.37). A suitable choice of water flux is the flux scale where l is a suitable macroscopic ice flow length scale.

At the sides, it is suitable to prescribe no transverse flux of water or fluvial sediment, and then where 2W is the width of the domain.

Finally, the front of the flow will typically be a land-based terminus or more likely a grounding line, but in either case we may put It is less obvious how to prescribe an outflow condition, but one which avoids the possibility of unnecessary boundary layers is to choose although other justifiable choices are possible.

Non-dimensionalization

We non-dimensionalize the equations by scaling the variables as where the subsidiary scales are defined by (2.39), and the last two being from Fowler [6].

The value of the bedload transport qb is estimated as follows. Values of qb/q0 for a turbulent Alpine proglacial stream are of order 10−4 [43]. A similar assumption under an ice sheet gives the estimate in table 1. As a typical bimodal material, there is no value for the grain size Ds. Tulaczyk et al. [44] give values of approx. 100 μm for the sandy till under Whillans Ice StreamQ1. Our value is comparable to this, but somewhat reduced so that the typical film water stress ∼ρighS is slightly greater than the critical Shields stress approx. 0.06ΔρswgDs necessary for transport to occur. Values of the other scales are given in table 2.

Table 1.

Assumed constants in the model.

symbol	meaning	typical value
B0	sediment transport coefficient	qb/τ0
di	ice depth	103 m
Ds	grain size of sediment	30 μm
g	gravity	9.8 m s−2
l	ice sheet length scale	500 km
l0	clast spacing	1.2 m
qb	bedload transport	0.3 m2 yr−1
Si	ice surface slope	∼10−3
u0	ice velocity	100 m yr−1
Γ	melt rate	3 mm yr−1
Δρsw	ρs−ρw	1.6×103 kg m−3
Δρwi	ρw−ρi	83 kg m−3
ηi	ice viscosity	1014 Pa s
ηw	water viscosity	1.8×10−3 Pa s
ηs	notional sediment viscosity	2.7×109 Pa s
μ	coefficient of friction	0.6
ρi	ice density	0.9×103 kg m−3
ρw	water density	103 kg m−3
ρs	sediment density	2.6×103 kg m−3
ϕ	sediment porosity	0.4

Table 2.

Variable definitions and associated scales if appropriate.

symbol	meaning	scale
A	till depth
b	sediment elevation
d	sediment topography	21.6 m
dT	till deformation depth	1.87 m
h,h0	water film thickness	3.9 mm
lD	drumlin length scale	625 m
l*,l0	distance between supporting clasts	1.2 m
N,N0	effective stress scale	1.76×104 Pa
q,q0	water flux	1.5×103 m2 yr−1
s	ice/water interface
t,t0	time	10.4 yr
zi	upper ice surface
Π	reduced pressure
τnn	deviatoric normal stress at bed
τb	basal shear stress	0.88×104 Pa
τw,τ0	water shear stress	3.4×10−2 Pa
ψ	hydraulic potential

This leads to the dimensionless set of equations where the parameters are defined by where lD is the drumlin length scale [42] The dimensionless deformable till thickness is, if we assume (2.32), If we take the kernel to have the form of (2.17), then its dimensionless expression is Typical values of the parameters are given in table 3, with values of variable scales given in table 2.

Table 3.

Approximate dimensionless parameter values.

symbol	typical value
a	2×103
r	0.9
α	6.2
β	0.96×10−5
γ	1.33
δ	1.8×10−4
ε	0.125
λ	1.1×10−7
Λ	0.64×106
ν	2×10−3
Π	2.0
σ	1.6×10−3
Ω	0.48×103

Simplification and reduction

The basic model in () can be reduced to consider the leading-order effects of interest for studying the water–sediment interaction in more detail. We define then, bearing in mind that Z′(x)=−S, we have providing six equations (along with (4.1a)) for h, s, N, Ψ, Φ and A.

We are interested in the possibility of a lateral, stream-forming instability in which downstream sediment transport is compensated by cross-stream sediment drift. A similar analysis was carried out by Fowler [6], using a model which is a slightly simpler version of (4.1), but his focus was on the interaction of water flow with ice deformation (and the formation of mega-scale glacial lineations (MSGLs)) rather than the formation of genuine finite depth streams. Indeed, his instability results were not pursued beyond the infinitesimal régime.

Our interest here is not with moulding of the subglacial topography, and so we will seek simplifications for the evolution equation (4.1e) for s. We might, for example, take the right-hand side to be prescribed.

There is a basic x-dependent uniform state in which whence also A=A(x), and the sediment surface b=s/δ−h is given by where N′=dN/dh. This causes a slow subsidence of the bed as the sediment is eroded. Since σ and β are small, this subsidence is essentially corresponding to an erosion rate of the order of a millimetre per year. Note that then st≈−δ(Au)x.

This basic state is unstable, analogously to the result of Fowler [6], and this can be due to the term −σ(BSh)x in (4.1c). In [6], the basic state was taken as uniform (h=const.) and the slope S was also constant, but allowance was made for the nonlinearity of B: instability occurs if B′>0. This is analogous to the Smith-Bretherton situation [45], but we can directly see that instability will still occur if B is constant, if also the slope increases down flow (S′>0), which it sometimes, but not always, does in practice [46,47]. Instability is also promoted by the till flux term Aui, if u increases with x as is likely, and A increases with h, i.e. A′(N)<0 (note this condition suppresses the ribbed moraine instability [42]).

Our specific concern is with the evolution of this instability as h grows to larger depths. When this occurs, it is clear from (4.1) that it is necessary to rescale y, and specifically, we write corresponding to a length scale of about 1500 m, somewhat larger than the MSGL width scale found by Fowler [6]. We substitute this into (4.1), and then neglect the small terms in ε. We use the fact that Ω≫1 to neglect the advective term in (4.1e). We also neglect the small bedload transport terms in σ. This is inadmissible if true separation of the water film occurs, since then A=N=0, so we initially restrict our attention to the case where this is not so. We then obtain the approximating system where and The values of the parameters α and γ are roughly O(1), as indicated in table 3. The reduction of the convolution integral in (4.1e) to the form in (4.6d) uses the fact that Φ varies relatively slowly in x compared with y, and then that the integral where K0 is the modified Bessel function of the second kind.

To be specific, we aim to solve the model (4.6) in the domain −Lwater flow through the sides: It is a consequence of this that although β/ν is small in (4.6b), integration across the catchment yields the conservation law for water flow, and we need to solve this in conjunction with (4.6c).

A reduced model

In this model, dependence on x is introduced only through the water flux (4.11), the advective till flux term in (4.6c) and the dependence of A on the slope in (4.8). The role of the x derivatives is to provide for change in mean water film depth (via (4.11)) and to provide an advection term for A (and thus N and thus h); neither of these has any essential effect on the mature downslope solution, and we therefore now focus on solutions which depend only on Y and t. It is possible to choose forms for S and u such that this is quantitatively justified (e.g. increasing downstream slope S=x and linearly increasing velocity u=u0+x), but we emphasize that such choices are only cosmetic. Actual surface slopes of ice streams are quite variable [46,47], but are generally more concave (thus decreasing slope) than the typical inland ice convex profile.

We can then take the approximate solution of (4.6b) to be Ψ≈0, in view of our earlier remark concerning the definition of ψ up to addition of an arbitrary function of time, so that the normal stress term and the model simplifies with ux=1 to

A simpler version of the evolution equation for s results from a formal assumption of large α; then the spiky nature of K0(α|Y |), together with the fact that suggests that we take K0(α|Y |)≈(π/α)δ(Y), and this leads to as a further simplification of (4.13b), and we shall use this as a basis for discussion. The implication of (4.15) or (4.13b) is that the basal ice surface and hence the bed relaxes to a state determined by the effective stress, itself determined by the water film depth; and where the film is thinnest (so N is high) the bed is most elevated.

The simplest discussion of the model (4.13), but using (4.15), assumes that s has relaxed to its equilibrium s=N(h), so that W=0 and (4.13) reduces to the system The precise choices for the functions N(h) and A(N) are immaterial for the finite amplitude solutions that we obtain, and for ease of exposition we will take A=|N′|=1, not the least because it is by no means clear that (4.16) even has a solution. The partial differential equation for h paradoxically describes sediment conservation, and the two terms on the right represent the downstream sediment flux dragged by an accelerating ice flow, and the squeezing of subglacial till due to lateral variations in the effective stress, respectively; the integral term represents conservation of the quantity of water in the film.

We thus consider the following reduced model: together with Because physically h≥0, the differential equation (4.17) should be interpreted as applying when h>0.

A different limitation occurs if N reaches zero, which will be the case if h reaches a critical value hc where proper flow separation occurs: the ice is afloat and the water film becomes a stream or cavity. In this case, A=W=0 also and the model breaks down. It is in this situation where the bedload terms in σ may become important. We initially ignore this possibility, formally supposing that N>0 everywhere.

It is obvious that the model (4.17) and (4.18) has a fundamental problem, and this is not caused by any of our simplifications, which are inessential. The solution of the differential equation with the boundary conditions is inconsistent with the integral constraint. The problem lies with the latent assumption that h>0. As discussed above, the model for h is more properly It is clear that the solution must immediately contain regions of finite support, and that there is a steady state consisting of a number of streams. The simplest solution consists of a single stream, for which the steady symmetric solution is where and we must have so that L>Ym. In the present situation, the meaning of the inequality is unclear, but its interpretation emerges in consideration of more general versions of (4.13a), where for example the choice W+A=h−1 shows that the inequality necessary for channel formation corresponding to (4.22) is precisely the condition for instability of the uniform steady state.

Numerical solutions

The only parameter in (4.19) is L, but it can be scaled out in the following way. Because the solutions develop finite support, the upper limit in the integral is irrelevant, and by rescaling the equation for h is unaffected, and the integral (assuming a symmetric solution) takes the form A canonical form is chosen by selecting and in the canonical form of the problem, the integral constraint is where Yf is the stream boundary. The steady-state value of Yf corresponding to (4.21) is

We now solve the canonical problem numerically. To do this, we follow a method devised by Fowler et al. [48] , where a very similar problem arises. First, we suppose the solution of (4.19) is symmetric, so that Next, we change variables from Y,t to X,T defined by The chain rule gives and thus also, applying these to the inverse function Y (X,T), Substituting these into (4.19), and writing T=t, we obtain the pair and We define the variable ξ is a strained coordinate which is useful in improving the analyticity of the solution. We then obtain an equation for w by differentiating (4.32a) and substituting for ht from (4.32b). This gives and the boundary conditions are the second of which implies ϕ=0 at ξ=ξf. Furthermore, this front position Ym, where h=0, is defined by To solve (4.34), we step forward h explicitly, and then calculate w by quadratureQ2. Although the equations appear degenerate, a local analysis of the behaviour near ξ=ξf shows that where wf=w(ξf,t).

There is some awkwardness in making the (finite difference) approximation near ξ=ξf. Define and The definition of ϕ in (4.34) is thus and is approximated in a simple finite difference scheme by where Δ is the step size and ξf=KΔ. The backward difference for the advective term is used if wi>0, otherwise the forward difference is used. In an explicit method, we only need to evaluate ϕ for steps 1 to K−1, but the quadrature for w also needs ϕ at i=0. From (4.40), we have , and therefore, using a 3-point approximation for h at ξ=0 (bearing in mind that h=0 there, so h0≈h1), we take

The prescription of ϕ from (4.41) is immediate, even if w<0, since we have hK=0. An alternative strategy uses the fact that ϕ=0 at ξ=ξf, together with (4.39), to show that we could then replace (4.41) by except that for i=K−1, we replace the approximation of J′ by

Figures 2 and 3 give two typical numerical solutions of (4.19) (in the canonical version with (4.26)). The solution in figure 2 is initialized with a small perturbation from steady state and that for figure 3 a much larger perturbation. It can be seen that both relax to the steady state as time progresses.

Figure 2.

Solution of (4.19), but using (4.26), small perturbation from steady state. (Online version in colour.)

Figure 3.

Solution of (4.19), but using (4.26), large perturbation from steady state. (Online version in colour.)

Discussion

The derivation of (4.19) from the reduced model (4.13) assumed that W=0, based on the equilibrium for s in (4.15). If instead we start from an ice flow which is out of equilibrium, for example, a uniform state with s=0, then (4.15) suggests and a rescaling of (4.13a) is then necessary. Taking the diffusion coefficient to be one, this is (bearing in mind the integral constraint and presuming that L∼O(1)) corresponding to new depth, breadth and time scales describing the size and evolution of streams of 1.3 cm, 38.5 m and 2.3 days. With this rescaling, the model takes the form (using (4.15)): The results of solving this are essential as before, except that the streams are a little deeper and less wide: on the short 2-day time scale, s≈0 is stationary and h tends to stationary finite width stream solutions of h+γN(h)=0.

In the final approach to the steady state, we have s≈N(h), and thus more accurately In that case, the model takes the form, with the original scales, together with the integral constraint, and the result is as before, albeit with the evolution on the much longer time scale of ice flow adjustment: the evolution of the ice roof acts as a huge brake on the development of the fluvial system.

Two complications merit discussion in our simple analysis of (5.5), or (4.13a), where we took A and |N′| as constant, and these are the consequences of taking more realistic forms for these functions.

The first concerns what happens as h→0. In reality, water is produced everywhere at the bed from the ice above and so the film thickness can never truly go to zero. This is manifested in the model through the closure relation (see (2.24)), which implies that N(h)∝(l*)−1, where l* is the supporting clast spacing. Effectively, this tends to zero as h→0, implying that , and this prevents h reaching zero. In addition, the deformable till thickness A in (3.6) becomes zero for N larger than a critical O(1) value, and so the bed becomes immobile, b is constant, and the water film thickness is determined geometrically from the solution of s=N(h)=h, assuming b=0.

Of apparently greater interest is what happens if N reaches zero, when also A=0. When this occurs, the ice surface separates properly from the bed and a stream or cavity forms. We consider briefly the dynamics of such streams. If N=0 the till flux and squeezing terms vanish, and the fluvial sediment transport terms become significant, on a longer time scale. If we bring back the fluvial terms in (4.1c) and ignore the till flux and squeeze terms, we find (again taking s=N and thus =0), that (5.5) is modified to where the integral term arises from the solution of (4.6b) correct to O(β/ν), and we have assumed a symmetric solution in which there is no water flow across the divide at Y =0; we have also taken for illustration (BS)′=1, although other choices are possible (in particular, we could have (BS)′<0 for a concave upper ice surface) [46,47]. The three terms on the right represent excavation of the bed by increasing downflow bedload transport, cross-stream bedload transport, and downbedslope sediment rolling, respectively.

Note that λ/β∼10−2, so that the natural length scale for the stabilizing bedslope term in (5.6) is corresponding now to a width of 166 m. We assume that h reaches a steady state, whose profile thus represents stream depth, and these steady solutions of (5.6) satisfy

Figure 4 shows numerical solutions of this equation, assuming a symmetric solution, and for various values of the centre-line dimensionless depth h0. A comment on the method of solution is necessary. We solve the equation by writing it in the form of a system with boundary conditions q(0)=w(0)=0, h(0)=h0.

Figure 4.

Solutions of (5.8) with hZ=0, h=h0 at Z=0, for h0=1.5,3. (Online version in colour.)

Of particular interest is whether significantly deep streams can be predicted from this theory: these would correspond to large numerical values of h0≫hc, where hc is the presumed critical value, where N reaches zero and separation occurs. Figure 4 indicates that as the maximum stream depth increases, the solution of (5.8) becomes approximately constant,[1] which suggests that large depth streams of constant depth are only possible providing hc≫1. However, this appears to require the presence of coarse clasts, which seems unlikely under the Ross ice streams. On the other hand, the prescription of a local clast spacing law as in (2.24) does not seem likely when streams form, and further discussion of the formation of genuine metre deep streams may require further elaboration of the film closure model.

Conclusion

In this paper, we have endeavoured to provide a mechanism for subglacial stream formation based on a first principles description of the interaction of ice, water and sediment. The results of our study suggest that for a Creyts–Schoof water film of thickness approximately 5 mm beneath an ice sheet, a uniform film is unstable, just as is the case for subaerial overland flow, and the model suggests that ‘streams’ of finite width will develop. While details will depend on the precise physical conditions, these streams appear to be more like swamps, having widths of hundreds of metres and depths of the order of a centimetre. Such a style of drainage is consistent with Alley's [10] notion of a patchy film beneath ice streams, and it is also consistent with Engelhardt & Kamb's [49] direct borehole observations beneath the Kamb Ice Stream of an ice/till gap of 1–2 cm in their borehole 5.

What we are less able to comment on is the observation by Engelhardt and Kamb of a 1.4-m deep water layer at borehole 9, indicative in their view of a subglacial channel. An alternative possibility is that this layer may represent a subglacial cavity formed by a developing drumlin field following the shutdown of the ice stream, but the issue remains as to the origin of the water. Concentration of an average 5 mm thick water film into a 1.4-m deep channel represents an amplification of 280 and such a stream of width, e.g. 5 m would drain all the water from subglacial melt from a surrounding width of 1400 m. While this seems feasible, our theory does not apparently provide the mechanism for such deep streams. On the other hand, the water-saturated till itself contains a large quantity of water: 8 m of till with porosity 0.4 has an effective film depth of 3.2 m! Following the Kamb Ice Stream's shutdown, the basal effective pressure would have risen, and this allows the water to be sucked out of the till. In that case, it would be possible to explain the presence of pockets of water, though not within the present theory, which ignores groundwater storage (cf. [50]). It is actually not difficult to include till storage in our model, allowing for infiltration to, or expulsion from the water film as the effective stress varies locally. However, the effect of the till is to act as a giant sponge, with the only dynamic effect being to provide a longer time scale for film thickness adjustment; essentially, this is because there is no net source associated with the till layer, it simply acts as a reservoir.

Even though the model is complicated, it does introduce simplifications that can be criticized. Examples of these are the assumptions of a constant viscosity for ice, a rheology for till deformation that assumes a yield stress based on solid friction and a prescribed effective viscosity, and the use of a simple linear bedload sediment transport law; however, in our view the most critical component of the model is the prescription of a closure relation between effective pressure and water film thickness, and it seems probably that more refined descriptions of this will be made in the future.

So what perspective does the present theory provide for present-day subglacial drainage? One of the things we know about subglacial water, at least beneath the Antarctic, is that there are numerous lakes, and that these form an ingredient of an active drainage system, in which lakes sporadically drain from one to the other [1,2,51]. It seems possible that basal water flows slowly through a patchwork of swamps, with occasional floods opening larger channels, although the dynamics of such floods remains to be elucidated. A further possibility, and one suggested by recent work of Christoffersen et al. [50], is that groundwater stored in subglacial sediment acts as a storage reservoir, which in conditions of stream formation allows expulsion of groundwater at topographic highs (where s≈N is large, h is small), but further examination of this possibility requires a specific inclusion of groundwater storage into the model.

The overriding motivation for studying the subglacial water system lies in a bid to explain feedbacks with the ice sheet flow. In this paper, we have concentrated on modelling on a small scale, with feedbacks from the water to the ice flow not included. Recent work by the authors has considered the feedbacks for a simplified version of the system presented here (not including sediment evolution) [41], but there is scope for a lot of work to be done in the future considering these feedbacks in more detail. Upscaling a model such as that presented here can be further used as a basis to develop subglacial water flow parametrizations for large-scale ice-sheet models.