Subglacial Hydrology for an Ice Sheet Resting on a Deformable Aquifer

Abstract

Subglacial hydrology is investigated for an ice sheet where the substrate consists of a deformable aquifer resting on an aquitard. If sliding velocities are low or absent, subglacial melt-water drainage is dominated by drainage through the aquifer to water channels. Drainage along the bed is negligible. Efficient melt-water drainage requires that a system of subglacial water channels exists; otherwise, pore-water pressures will exceed the overburden pressure. In general, aquifer deformation near (away from) the terminus is most likely to occur during the winter (summer). The effect of short-term high channel pressures is, in general, not critical to aquifer deformation because the pressure pulse does not propagate far into the aquifer. (For aquifers of high permeability, short periods of high channel pressures constitute the most critical condition.) Aquifer deformation at the terminus is very likely to occur if the terminus ice slope exceeds tan ϕ, where ϕ is the Coulomb friction angle of the aquifer material. Upwelling of basal melt water near the terminus will normally cause soil dilation if the aquifer has a low permeability (e.g. till). Maximal profiles are computed corresponding to various aquifer materials using channel spacings which provide efficient drainage. (A maximal profile is the highest ice profile which the aquifer can sustain without deformation.) In general, maximal profiles lie well above observed profiles (such as h(x) = 3x^1/2 (m)) except near the terminus. However, if channel spacings are sufficiently large, pore-water pressures are increased and maximal profiles can lie well below h(x) = 3x^1/2.

Résumé

On étudie l'hydrologie sous glaciaire d'une nappe de glace dont le substratum est un aquifère déformable reposant sur une couche imperméable. Si la vitesse de glissement est faible ou nulle le drainage sous glaciaire des eaux de fusion se fait, à travers l'aquifère, vers les chenaux. Le drainage le long du lit est négligeable. Un drainage efficace des eaux de fusion nécessite l'existence d'un réseau de chenaux sous glaciaires, sinon la pression intersticielle dépasserait la pression lithostatique. En général la déformation de l'aquifère près (loin) du front a plus de chances de se produire en hiver (été). Les épisodes brefs de haute pression dans les chenaux n'ont généralement pas de conséquences critiques sur la déformation des aquifères parte que l'onde de pression ne se propage pas trts loin dans l'aquifère. (Par contre pour des aquifères très perméables les courtes périodes de haute pression dans les chenaux représentent les conditions les plus critiques.) La déformation des aquifères au front du glacier a le plus de chances de se produire lorsque la pente de la surface de la glace dépasse tg ϕ où ϕ est l'angle de frottement de Coulomb du matériau constituant l'aquifère. La remontée des eaux près du front doit provoquer une dilatance du sol si l'aquifère a une faible perméabilité (ex: argile à blocaux). On calcule les profils maxima correspondant à différentes compositions de l'aquifère et une densité de chenaux permettant un drainage suffisant. (Un profil maximum est le profil le plus épais que puisse supporter un aquifère sans se déformer.) En général les profits maxima sont bien au-desus des profils observés (comme sauf au voisinage du front. Néanmoins si les chenaux sont suffisamment espacés les pressions intersticielies augmentent et les profits maxima peuvent-être bien en dessous de

Zusammenfassung

. Für ein Eisschild, dessen Untergrund aus einer deformierbaren wasserführenden Schicht, aufgelagert einem Stauhorizont, besteht, wird die subglaziale Hydrologie untersucht. Wenn die Gleitgeschwindigkeiten gering sind oder ganz fehlen, wird das Schmelzwasser hauptsächlich über die wasserführende Schicht in Wasserkanäle abgeleitet. Der Abfluss über das Bett ist vernachlässigbar. Wirksamer Schmelzwasserabfluss erfordert das Vorhandensein eines Systems subglazialer Kanäle; sonst überschreitet der Porenwasserdruck den Auflagedruck. Im allgemeinen deformiert sich die wasserführende Schicht nahe bzw. fern dem Gletscherrand meist im Winter bzw. Sommer. Die Auswirkung kurzzeitiger hoher Drucke in den Kanälen ist im allgemeinen unkritisch für die Deformation, weil der Druckimpuls nicht welt in die Schicht eindringt. (Für wasserführende Schichten mit hoher Durchlässigkeit stellen kurze Perioden hohen Kanaldrucks die kritischsten Zustände dar.) Deformation am Eisrand tritt mit grosser Wahrscheinlichkeit dann ein, wenn die Neigung des Randeises tan ϕ überschreitet, wobei ϕ der Coulomb'sche Reibungswinkel des wasserführenden Materials ist. Vom Untergrund aufquellendes Schmelzwasser nahe dem Eisrand wird gewöhnlich zur Ausdehnung des Bodens führen, wenn die wasserführende Schicht nur wenig durchlässig ist (z.B. Geschiebe). Mit Kanalabständen, die einen wirksamen Abfluss ermöglichen, werden Maximal-Profile für verschiedenes wasserführendes Material berechnet. (Ein Maximal-Profil ist das höchste Eisprofil, das die wasserführende Schicht ohne Deformation tragen kann.) Im allgemeinen liegen die Maximal-Profile weit über beobachteten Profilen (wie etwa h(x) = 3x^1/2(m)), ausser in der Nähe des Eisrandes. Überschreiten jedoch die Kanalabstände ein gewisses Mass, dann steigt der Druck des Porenwassers, und die Maximal-Profile können durchaus unter h(x) = 3x^1/2 liegen.

Introduction

Observations by Boulton and Jones (1979) show that a large proportion of the forward motion of a temperate glacier lying on a deformable bed (e.g. glacial till or sand) can be contributed by the deformation of the bed. In modelling bed deformation they assumed that the evacuation of basal melt water was by horizontal seepage through the permeable bed which is assumed to rest upon an impermeable substrate. It follows from Darcy's law that, for horizontal seepage, pore-water pressure differences depend linearly on (i) the horizontal distance travelled, and (ii) the reciprocal of the thickness of the permeable layer. Thus, if calculations are made for permeable beds a few meters thick, pore-water pressures higher than the pressure of the over-riding ice are developed over a distance of order 50 km (Boulton and others, 1974). (Boulton and Jones did not specify bed thickness in their calculations.)

