Channelized subglacial drainage over a deformable bed

Abstract

We develop theoretically a description of a possible subglacial drainage mechanism for glaciers and ice sheets moving over saturated, deformable till. The model is based on the plausible assumption that flow of water in a thin film at the ice-till interface is unstable to the formation of a channelized drainage system, and is restricted to the case in which meltwater cannot escape through the till to an underlying aquifer. In describing the physics of such channelized drainage, we have generalized and extended Röthlisberger’s model of channels cut into basal ice to include “canals” cut into the till, paying particular attention to the role of sediment properties and the mechanics of sediment transport. We show that sediment-floored Röthlisberger (R) channels can exist for high effective pressures, and wide, shallow, ice-roofed canals cut into the till for low effective pressures. Canals should form a distributed, non-arborescent system, unlike R channels. For steep slopes typical of alpine glaciers, both drainage systems can exist, but with the water pressure lower in the R channels than in the canals; the canal drainage should therefore be unstable in the presence of channels. For small slopes typical of ice sheets, only canals can exist and we therefore predict that, if channelized meltwater flow occurs under ice sheets moving over deformable till, it takes the form of shallow, distributed canals at low effective pressure, similar to that measured at Ice Stream B in West Antarctica. Geologic evidence derived from land forms and deposits left by the Pleistocene ice sheets in North America and Europe is also consistent with predictions of the model.

Mathematical Symbols

• α Exponent in till flow law

• A _i Ice flow-law parameter

• A _s Till flow-law parameter

• b Exponent in till flow law

• b ₁,b ₂,b ₃ Numerical constants

• c Average suspended-sediment concentration as mass fraction

• C_p Specific heat of water

• C _t Rate of change of melting temperature with pressure

• d Ice thickness

• D _B “Characteristic” size of suspended sediment

• D ₅₀ Median grain-size of sediment

• D Deposition rate of suspended sediment

• Ē Rate of suspended-sediment erosion and entrainment

• f _R Dimensionless stream-bed roughness factor

• g Acceleration due to gravity

• h Stream depth

• [h] Magnitude of stream depth in metres

• h _T Heat-transfer coefficient

• H Conduit spacing transverse to ice-flow direction

• H _t Thickness of till

• k _w Thermal conductivity of water

• K Dimensionless constant in bed load-transport equation

• K _i,K _S Dimensionless constants related to conduit closure

• K ₁,K ₂ Dimensionless constants related to suspended-sediment transport

• l Perimeter of subglacial drainage conduit

• l _f Characteristic distance along glacier

• l _i Part of l bounded by ice

• l _s Part of l bounded by till

• L Heat of fusion of ice

• 
  Rate of meltwater input to basal till layer

• 
  Rate of melting of ice at tunnel wall

• 
  Rate of sediment erosion on conduit floor

• 
  Rate of sediment input from tributaries

• 
  Rate of meltwater input from tributaries

• n Exponent in flow law for ice

• Ν P_i - P_w

• N _c P_i - P_c

• ΔN _C Variation in N _c

• [N _c] Magnitude of N _c in bar

• Nu Nusselt number

• P _c Water pressure in subglacial tunnel

• P _i Ice-overburden pressure

• p _w Water pressure in till far from a conduit

• Δp _w P _w - P _c

• 
  Characteristic value of N _c separating channel and canal solutions

• Q Water flux

• [Q] Magnitude of Q in m³ s ⁻¹

• Q _b Bed-load flux

• Q _S Suspended-load flux

• R Conduit radius for idealized semi-circular geometry

• Re Reynolds number for flow in subglacial conduit

• R _f Dimensionless constant related to suspended-sediment transport

• s Curvilinear co-ordinate measuring distance along conduit

• S Conduit cross-sectional area, equal to S _i+S _s

• S _i Conduit area bounded by ice

• S _s Conduit area bounded by till

• t Time

• T Celsius temperature

• ΔT Temperature drop across boundary layer between water and ice at wall of conduit

• u Mean speed of water flow in conduit

• u _b, Glacier-sliding speed

• v _S Particle-settling speed

• V _i Rate at which ice flows into conduit

• V _S Rate at which sediment flows into conduit

• α Slope of glacier surface

• α_B Suspended-sediment concentration as volume fraction

• γ ρ_wC_pC_t = 0.32

• έ Shear strain rate

• θ Slope of glacier bed or stream bed

• θ _i,θ_s Angles characterizing conduit shape

• K _W Thermal diffusivity of water

• ^ Dimensionless measure of till permeability

• 
  Dimensionless parameter related to critical stress for bed-load transport

• ρ _i Density of ice

• ρ _s Density of sediment

• ρ _w Density of water

• Δρ ρ_s - ρ_w

• T Average shear stress exerted by stream on boundary of conduit

• T_b hear stress exerted by glacier on its bed

• T_c Critical value of T to mobilize bed-load

• 〈〉 Designates average over tunnel cross-section

1. Introduction

When seeking to understand the dynamics of ice sheets and glaciers, it is necessary to pose a sliding condition for the basal velocity when the basal ice reaches the pressure-melting point. This relates the velocity u_b to the basal stress T _b via a sliding law. It is now well recognized that this sliding law, which gives u _b as a function of T _b, also depends on the basal water pressure through the effective pressure N, defined as the difference between the overburden pressure p _i and the water pressure p _w

. (1.1)

This dependence has been elucidated both theoretically (Lliboutry, 1968) and observationally (Iken and Bindschadler, 1986). On the assumption that such a sliding law is available, it is then necessary to provide a prescription for N and this requires a theoretical model for the subglacial drainage pressure.

