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⁻¹, 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⁴ (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

(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 (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.

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

(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 (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

(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 (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⁻¹² to 10⁻⁶ m s⁻¹ (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⁻⁶ 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 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

(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 (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

(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 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

(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 Reference Shoemaker and Morland(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 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⁻¹. 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 (Reference Shoemaker and MorlandShoemaker 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 Reference RöthlisbergerRö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 (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)³. (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 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) < 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. (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 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 (Reference WeertmanWeertman, 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.