Even for highly permeable beds similar results are obtained. For example, a 2 m thick gravel bed evacuating basal melt water, with a melt rate of 3 cm a⁻¹, will develop a pore pressure of at least 500 m (water equivalent) within 150 km of the divide. In this calculation the evacuation of surface melt water which reaches the base was not considered. This could increase the subglacial flow rates and pressure gradients by several orders of magnitude. Such results suggest that, at least for ice sheets, basal melt water cannot be evacuated exclusively through the permeable bed. The approach taken here is that water channels participate in the drainage of basal melt water. (Sheet flow will later be ruled out.)

Normally permeable beds exist in layers. There are two situations of particular interest: an upper layer of relatively high permeability (aquifer) resting upon a lower layer of much lower permeability (aquitard), and the reverse. The modelling and results of these two cases are so different that we will consider only the first case.

Examples of aquifers upon aquitards and representative ratios of permeabilities are: medium sand upon till, ∼10⁵; gravel upon medium sand, ∼10³; medium sand upon clay, ∼10⁵; gravel upon fractured igneous and metamorphic rock, ∼10⁴ (Freeze and Cherry, 1979). For any such combinations, if the aquitard is continuous and sufficiently thick, we can assume that the aquitard does not participate in the drainage and therefore represents an impermeable barrier to water flow. This was the situation considered by Boulton and Jones.

It will be shown that pore-water pressures resulting from the channel model can be much lower than those predicted by the Boulton–Jones model. Moreover, provided the distance between channels is not too large, pore-water pressures are close to channel pressures. (Beds of low permeability such as clay or till, which are really aquitards, represent special cases of difficult drainage.) Under constant-state conditions, channel pressures are higher in the winter than in the ablation season. Thus, the critical constant-state condition where the bed is most likely to deform can be expected to occur in the winter. However, the diurnal fluctuation in channel flow rate during the ablation season can change this conclusion. This causes an increase in mean channel pressure away from the terminus. A critical transient condition occurs in the spring when water channels of relatively small diameter are suddenly forced to evacuate surface melt water. This complex problem is not analyzed extensively but qualitative remarks are made.

Ice-sheet profiles are frequently closely approximated by the parabola

(1)

where h is height above the horizontal base and x the distance up-stream from the terminus. Coefficient A has a value of 4.7 m^½ for both the Antarctic ice sheet inland from Mirniy (Hollin, 1962) and the West Greenland ice sheet (Mathews, 1974). These profiles can be produced theoretically by assuming no basal sliding and a basal shear stress of 1 bar (Orowan, 1949; Nye, 1952). It has been estimated that several examples of late Pleistocene ice-sheet lobes within mountainous terrain of North America and New Zealand had A values ranging from 2.9 m^½ to 4.1 m^½; several ice-sheet lobes in the south-western part of the late Pleistocene Laurentide ice sheet had A values estimated from about 0.3 to 1.0 m^½ (Mathews, 1974). These latter values correspond to a basal shear-stress range of 0.004–0.045 bar, which is extraordinarily low.

It is generally believed that low basal shear-stress values imply that the deformable bed is the weak link in the bed-ice system. One focus here will be an assessment of the importance of various bed and channel parameters in determining maximal ice-sheet profiles which the beds will sustain without deforming. These maximal profiles will be compared with profiles satisfying Equation (1).

Clarification should be made of the basal temperature of continental ice sheets. On theoretical grounds as well as from field data, it is known that the combination of the insulation property of ice and geothermal heating from below can produce temperate basal conditions in an otherwise cold glacier, even in central regions of the present Antarctic and Greenland ice sheets (Paterson, 1981). The presence of extensive areas of subglacially deformed sediment in regions occupied by Quaternary ice sheets is one indication that the beds of these ice sheets were at the pressure melting-point during long periods, e.g. Slater, 1926; Kupsch, 1962; Woldstedt, 1965; Lavrushin, 1971. The near-terminus regions of ice sheets located far inland would likely have had cold basal conditions (Wright, 1973; Moran and others, 1980). We shall restrict our consideration to ice sheets with temperate bases throughout. This should have been the condition of the many marine ice-sheet lobes such as the Puget lobe of western Washington State, U.S.A. There, the mean annual temperature was positive during the maximum of the last glaciation (Porter, 1977).

Sliding or Deformation?

There is observational evidence of ongoing shearing in subglacial deformable beds; Boulton and Jones (1979) and evidence of shear structures in lodgement till abounds and will not be referenced. Thus, subglacial unlithified beds do deform. What of sliding? Our interest in sliding is in the effect this has on subglacial hydraulic conditions, particularly pore-water pressure. As sliding speeds increase, bed separation and bed hydraulic connectivity increases. We define three basal hydraulic regimes.

Type I bed: basal melt-water drainage is along the bed

This is normally the model assumed for valley glaciers with bedrock beds and sparse entrained debris, sliding at moderate to high speeds. Where the ice is in contact with bedrock, basal melt water drains locally through a water film (Hallet, 1979; Vivian, 1980). A non-arborescent but connected system of N-channels and pockets of separated flow provides connections between major N-channels and R-channels. (Walder and Hallet (1979) and Hallet and Anderson (1980) examined the bedrock in front of two retreating glaciers and showed that regions of separated flow plus N-channels comprised 20–30% of the bed.) Major cavities may form where the bedrock slope is large and sliding velocities are the order of 1 m d⁻¹ (Vivian, 1980). Sheet flow (Weertman, 1969, 1972) may also be present but is not considered likely by this author.

The above outlines a system of drainage along the bed. Drainage through the substrate occurs but is negligible, provided the substrate permeability is sufficiently low. Thus, the realization of a type 1 bed requires that bed resistivity be negligible compared to substrate resistivity. For bedrock beds and moderate to high sliding velocities, unless the bedrock is highly fractured and contains faults, it is probably safe to neglect substrate drainage. Bed hydrology is now a partially decoupled problem. The bed pressure distribution is independent of the ground-water flow. The ground-water pressure field has its boundary conditions (pressure distribution at the bed) determined by bed hydrology. However, as we still see, these boundary conditions may be impossible to determine for a type 1 bed.