The classical theory of subglacial drainage is due to Röthlisberger (1972) and Shreve (1972), and envisages channels cut into the ice, and maintained open (with positive N, which promotes channel closure through ice creep) by melt-back due to the frictional heating by the turbulent water flow. This theory applies to hard beds, that is to say, the substrate over which the ice moves is considered rigid and impermeable. The same assumption is made in the classical theory of sliding. In the last several years, however, a radically different viewpoint has arisen. Observations of subglacial till deformation by Boulton and Hindmarsh (1987), Blake and Clarke (1989), Kohier and Proksch (1991), Blake (1992) and Humphrey and others (1993), along with the recognition that the rapid movement of Ice Stream B in West Antarctica may be due to shearing of weak subglacial sediment (Blankenship and others, 1986, 1987; Alley and others, 1986; Engelhardt and others, 1990), have caused a realization that the traditional theory of glacier sliding over an undeformable bed is of questionable relevance in some circumstances.

In this paper, we develop a model for subglacial drainage for the case in which the ice is underlain by a layer of sediment, which for convenience we will term till. We modify the Röthlisberger analysis to show that, in addition to channels carved into ice, “canals” can exist which are cut into the till, and which are maintained open by fluvial removal of inwardly creeping till. In the following section, we discuss the physical assumptions and ingredients which such a model requires. In section 3, we present a mathematical model for channelized flow, then show how it may be approximately “solved” in section 4. Our conclusions are presented in section 5. A list of mathematical symbols used is also provided at the beginning of this paper.

2. Physics of Drainage Through Till

Before addressing a quantitative model of subglacial drainage, we must consider the form that this drainage takes. A number of different possibilities have been considered in the literature. Most simply, drainage could take place by evacuation through a subglacial aquifer or by channels incised into the ice or the bed (here the till), through a system of linked cavities, or through a water film of variable thickness.

In the present analysis, we assume that there is no aquifer beneath the till and that all drainage is through the till or along the bed. The case in which there is a sub-till aquifer has been considered by Boulton and Hindmarsh (1987), Clarke (1987) and Shoemaker and Leung (1987).

A number of earlier workers (e.g. Boulton and Jones, 1979; Alley and others, 1986; Brown and others, 1987) have shown that for plausible permeabilities and thicknesses of subglacial till, Darcian flow through the pore space of a subglacial till layer israrely sufficient to evacuate all the meltwater flowing at the glacier bed. In certain special cases, meltwater may flow englacially (Hooke and others, 1988) or, if the subglacial till is being sheared rapidly enough, be simply adverted along with the till (Alley and others, 1986, 1987); in general, however, meltwater must escape via some sort of localized drainage at the glacier sole.

We consider that there are three main likely types of drainage system for sediment-floored glaciers: conduits incised in the ice (R channels); conduits incised into the till, which we henceforth call “canals”; and a water film of variable thickness (Weertman, 1972; Alley, 1989).

We will not offer the last word on films versus conduits in this paper but a few comments are in order. It is possible to envisage a water film of variable thickness at the ice/till interface, and unpublished data of H. Engelhardt and B. Kamb (cited by Kamb (1991, p. 16,589)) are consistent with (but do not require) existence of such a film beneath Ice Stream B in West Antarctica. However, the existence of a variable-thickness film is threatened by an instability due to viscous dissipation in the water. Walder (1982) demonstrated that for a rigid bed this instability leads to incipient initiation of R channels. There is no reason to suppose the same instability would not exist in the case of deformable bed. Various heuristic arguments can be advanced to suggest that such an instability will not in fact be manifested until the film thickness exceeds several millimetres; the point is debatable. Weertman and Birchfield (1983) argued that, despite the film instability pointed out by Walder (1982), channelized flow itself is unstable if all meltwater is subglacially derived, due to the supposed inability of channels to collect water from large lateral distances. We are skeptical of this conclusion, because it depends on the assumption of straight and parallel channels, and it is more natural to suppose that a network channel will not be parallel. Moreover, Weertman and Birchfield’s argument, which relies on peculiarities of the stress field around a channel, applies only for rigid, impermeable beds. If the bed consisted of permeable sediment, there would be less impediment to lateral water movement to channels.

For a sediment bed, another instability should come into play: the tendency for rill formation, just as in subaerial overland flow (Smith and Bretherton, 1972; Loewenherz, 1991). Thicker parts of a film will exert a greater shear stress on the sediment bed than thinner parts and will thus be relatively more efficient at eroding the bed and transporting sediment. Again, flow localization occurs. Thus, there are two instabilities tending to produce channelized flow over a sediment bed. Accordingly,we focus in this paper on the characteristics of a channelized drainage system over a deformable bed.

3. A Mathematical Model

We consider the situation shown in Figure 1, where a conduit exists at the ice/till interface, in general with parts cut into both the ice and the till. The ice-overburden pressure (and thus also the far-field till pressure) is p_i, while the water pressure in the channel is p_c. In general, we expect p_i > p _c, and our aim is to find an expression relating the effective channel pressure N _c = p _i – p _c to the water flux Q through the channel, as well as to the effective pressure N = p _i – p _w, where P _w is now understood as the pore-water pressure in the till far from the conduit.

Fig. 1. Idealized subglacial conduit at the ice-till interface. Ice and till move at speeds V_i and V_s, respectively.

