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 
Introduction
Observations by Reference Boulton and JonesBoulton 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 (Reference Boulton, Boulton, Dent and MorrisBoulton 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−1, 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, ∼105; gravel upon medium sand, ∼103; medium sand upon clay, ∼105; gravel upon fractured igneous and metamorphic rock, ∼104 (Reference Freeze and CherryFreeze 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
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 (Reference HollinHollin, 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 (Reference OrowanOrowan, 1949; Reference NyeNye, 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½ (Reference MathewsMathews, 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 (Reference PatersonPaterson, 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. Reference SlaterSlater, 1926; Reference KupschKupsch, 1962; Reference WoldstedtWoldstedt, 1965; Reference LavrushinLavrushin, 1971. The near-terminus regions of ice sheets located far inland would likely have had cold basal conditions (Reference WrightWright, 1973; Reference Moran, Moran, Clayton, Hooke, Fenton and AndriashekMoran 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 (Reference PorterPorter, 1977).
Sliding or Deformation?
There is observational evidence of ongoing shearing in subglacial deformable beds; Reference Boulton and JonesBoulton 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; Reference VivianVivian, 1980). A non-arborescent but connected system of N-channels and pockets of separated flow provides connections between major N-channels and R-channels. (Reference VivianWalder and Reference HalletHallet (1979) and Reference Hallet and AndersonHallet 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−1 (Reference VivianVivian, 1980). Sheet flow (Weertman, 1969, Reference Weertman1972) 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 0, 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 0 at water channels and regions of large separation, pressures must be greater than P 0 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 0.
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 0. 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 0 and, presumably, P GI ⟶ P 0; 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 0, 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.
Reference McCallMcCall (1952) investigated a small cirque glacier with sliding velocities of order 0.005 m d−1. 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.
Reference HodgeHodge (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 (Reference Kamb and LaChapelleKamb 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 (Reference Engelhardt, Engelhardt, Harrison and KambEngelhardt 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. Reference Wold and ØstremWold and Østrem (1979) were unable to distinguish sliding from deformation at the sole of Bondhusbreen, Norway. Reference Boulton and JonesBoulton and Jones (1979) claimed to have measured a sliding velocity of about 0.6 cm d−1 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 (Reference Boulton, Boulton, Dent and MorrisBoulton 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 0, the penetration distance is zero. If PGI = 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−1. 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 (Reference RöthlisbergerRothlisberger, 1968, Reference Röthlisberger1972) or N-channels incised into bedrock (Reference NyeNye, 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
where m is the uniform basal melt rate and d/2 the estimated length of CC′. The hydraulic gradient is
where K is the hydraulic conductivity of the permeable material. The critical hydraulic gradient is about i c = 1.2 (Reference Terzaghi and PeckTerzaghi 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
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 (Reference Terzaghi and PeckTerzaghi 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−12 to 10−6 m s−1 (Reference Freeze and CherryFreeze and Cherry, 1979). This encompasses the range from hardpan to dilatant till. Reference Boulton, Boulton, Dent and MorrisBoulton and others (1974) measured K values in the range 1.1 × 10−6 to 2.25 × 10−6 m s−1 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−1, t = 5 m, and i c = 1.2. We choose K 10−7 m s −1 (silty clay), K = 10−6 m s−1 (dilatant till and silty sand), K = 10−4 m s−1 (medium sand), and K = 10−4 m s−1 (medium gravel). The corresponding values of d c from Equation (4) are: 4 × 103 m (silty clay), 4 × 104 m (dilatant till and silty sand), 4 × 106 m (medium sand), and 4 × 108 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−12 to 10−8 m s−1, 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 Reference WeertmanWeertman (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
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 (Reference WeertmanWeertman, 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
is balanced by the melting due to water flow
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 Reference HookeHooke, 1984).) B is a constant different from B*. The infusion of heat from springs has not been included (Reference LliboutryLliboutry, 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 Reference NyeNye'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
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 Reference Shoemaker and Morland(Shoemaker and Morland, 1984), so that the classical equilibrium equation
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
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
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
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
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 × 104 m, depending upon the amount of surface melt water being drained through the aquifer at the terminus. In Equation (13), taking d = 103 m, K = 10−6 m s−1 (silty sand) and m = 1 cm a−1, we obtain ΔP p = 0.8 bar, a reasonable value. If we take d = 5 × 103 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−1.
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−7 to 10−2, 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 0 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 × 103 m (medium sand), and 11.4 × 104 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−9 ms−1 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 = 105 Pa. Values of d are computed from Equation (13). Data from Reference BjerrunBjerrun (1967) and Reference HeadHead (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−1. 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
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
where the approximation p = p* has been made. Define
Equation (15) is then written as
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 (Reference Shoemaker and MorlandShoemaker and Morland, 1984).
The more general equilibrium equation is
Near the terminus, Equation (15) must be replaced by
and Equation (15′) by
where
Equation (19) may be used to solve for txx 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
(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
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
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 103 m, and 5 m a−1, respectively.
The initial conditions for the maximal profile problem are
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:
-
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.
-
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 Reference RöthlisbergerRöthlisberger, 1972).
-
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 = 103 km and ΔPp = 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 (Reference Kamb and KambKamb 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)3. (If this factor was (ρgh − Pc), the average u c would be unchanged.) Of course, um 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 Reference EllistonElliston'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) < H0. In this lower region the mean flow rate Q = (Q M + Q m)/2 will be assumed constant. Thus
where L* and M are consistent with Equation (24).
The Darcy–Weisbach formula states that dP c/dx is proportional to Q 2(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 = k2 P m. (Reference Boulton and VivianBoulton 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 Pm and P
M m terms of average pressure 
Thus Equation (6) is replaced by
where 
From Equation (7) the melt rate is proportional to QdP c/dx which, by the previous assumptions, is proportional to
since 
and 
Equation (7) is replaced by
An equation analogous to Equation (5) is now obtained by setting u c = u m and eliminating a by using the equation
where p * = ½ and q * = 2/3 (Reference WeertmanWeertman, 1972). Here R is a roughness factor which is taken as 0.58 m11/6 N−½ s−1, corresponding to a concrete surface. the result is
where
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 Reference Freeze and CherryFreeze 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 Pc(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
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
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
which upon integration gives
The initial channel radius is calculated from the winter flow rates using Equation (33). After substituting, this results in
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−6 ms−1) 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−7ms−1), 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−3 m s−1), 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.
Acknowledgements
The research was partially supported by NSERC of Canada. The numerical calculations were skillfully performed by H. Leung. The author thanks C. MacDonald for the use of his library card.