Note that, if there is no debris entrained in ice at the bed and in bed contact, the average water pressure at the bed must be P ₀, the overburden pressure. This is because ice is everywhere separated from the substrate by water. (Rock-rock contact at the interface has been ruled out by assumption.) The water film will be very thin where the ice moves over obstacles (and pressures correspondingly large). Because bed pressures are less than P ₀ at water channels and regions of large separation, pressures must be greater than P ₀ in regions of close bed contact.

If there are entrained debris fragments at the bed, these must be on average, pushed into the ice by substrate contact at speed μ_m, the basal melt rate (Shoemaker, in press). The existence of rock-rock contact forces imply a reduction in the average interface water pressure. However, unless there is a favorable combination of very high melt rate, dense debris, or thin ice, this pressure reduction is negligible compared to P ₀.

The ground-water pressure, P _GI, at the ice-rock interface is determined by the pressure of water entering (or leaving) the aquifer at aquifer cracks. This implies that in regions of close ice–rock contact P _GI reflects bed pressure where flux is a local maximum and bed pressure a local minimum, that is, at gaps between points of nominal ice–substrate contact. (The flux at points of nominal contact is negligible and does not directly affect the aquifer pressure.) Therefore, it may be that P _GI < P ₀. However, because P _GI depends upon sliding velocity and other factors, the problem of determining P _GI for a type I bed appears to be intractable. Note that as channel pressures, P _c, increase bed pressures become more uniform because separation distances increase. A situation of sheet flow is approached as P _c ⟶ P ₀ and, presumably, P _GI ⟶ P ₀; the glacier floats.

The significance of the last paragraph will become apparent when we consider a deformable substrate. If we required, incorrectly, that P _GI = P ₀, the effective pressure in the soil would be zero and the substrate would (always) have zero strength, a physical absurdity. For the case to be considered (type III bed, below) we will be able to determine P _GI.

Type II bed: bed and substrate hydrology fully coupled

As sliding velocities and/or bed slopes decrease, the ice comes into closer contact with the bed. Separation distances decrease, except at water channels. Drainage along the bed becomes more difficult. For substrates of moderate to high permeability (an aquifer), there should exist a sliding velocity where basal melt water drains equally along the bed and through the aquifer. It follows that there is a range of sliding velocities where basal drainage is a fully coupled problem involving both bed and ground-water drainage.

Note that ground-water flow is controlled by Darcy's law and the flux is proportional to the aquifer pressure gradient. Therefore, ground-water flux will tend to increase linearly with pressure P _GI However, because bed resistance changes with water pressure (bed separation changes), drainage along the bed can be expected to change more rapidly than aquifer drainage. Thus, a modest increase in bed pressure could change the bed hydrology from type I to type II. The type I hydrology is really a special case of a type II hydrology where we are allowed the luxury of neglecting the contribution of ground-water flow to bed drainage. Both types involve basal sliding.

McCall (1952) investigated a small cirque glacier with sliding velocities of order 0.005 m d⁻¹. The sole was in close bed contact except on the lee side of obstacles where gaps of order 1 cm existed. No mention was made of a water film or of channelized bed flow existing locally. A bed which is hydraulically inactive, except for channels connected to the moulin–crevasse system, is an indication that basal melt water drains through the substrate.

Hodge (1979) drilled 24 holes to the base of South Cascade Glacier, Washington State, U.S.A. Only two achieved any connection with the basal water system. Hodge concluded that possibly as much as 90% of the bed is hydraulically inactive.

We note that agiven glacier may contain type I and type II sub-regions. For example, Blue Glacier, Washington State, U.S.A., has been observed subglacially in a region just above an ice fall where the ice is approximately 45 m thick and the basal ice relatively clean (Kamb and LaChapelle, 1964). Below this region the sole separates from the bed as the bed slope increases. The sliding velocity is approximately 90% of the surface velocity. In another region of the glacier where the ice is about 120 m thick the bed was investigated by bore-hole photography (Engelhardt and others, 1978). The bed consists of a debris layer and apparent sliding velocities are less than half the surface velocity. (In some holes no sliding was observed.) The bed is, in general, hydraulically inactive as evidenced by the long time periods, 10–80 d, required for the pressurized water in the holes to establish contact with the bed channel system. In these two regions the beds are probably types I and II, respectively.

Type III bed: basal melt-water drainage is dominated by ground-water flow

As basal debris content increases and bed slopes decrease, sliding virtually ceases. There are few observations corresponding to this situation. Wold and Østrem (1979) were unable to distinguish sliding from deformation at the sole of Bondhusbreen, Norway. Boulton and Jones (1979) claimed to have measured a sliding velocity of about 0.6 cm d⁻¹ beneath Breidamerkurjekull, Iceland; the upper substrate is till. From the data given, their figure 1, Ibelieve that one cannot conclude that sliding was taking place although cavities were observed (Boulton and others, 1974). Cavities can also be formed by highly localized deformation induced by high rock–rock frictional forces.

We shall assume that type III beds exist and shall be concerned exclusively with this case. Low or zero sliding velocities are most likely to occur near the terminus of valley glaciers or beneath ice sheets. If basal melt water drains exclusively through the substrate, pressure P _GI may be determined. This will be demonstrated later. Note, however, that with the fluid velocity known along an aquifer stream line and the terminal (channel) pressure known, the initial pressure (P _GI) is determined through an integration of Darcy's law.

There is an important distinction between the situation of slow sliding and no sliding. If sliding ceases, regelation ice will penetrate the soil to form an ice-debris layer. The constant-state penetration distance will depend upon bed permeability, overburden pressure, basal melt rate, and pore pressure P _GI. Here, P _GI is the pore pressure at the soil–ice debris interface.) For example, if P _GI = P ₀, the penetration distance is zero. If P_GI = 0, b the penetration distance is maximized and is, by calculation on the order of 1 m for till, with a 100 m ice overburden, and a melt rate of 1 cm a⁻¹. Normally, because of higher melt rates and high P _GI we would expect the ice-debris layer to be only centimetres thick.