When N _c < 0, the conduit tends to close by creep but is maintained open by melt-back due to viscous dissipation. In a similar way, it is recognized that wet till also creeps, with a yield stress on the order of 50 kPa. In many circumstances, the shear stress at the bed will far exceed this, so that the till also will tend to creep into the conduit. In this case, the till floor may be maintained by erosion and fluvial transport of sediment, either as suspended load or bed load. Our model is predicated on the effectiveness of these processes in maintaining a steady channel such as that in Figure 1. We now describe the several constituents of our model.

Creep of ice

The creep-closure rate of a tunnel in ice is given by

(3.2)

where R is the tunnel radius, and A _i and n are the constants in Glen’s flow for ice

where

and T are the second invariants of the strain-rate tensor

and the stress deviator tensor T_ij). This result applies directly to the closure of a cylindrical borehole (Nye, 1953) and is reasonably presumed to apply approximately in other geometric configurations.

Creep of till

Boulton and Hindmarsh (1987) proposed a flow law for till of the form

(3.3)

where

is the shear rate, T is the shear stress, N is the effective pressure P _i – p _w and A _s,a,b are positive constants. This law, based on a fit to seven data points, has been criticized by other authors (e.g. Kamb, 1991; Blake, 1992) but, in the present context, we consider it the simplest law which captures the idea that Equation 10 should increase with N, although the singularity at N = 0 is unphysical.

Boulton and Hindmarsh (1987) conjectured, following Nye’s (1953) solution for closure of a borehole in ice, that the consequent closure rate of a canal due to creep of till would be

where implicitly N_c = N, i.e. P _c = p _w With a = 1.33<b = 1.8 (as measured), sediment creep is then most effective at low values of N _c, which has an important bearing on the stability of the canal system, as we shall discuss later.

More recendy, Fowler and Walder (1993) have solved the Nye borehole-closure problem for the case of saturated, deformable till where in general p _c ≠ P _w. They showed that the closure rate depends on the parameter ^ a dimensionless measure of permeability. Estimates of ^ are in the range 10^-;6 < A < 1, for a permeability range 10^-;19–10^-;13m² (Freeze and Cherry, 1979, p. 29) and an apparent till viscosity of 10⁹Pas (cf. Blake and Clarke, 1989; Blake, 1992; Humphrey and others, 1993). For most tills, then, the asymptotic case

applies and the correct generalization of Equation (3.2). is as long as

, (3.4)

as long as

. (3.5)

If the condition Δp_w is not satisfied, the analysis predicts the existence adjacent to the hole of a “failed” zone in which piping probably occurs and permeability is locally enhanced. An analogous situation arises subaerially where streams are recharged by ground water: as stage falls following a flood, the hydraulic gradient driving seepage through the stream banks may become great enough that piping and bank failure occur (e.g. Bradford and Piest, 1980). Fowler and Walder (1993) considered how the solution to the borehole-closure problem must be modified to account for this effect; they showed the “corrected” closure rate is not greatly different from Equation (3.4)., so for simplicity we assume Equation (3.4). is generally valid in the rest of our analysis.

Pore-water pressure

In general, we expect p _w > p _c: meltwater generated at the ice/till interface makes its way towards conduits (although the opposite situation could occur transiently if, for example, flooding occurs via an englacial waterway). If the till is very permeable, or the melt supply is low, then the water will be evacuated towards conduits through the till. If these are spaced a distance H apart and the till is of thickness H_t, then the excess pore-water pressure necessary to drain water supplied at a rate

per unit area of ice-till contact is

where k is permeability and μ_w is viscosity (Shoemaker, 1986). Taking H _t = 10 m,

, μ_w = 2 x 10^-;3Pas, k = 10^-;16m² (a median value), then μp_w = 1 bar; for H = 100 m, μp_w = 100 bar. In general, μp_w/N may be significant; for example, under Ice Stream B, N is only about 1 bar (Engelhardt and others, 1990).

Various criteria may be offered for the initiation of conduits over a sediment bed (cf. Shoemaker, 1986). A physically reasonable criterion is that conduit spacing adjusts itself so that sediments at the edge of conduits are marginally stable against piping; from Equation (3.5)., this means N_c = aN. If conduits were closer, piping would occur, increasing the permeability of sediment near the conduit and thereby drawing down pore pressure p_w at the drainage divides until the condition N_c = aN was satisfied. The sediment creep law would then become

. (3.6)

In subsequent analysis, we will adopt the form of Equation (3.6)..

Sediment erosion and transport

The inward creep of till into a conduit is balanced by erosion of the sediment and its subsequent transport along the channel. We base our discussion on that of alluvial channels in non-cohesive sediment and give only the basic details: in particular, the formulae we use are very rough approximations.

Erosion occurs only if the shear stress T exerted by the water flow is larger than some critical stress T_c, which for sorted sediment depends on grain-size D_s. Beyond this critical stress, particles can be dislodged and they are transported downstream either as bed load if they are large enough (D_s > 0.1 mm) or as suspended sediment.

Various empirical expressions are available for bed-load transport; we take, for simplicity, the common Meyer-Peter and Müller (1948) equation, which we write as

, (3.7)

where (a;)_(x)+ = max(x, 0), Q_b is the bed-load-transport rate (mass per unit time), K is a constant approximately equal to 10, l_s is the width of the stream bed, g is acceleration due to gravity, ρ_w is the density of water and μρ = ρ_s - ρ_w.

Theoretical justification for expressions such as Equation (3.7). has been developed by Wiberg and Smith (1989). For stream beds with a mixture of sand-and gravel-sized clasts — probably a reasonable description for many subglacial streams (e.g. Engelhardt and others, 1990; Humphrey and others, 1993) — a condition of “equal mobility” holds, wherein all grain-sizes begin to move at about the same stress, namely,

(3.8)

where

is an empirically found, dimensionless constant, and D ₅₀ is the median grain-size (Parker and others, 1982).

Expressions for the rate of erosion of suspendable sediment are difficult to come by. Parker (1978a) emphasized that this is pardy because the usual theory of suspended-sediment concentration (e.g. Richards, 1982, p. 101) is a local, steady-state theory in which erosion rate Ė and deposition rate Ḋ balance exactly. In principle, however, Ė may be prescribed as boundary conditions on the concentration profile. Parker (1978a) proposed that the erosion rate is given by

(3.9)

where K ₁ ≈ 0.1 and v _s is the particle-settling velocity, and that the deposition rate is

(3.10)

where c is the depth-averaged suspended-sediment concentration and K ₂ ≈ 6. The net-erosion rate of suspendable sediment is then

. (3.11)

An important feature of subaerial stream beds is that they evolve their own shape through the processes discussed above. We can expect this to be true also for subglacial conduits cut into sediment. For non-cohesive sand/silt stream beds, Parker (1978a) argued that a stable bed shape could be maintained by a secondary lateral circulation of sediment, as follows. In the central parts of the stream, the shear stress acting on the bed is greatest and consequently greater erosion there leads to a larger suspended load than at the banks. This causes a lateral eddy diffusion of sediment towards the banks, which balances the bed-load transport that tends to bring sediment back down-slope toward the centre. Thus, even though the banks are being actively eroded, the river can maintain a constant cross-section via this dynamic equilibrium. In alluvial rivers, this leads to channel cross-sections with high aspect ratio (i.e. width l much greater than depth h). In gravel-bed rivers, the adjustment mechanism is more complicated. Since sediment is not being carried in suspension toward the banks, a stable channel cannot exist unless no erosion (and thus no bed-load transport) occurs on the banks. Accordingly, the shear stress must be great enough near the channel centre to transport the imposed sediment load but below critical on the banks. Parker (1978b) resolved this by showing that lateral turbulent flux of downstream momentum can produce the required shear-stress distribution. However, it is generally the case that T ≈ T_c for alluvial gravel-bed rivers. Again, large aspect ratios are typical.

Channels eroded into deforming subglacial till will differ from alluvial channels in one important aspect: whereas for an alluvial channel a steady-state cross-section exists only if there is no net erosion, for the subglacial channel a steady state will exist only if the net erosion balances the creep of till toward the channel. Thus, for a sand/silt-floored subglacial channel, there will not be a macroscopic balance between Ė and Ḋ, but rather

(3.12)

where the brackets

denote cross-sectional averages. For a gravel-floored subglacial channel, t will not drop below T_c on the banks, but rather T will exceed T_c by just enough to erode the banks and transport the inwardly creeping till at the appropriate rate. Accordingly, (Parker’s 1978a, b) analyses of the shapes of equilibrium stream channels cannot strictly apply to the subglacial case; in fact, his analyses should give lower bounds on the magnitude of bed shear stress. Qualitatively, however, the processes described by Parker for subaerial streams should apply to subglacial ones,hence we fully expect subglacial streams over non-cohesive sediment to tend towards high-aspect ratio cross-sections.

Governing equations

Having described the constituent ingredients of the process, we now write down governing equations for the flow. We consider a conduit at the ice–till interface (Fig.1). In general, we expect the conduit to extend both upwards into the basal ice and downwards into the till. We can thus define cross-sectional areas S _i and S _B of ice-or till-bounded conduits, respectively, and the total cross-sectional area S by

. (3.13)

The differential equations governing the physics are straightforward and are based on those elucidated by Röthlisberger (1972), Nye (1976), Spring and Hutter (1981) and Lliboutry (1983) but, here augmented by relations describing sediment creep, erosion and trans port. They must properly be framed in the context of a two-phase flow and we assume the simplest such model, that of homogeneous flow; in other words, we assume the sediment moves with the same velocity as the fluid. Let s be the downstream longitudinal coordinate and t be time. Equations of conservation of mass are, for the water,

(3.14)

and for the sediment

(3.15)

where

is the ice-erosion (i.e. melting) rate per unit downstream distance;

and

are, respectively, the water- and sediment-supply rates due to flow from tributaries, moulins or veins;

is the volume fraction of suspended sediment and Q_a is the suspended-sediment flux, defined by

(3.16)

where u is the mean flow rate; other symbols are as defined previously. We include the bed-load transport rate Q _b in Equation (3.15)., because we envisage the bed load as existing as a thin skin at the bed, occupying negligible volume.

Using Equations (3.2) and (3.6), the kinematics of conduit closure are described by

, (3.17)

and

(3.18)

where l _i and l _s are defined in Figure 1 and K _i and K _s and O(1) shape factors, partly representing the geometric difference between the conduit in Figure 1 and the idealized cylindrical shape assumed in the closure-rate calculations of Nye (1953) and Fowler and Walder (1993). Equations (3.17) and (3.18) are appropriate generalizations for non-cylindrical conduits. If the conduit is actually cylindrical, then

and the creep-closure rate

as expected. For a wide, flat channel in ice (cf. Hooke and others, 1990) with height h and width

(thus l≈l _i), we expect

thus

, because the average closure rate depends only on the width of the “hole” into which ice is flowing. A similar argument applies to a wide, shallow conduit incised into till.

Neglecting inertial effects, the momentum balance for the flow can be replaced by a force balance between pressure forces and drag on the channel walls,

(3.19)