The presence of an ice-debris layer does several things. It strengthens the mechanical coupling between the ice and substrate, thereby reducing the possibility of sliding. In order for sliding to occur, either the ice-debris layer must be melted out or fracture ofthe ice or ice-debris layer must occur in order to create a velocity discontinuity. The latter possibility would be unusual because ice or ice debris should easily withstand a nominal shear stress the order of 1 bar. (Note that the ice pressure decreases rapidly through the ice-debris layer. Therefore, debris inter-particle forces increase towards the bottom of the ice-debris layer. As a consequence, the ice-debris layer will normally be weakest at the top.)

The ice-debris layer also serves to reduce the hydraulic coupling between the bed and substrate. The bed is effectively sealed off from the pore water. Basal melt water is generated at the base of the ice-debris layer, not at the bed. Finally, the presence of an ice-debris layer rules out drainage by water film or sheet flow. These do not exist since there is no sliding interface.

Subglacial Water Channels

Consider a type III bed where the substrate is an aquifer resting upon an aquitard. A system of subglacial water channels is assumed to exist. We shall consider the likely channel geometry and the mechanics of channel formation. In the case of temperate glaciers resting on bedrock, the widely held view is that the final drainage of surface and basal melt water is by subglacial channels, either by R-channels cut upwards into the ice (Rothlisberger, 1968, 1972) or N-channels incised into bedrock (Nye, 1973). A model for the creation of an R- or N-channel has never been presented. Lliboutry (1983) believed that fast sliding and ice–bedrock separation are a necessity. If, indeed, fast sliding is a prequisite, this presents a problem for a type IIIbed. How can a channel form if the ice is in close bed contact and, in the case where sliding vanishes, there is an ice-debris layer?

A model for the creation of a subglacial channel is suggested by Figure 1 for the case of an advancing temperate ice sheet resting upon soil. (The argument is not directly applicable to a sub-polar ice sheet with a cold terminus.) AA′ and BB′ are established channels, distance d apart, and the ice sheet is advancing to the right in radial flow. It is reasonable to assume that as the ice sheet advances the established channels will simply be extended in the general direction of A′A″ and B′″, most likely along an already existing surface stream bed if this is not being obliterated by a moraine. In any case the erosive effect of the water is certain to provide a channel continuation along the path of least resistance.

Fig. 1. (a) Plan view of the terminus of an ice sheet expanding in radial flow. (b) Simplified aquifer hydrology at the terminus midway between existing water channels.

Now consider the process involved in the formation of a new channel C′ C”. Near the terminus the basal melt water drains through the aquifer across the terminus, not through the channels. Thus, melt water from region A′ C′ B′ CA′ drains across A′ B′. Assuming symmetry about CC′ and uniform basal melting, it follows from Darcy’s law that the maximum hydraulic gradient in region it A′ C′ B′ CA′ occurs at C′.

In Figure 1b we consider the flow of water melted from CC′. In estimating the hydraulic gradient at C′, it is appropriate to consider uniform flow in a channel CC′ of depth t of the aquifer. The average water velocity across C′ D is then

(2)

where m is the uniform basal melt rate and d/2 the estimated length of CC′. The hydraulic gradient is

(3)

where K is the hydraulic conductivity of the permeable material. The critical hydraulic gradient is about i _c = 1.2 (Terzaghi and Peck, 1948)

For the radial flow situation of Figure 1, d increases as the ice sheet expands. A new channel will tend to form when i = i _c or when d increases to a critical distance d _c found from Equations (2) and (3) as

(4)

The consequence of reaching or exceeding i _c is that soil boiling or heaving will occur along C′D due to upward flow of drainage water. This is accompanied by soil dilatation. Heaving initiates a sequence of events. (i) Debris erosion begins along C′D. This is frequently initiated by a spring. (ii) As erosion proceeds, the neighbourhood of C′D becomes a more efficient water collector which increases the erosion rate. (iii) Backward sub-surface erosion (scour) or piping failure (channel erosion) takes place along EF probably at the base of the debris layer since bedding planes are favored (Terzaghi and Peck, 1948). (iv) Once a backward channel is created, forward surface erosion will proceed along C′C″ as it does along A′ A″ and B′ B″.

It is important to note that soil dilatation increases porosity and hydraulic conductivity. Thus, the range of hydraulic conductivity for till, as measured in situ in a de-glaciated environment, is 10⁻¹² to 10⁻⁶ m s⁻¹ (Freeze and Cherry, 1979). This encompasses the range from hardpan to dilatant till. Boulton and others (1974) measured K values in the range 1.1 × 10⁻⁶ to 2.25 × 10⁻⁶ m s⁻¹ in saturated till near the terminus of Breidamerkurjokull, Iceland. They concluded that this till had dilated. (Very high till permeabili-ties have also been measured by G.K.C. Clarke (personal communication) near the terminus of Trapridge Glacier, Canada.) Thus, there is evidence that near-terminus upwell-ing does produce dilatation which is the first stage in piping failure and (by our model) channel formation.

In applying Equation (4) we must keep in mind that, not only is there a broad range of K values for any material, but the range of K values which applies beneath a glacier are unknown. In addition, K depends upon effective pressure which varies with position. We will calculate d _c values for various aquifer materials, taking m = 1 cm a⁻¹, t = 5 m, and i _c = 1.2. We choose K 10⁻⁷ m s ⁻¹ (silty clay), K = 10⁻⁶ m s⁻¹ (dilatant till and silty sand), K = 10⁻⁴ m s⁻¹ (medium sand), and K = 10⁻⁴ m s⁻¹ (medium gravel). The corresponding values of d _c from Equation (4) are: 4 × 10³ m (silty clay), 4 × 10⁴ m (dilatant till and silty sand), 4 × 10⁶ m (medium sand), and 4 × 10⁸ m (gravel). These d _c values must be reduced to reflect the fact that during any melting some of the surface melt water will reach the bed through crevasses near the terminus. This could easily increase m (reduce d _c) by two orders of magnitude; a reduction by a factor of 100 renders reasonably looking d _c values for silty sand and medium sand. Note that similar calculation made for non-dilatant till, using K values in the range 10⁻¹² to 10⁻⁸ m s⁻¹, would give a range of d _c values of 0.04 to 400 m, before any allowance for drainage of surface melt water. This results in unacceptably small channel spacing and suggests that a dilatant till bed may be normal near the terminus.