where sin θ is the bed slope and l is the channel perimeter, defined by

. (3.20)

Assuming turbulent flow, which is valid except for very small tunnels (cf. Weertman, 1972), the average wall shear stress is given by

(3.21)

where the dimensionless friction factor fR is weakly dependent on the Reynolds number and is typically ~0.1 for subaerial streams. In general T_w (the average shear stress over the sediment bed only) and T may differ but we assume for simplicity that they are the same. This seems reasonable for pressurized flow but may be less realistic for open-channel flow.

We write the energy equation in the form

(3.22)

where Cp is specific heat, T is temperature, ΔT is the temperature drop across a thermal-boundary layer adjoining the ice wall (but not the stream bed, the temperature of which is unconstrained) and h _T is a heat-transfer coefficient. We neglect work done in eroding and transporting sediment (Bagnold, 1961). The temperature T will be related to the temperature of the ice, which we expect to be at the pressure-melting point, by

(3.23)

where c _t is the rate of change of melting temperature with pressure and ΔT is the temperature step across the thermal boundary layer at the ice/water interface. The rate of heat transfer across the thermal boundary layer is given by

(3.24)

where L is the latent heat and we neglect any small temperature gradient in the ice (cf. Weertman, 1972).

The Equations (3.14), (3.15), (3.17)–(3.19) and (3.22) are six differential equations for the variables

and T; in addition, the seven algebraic relations in Equations (3.7), (3.10), (3.13), (3.16), (3.21), (3.23) and (3.24) determine the subsidiary variables

. Conduit perimeter I is determined by Equation (3.20)., so provided that l _i and l _s are prescribed functions of S _i and S _s respectively, the model is complete.

4. Analysis

Our aim is to solve the model posed in section 3 to determine the effective channel pressure, N _c = p _i – p _c. This depends on the inputs and most notably on the water supply. Specifically, we can, in principle, determine the water flux Q by integrating Equation (3.11). and then N_c is determined in terms of Q. In what follows, we consider only steady-state solutions. We hope this provides some average description of the flow, although water pressure often fluctuates very irregularly, at least beneath alpine glaciers.

Critical effective pressure

From (3.17) and (3.18), we have at steady state

, (4.1)

(4.2)

whence the relative rate of ice and till creep is given by

(4.3)

where we define the critical effective pressure to be

(4.4)

Taking K _i ≈ K _s, l _i ≈ l _s, and with ρ_i = 9 x 10² kg m^-;3, ρ_s = 2.65 x 10³ kg m^-;3, A _i = 7.36 × 10_-;24 Pa^-;3 s^-;1 (cf. Lliboutry, 1983), n = 3 and A _s = 3 × 10^-;5 Pa^b-as^-;1 (after Boulton and Hindmarsh, 1987), a = 1.33, b = 1.8 we find

. (4.5)

Thus, for

, till creep is larger, while for

, ice creep is. This leads us to consider two end cases, when either one or the other creep process dominates conduit closure.

Conduit geometry

In order to proceed, we need to describe the “constitution” of the conduit. Since the conduit in Figure 1 has a boundary partly of ice and partly of till, the respective closure rates will generally be different, i.e.

. Since both till and ice are viscous, a boundary condition applies at their interface. In what follows, we assume a no-slip condition. In this case, the ice and sediment speeds towards the conduit at the interface, V _i and V _s, must be equal. In a slightly crude way, we identify

(4.6)

where k = i or s and θ(θ_s) is the angle of the ice roof (till flow) to the horizontal. It then follows that

(4.7)

and thus the two extreme cases

correspond to values

, i.e. level-floored ice channels and flat-roofed till canals, respectively. That is to say, if

, the conduit is essentially a Röthlisberger channel in ice, floored by a relatively stiff till, while if

it is essentially a canal cut into soft till and roofed by relatively stiff ice.

Note that the determination of

depends critically on the assumed till rheology. However, the identification of two asymptotically distinct conduit shapes is robust.

Determination of water pressure–flux relationship

The way to solve the equations is now as follows. First, we recognize that, based on observations c

so

, and the volume flux Q = uS is determined from Equation(3.14).. From Equation (4.1)., it is easy to calculate that melting due to frictional dissipation is likely to be eligible in most circumstances, and therefore the water flux Q(s) in a tunnel is essentially determined by the input from surface-water and basal tributaries. In effect, we can consider Q(s) to be a prescribed function. Equation(3.15). determines the suspended-sediment load but this is not needed for the solution of the rest of the problem.

Using Equations (3.23) and (3.24), the energy Equation (3.22). can be rewritten in the form

(4.8)

where

(Röthlisberger, 1972), K _w = k _w/p_wC_p is the thermal conductivity. The Nusselt number Nu, a dimensionless measure of heat transfer at the ice–water interface, is defined by

. Equation (4.8). determines

in terms of T, u and p _c, and in addition we have Equations (4.1), (3.19), (3.21) and

(4.9)

to determine the five unknowns S, u, T, p_c and m, (given l _i and l_s ). It is straight-forward to combine these to form a single, non-linear, second-order differential equation for Pc but it is probably more illuminating, and certainly easier, to simplify the system using two judicious approximations. For the first of these, note that, since N_c = p _i – p _c and P _i = ρigd, where d is the ice depth, then

(4.10)

where α is the ice-surface slope (Weertman, 1972). For bedrock, that varies sufficiently slowly, the term

will be negligible in this expression. For example, if the bedrock slope changes significantly over a length l _f along the flowline, then

(4.11)

where ΔN _C is the variation of N _c. For an alpine glacier, we can expect ΔN _C <