The present criterion for determining channel spacing results in a very large spread of values over the various soil types because d _c ~ K. For this reason, our later choice of d values will be based upon another criterion.

Water channels, particularly connecting channels, might also form wherever surface melt water finds a new path to the bed through a moulin-crevasse. This should result in localized melting of the ice-debris layer and eventual connection to an established channel.

A cartoon cross-section of a water channel eroded into the top of the debris layer is shown in Figure 2. Will the channel remain at the top; what of downward erosion? Certainly, if the debris layer is sand or silt, rapid downward erosion will occur. In the case of till, after the fines have been washed out a stable gravel bed should remain. The possibility of a channel existing at the bottom or interior of the debris layer must, however, be ruled out. Such channels would be highly prone to collapse and to piping failure. Note that water flowing downward through debris into a bottom channel is a much more critical condition than water flowing upward into a top channel. Also note that the intruding ice tongue is a stabilizing factor in the case of a top channel. With the exception that the channels are cut into debris instead of bedrock and that basal melt water flows towards the channels through the debris layer rather than along the bed, the picture is similar to a network of N-channels. However, there is an important difference. In the case of N-channels, if the channel is not in the direction of the ice flow there is a component of sliding velocity normal to the channel which serves to diminish the depth of ice penetration into the channel. In the present case, sliding may be non-existent, except possibly in a surge stage where the ice-debris layer is absent, so that the channels behave as R-channels as far as ice penetration is concerned. If the debris layer deforms the channels move with the debris.

Fig. 2. Cross-section of a water channel at the top or center of a debris layer where the channel bed is stable to erosion because the fines have been washed away.

The Debris-Layer Deformation - Channel-Flow Dependency

We consider initially the situation of constant-state flow in the debris layer–channel system. Because it can take on the order of a month or longer for channels to adjust their size to a change in flow rate, the constant-state assumption has possible application to only two situations: (i) low constant flow rates during the winter, and (ii) high constant flow rates during the ablation season. We are really considering an average flow rate during the ablation season, ignoring, for example, the diurnal fluctuation. The effect of fluctuations in flow rates will be addressed later. In addition, we assume plane flow with a system of identical parallel water channels which are aligned in the direction of the gradient to the surface of the ice sheet.

As shown by Rothlisberger (1972) and Weertman (1972), the turbulent constant-state flow of water in a horizontal subglacial R-channel, whose cross-sectional area is controlled by the creep of intruding ice, satisfies

(5)

Here, x is measured up-stream from the terminus, ρgh is the overburden pressure, P _c the channel water pressure, B ^* a constant dependent upon the heat of fusion of ice and certain constants describing the creep properties of ice. Q is the channel flow rate and dP _w/dx the generalized pressure gradient. Corresponding to an n = 3 Glen creep model, p = 11/24 and q = 1/12 (Weertman, 1972), Lliboutry (1983) has shown that dP _w/dx can be replaced by dP _c/dx because the velocity gradient is negligible. (This approximation may be in error near the terminus.)

Equation (5) expresses the steady-state situation whereby the closure rate of the channel by intruding ice

(6)

is balanced by the melting due to water flow

(7)

Here, F is the heat of fusion of ice modified by assuming that only two-thirds of the heat released goes into melting the ice (Röthlisberger, 1972). The conduit is assumed to be circular with radius a and, because the ice sheet rests on a horizontal bed, the channels are filled with water. (The channels are nowhere at atmospheric pressure which is an important distinction between ice sheets and valley glaciers (see Hooke, 1984).) B is a constant different from B*. The infusion of heat from springs has not been included (Lliboutry, 1983). Equation (6) also neglects the drag on the intruding ice exerted by the sides of the water channels. This effect would obviously depend on the shape of the channels. It is felt that the geometry is not radically different from that of Nye channels and that Nye's (1953) analysis, upon which Equations (5) and (6) are based, is appropriate.

The deformation, or failure, of most soils is adequately modelled by the Coulomb inequality. In the context of the present problem this takes the form

(8)

where P _p is the pore-water pressure and τ_b the basal shear stress. Both c, soil cohesion, and ϕ, Coulomb friction angle, will be taken constant here. Equality (8) models the failure of soil at the bottom of the ice-debris layer. Assuming horizontal pore-water flow, the soil below this level will not fail because of additional strength gained by the soil weight.

Under constant-state conditions, the pseudo-hydrostatic model (neglect of deviatoric stresses) will be assumed to apply (Shoemaker and Morland, 1984), so that the classical equilibrium equation

(9)

holds. It is to be understood that h is smoothed over a distance of order h.

Pore-water pressure P _p is governed by Darcy's law

(10)

where ρ* is the density of water, q is pore-water flow rate, K hydraulic conductivity, and y is measured towards a channel from a point midway between channels. The vertical pressure drop through the aquifer is neglected and the debris layer is assumed to have a uniform permeability. We are assuming that pore water moves perpendicular to the channels. We shall see that this last restriction does not affect the conclusions.

Assuming a constant basal melt rate m, pore-water flow rate q is given by

(11)

Here, m may be modified to include a contribution due to surface melt water if it is assumed to reach the bed uniformly. Most of the surface melt water is believed to drain directly into the channels for a type III bed. Equations (10) and (11)) describe the motion of basal melt water as it moves through the aquifer towards the channels.

Substituting Equation (11) into Equation (10) and integrating gives

(12)

where the boundary condition P _p(d/2) = P _c has been applied. From Equation (12) the pressure drop ΔP _p along the pore-water path is constant and given by

(13)

This result will be used to evaluate the earlier calculation concerned with the formation of new channels. There, we saw that the distance d between channels, for example for silty sand and t = 5 m, should be less than 4 × 10⁴ m, depending upon the amount of surface melt water being drained through the aquifer at the terminus. In Equation (13), taking d = 10³ m, K = 10⁻⁶ m s⁻¹ (silty sand) and m = 1 cm a⁻¹, we obtain ΔP _p = 0.8 bar, a reasonable value. If we take d = 5 × 10³ we obtain ΔP _p = 19 bar which is very high. This indicates that the d values should be perhaps two orders of magnitude less than the values given by Equation (4) using m = 1 cm a⁻¹.

Several additional comments should be made concerning Equation (13). First, we could determine d values such that ΔP _p is a universal constant, say I bar. The resulting values of d are proportional to √k and therefore the range of d values, as K varies from 10⁻⁷ to 10⁻², is much narrower (and appears more reasonable) than the range which results from applying Equation (4). From the latter, d is proportional to K and pore pressures greater than P ₀ are predicted corresponding to small K values. For example, taking t = 5 m and ΔP _p = 1 bar in Equation (13) and the previousK values, the corresponding d values are computed as 360 m (silty clay), 1140 m (silty sand), 11.4 × 10³ m (medium sand), and 11.4 × 10⁴ m (gravel). These values are given in Table I. It is still true that for low permeabilities d values are unreasonably low. For example, taking K = 10⁻⁹ ms⁻¹ for till gives d = 36 m. Based on these results, it is convenient to adopt a uniform ∆P _p criterion for channel spacing in future calculations. Note, however, that we have offered no physical model to explain why channel spacing should be governed by a uniform ΔP _p criterion.

Table I. Standard parameter values for four soils used in calculations

(e values are computed from Equation (16) using ΔP _p = 10⁵ Pa. Values of d are computed from Equation (13). Data from Bjerrun (1967) and Head (1980))

Secondly, from previous arguments concerning a type III bed, most of the surface melt water which drains sub-glacially will do so through the main channel system. A newly opened moulin-crevasse may require temporary local aquifer drainage until a feeder channel is formed. There is another contribution to basal melting due to losses from ground-water flow. This source of losses is normally small. Therefore, under most of the ice sheet, effective basal melt rates should not be much greater than 1 cm a⁻¹. The terminus region may be an exception.

We have assumed that aquifer flux is much higher near the (active) terminus than elsewhere. Thicker ice implies a thin ice-debris layer and more possibility for crevasses to reach the aquifer. In addition, surface melt-water flux is much higher near the terminus. A third factor is that, as we shall see, substrate deformation can be expected to occur at the terminus, particularly where ice with a high frontal relief is advancing. The substrate deformation energy release would add to m and further reduce the ice-debris layer thickness. Thus, it is consistent to prescribe a large effective m value in Equation (2) and a small m value in Equation (13) away from the terminus.

Maximal profiles

We shall now consider the problem of obtaining maximal profiles, h(x), determined by the condition that the top surface of the deformable aquifer is on the verge of failure everywhere, i.e. that Equation (8) is satisfied. We assume that d is determined by the uniform condition ΔP _p = 1 bar (the values of Table 1). (We will show that the conclusions are not affected by the value of ΔP _p, provided it is not large compared with 1 bar.) It is convenient to employ an average P _p at section x. Thus, Equation (12) is replaced by

(12’)

This amounts to the tacit assumption that the bed will either deform uniformly or not deform at all at a section x = constant.

Consider Equations (5), (8), (9), and (12′) in variables h, τ_b, P _p, and P _c. Elimination of P _p and τ_b results in

(14)

(15)

where the approximation p = p* has been made. Define

(16)

Equation (15) is then written as

(15a)

Table I gives typical parameter values for various soils. Since, in general, e < 0, it follows from Equation (15′) that there is no solution in the neighbourhood ofx = 0 satisfying P _c(0) = h(0) = 0. At the terminus the overburden pressure vanishes and with P _p < 0. Equation (8) implies that if P _p tan ϕ ≥ c no shear stress can be sustained. Yet, at the terminus we expect that h′ = 0(1) and Equation (9) insists that the shear stress be positive. The culprit is the equilibrium Equation (9). Near the terminus, longitudinal deviatoric stresses t _xx must be present (Shoemaker and Morland, 1984).

The more general equilibrium equation is

(17)

Near the terminus, Equation (15) must be replaced by

(18)

and Equation (15′) by

(18)

where

(19)

Equation (19) may be used to solve for t_xx near the terminus once h is known.

In order to achieve a simple regular solution at the origin, we prescribe e*(0), h′(0), and t _xx(0) as

(20)

(21)

(22)

(Note that, if e*(0) > 0, it would be possible to have singular solutions satisfying h′(0) = W. Our aim is to consider parameter values which yield low maximal profiles.)

The prescription of e*(x) near the terminus introduces arbitrariness into the solutions but, as we shall see, the qualitative results are not affected provided that Equations (20), (21), and (22) are satisfied. Function e^*(x) is arbitrarily defined as

(23)

Unless otherwise stated, all results correspond to l ^* = 2 km. Note that changing l ^* by a factor of 20 leaves the results unchanged to four places outside the near-terminus region. Thus, the arbitrariness of the solution near the terminus has little effect upon the global solution.

Function Q(x) in Equation (14) has a contribution from basal melt water and, in general, from surface melt water. The function is written as

(24)

The second term in the summer formula provides for a surface-melt contribution which varies linearly with elevation with a cut-off at h(x) = H; Here H = h(L ^*) defines L ^*. The formulation is implicit but the winter H(x) can be used as the initial function in an iteration scheme. The standard values ofH and M were taken as 10³ m, and 5 m a⁻¹, respectively.

The initial conditions for the maximal profile problem are

(25)

Equations (14) and (15′) are solved on 0 < x < L. Results will be given for L = 100 and 1000 km. Using the e values of Table 1 it is clear that t _xx is less than 1 bar; quantitative results for t _xx will not be given.

Figure 3 illustrates various maximal profiles corresponding to the data of Table I for L = 1000 km, ΔP _p = 1 bar, assuming constant-state winter hydrology conditions. Several conclusions follow:

1. All curves lie above 3√x except at the terminus and curve A (silty clay) near the divide. For smaller ice sheets or summer conditions, this conclusion is strengthened. Thus, the occurrence of very low-relief ice sheets cannot be attributed to substrate deformation under constant-state conditions for the case of an aquifer resting upon an aquitard if efficient melt-water drainage is provided by sufficiently small d values.

2. The profiles become very flat well before the divide where h′(L) = 0. This is because P _c/∂gh increases very rapidly and approaches unity (see Röthlisberger, 1972).

3. In moving from silty clay to medium gravel (increasing K) h increases. This is primarily a reflection of increasing ϕ.