is negligible if

, typically about 10 ice thicknesses. On an ice sheet, we might expect

, but for ice streams, at least, geophysical evidence from (Blankenship and others 1986, 1987) suggests strongly that

. In either case, it is therefore useful as a rough estimate to neglect

in Equation (4.10).. A similarly small error is made if we then approximate the hydraulic gradient using the surface slope, thus

. (4.12)

This approximation breaks down near the glacier terminus, where variations in p _c are substantial (Fowler, 1987).

The second simplification results from Equation (4.8).. We note that Q/k _wNu is dimensionally a length, and is in fact the length over which

relaxes to its equilibrium value. Taking values Q = 1 m³ s^-;1, K _w = 10^-;6 m² s^-;1, Nu ≈ 0.1 Re⁰⁸ (Spring and Hutter, 1981), Re ≈ 10⁶ (based on u = l ms^-;1, l = 1 m), so that Nu ≈ 10⁴, we have Q/k _wNu ≈ 10² m, so the derivative of

in Equation (4.8). is negligible almost everywhere. The other derivative term in Equation (4.8). is

. (4.13)

It is consistent with our previous approximation that

and thus using p _c = p _i – N _c,

. (4.14)

Finally, then Equation(4.8). can be approximated by

(4.15)

and we have only to solve algebraic equations to find N_c, since p _i and Q are prescribed functions of s. The square-bracketed term in Equation (4.16). represents that part of the viscously dissipated energy that does not go into melting the ice roof. Röthlisberger (1972) showed that for the simplest case of a horizontal tunnel and constant discharge Q, this term would be a fraction γ of the total viscous dissipation. More generally, simple estimates of the square-bracketed term in Equation (4.16). suggest it is typically about 0.1–0.3 times the term Tul/L, so for algebraic convenience we neglect it.

Model summary

We now simply have five algebraic equations:

(4.16)

we take l _i ≈ l and

to find

(4.17)

where

(4.18)

Thus, in order finally to determine N_c in terms of Q, we must prescribe l in terms of S. It is at this point that the conduit shape becomes critical. Specifically, we now consider the two distinct cases

(channel) and

(canal), and seek to prescribe l(5) appropriately in each case.

Channel drainage

In the case

, the tunnel protrudes into the ice, forming a Röthlisberger-like channel. We then expect a more or less semi-circular cross-section, and accordingly we choose

. (4.19)

Equation (4.17). then yields

, (4.20)

which is essentially Lliboutry’s (1983) form of the standard Röthlisberger result. Note that as Q increases, N_c increases and thus p_c decreases, which is consistent with formation of an arborescent network of channels, as we discuss below. Because N_c depends weakly on Q, we can reasonably take Q ≈ constant to evaluate N_c.

Canal drainage

Suppose now that

. In this case, drainage is through a canal incised into the sediment. The mechanics of sediment erosion and transport now become very important in determining the cross-sectional shape. In particular, canals are typically wide and shallow, in which case Equation (4.19). no longer applies. If h is the canal depth, then

. (4.21)

and the determination of h depends on the precise nature of the stream bed, gravel, sand/silt and clay beds having different prescriptions. We discuss these subsequently. For now, suppose that h is prescribed in terms of the bed composition. Then, using Equation (4.17)., we have

, (4.22)

Note the crucial result that, in contrast to Röthlisberger channels, canals incised into sediment lead to a water pressure that increases with increasing water flux. This is consistent with a distributed rather than an arborescent drainage network as for flow in cavities (Walder, 1986; Fowler, 1987; Kamb, 1987). We emphasize that the physics leading to a non-arborescent subglacial canal system is distinct from what is involved in the instability causing braiding of subaerial streams. That instability involves details of the sediment-transport process, includ ing non-steady phenomena such as growth and destruction of sediment bars and passage of bed-load pulses (e.g. Parker, 1976; Ashmore, 1991). The predicted non-arborescent nature of a subglacial canal system, on the other hand, is simply the result of the relation between Q and N_c.

Numerical estimates

The coefficients b ₁, b ₂ and b ₃ depend only on sin α, and taking n = 3 they are given by

(4.23)

Assuming

, Röthlisberger-like channels exist, and the relation in Equation (4.20). between N_c and Q may be expressed numerically as

(4.24)

where Q = [Q] m³ s^-;1 and N _c = [N _C] bar. Inverting Equation (4.24). gives

. (4.25)

Therefore, for sin α = 0.1, [Q] = 1, we have [N_c] ≈ 38; for sin α = 10^-;3, [N_c] ≈ 4. Thus, we can expect that in fact

for valley glaciers, although open-channel flow may be common, as pointed out by Hooke (1984). The calculated value of N_c is, however, probably less than

for very small slopes, contrary to assumption, so R channels are unlikely to exist beneath ice streams and ice sheets, where these are underlain by deforming sedimentsor till.

If, on the other hand, we assume

, then the expression in Equation (4.22). for water flux is, numerically,

(4.26)

where h = [h] m. For [Q] = 1, we find [N _c] = 85 [h] for sin α = 0.1, [N _C] = 4 [h] for sin α = 0.001. Thus, for steep valley glaciers, canals are only viable if h < 10 cm, while for ice sheets, canals can exist for stream depths up to about 2 m. In order to estimate what h is, we must consider further the mechanisms which determine it.

Determination of channel depth

Gravel

As pointed out previously, alluvial gravel-bedded streams are wide and shallow, and adjust their shapes so that the average shear stress T exerted on the stream bed is approximately equal to the critical stress T_c for bed-load transport, given by Equation (3.8). for mixed grain-size beds. Thus, from Equations (3.8), (4.17)₂ and (4.21), we have

(4.27)

for a subglacial gravel-bedded stream. The net input of sediment by creep actually requires a slightly higher average shear stress, as argued in section 3. A crude estimate using the Meyer-Peter and Müller bed-load transport Equation (3.7)., together with measured values of

on alpine glaciers (Gurnell, 1987), suggests that expression (4.27) need not be significantly modified. Thus, we have h/D₅₀ = O(1) for alpine glaciers (sin α~0.1)— an estimate consistent with observations on steep alluvial streams — and h/D50~O(10²) for ice streams (sin α~10^-3).

Sand/silt

Here the sediment is non-cohesive, but suspended-sediment transport is important. We again expect broad, shallow channels with S ≈ lh, much as in the alluvial case. An estimate of the flow depth comes from Parker’s (1978a) analysis of sand/silt bedded streams:

(4.28)

where D _s is, after Parker (1978a), a “characteristic” size of suspended grains. (In his comparison of theory to field data, Parker chose D_s = D ¹⁵, i.e. the grain-size such that 15% of all grains were smaller.) Rf is a dimensionless parameter that depends on grain-size and is typically O(10^-;1) for fine sand and coarse silt; hence

for

for sin α = 10^-;3.

Clay

In contrast to a non-cohesive sediment, canals cut into a cohesive bed can reasonably be expected to maintain steep banks. In this case, the cross-section may be more or less semi-circular, and we obtain the Röthlisberger result (S ≈ l²) once again. However, if the sediment layer is thin and the canal erodes through it to bedrock then h is given by the thickness of the till layer and we would again have S≈lh.

Subglacial till

Till typically has a bimodal grain-size distribution (e.g. Drewry, 1986) with a matrix rich in clay and silt-sized grains, presumably derived from comminution processes, supporting coarse clasts. The cohesive nature of the matrix may allow steep banks to occur, but since till creep and erosion also occur, we can plausibly expect that a lag will tend to form at the stream bed, comprising sandy and gravelly material. The presence of such a layer will most likely cause the canal to evolve towards a broad, shallow shape.

At the risk of generalization in the absence of much observational verification, we therefore suggest that, on the whole, canals incised into subglacial till will be broad and shallow, with a prescribed depth h that is a function of slope and grain-size D _s (and, probably, discharge), with

. Although we have applied the best present understanding of stream mechanics and sediment transport, our estimates of h are poorly constrained, with likely values in the range 10 mm-1 m, the value increasing as sin α decreases.Specifically, for gravelly till with D _s~10 mm, 1 cm < h < 1 m, while for sandy/silty till with D _s~0.1 mm, 1 cm < h < 10 cm.

Canal effective pressures

Since h will depend on sin a as well as D_s and bed type, it is impossible adequately to summarize the resulting variations in N _c. Nevertheless, the gist of the above discussion can be summarized in Table 1, where we take Q = 1 m₃s^-;1, We see that

for both ice sheets and glaciers, so that we could expect a canal system beneath ice sheets, whereas for valley glaciers it appears at first sight that either R channels or canals might exist. We will have more to say on this matter below.

Table 1.

5. Discussion

The results of our analysis so far are summarized in Figure 2, which indicates the two limiting cases of canals

or channels

. In particular, both types of solution can exist for large enough slopes (α~10^-;1), while only the canal solution will exist for very small slopes (α~10^-;3). It is particularly important to under stands how the rest of Figure 2 can be completed, since for given Q there are apparently two stable solutions when α is large.

Fig. 2. Plots of effective pressure vs discharge for Röthlisberger-like channels and canals, for several values of surface slope (sin α = 10^-;3, 10^-;2, 10^-;1). The curves for R channels intersect the line

at lower values of discharge than do curves for canals.

Stability

In discussing stability, we distinguish between the stability of an individual solution branch (here used in a mathematical sense), i.e. the stability of a particular conduit to perturbations in that conduit, and what we will term collective stability, which is the stability of a particular drainage system against perturbation to another system, e.g. two (or more) channels, or two or more canals. First we discuss branch stability.

We suppose that each end-member solution (channel or canal) derived above is stable; we have no reason to suppose otherwise. Now consider what happens as Q increases from zero. For very small Q, Darcy flow through the till will suffice to drain all melt water. As Q increases, the till pore-water pressure will increase until N= 0, at which point flotation is initiated and a water film forms at the ice/till interface. As shown by Walder (1982), this sheet flow is always unstable, and we therefore expect that a bifurcating branch representing interfacial conduits exists with N_c > 0. This situation is portrayed in Figure 3.

Fig. 3. Probable conditions for very low values of discharge. Darcy flow can drain all meltwater for sufficiently low discharge. For Q = Q _c, sheet-water flow commences at zero effective pressure (broken line) but this drainage configuration is always unstable relative to channelized drainage.

At larger Q, if canals and channels can exist, then we surmise that these solutions arise as bifurcating branches from an (unstable) solution in which

, with conduits more or less equally incised into ice and till. This branch will be the continuation of the bifurcating branch at low Q. The situation is portrayed in Figure 4. We have explored other possibilities, but this one seems the likeliest. There are certain immediate implications. In the case that channels and canals can coexist (as in Figure 4 ), then because the channel-solution curve intersects the line