Fig. 3. Maximal profiles corresponding to data of Table I for L = 10³ km and ΔP_p = 1 bar. A. Silty clay; B. Silty sand; C. Medium sand: D. Gravel; E.

As e increases or ΔP _p increases, h decreases. Such an increase could result froma decrease in K beyond the values assumed in Table I. Table II illustrates the effect of a ten-times increase in e upon maximal profiles. Channel spacings remain fixed, so Q is unchanged. It is clear that the effect is minor. Therefore, the arbitrariness in fixing ΔP _p and in defining e(x) does not necessarily jeopardize the qualitative results.

Table II. Ratios of: Winter h (Standard e Values)/h (10 × Standard e Values) at Mid-Span and Divide for 100 km and 1000 km Ice Sheets. Channel Spacings Unchanged

Increased channel flux during the ablation season causes a reduction in P _c and a corresponding increase in bed strength and maximal profiles. However, the effect is not large. If the Q(x) functions of Equations (24) are used, the increase in h(x) never exceeds 10% for any of the four standard soils.

Table III illustrates the effect of an increase in d by a factor √10 with e and ΔP _p undergoing a corresponding increase by a factor of 10 (see Equations (13) and (16)). Here Q increases since d increases, so that the weakening effect of an increase in ΔP _p tends to be offset by the strengthening effect of a decrease in P _c. It is clear that the two effects almost cancel.

Table III. Ratios of: Winter h (Standard e Values)/h (10 × e Values) at Mid-Span and Divide for 100 km and 1000 km Ice Sheets. Channel Spacings Satisfy e ~ √d. Silty Clay e Value Remains at Zero

On the basis of Tables Table II and Table III, the general conclusions given previously hold unless ΔP _p values are much higher than 10 bar.

Experience with alpine glaciers with presumed type I beds indicates that there is usually more than one main channel across a valley of width 2 km and less, e.g. Kamb and others (1985). The impression is that the bed–channel system provides very efficient drainage. The response time for dye-tracing experiments on fast-sliding glaciers provides some evidence for this assertion (Kamb and others, 1985). It does not follow that drainage will be this efficient or channel spacings this small for type III beds. One reason is the increased difficulty of channel formation. Thus, until evidence for channel spacing is found, no definite statements can be made about values for ΔP _p and e, or about maximal profiles.

Temporal Effects

We shall first consider the effect of a diurnal fluctuation in flow rate Q during the ablation season. An examination of Equations (14) and (15) reveals that pore water pressures are increased if the closure rate of the channels is increased (or melt rate decreased), because u _c and u _m in Equations (6) and (7) will find equality at a reduced diameter. A smaller-diameter channel causes an increase in P _c and P _p. If Q, and hence P _c, fluctuates periodically about a mean value, the average closure rate will be increased because u _c in Equation (6) is proportional to the time integral of (ρgh − P _c)³. (If this factor was (ρgh − P_c), the average u _c would be unchanged.) Of course, u_m in Equation (7) will also be affected by fluctuation inQ and, as will be shown, the two effects tend to cancel each other.

It is possible to analyze diurnal fluctuations in Q by amending the previous constant-state analysis. Assume that during the ablation season Q(x,t) is periodic in time, 12 h at a maximum flow rate Q _M(x) and 12 h at a minimum flow rate Q _m(x). Let Q _M(x) = kQ _m(x). Coefficient k could be as large as two. Such an upper bound would be consistent with Elliston's (1973) measurements of the diurnal fluctuations in the flow rates of glacial streams. The temporal variation in Q(x) will be confined to a lower elevation 0 < h(x) < H₀. In this lower region the mean flow rate Q = (Q _M + Q _m)/2 will be assumed constant. Thus

(26)

where L^* and M are consistent with Equation (24).

The Darcy–Weisbach formula states that dP _c/dx is proportional to Q ²(X). Thus, we can assume that for x < L ^* the channel pressure fluctuates periodically between P _m(x) and P _M(x) where P _M = k² P _m. (Boulton and Vivian’s (1973) measurements of diurnal pressure fluctuations indicate k ≃ 1.3.) The pressure also fluctuates in the upper channel, x > L _*, but we will not consider this.

We assume that the channel diameter does not change significantly in 12 h and therefore is constant in time. From Equation (6), the closure rate is proportional to (ρgh – P _m)3. The average channel-closure rate is therefore given by ½aB[(ρgh − P _M)3 + (ρgh − P _M)3]. It is convenient to express P_m and P _M m terms of average pressure Thus Equation (6) is replaced by

(27)

where

From Equation (7) the melt rate is proportional to QdP _c/dx which, by the previous assumptions, is proportional to

(N1))

since and Equation (7) is replaced by

(28)

An equation analogous to Equation (5) is now obtained by setting u _c = u _m and eliminating a by using the equation

(29)

where p ^* = ½ and q ^* = 2/3 (Weertman, 1972). Here R is a roughness factor which is taken as 0.58 m^11/6 N^−½ s⁻¹, corresponding to a concrete surface. the result is

(30)

where

(31)

It is crucial to recognize that for most soils the diurnal channel-pressure fluctuation will affect only a small region of the aquifer near the channel. Table IV gives the estimated distance at which a diurnal sinusoidal pressure wave will attenuate by a factor of 1/e. Gravel is a possible exception to the above conclusion.

Table IV also gives distances at which annual pressure fluctuations will attenuate by a factor of 1/e. Even here, with the exception of gravel, the majority of the aquifer will be unaffected by annual pressure variations. Therefore, increased channel pressures which probably occur during the onset of the ablation season will, at most, cause local aquifer deformation in the vicinity of channels. Overall aquifer failure can be adequately modeled, as in Equations (29) and (30), by prescribing a time-averaged and in a constant-state analysis. We have already seen that there is little difference between summer and winter constant-state solution.(because Q is raised to the one-twelfth power in Equation (14)). We caution that these conclusions are based upon the channel spacings given in Table I and that the data, particularly for hydraulic conductivity, can have wide ranges

Table IV. Hydraulic Conductivities K, Aquifer Compressibilities α, Porosities n, and Computed Dist Ances x^* AND x^** at which the Amplitude of Daily and Annual Channel-Pressure Fluctuations Attenuate By 1/e. (Data From Freeze and Cherry, 1979)