at a lower value of Q than the canal curve, channels will be the preferred solution as Q is increased. Moreover, if canals and channels were to coexist at larger Q (each system being in itself stable), then because the effective pressure is higher for channels (and hence the water pressure is lower), leakage through the till will drain water from canals towards channels. (This is an example of collective stability of channels.) The net effect of these factors is that if channels can exist, they will be the preferred stable drainage style, and this is likely to be the case for valley glaciers. If the surface slope is sufficiently small, as for ice sheets or ice streams, canals will form the stable drainage network.

Fig. 4. Qualitative picture of the entire pressure vs discharge plot. The line

probably coincides (for large enough values of Q) with an unstable solution branch representing conduits inched equally into ice and till, and joins with the incipient-channel solution branch beginning at Q = Q_c. For small slopes (sin α~10^-;3), the channel solution does not exist, whereas for large slopes (sin α~10^-;1), the canal solution does not exist. At intermediate slopes (sin α~10^-;2), both solutions exist but canals are always unstable relative to channels.

Hydrological and geomorphological implications

Our analysis predicts that drainage at the base of ice sheets and ice streams should take place via canals of depth ~10 cm (plus, presumably, Darcy flow) at effective pressures of only a fraction of a bar, for fine-grained sediments. This seems consistent with basal water-pressure values at Ice Stream B, West Antarctica, inferred on the basis of geophysical measurements (Blankenship and others, 1986, 1987) and subsequently confirmed by borehole measurements (Engelhardt and others, 1990), although we should emphasize that actual data on the geometry of the basal drainage system is practically nil. Our qualitative prediction is independent of the exact rheological properties of till; as long as there is some value of effective pressure below which till can creep more quickly into conduits than can ice, distributed canals as we have described will be the stable drainage type. Strictly speaking, this requires

. The results of Fowler and Walder (1993) give

, but this particular value is by no means critical.

There is indirect sedimentological evidence for the existence of broad, shallow canals beneath ice sheets. Brown and others (1987) studied deposits left by the Puget lobe of the Pleistocene Cordilleran ice sheet. Their detailed stratigraphie work in the Vashon Till of the Puget Sound region, Washington State, U.S.A., revealed the existence of lenses of water-laid, sorted sediment within the till. Brown and others interpreted these as having been deposited in subglacial conduits and estimate the width of such conduits as ranging from 0.1 to 8 m. They further inferred, from several lines of evidence, that water pressure beneath the Puget lobe was very near the ice-overburden pressure. Such low effective pressures in channelized flow would be consistent with canals but not with R channels in the ice. We further note that the discontinuous character of the water-laid subglacial deposits described by Brown and others (1987) is consistent with a picture of broad, braided canals (similar to alluvial gravel-bedded streams), reflecting the constantly shifting geometry of the braids.

In northeastern England, an area of fine-grained till, Eyles and others (1982) described “shoestring” deposits, which they interpreted as channel fills, and which may also have been deposited within a canal system. These deposits are elongated sub-parallel to the former direction of ice flow. They have nearly flat upper surfaces and irregular, concave-up lower surfaces, and are much wider than thick. Gravel lags at the deposit base are overlain mostly by sand but also include laminated clay. Eyles and others (1982, p. 323) proposed that these deposits represent cut-and-fill structures formed by subglacial drainage, with grain-size variations within the deposit indicating variable discharge, including “ponding episodes possible through closure of the channel”. Diapiric deformation at the channel bases indicates low effective pressure in the till at the time of channel formation. Lenses of till that occur within some of these deposits may represent sediment masses that dropped into the channels from overlying ice.

Another line of geologic evidence supporting our picture of drainage beneath ice sheets comes from the distribution of eskers over the region of North America formerly covered by the Laurentide ice sheet. Eskers are very rare in areas where the ice flowed over thick till but ubiquitous where the ice-sheet bed comprised bedrock with only thin, discontinuous till cover. We would predict that, in the former instance, the stable drainage system would have comprised widespread, shallow canals incised downward into the till, whereas in the latter instance, with an essentially rigid bed, we would expect an arborescent R channel network, as in the classical view of subglacial drainage. Based on observed morphological and sedimentological aspects of eskers, we believe their formation would be very unlikely in the distributed-canal system, but definitely possible with the largest channels of an arborescent network, in accord with the observed distribution. This argument has been elaborated by Clark and Walder (in press).

Shoemaker (1986) and (Alley 1989, 1992) have previously presented analyses of channelized drainage for sediment-floored glaciers but with different results. Shoemaker (1986) never considered creep of till into sediment-floored tunnels; furthermore, he did not consider the way in which sediment-transport mechanics affect conduit shape. Consistent with the neglect of these factors, Shoemaker further assumed that the relation between Q and N _c for sediment-floored conduits is simply Röthlisberger’s (1972) result. We have demonstrated that this is not necessarily correct. (Alley 1989, 1992) dealt more carefully with the issues of till creep and sediment transport, but did not take into account the fact that sediment-transport mechanics may strongly affect conduit shape.

Finally, we remark that the present analysis does not point at any sort of glacier-surge mechanism associated with a switch between R channel and canal drainage. Unlike the case of R channels and cavities over a rigid bed (Fowler, 1987), we do not expect the two different types of drainage systems to coexist stably. Additionally, nothing in the present analysis indicates a dependence of the drainage-system characteristics on the ice-sliding speed (as for cavities). A surge mechanism that involves a switch between R channels and cavities is still conceivable for a glacier resting on deformable sediment. The glacier bed will not be flat and cavities can still exist to the lee of relatively large rock clasts that protrude into the ice. Thus, the potential for a linked-cavity system still exists for a deformable sediment bed, as pointed out by Kamb (1987).