Figure 4 graphs the ratio where is computed from Equation (30) and P _c(x) from Equation (14) (or Equation (30) with k = 1) using h(x) = 3√x Curves. A and B illustrate that diurnal mean pressures are much lower near the terminus and moderately higher elsewhere as compared to winter constant-state pressures. Curves c and D illustrate that this effect is smaller when the same channel flux is used to compute and P _c. Figure 4 assumes that channel spacings correspond to a medium sand aquifer. The results are qualitatively similar for other channel spacings. In general, the critical substrate-pressure conditions are encountered during the winter near the terminus and during the summer near the divide. However, if the diurnal flux variation occurs in only a short channel length L ^* ≪ L, the winter pressure may exceed the mean diurnal pressure almost everywhere.

Fig. 4. Ratio of (Equation (30)) to P_c(x) (Equation (14)). Curves C and D use from Equation (26). Curves A and B use from Equation (24) (winter). κ values are 1.2, curves A and C, and 2.0, curves B and D.

Maximal profiles corresponding to the diurnal pressure condition can also be obtained. Equation (15′) is retained in form with and h → ĥ The parameters are unchanged as are the e ^* functions. The results are ambiguous. In the case of gravel and medium sand, the constant-state maximal profiles are so elevated that any reasonable H value, say H = 1 km, implies that only a very short channel section will be affected by diurnal pressure fluctuations; there is a negligible effect upon pore pressures. For the case of silty clay, there is no maximal profile solution for k > 1.2. This does not mean that maximal profiles do not exist but only that they do not exist within our simplified model. The case of silty sand and k = 2 produces a 2 produces a 9.5% decrease in the maximal profile at L ^* = 4.4 km.

Finally, we consider the transient effect of an increase in melt rate during the transition from winter to the ablation season. Because the channels cannot react instantaneously by adjusting their diameters to accommodate an increase in melt water, the moulin–revasse system will initially store water. This produces an increase in P _c and Q. Thus, u _m(u _c) increases (decreases) from Equations (7) and (6), respectively. We are concerned with estimating the time interval ∆t required for the system to again achieve steady state under an increased Q(x). A rigorous solution of the problem will not be attempted here.

During the Δt period of channel enlargement P _c(x) is evevated. We will assume that

(32)

which is the maximum realistic value that can be assigned. This choice for P _c implies, from Equation (6), that u _c = 0. This affords a major simplification to the problem of determining Δt but it must be recognized that Δt values based upon Equation (29) will underestimate the true values.

The flux Q in a circular pipe is expressed by the Manning formula for turbulent flow. Thus

(33)

where R is a roughness factor. By our assumption, dP _c/dx = ρgh′ during the transition period. Near the terminus Q increases by a factor of perhaps 500 between the two constant-state conditions corresponding to winter and summer. From Equation (33), radius a therefore increases by a factor of about ten. We will approximate the time Δt required for channel enlargement by a factor of ten.

From Equation (7), u _m = da/dt with dP _c/dx = ρgh′. Substituting for Q from Equation (33) gives

(34)

which upon integration gives

(35)

The initial channel radius is calculated from the winter flow rates using Equation (33). After substituting, this results in

(36)

Here Q _i is the winter flux from Equation (24).

Table V illustrates Δt values for the four prototype winter maximal profiles as well as the more realistic 3√x profile, all for L = 100 km. It is assumed that channel enlargement occurs in the lower elevation h(x) < H = 2 km.The corresponding L ^* values are noted.

It is clear from Equation (36) that, with dP _c/dx = ρgh′, the critical variable in determining At is the surface slope. As h′ increases. At decreases rapidly. Thus, near the terminus of any active glacier or ice sheet, channels can even adjust to diurnal pressure fluctuations

Channels in a given ice sheet respond on at least three different time-scales. Adjoining the near-terminus region, the response is rapid enough that channel diameters fully adjust to seasonal fluctuations in Q. In a central section the response is so slow that channels never adjust fully to constant-state winter and summer conditions. (Note, however, that if this central region is at an elevation above H there is no seasonal flux variation.) In a third region adjoining the divide, the response time is so long that, even if there are seasonal flux variations, channel diameters may be assumed to be constant.

Table V. Δt Values in Days at Various x/L^* Values for Constant-State Winter Maximal Profiles and 3√x. L^* Values are: Silty Clay, 100 km; Silty Sand, 38.3 km; Medium Sand, 10.9 km; Medium Gravel, 2.6 km; 3√x Profile, 100 km. Channel Spacing is 400 m for last case

Conclusions

• We have shown that it is possible for a system of subglacial water channels to evacuate basal melt water from beneath an ice sheet with such efficiency that for the case of a substrate consisting of an aquifer resting upon an aquitard, if the aquifer has moderate to high permeability (K > 10⁻⁶ ms⁻¹) the aquifer will, in general, not deform beneath ice sheets of normal profiles such as h(x) = 3√x. There are exceptions to this rule. First, substrate deformation can be expected near the terminus where the terminus slope exceeds tan ϕ. Secondly, local deformation near channels may occur during periods of high channel pressures approaching ρgh. This local deformation may become global in the case of highly permeable aquifers such as gravel.

• The analysis neglects longitudinal drainage through the aquifer. The neglected term only adds to the drainage efficiency. It becomes important for very large channel spacings.

• The critical unknown has been shown to be channel spacing. If channel spacing becomes so large that ΔP _p values exceed 10 bar, the channel system becomes relatively inefficient. High pore pressures result and substrate deformation becomes much more likely.

• Aquifers of low hydraulic conductivity (K < 10⁻⁷ms⁻¹), such as silty clay and till, require very narrow channel spacings in order to achieve efficient meltwater evaluation.

• It is likely that soil dilatation adjacent to the terminus is a natural consequence of the critical hydraulic gradient being exceeded by the upwelling of melt water being evacuated across the terminus. It is dangerous to infer from this that aquifer dilatation is a general condition beneath a glacier or ice sheet.

• The aquifer–channel hydraulic system, for the case of moderate to low hydraulic conductivity (K < 10⁻³ m s⁻¹), appears to experience its critical condition during the ablation season. This condition may be modelled by extending the diurnal model to elevations h(x) > H.