Ice-stream surface texture, sticky spots, waves and breathers: the coupled flow of ice, till and water

Abstract

This paper addresses the coupling of flows of ice, till and water, and the issue of whether such coupling provides mechanisms for meso-scale (kilometres to tens of kilometres) variability in ice-sheet flow and texture. The question of whether effective pressures at the ice-bed interface are statically or hydraulically controlled is examined in this paper. The answer is scale dependent, and has a significant effect on the relationship between ice surface and basal topography.

The consequences of these considerations on till flow, coupled ice–till flow and coupled ice, till and water flow are examined. An analysis of till-flow kinematics and shock formation is carried out. The linear stability of coupled long-wavelength ice-till flow is analysed, and regions in parameter space where this flow is unstable, with rather small rate constants are found. Upstream-moving ice surface waves are predicted. The linear stablity of coupled ice–till–water flow is examined, where water flow is modelled using a basal flow system with effective-pressure-dependent properties. Again, regions in parameter space where the system is linearly unstable are found, this time with relatively rapid rate constants. The water pressure exhibits “breather” modes.

These analyses assume that there is a substantial basal traction. A problem with models of ice streams wholly restrained at the side is identified: they seem to predict erosion rates which are unfeasibly large.

There appears to be sufficient variability in the ice–till– water system to potentially explain texture in ice-stream surfaces, variations in ice-stream thickness of tens of metres not directly relatable to topography, and waves moving upstream or downstream. Most importantly, the ice-stream–bed system is shown to exhibit meso-scale variability simply by coupling ice flow according to the shallow-ice approximation, till flow according to the hydrostatic thin-till approximation and water flow according to an effective-pressure-dependent hydraulics.

Notation

1. Introduction

Attention is being increasingly focused upon meso-scale variability in ice sheets and ice streams; the switching of ice-streams (Rose, 1979; Relzlaff and Bentley; 1993), thermal controls on ice-sheet and ice-stream behaviour (MacAyeal, 1992b, 1993; Payne, 1995; Greve and MacAyeal, 1996; Payne and Dongelmans, 1997; Marshall and Clarke, 1997), medium-scale variability in ice-stream flow (Bindschadler, 1997), ice-stream texture (Bindschadler and Vornberger, 1990; Bindschadler and others, 1996), variations in ice-stream velocity and basal traction (Whillans and others, 1989; Alley, 1993; Whillans and Van der Veen, 1993; MacAyeal and others, 1995), the role of ice-stream margins (Echclmeyer and others, 1994; Van der Veen and Whillans, 1996; Jackson and Kamb, 1997), migration of inter-stream divides (Nereson and others, 1998), migration of smaller scale wave-like forms (Hulbe and Whillans, 1997) and putative instability mechanisms (Kamb, 1991; Fowler and Johnson, 1995, 1996; Anandakrishnan and Alley, 1997a).

It seems that ice-stream texture is something of a puzzle; for example, there is not always a clear relationship between surface and basal topography (Whillans and Van der Veen, 1993; Bindschadler and others, 1996). Figure 1 shows an image of the Siple Coast ice streams, and shows in particular the small-to-medium scale (tens of kilometres) topography that gives the fast-flowing ice streams surface texture. If this texture is not due directly to form drag induced by flow around obstacles, as Bindschadler and others (1996) state, then it can only be due to variations in friction at the base (Alley’s (1993) sticky spots) or variations in the properties of the ice (e.g. Whillans and Van der Veen, 1993; Hulbe and Whillans, 1997). In either case, we expect dynamism in these forms, and there is some suggestion that this is true; Hulbe and Whillans (1997) in particular identify surface forms moving upstream, explaining this in terms of variations within the ice fabric.

The implication of this is that ice-stream texture arises from the same mechanisms which are believed to control the large-scale flow of ice streams. This paper concentrates on basal mechanisms which give rise to surface topography not strongly related to basal topography and basal mechanisms which might give rise to spatial and temporal variability in ice-stream flow. It shows that simple physical models of ice, till and water flow give rise to a wide variety of possible behaviour; this is a desirable feature of models when the message from observations at the moment is that rather different types of behaviour are possible. A recent related study, by Gudmundsson and others (1998), has found a strong dependence of surface topography on basal-friction variations. There is no intention to argue that variations in ice rheology, due in particular crystal fabric, cannot produce the same sorts of effects.

A primary concern of this paper is the effect effective-pressure distributions at the ice–till interface have on ice dynamics. In particular, it is concerned with the coupling between ice-sheet dynamics and long-wavelength (i.e. longer than the ice is thick) variations in the interfacial effective-pressure distribution arising from both static and hydraulic effects. This is a complex and controversial subject, and has implications for, and ramifications on, variability in ice-stream dynamics, ice-surface texture, water routings beneath an ice sheet and subglacial-sediment deformation. At its heart is the interplay between static pressure gradients in ice, water and sediment and the potential gradients induced by the need to discharge subglacial melt, and the effect this has on the coupled flow of ice, till and water.

A significant difference in the hydraulic theories of glaciers overlying porous media and those overlying impervious media is that, in the latter case, water is generally believed to drain through a thin film or distributed system which is more or less at the local ice pressure (Paterson, 1994; Hooke. 1998). Thus, the effective pressure is, on this level of theorizing, independent of elevation. The implication of this is that hydraulic potential gradients are important compared with static pressure gradients over short length-scales. In contrast, for glaciers overlying porous media, hydraulic gradients are small compared with static pressure gradients over sufficiently short horizontal length-scales. Thus, for example, as one moves up along the surface of a porous subglacial hummock, the interfacial effective pressure increases owing to the difference in density between ice and water.

Fig. 1. AVHRR image of the Siple Coast. Upstream areas of ice streams marked with letters. Flow is separated into slow domes, with smooth texture and dear ridges, fast-flowing ice streams with lots of texture, and the smoother, stagnant Ice Stream C. After Bindschadler and Vornberger (1990) to whom thanks are expressed for providing this image. The vertical dimension is about 500 km.

A message of this paper is that contrasting the theories in such a manner rather oversimplifies matters. Over short wavelengths (how short is discussed extensively here), interracial effective pressures are significantly affected or controlled by hydrostatic and lithostatic effects. We investigate length-scales over which drainage by ground water alone is possible, and consider the effects of basal and ice topography on interfacial effective pressure and drainage routings. If all water could be drained by ground-water flow, the matter would be simple; but it cannot in general, and there is in consequence a need, especially for ice sheets advancing over the continental shelf, for theory to explain drainage between ice and till, for example that provided by Walder and Fowler (1994). A more general theory from a dynamical point of view, which includes the Walder–Fowler theory (but does not make the same specific predictions about sub-Glaciol sedimentology associated with drainage) and is also related to that presented by Alley (1996), is used in this paper to examine the interactions between drainage, ice cover and static effective-pressure gradients.

While theories of subglacial water flow which assume that the water pressure is of the ice-overburden form are well known, the lesser used theory of static effective-pressure gradients has been exploited (Lliboutry, 1983; Boulton and Hindmarsh, 1987, 1996; Hart and others, 1990; Boulton 1996a, b; Hooke, 1998) to construct theories which try to explain a variety of Glaciol-geological features; certain properties of eskers, tunnel valleys, till thicknesses and observations of depth of tectonization and drumlinization. It is fair to say that in each of the cases the static effective-pressure-gradient theories have not been ruled out as significant mechanisms.

An apparent, implication is that at increased elevation subglacial relief draped by or composed of porous sediments will be more resistant to glacier flow owing to its increased strength. This effect is quite significant; where hydraulic potential gradients are small, a 100 m hillock will produce a 80 kPa difference in effective pressure between the top and the bottom. Since such effective-pressures differences are comparable with typical basal shear stresses, and it is widely believed that for a bed to be deforming, effective pressures must be comparable with the applied shear stress, there is a strong suggestion that even relatively gentle relief can significantly affect ice-stream dynamics and may be one cause of the sticky spots (areas of slower ice velocity) discussed by Alley (1993). The basal relief of Ice Streams D and E is of the order of hundreds of metres (Bindschadler and others, 1996), as is that of the Rutford Ice Stream (Smith, 1997). However, if interfacial pressures are statically controlled, a increase in ice thickness of 10 m can also increase effective pressures by 10⁵ Pa and could amplify or damp effects induced by variations in basal topography. The relationship between ice and surface topography where basal flow is effective-pressure dependent is investigated in this paper.

The theory of interfacial pressure distributions is reviewed in section 2. Here it is shown that over relatively short length-scales, interfacial effective pressures are statically controlled, while over long length-scales, interfacial effective pressures are hydraulically controlled. What is meant by short depends on the hydrogeological properties of the aquifer beneath the glacier and the recharge rate. Above sedimentary basins, static gradients are expected to occur over length-scales of hundreds of metres to kilometres. This analysis is extended by the introduction of the so-called “linked-cavity” type of drainage system (Walder, 1986; Kamb, 1987; Paterson, 1994; Walder and Fowler, 1994; Hooke, 1998), where storage and transmissibility decrease with effective pressure; this is not a claim that linked cavities are the main drainage routing but is a claim that R-channels are not the preferred water routing under ice streams (e.g Alley, 1989b, Walder and Fowler, 1994). There have been some important recent borehole observations of subglacial hydraulic system by Murray and Clarke (1995) of Trapridge Glacier and Engelhardt and Kamb (1997) of Ice Stream B. They describe the small-scale variability of basal water-pressure systems and Engelhard and Kamb contrast these with the much less variable large-scale flows. Hindmarsh (1997a) and Engelhard and Kamb (1997) discuss the causes and consequences of this scale-dependence of variability; however, this remains an important unsolved problem.

In between the two poles of statically controlled and constant interfacial effective pressures are cases where both hydraulic and gravitational factors are significant. The interaction between bed and ice topography, effective pressure and drainage routings are examined in a plane-flow model (section 2.5). This shows that coupling between static and hydraulic gradients reduces, but does not eliminate, interfacial effective-pressure gradients over short length-scales.

The above analyses do not consider the additional coupling induced by the flow of ice. This is introduced by assuming that the ice motion is entirely due to till deformation. The ideas of till deformation in this paper rest upon the notion of a viscous rheology for till, which is becoming increasingly controversial; laboratory observations and point measurements under glaciers indicate plastic-type behaviour (Kamb, 1991; Iverson and others, 1995; Murray and Clarke, 1995). Hindmarsh (1997a) has reviewed this controversy, arguing that while on the small scale till behaves plastically on the larger scale the net effect of failure events is that of viscous behaviour. One way of testing that till is viscous on the large scale is by looking at tine large-scale consequences of the idea. To this end, till kinematics and shock formation (jumps in till thickness) are considered (section 3). This uses a model of subglacial till deformation, which has a constant shear stress with depth and where effective pressures are controlled by static gradients. We shall call this an HTTA theory, for hydrostatic-thin-till approximation. It was developed in the late 1980s by several glaciologists (Boulton and Hindmarsh, 1987; Clarke 1987; Alley, 1989a).

The coupling of ice and till flow and topography are examined in section 4. The consequences of these effects on sticky spots (Alley, 1993; Whillans and Van der Veen, 1993) are considered, as well as the possibility that zones of drumlin formation may have dynamics; that is, subglacial conditions may change so as to cause zones of drumlin formation to move. It is found that small variations in till topography can induce much larger variations in the ice-surface topography. In this section, the problem of how till and ice topography change through time is not considered.

Some aspects of the dynamic problem are considered by conducting a linear stability analysis (section 5.1). Hindmarsh (1998b) has found that the HTTA model coupled with Newtonian ice allows unstable growth of till relief for short wavelengths (i.e. less than the thickness of ice). In this paper, we carry out a corresponding linearized analysis for long wavelengths. Regions of parameter space where till flow is linearly unstable are found. Shock formation is suggested as a mechanism for non-linear quenching of these instabilities. It is known that till kinematic waves can move backwards (Hindmarsh, 1996) and it may be that observations of ice features moving upstream (Hulbe and Whillans, 1997) are a consequence of this, a possibility that is considered here.

This model is extended to include water pressure as a prognostic quantity (section 5.31), its evolution being determined by a non-linear diffusion equation, analogous to the Darcy equation but with effective-pressure-dependent parameters which represent a distributed system (Alley, 1996). Again, a highly unstable system is found, with unstable modes with very high growth-rate constants. “Breather” modes in the water pressure occur. It is suggested that these instabilities are related to observations of meso-scale variability in ice streams (Rose, 1979; Retzlaff and Bentley, 1993; Bindschadler, 1997). These analyses differ from those of Kamb (1991) because they are not coupled to the heat equation.

These linear-stability analyses are based on a model of ice-stream mechanics which ignores lateral shear. There have been some fairly strong assertions in the literature (e.g. Whillans and Van der Veen, 1997) that basal resistance in some ice streams is negligible. This conclusion is obtained on the basis of force-balance measurements and assumptions about ice rheology in the marginal zone. However, it is shown in the Appendix (A.4) that ice streams resting on a till layer less viscous than that suggested by Boulton and Hindmarsh (1987) or Anandakrishnan and Alley (1997b) will be subjected to unfeasibly large erosion rates. An implication of this paper is that ice-stream texture can potentially be explained by basal interactions, and it is the belief of this author that the debate over the partitioning of stresses between base and flanks will be settled by explanations of the smaller scale features within ice streams.

2. Drainage Theories and the Datum Effective Pressure

Much theory in soil mechanics is based on the idea of vertical effective-pressure gradients determined by static conditions. This assumption is valid over a wide range of conditions for vertical gradients within the soil but the situation is more complex along the ice–water interface. It seems that the interfacial effective-pressure gradient depends upon (i) whether ground-water flow is able to discharge Glaciol melt and (ii) whether any interfacial drainage system (e.g., channels cut into ice or till, or sheet flow or some other sort of distributed system) exists. In this section, we analyse interactions of static effective-pressure gradients and potential gradients necessary to discharge meltwater.

2.1. Pressure and stress fields: static interfacial gradients

We let (x, y) be the horizontal coordinates, z the vertical coordinate and t represent time. The effective pressure is defined by p _c = p – p _w where p = φp _s + (1 – φ)p _w is termed the bulk pressure of the sediment, φ is its porosity, p _w is the water pressure and p _s is the pressure in the sediment grains. Under static conditions the bulk stress p = p _i(D) + ρg (D – z) where ρ is the bulk density of the sediment, g is the acceleration due to gravity and D represents the upper surface and thickness of the till body. The water pressure is given by p _w = p ⁱ(D) – p _f + p _w(D – z) where ρ _w is the density of water and p _f = p _e(D) is the effective pressure at the ice–till interface. The effective pressure within the body is given by p _e = p _f + (1 – φ) (ρ _s – ρ _w) g (D – z) and at the base of the body by

(1)

(2)

The horizontal pressure gradients along the upper surface of the till are

which when combined with the definition of effective pressure yields

(3)

which integrates to

(4)

(5)

that is elevation causes an increase in the interfacial effective pressure. The quantity p _c is the datum ice–bed interface effective pressure. This is an important idea, representing the means by which hydraulic theories interact with the present theory. The quantity r is a drainage-model switch, taking on the values of 1 (interfacial effective pressures statically determined) or 0 (interfacial effective pressures constant, sec section 2.2).

The integration of Equation (3) can be carried out along an arbitrary line on the (x, y) plane; the theory is not restricted to one dimension. Over relatively long length-scales (depending on the hydrogeology) hydraulic potential gradients are comparable with static gradients and in such cases dp _e/dx| _z=D is given by a hydraulic theory such as the one developed by Walder and Fowler (1994). This is discussed further below.

The effective pressure at the base of the till body (i.e. at z = 0) is given by combining Equation (4) with Equation (1) to obtain

(6)

where

(7)

We also assume that the shear stress τ _b within the ice can be transmitted to the till body and does not vary greatly with depth. This is the “thin-till” approximation (Alley, 1989b). This assumption is mechanically dubious under some of the circumstances investigated in this paper and requires further investigation. It does represent the simplest mechanical configuration able to produce discharge of till and, as we shall see, produces quite varied results.

2.2. Pressure and stress fields: no interfacial gradients

The classical glaciological hard-bed assumption is that the effective pressure is constant (actually, zero) over the whole bed. Having a zero interfacial effective-pressure gradient gives rise to the usual glaciological calculations of water-pressure gradients which are set by ice-surface gradients (Paterson, 1994). We do still expect static gradients within the sediment and the effective-pressure formula is thus

Horizontal water-pressure gradients are given by the usual formula

where H is the thickness of the ice. This configuration applies (i) to the ease where the drainage system is very sensitive to small changes in the effective pressure, and can maintain an effective pressure which is constant, or (ii) the water pressure equals the overburden pressure and drainage is through a thin film. These ideas have been discussed by Alley (1996).

The same equations for surface and basal effective pressures (Equations (4) and (6)) can be obtained by setting r = 0, ensuring that the interfacial effective pressure remains constant.

2.3. Combined theory

In fact, we expect the situation to be more complex and for effective pressure to be determined by a mixture of static and hydraulic effects (e.g. Boulton and Hindmarsh, 1987). In particular, for any combination of aquifer and basal topography with relief magnitude [R], there exists a horizontal length-scale [L] such that the potential variations needed to drive ground-water flow are comparable with the static interfacial effective-pressure gradients. Over these length-scales we expect topographic variation to be reflected in the interfacial effective pressure.

There is some controversy over the spatial scale of drainage structures beneath ice sheets. If water flow is fully channelled then it is reasonable to assume that the channel spacing will be determined by [L] (Boulton and Hindmarsh, 1987; Walder and Fowler, 1994). In these cases, the hydraulic theories set the effective pressure in topographic depressions, where the effective pressure is lowest and drainage channels most likely to form.

Theories have also been proposed where the drainage structures have small spatial scales (i.e. less than the horizontal length-scale of the relief) (e.g. Alley, 1993, 1996; Walder and Fowler, 1994). These drainage structures are considered in the following analysis; low-pressure R-channels are specifically excluded. Closely spaced channels can co-exist peacefully with the large-scale drainage theories, as they may describe drainage of the watershed into the larger channels. It is expected that the interfacial drainage systems will have effective-pressure-dependent transmissibility; this is important because it can affect interfacial effective-pressure gradients.

We consider these situation in more detail by first showing how to compute [L] and then by considering the possibility of effective-pressure-dependent interfacial drainage.

2.4. Length-scales over which discharge by ground-water flow is possible

Elevation differences large enough to be of significance require relatively large bed slopes and thus relatively large static pressure gradients as one moves along the ice–bed interface (see Equation (3)). Over longer horizontal distances, where average slopes are less, the hydraulic gradients necessary to discharge subglacial melt become comparable with static water-pressure gradients, and the interfacial effective pressure is affected by hydraulic gradients.

It is a straightforward matter to investigate length-scales over which hydraulic gradients become significant and the following calculation is in essence the same as that presented by Boulton and Hindmarsh (1987), who hypothesised that tunnel valleys were spaced over length-scales where hydraulic gradients became large.

We first seek those length-scales where drainage is possible by ground-water flow alone. A necessary condition for this to be useful is that these length-scales are greater than the fluctuation length-scale. For length-scales shorter than those that appear in the following analysis, we expect interfacial effective pressure to depend primarily on elevation, while for longer length-scales, we expect to find the water pressure close to the ice overburden and regulated by an interfacial drainage system.

One constructs a vertical scale [R], a horizontal length-scale [L], and finds that the static interfacial-pressure-difference magnitude is (ρ – ρ _w) g [R]. Under the Dupuit assumptions (Bear, 1972, p. 361–366) the Darcy flux q _w over the depth of the underlying porous medium is given by

(8)

where ψ is the hydraulic potential, μ is the viscosity of water. m is the basal melt rate, Δ is the thickness of the porous medium and k is its permeability. The potential difference φ arising over a distance [L] is given by

At the length-scale where hydraulic and static gradients become equal we have

and can solve for the length-scale

Let us take μ = 0.001 Pa s, m = 10^-10 m s^-1, (ρ _w – ρ _i)g = 800 Pa m^-1, k△ = (10^-16 → 10^-12) m³, (see, e.g., Bear, 1972, table 5.5.1). The corresponding length-scale [L] is computed and shown in Table 1.

The lowest length-scales correspond to the most impermeable strata, and the lowest permeabilities chosen for illustration are very low indeed. Over good aquifers (e.g. k = 10^-15 m², D ≳ 10 m) the length-scales are comparable with those of drumlins. Such aquifer transmissibilties could be provided by a till deposit alone.

2.5. Drainage at the interface

We now consider the effect a distributed drainage system might have on the effective-pressure distribution. We use a heuristic function κ to describe effective-pressure-dependent flow at the ice–bed interface

(9)

We model κ as

(10)

where κ _c is a constant. This is, in functional form, the same as the Walder–Fowler model if we set λ = n, the Glen exponent, but we do not wish to restrict consideration to a particular drainage mechanism. This model causes the transmissibility of interfacial drainage routes to increase as the water pressure reaches the ice pressure and has been examined by Alley (1996). The basic idea is that as the water pressure increases, more and more roules become available to drain the water as the ice and till decouple.

Table. 1. Dependence of length-scale [L] of hydraulic/static equality on the hydrogeology k△ and the relief [R]

The value of κ _c determines whether interfacial routes play a role in subglacial drainage, while the exponent λ determines the range of effective pressure over which interfacial routes play a significant role; the greater λ, the narrower the range of effective pressures where interfacial routes drain significant quantities of water. The steady conservation equation

is then solved numerically using finite differences.

Firstly, it is solved in one dimension, across an infinitely long valley symmetric about the talweg with constant slope across the valley. Various cases are solved with k△/m³ ε {10^-17, 10^-15, 10^-13, 10^-11, 10^-9}. A no-flow condition is applied at the hill crest, while the effective pressure is set along the talweg, simulating the presence of a subglacial river. In all cases, this effective pressure was set al 10⁴ Pa.

We vary the constant κ ₀ as part of the experiment. A magnitude for the constant κ ₀ was specified as follows. We write κ ₀ = νκ _s. We suppose that the hydraulic potential is roughly of the same order as the channel potential. We then define κ _s as that constant which permits water to be discharged at around the channel effective pressure. Thus, to scale

whence

We then vary ν as a parameter; if ν < 1, then interfacial drainage can only occur at very low effective pressures and in consequence we expect effective pressures to be very low, and water pressures to track the ice pressure, while if ν > 1, interfacial drainage can occur at relatively high effective pressures.

Some typical calculation results are shown in Figure 2. Figure 2a shows cross-valley effective pressures for various transmissibilities; for each case the multiplier ν ε {0.01, 1, 100}. The results show that for low aquifer transmissibilities k△, the value of the multiplier ν determines the drainage style, with effective pressures being somewhat less than the channel effective pressure (ν = 0.01), of roughly the same magnitude (ν = 1) and statically determined (ν = 100). Figures 2b, c and d show the hydrogeology of a valley with a bump on its side (Figure 2b), the computed effective pressure (Figure 2c) and the transmissibility enhancement κ/k△ (Figure 2d). The transmissibility k△ = 10^-13 m³ and ν = 0.01. There is an obvious increase of effective pressure around the bump but the average vertical gradient following the surface is about 25% of the static value. That is, effective pressure does increase with bump elevation but not as rapidly as if effective pressures were statically determined. Very approximately, we can view the parameter r lying between zero and one. The transmissibility enhancement κ/k△ is shown in Figure 2d. The area of the bump has lower transmissibility, while its flanks have higher transmissibility and water from upstream is thus diverted around the bump. The same experiment, but for ν = 1, is shown in Figures 2e and f. In this case, detailed inspection of the results shows that while the effective pressure along the valley flank is not statically determined, the effective-pressure gradient following the bump is statically determined.

In short, we do not know enough about interfacial drainage in ice sheets to determine how interfacial effective pressure varies with elevation for bumps of different horizontal dimension but it is certainly possible that interfacial effective pressure increases with elevation, possibly at nearly static rates, and that this becomes more likely if the span of the bump decreases.

2.6. Effect of ice topography

An interesting feature of the results above is that they can be applied in a straightforward way to understand the effect of variations in the ice thickness. From a hydrogeological point of view the bump can also be considered to be a bump in the ice thickness reduced by a factor of α/ρ_i g. In this case the effective pressure is increased by the additional weight of ice. A typical value of α is about 0.1. Thus, small variations in the ice-surface topography can have marked effects on basal water pressures. By choking the drainage system, the increased effective pressure causes water pressures to rise under the ice mound. This reduces the increase in effective pressure and allows more drainage to occur.

2.7. Coupling of ice, till and water flows

Such mounds of ice or sediment arise in response to the coupled dynamics of the ice and the bed, and anything which affects the water pressure in this way will also affect the flow of ice and of till. This means that we need to consider the coupled flow of ice, till and water. We shall do this by looking at the kinematics of sediment and shock formation, an important topic for drumlin formation, and then we shall use the analysis in linear-stability analyses of coupled ice, till and the water-flow system. For this purpose it is now convenient to develop a linearized ground-water model.

We model potential and effective-pressure distributions, using the heuristic function κ (see Equations (9) and (10)) to describe effective-pressure-dependent flow at the ice–bed interface. Since we are now considering a dynamic theory, the water storage has to be considered, and we assume it to be given by (cf. Alley, 1996)

where the first term on the righthand side represents the normal aquifer storage and the second term represents storage in ponds, etc., at the interface which we also model according to some sort of power law

This is motivated by similar considerations that led to the transmissibility-effective-pressure relationship (Equation (10)); as water pressure increases, the ice and till decouple, leading to the creation of storage volume for water at the interface. Linearizing, we find that

where

and the first-order storage is given by

and we can compute

to find the first-order water flux

Under the Dupuit assumptions (e.g. Bear, 1972, p. 361–366), potential lines in aquifers are independent of depth. This means that we can identify the first-order water pressure (which is defined at the glacier base) with the first-order potential. Conservation then gives us the first-order water-pressure-evolution equation

(11)

where we have used

and defined a diffusion coefficient

(12)

We shall use this theory later.

Fig. 2. The relationship between effective pressure, hydrogeology and topography. The quantity ν determines at what effective pressure interfacial drainage becomes significant; high ν means it becomes significant at a relatively high effective pressure. (a) Dependence of effective pressure across a valley of uniform slope with varied aquifer transmissibility k△ and interfacial drainage coefficient ν. For each transmissibility the effective pressure is plotted as a function of distance x from the valley axis. Each transmissibility has three cases corresponding to ν ε {0.01, 1, 100}; the relative position of the lines is the same for each case. For high ν or high transmissibility, effective pressures are statically controlled (i.e. increase with elevation). (b–f) Study of the effect of a bump in a valley side on pressure and water flow. Basal axes are positions x and y. (b) Valley topography with 30 m bump. Case (a) corresponds to the same valley with no bump. (c) Variation of effective pressure p_e over valley side, ν = 0.01. (d) Variation in space of transmissibility enhancement arising from interfacial drainage κ/k△ relative to aquifer transmissibilty; ν = 0.01. (e) Variation of effective pressure p_e over valley side; ν = 1. (f) Variation of κ/k△ relative to aquifer transmissibilty; ν = 1. Note constrast in effective pressures between (c) and (e).

3. Sediment-Flow Kinematics

In order to couple ice and till flow we need to examine sediment kinematics. Our next step in investigating till flow will be to produce a kinematic-wave theory (e.g. Whitham, 1974) for which we need to know the sediment discharge-thickness relationships. We also consider some aspects of shock formation, which is postulated to be an important part of drumlinization and is one of the non-linear mechanisms which acts to control the growth of sediment-thickness instabilities.

3.1. Internal deformation flux relationship

The flux contribution arising from internal deformation can be computed from a postulated viscous relationship for till (Boulton and Hindmarsh, 1987)

(13)

where the subscript d refers to internal deformation. This flux computation is discussed in detail by Alley (1989b) who treats the special cases b = 1.2. Since the effective pressure varies linearly with the depth, we can compute the integral by the change of variables

whence

(14)

(15)

For b = 2 a special form exists

A special form also exists for b = 1 but we shall not consider this case. Both forms were derived by Alley (1989b).

3.2. Sliding

Two till-sliding laws have recently been proposed, a quasi-plastic (Cuffey and Alley, 1996) and a viscous one (Hindmarsh, 1996), which is the form used here. They need not be mutually exclusive; plastic behaviour may be an appropriate description for small length- and time-scales, while the viscous one is more appropriate for larger scales.

Hindmarsh suggested that the till-sliding velocity is given by

(16)

The flux contribution from sliding is Du _b , and the total till flux q is given by the formula

(17)

3.3. Scaling

There are a large number of free parameters in the flux expressions but these can be reduced by scaling.

The scale of the Glaciolly applied shear stress [τ] is regarded as an externally determined parameter. There are a number of reasonable choices in the way we scale the depth. It can be done by setting

(18)

i.e. setting the depth-scale equal to the scale of the depth of the base of deformation. It might seem logical to set the depth-scale according to the initial condition, but owing to the fact that there is a non-monotonic dependence of the flux on the thickness in a large volume of parameter space, this can result in the flux-scale computed below being completely unrepresentative. It is possible that in these cases the most suitable depth-scale is that depth which produces the greatest flux, but this scale is not always bounded, and in any case this depth is only definable implicitly (sec section 3.5). We thus use Equation 18 to set the depth-scale. Henceforth in this paper, apart from presentation of results of computations in sections 5.2 and 5.4, discussion is presented in dimensionless units.

The effective pressures are scaled

(19)

whence in dimensionless form

where

(20)

The messy denominator (2 – b) (1 – b) in the flux expression is essential; this ensures that q is an O(1) function for reasonable values of b. It is easy to show that in scaled form

and after defining

(21)

we can write the flux relation in scaled form as

(22)

The parameter δ is given by

(23)

We are free to choose [q], and the natural choice is to set it so that one of the two coefficients A _d, A _s is unity, with the other coefficient being the lesser. Thus, we compute

(24a)

(24b)

and set

(25)

and use this flux scale to compute the parameters A _s, A _d.

The parameters δ and Ψ depend upon densities and the acceleration due to gravity, which are well known, and the porosity φ, which can reasonably be expected to vary between 0.2 and 0.4, which is a small variation when compared with that conceivable in the other parameters. With ρ _s = 2700 kg m^-3, ρ _w = 1000 kg m^-3, ρ _i = 917 kg m^-3, φ = (0.2 → 0.4) gives δ = (0.06 → 0.08) and we shall take δ = 0.07 in this paper, which roughly corresponds to φ = 0.3. It is significant to have δ non-zero, as this (i) leads to reverse shock motion when internal deformation only is occurring and the till thickness is sufficiently large, and (ii) leads to the possibility of relief amplification when internal deformation only is occurring and the till thickness is sufficiently large. These processes can occur when the till is sliding even when δ = 0.

We see now that the parameters are the exponents b, d, the datum effective pressure p _c and the rate-factor ratio A _d/A _s. In this paper we shall not vary b and d independently. The fourth free parameter is B, the initial thickness of the sediment body.

3.4. Analysis of the flux expressions

For the purposes of illustrating the dependence of flux upon thickness, we consider the contributions from internal deformation and sliding separately. Defining P = D/p _c, we find from Equation 22 that we can write the flux relationships in an alternative form as

(26)

(27)

The monotone direct dependence of Q on p _c, expressed by the middle relationships, implies that whether it decreases with p _c for a given P simply depends on whether b > 2 and d > 1 for internal deformation and sliding respectively. Thus, relations in Equations (26) and (27) show how p _c enters as a rate factor. Since we do not know A _s, A _d, this is not as significant as the influence p _c has through its appearance in the variable P = D/p _c.

We illustrate Q(P) for both flow mechanism for varions b, d (Fig. 3). The constructions in Equations (26) and 27 show that the datum effective pressure has two clearly defined roles; as a rate factor, through its appearance multiplying the coefficients A _d, A _s, and as a depth-scaler, through its appearance in the term P = D/p _c. The thickness P determines whether a sediment body is thick enough to experience kinematic-wave-velocity maxima and zeros which are crucial to shock-generation properties.

Figure 3 shows that for small thickness, sediment flux increases with thickness in both cases. For sliding, a maximum flux is reached for all cases of d considered, and thereafter the flux declines with thickness, and asymptotes to zero. The flux increases because more sediment is being transported in a plug flow, despite the decrease in the sliding velocity. Eventually, the decrease in sliding velocity with thickness (effective pressure at the base) becomes a significant factor. A maximum flux is reached and thereafter the flux declines asymptotically to zero if there is no internal deformation present.

For internal deformation, the situation is more complex and is analysed further in section 3.5. The flux increases but the rate of increase falls as sediment thickness increases. This is because the average viscosity increases with sediment thicknesses since the effective pressure and till viscosity are high at depth. Moreover, as the sediment-body elevation increases, the effective pressure at the interface increases owing to the density difference between ice and water. This eventually causes a decrease in the flux for b > 2. This is also analysed in section 3.5.

It is convenient from the point of view of presentation of the flux relations to consider sliding and internal deformation separately, and to consider the dependence Q upon P rather than q upon D. Where sliding and internal deformation are operating together, this convenience is lost. In much of the rest of the paper, however, we consider B, the initial maximum thickness (which determines the subsequent evolution), and p _c as being separate parameters. This is essentially a matter of taste.

3.5. Kinematic waves and shocks

In Hindmarsh (1998a) drumlinization as a shock-formation process is modelled using numerical codes. Shock formation occurs when characteristics of the hyperbolic conservation equation

cross. Kinematic waves move along these characteristics and in this section we consider kinematic waves and shock formation in rather more detail than in Hindmarsh (1998a).

Standard kinematic-wave theory (Lax, 1973; Whitham, 1974) yields the following expression for the kinematic-wave velocity

We can construct

(28)

thus, W is proportional to the wave velocity ν. This fact is used in Figure 3 which illustrates the dependence of W and thus ν on P = D/p _c. i.e. on a combination of the thickness and the effective pressure. It should be remembered that the relationships in Equations (26) and (27), which define the ratio Q/q, imply a further monotone dependence of W upon p _c which does not affect the kinematic-wave analysis which follows.

Fig. 3. Till kinematics showing the dependence of Q(P) ∝ flux q and W (P) ∝ kinematic-wave velocity ν on thickness for internal deformation (a–c) and sliding (d–f). Parameters are b (internal deformation), d (sliding) and P = D/p_c. (a, d) Graphs of Q (never negative) and W for internal deformation (a) and sliding (d). (b, d) = 3. (b, e) Three-dimensional plots of Q on P and b or d for internal deformation (b) and sliding (e). Note viewing angles differ. (c, f) Three-dimensional plots of W on P and b or d for internal deformation (b) and sliding (e). Note viewing angles differ. Note that kinematic-wave velocities are negative in some regions of the parameter space.

3.6. Kinematic-wave expressions: sliding

Let us consider the dependence of kinematic-wave velocity on thickness. Consider first the case of sliding only.

(29)

(30)

and from Equation (28) we obtain

(31)

The ratio D/(D + p _c) < 1 but increases with D while the quantity d is expected to be greater than one; plasticity theory gives it as ∞. If d > 1, this expression permits kinematic waves which move both forwards and backwards, with velocities becoming more negative as the till becomes thicker. This can be seen in Figure 3. The implication of this is that shocks will form on the upstream side of the bodies, which will therefore have blunt upstream ends, as do drumlins. The thickness D _s at which the kinematic-wave velocity becomes zero is given by

3.7. Kinematic-wave expressions: internal deformation

For internal deformation alone, the wave velocity, its large thickness asymptote and its derivative are given by

(32a)

(32b)

(32c)

The corresponding expression for W (P) is

(33)

which is also illustrated in Figure 3. This shows that kinematic-wave velocities increase with P (thickness) up to a point, whereafter they decline, sometimes reaching negative values.

3.8. Further analysis

3.8.1. Conditions for negative wave velocities, internal déformation

The asymptotic formula for large D (Equation (32b)) shows that cases exist, depending on b and δ (recall Ψ = (1 – δ) (b – 2) + 1), where the kinematic velocity is positive or negative for sufficiently large D. We can write the condition for ν(D → ∞) → 0 (assuming b > 1 ) is

Since by construction 0 ≤ δ ≤ 1, we can easily show that both sides of the equation are negative (positive) as b is greater than (less than) 2, with equality occurring for b = 2. Thus, negative wave velocities are expected for some thickness when b > 2. To be precise, if δ = 0, (i.e. horizontal hydraulic gradients are significant) the wave velocity asymptotes to zero for large thicknesses but never becomes negative. This is the only significant qualitative difference between δ = 0 and δ > 0.

We are particularly interested in whether negative wave velocities can occur for realistic drumlin thicknesses. The length-scale [D] is maximally of the order of ten metres, while drumlins a hundred metres thick represent a very approximate observational upper limit. We take the upper limit of scaled drumlin thickness to be 40 and, using Equation (32a), look for those parameter ranges which give a flux maximum, and thus a kinematic wave of zero, and compute the (unique) corresponding thickness D _q if it exists in this range. These values are shown in Table 2 and demonstrate that the maximum thickness decreases with p _c and decreases as b increases. There is an increasing expectation of the thickness of maximum discharge being greater than the thickness of the original sediment body as b decreases and as p _c increases.

Table. 2. Dependence of thickness of maximum flux D_q on the parameters b and p_c when a sediment body is deforming internally. The maximum was sought for in the range 0D40; where it was not found in this range, a NaN (i.e. the IEEE undefined number) is specified. The search range is a plausible upper limit to the range of scaled drumlin

Of course, there is a linear dependence of D _q upon p _c, which arises from the symmetry of terms in D and p _c in the righthand side of Equation (32a) which we have previously expressed through the ratio P = D/p _c. Inspection of the results in Table 2 shows that

which is analogous to the expression for the case of sliding only (Equation (29)), but this is only approximate and not deducible from the wave-velocity expression; it may be an asymptotic result.

3.8.2. Thickness of maximum kinematic-wave velocity

In a similar way, the thickness D _v corresponding to the maximum kinematic-wave velocity may be computed. Where a maximum exists, again we expect reverse-facing shocks, as kinematic waves from thinner, faster regions catch up with those from thicker, slower regions. Thicknesses corresponding to maximum kinematic-wave velocities are shown in Table 3. The same symmetries ensure that D _v ∝ p _c, and inspection of Table 3 shows that there is an inverse relationship between D _v and b. The significance of D _v is that any sediment body which has an initial thickness less than D _v will not form upstream-facing shocks. It can be seen that this is favoured by low effective pressures and high b. The maximum kinematic-wave velocity for sliding occurs when D = 0. Upstream-facing shocks are therefore inevitable for sliding.

3.8.3. Downstream-edge shocks

Λ further point of importance is that it can easily be shown that for flow by internal deformation

which means that at zero the function is concave and up to a certain thickness we expect shocks to form at the downstream edge whatever the rheological index might be.

In the case of sliding, it is easy to see from Equation (29) that at zero thickness the kinematic-wave velocity is the sliding velocity; and that dν/dD = , meaning that wave speeds become faster as the sediment thins. No shocks are expected at the downstream edge.

Table. 3. Dependence of thickness D_v at which the maximum kinematic-wave velocity occurs on the parameters b and p_c. It depends on p_c, b and δ and no other parameters

3.8.4. Equality of kinematic-wave speed and shock speeds

It is of interest to know if there is any thickness where the kinematic-wave velocity is equal to the shock speed, where the till thickness on the other side of the shock vanishes. It is easy to show for sliding that this does not occur at any finite thickness.

For internal deformation, equality occurs when

which we can write in terms of P = D/p _c

For 1.5 ≤ b ≤ 4, the P at which equality occurs declines from about 4 to a bit less than 1. For example, if p _c were 0.1, P would range from 0.4 units to 0.1 units: roughly 4 m to 1 m. For values of b less than about 1.5, P increases very rapidly, above likely original thicknesses.

3.8.5. Effect of shocks on total relief

The inevitable effect of shock formation is to reduce the rate at which relief grows or to cause it to become negative. Shocks move through the flowing material at different rates to the kinematic waves, meaning that a shock will eventually move to an elevation maximum. At this point, the elevation of the higher shock edge is that of the local maximum, and as the shock continues to move, the elevation of the higher shock edge will decrease. In this way, the rate of change of relief is made (more) negative.

Where instabilities exist to cause sediment relief to increase, we expect shock formation to occur. It is likely that this could be sufficient to control the instability in an average sense.

3.8.6. Smoothing of shocks

A shock is an ideal state where the thickness of the sediment jumps. In practice, the jump is smeared by diffusive processes. These processes are likely to express themselves over length-scales equal to the thickness of sediment deformation. Firstly, at this length-scale, the weight of the till affects the shear-stress distribution. Including this term will cause a diffusive term in the sediment thickness to appear, which is of significance at the length-scale of the deforming-layer thickness. Secondly, over this length-scale, it is not at all clear whether till behaves as a viscous fluid anyhow (Murray and Clarke, 1995; Hindmarsh, 1997a). The effect of this is unknown but is unlikely to operate over length-scales much longer than the deforming-layer thickness. This length-scale is of the order of a few metres and we thus expect the expression of shocks to be blunt, but not vertical faces with horizontal expression of at the very most a few tens of metres. Diffusion also leads to the decay of shocks, but this also occurs in the absence of diffusion, so is not a serious objection to the drumlinization as shock formation discussed in Hindmarsh (1998a).

4. Sediment-Thickness Variations and Ice-Surface Texture

The increasing spatial coverage of ice-stream velocity fields has led to the observation that there are spatial variations in ice-sheet velocity over scales somewhat smaller than the ice streams (MacAyeal, 1992a; Alley, 1993; Whillans and Van der Veen, 1993). Of interest are recent observations by Bindschadler and others (1996). They find that there is increased roughness in the surface topography of Ice Streams D and E associated with areas of increased bed friction, and conclude that this variation in surface topography is due to variations in bed friction rather than basal topography. We ask whether this might be explained by coupling of the HTTA theory with ice-flow theory. Two tractable situations are long-wavelength variations in basal topography (i.e. wavelength much greater than ice-sheet thickness) and short-wavelength variations in basal topography (wavelength much less than ice-sheet thickness). In the latter case, it is known that under certain conditions, sheet flow of till is unstable (Hindmarsh, 1998b) and it is hypothesized that this might lead to zones of drumlin formation, which increase basal drag (Smalley and Unwin, 1968; Boulton, 1987; Menzies, 1989). We look at long-wavelength perturbations below.

4.1. Long-wavelength bed variation and ice-stream surface texture

We investigate relationships between basal and long-wavelength ice-sheet topography using a perturbation method. This is a long-wavelength theory and can determine the necessary basal relief to create sufficient frictional variation to induce observed ice-sheet texture.

We first need to establish the bottom-velocity relationship for ice. We define scales for basal-ice velocity by

where subscripts d, s refer to deformation and sliding respectively. This scaling is equivalent to the flux scaling.

We need to consider the scaling of the applied shear stress τ a little more carefully. Using the variant of the Hutter–Morland–Fowler scaling (Hutter, 1983; Morland, 1984; Fowler, 1992) described by Hindmarsh (1997b) we can write the shear-stress scale in terms of the ice-sheet-thickness scale [H], the ice-sheet-aspect ratio or bed slope ε and the density of ice

and use the same scalings for the shear stress, effective pressure and vertical distance. We also define a characteristic length [X] = [H]/ε. We now need to introduce the effect of varying ice thickness on the interfacial effective pressure p _f.

At the point where we define the datum effective pressure, we define a datum thickness H _D, and we then construct

The hydrostatic-pressure scale is ρ _i g[H]. which implies that a term A _*/ε accounts for the contribution to effective pressure from the varying thickness of ice. In scaled form therefore, the additional thickness enters divided by a small parameter, and is thus likely to play a very important role. The surface effective pressure is given by p _f = p _c + rH _*/ε + δD, and the basal effective pressure by p _b = p _c + rH _*/ε + D. We know from the shallow-ice approximation that τ = –H∂ _x (H + μD) where

is the conversion factor from “till-thickness units” to “ice-sheet-thickness units”. We can also write

whence

Again, the parameter r enters; if r is zero, then the inter-facial effective pressure is given by p _c (recall δ is a product of r and some other terms), and this configuration represents the usual glaciological hard-bed case; if r = 1, interfacial effective pressures are statically determined.

We find in scaled form that

(34)

where

(35)

where

We have chosen the length-scale [X] to be the characteristic span of the ice sheet, so that ε = [H]/[X]. This choice of scales ensures that ice- and till-velocity units are scaled similarly. We need to ensure that we choose the till-flux scale [q] (see relationships in Equations (24a) and (24b)) such that [q] = [u][D]. We do expect the scaled velocities to be order one, otherwise our assumption that all the deformation is occurring by till internal deformation cannot be correct, but we can no longer expect the till flux in be order one in general.

Let us write the velocity functions from Equations (34) and (35) as

and consider sensitivity coefficients

(36)

Finally, we introduce a length-scale for the bed forcing Ξ < 1 (in dimensionless units) and write a further scaling

The perturbed ice flux is given by

In steady state q ₁ = 0. Consider first, for the sake of argument, the case where, f ₁, D ₁ are prescribed. We expect the coefficients R _t, R _n, R _D to be O(1), and R _f = O(rδ).

Consider first the case where Ξ > O(μ, ε). If r is O(1), it is clear the dominant forcing terms (d) and (e) will be balanced by the (a) term on the righthand side, and that we can expect H ₁ = O(εD ₁, εδf ₁). If, on the other hand, r = 0, terms (c) and (e) are zero and we expect the (d) term to be balanced by the (a) term and we thus expect to find that H ₁ = O(ΞD ₁). We can illustrate by taking [H] = 1000 m, [D] = 10 m, ε = 0.01. Then, if r = 1, δ = 0.07 and f ₁ = 0.1 (i.e. 1 m), since H ₁ = O(εδf ₁) this should force a change of 70 mm in the elevation of H. If on the other hand it were sediment relief of 1 m (i.e. D = 0.1), we should expect H ₁ = O(εD ₁), a perturbation of 1 m. Now let us consider the case r = 0. Since H ₁ = O(ΞD ₁) this could force quite a large change in H ₁ for longer-wavelength variations in the till thickness; the ratio of the ice-thickness change to the till-thickness change in physical units is Ξ/μ. Thus, if Ξ = 0.1 and μ = 0.01, the ratio would be 10, and a 10 m till mound would produce a 100 m ice mound. However, we do not expert a till mound to persist in one position but to exist as a travelling wave, with the implication that substantial ice relief due to till mounds will also be a travelling wave. This is investigated in the next section. Finally, consider the case r = 0, Ξ = O(μ, ε). Then, we find that H ₁ = μD ₁, which in physical terms means that a 1 m change in D ₁ or f ₁ produced a change of the same order of magnitude in the ice thickness.

Note that the dependence on bedrock topography arises for different reasons (the effective-pressure variations) than the classic reasons discussed by Budd (1970) and Jóhannesson (1992) and that this topography is measured in till units (i.e. orders of a few metres) rather than ice-sheet units (order of a kilometre). In both cases we do not expect to see a very obvious relationship between ice topography and basal topography; when r = 1, a small additional ice weight (the response from term (c)) can suppress topographic influence, while when r = 0, the influence of basal topography on ice-surface topography is primarily controlled by variations in till thickness and the length-scale over which this happens (response from term (a)). Amplification (in the sense that the ice-surface relief is greater than the inducing bed-surface relief) can occur.

How does this analysis accord with such observations of ice-surface rugosity? Typical normalized short-wave fluctuations in the Ice Stream D and E profiles are maximally 0.05 (Bindschadler and others, 1996), which can be obtainable from small variations in till thickness in the appropriate part of parameter space. We can say that the HTTA theory suggests correlation between bed and surface topography should not be strong, as the influence of till thickness and effective-pressure variations will dominate, in general agreement with Bindschadler and others’ results, but there is no specific validation. The HTTA scale theory does not provide a simple explanation of why high traction and ice texture should be associated; we shall see in the next section that more sophisticated calculations do suggest the HTTA theory can explain this phenomenon.

It is not really possible to produce a sticky spot in a one-dimensional analysis such as this. If the ice thickness is perturbed by 10%, then continuity implies that the steady velocity will also be perturbed by 10%. An investigation of sticky spots really requires a two-horizontal-dimension study to determine whether basal perturbations such as those described in this section can produce sticky spots.

4.2. Sticky-spot migration and the Smalley–Unwin bifurcation

Smalley and Unwin (1968) proposed dial drumlin fields arose as a result of spatial variation of stability properties. I do not agree with their proposed method of relief amplification but the idea that there is a spatially varying stability parameter which determines whether drumlins form or not is of great significance. The possibility of some counter-intuitive effects occurring exists; for example, Hindmarsh (1998b) shows that for the short-wavelength instability, an increased effective pressure tends to stabilise the sheet flow. An implication of this is that an increased effective pressure would cause a drumlin field to disappear. If the drumlins, by protruding into the ice are creating an obstruction to flow and increasing the basal drag, their disappearance implies a reduction in the form drag. To see whether the disappearance of drumlins can cause overall bed roughness to decrease with increasing effective pressure, we need to look at the drag contribution from drumlins. The following analysis parallels that presented by Alley (1993) with some emphasis on ‘‘skin friction” (the drag at the interface) as compared with “form drag” (increased stresses generated within the ice).

We consider flow around cuboid shapes of dimension (L _x , l _y , l _z ), and consider the increase in skin friction only along the upper surface. We configure the problem such that the drumlin has a blunt upstream face of width l _y , elevation l _z and upstream slope l _z /l _x , and overall length L _x . The strain-rate is u/l _x , the stress and the force . Here, . The extra skin friction is given by rαl _z l _y L _x . If the drumlin density in terms of drumlins per unit area is denoted Φ, we have Φ = l _y L _x / Ξ and Ξ = l _y L _x /Φ, which means that the increase in shear stress is given by

Table. 4. Contribution to basal shear stress from drumlins assumed 30 m high, occupying 20% of the bed. The rows represent drumlin length, the columns the horizontal length of the drumlin blunt end. Stress in Pa. The contribution of the skin friction is around 5 ×10³ Pa, independent of the plan dimension. Stress contributions of both components are proportional to drumlin elevation and to the ground covered by drumlins. Drumlins can provide a substantial proportion of the driving stress of a typical ice stream

For cases where the drumlins make a significant contribution to typical basal shear stresses (in ice streams they are typically 20–30 kPa), form drag is the dominating component (Table 4). If till-sheet flow is unstable to sediment-thickness perturbations, drumlin formation can occur very much more quickly and could play a role in ice-sheet dynamics, and this would only occur for soft-bedded ice sheets. Migrating areas of drumlin formation or annihilation could play a significant role in ice-sheet dynamics.

5. Linear-Stability Analysis of A Coupled ICE–till–Water Flow System

We now introduce dynamic coupling between ice, till and water flows though considerations of the linearized equations. Linear-stability analyses of ice sheets date back to the 1950s (Bodvarsson, 1955; Nye, 1959). For an ice-sheet body with a finite span, the linearized eigenvalue problem is of Sturm–Liouville form (see also Hindmarsh, 1997b) which implies that the basic response is relaxation by diffusion. For the infinite plane (e.g. Paterson, 1994) perturbations propagate as well as relax. In either case, perturbations disappear very quickly by a basically diffusive mechanism. Further modes of behaviour are introduced by the flow of till and water, but the diffusive mode always seems to exist. We call this mode the Nye diffusion mode.

5.1. Model formulation

Firstly, we consider the coupling of ice flow and till flow. The model configuration is an infinite plane with slope ε, ice of thickness H ₀ overlying a deforming bed of thickness D ₀. The till rate factor and constant effective pressure p _e0 are chosen such that the velocity on the top of the till is μ ₀. By assumption, internal deformation within the ice is negligible, so the ice is moving as a plug at that speed.

In this section we analyse the linear stability of a system of coupled ice and till flow, with and without hydrological modelling. The mechanical model of the ice is the shallow-ice approximation, meaning that this is a long-wavelength theory (wavelength longer than the ice-sheet is thick), in contrast to the short-wavelength theory discussed by Hindmarsh (1998b) where unstable growth was found to occur. The perturbations are periodic along the plane, and we consider perturbations for wavenumbers k.

First of all we consider cases where the water pressure is uncoupled. To first order, we have

where changes in the interfacial effective pressure due to sediment-thickness changes are included in the R _D term. Hence we can write

where the R _i = {R _t, R _n, R _D, R _f} are defined in Equation (36) and section A.1, and continuity for ice gives us

(37)

Similarly, we can write down the till-continuity equation as

where

(38)

to obtain

(39)

Expanding the perturbations in terms of trigonometric functions

where k is a wavenumber, and substituting in the evolution Equations (37) and (39) gives us mode-evolution equations for ice thickness

(40a)

(40b)

and till thickness

(41a)

(41b)

where

Stability is determined by the eigenvalues of this system.

An important point is that the coefficients in this equation system are all proportional to the till deformation rate factor. This means that increasing the rate factor does not effect the location of bifurcations in the parameter space. However, an increase in rate factor will lead to a proportional increase in the basal velocity and in the rate constants (i.e. the eigenvalues of the system). Having more than one order of spatial derivatives means that scale dependencies are not simple.

5.2. Computations for the linearized ice–till system

In this section results are given in dimensional units. Parameter space for the coupled ice–till system (Equations (40) and (41)) has been explored by sampling at approximately 60 000 points for r = 0 and r = 1 (constant interfacial effective pressure and statically-determined interfacial effective pressure respectively). The set of sample points S _wt comprised the direct product of the sets D ₀/m € {1, 2, 5, 10, 20, 50}, H ₀/m € {100, 200, 500, 1000, 2000}, L/m € {200, 500, 1000, 2000, 5000, 10 000, 20 000}, p _e0/10⁵ Pa € {0.05, 0.1, 0.2, 0.5, 1}, τ₀/10⁵ Pa € {0.2, 0.5, 1, 2}, b € {1.02, 2.03, 5.10}, a € {1, 2, 5, 10}, where L = 2π/k. Excluded from S _wt are the points in parameter space where L < H ₀, as in this region the shallow-ice approximation is not valid. The basal velocity was set to be 100 m a^-1, by adjusting the rate factor. This study was carried out for internal deformation only. For r = 1, only 1% of cases were unstable (had positive growth rates), while 50% of parameter space contained instabilities for r = 0. This is the opposite dependence from the shorter wavelengths considered by Hindmarsh (1998b). Figure 4 shows the proportion of cases with positive-growth-rate binning along each parameter in turn for r = 0. This is a rough indication of how the parameter affects stability. These include very strong dependencies on the thickness and the effective-pressure index b. The system is most unstable for rather thin tills; this is the case the least likely to develop into a surface of large relief.

Closer inspection of the results shows that there is a complex conjugate pair of eigenvalues with large negative parts (relaxation times of less than a year), and these can be shown to be approximately what would be expected if we were simply dealing with short-wavelength ice-sheet diffusion (Nye, 1959; Paterson, 1991). We shall call this rapidly decaying mode the Nye diffusion mode. The robustness of this mode stems from the fact that it produces diagonally dominant elements on the matrix corresponding to the Equation system (40) and (41) which dominate one of the pairs of complex conjugate eigenvalues.

The other two eigenvalues form a complex conjugate pair with a much longer time constant, between 100 years and 100 000 years or longer, and sometimes with a positive growth rate. The eigenvectors, and more specifically, the orthogonal matrix associated with the Schur decomposition, show that these modes are associated with till-profile evolution, as opposed to the fast stable modes, which are strongly associated with ice-profile evolution. Experimenting with model perturbations showed that these eigenvalues were nearly associated with neutral stability. A consequence of having a complex conjugate pair of eigenvalues is that the solutions have a travelling-wave component. However, with the Nye diffusion modes, the decay is so fast that the wave-like aspect of the behaviour is barely discernible. In contrast, the pair of eigenvalues associated with the nearly neutral modes can have a discernible travelling-wave component. Figure 5 shows the proportion of cases with negative wave-velocity binning, the whole result set against each parameter in turn, for various wave velocities. Negative wave velocities only occur for r = 1 (for internal deformation) which is consistent with the kinematic-wave analysis presented in section 3. The figure shows a complex set of dependencies, with the proportion of waves moving backward increasing with sediment thickness (consistent with the kinematic-wave analysis) and decreasing with length-scale (not covered by the kinematic-wave analysis).

Fig. 4. The coupled flow of ice and till. Parameter study of bed-mode growth rates on the following parameters: till thickness D₀/m, ice thickness H ₀/m, wavelength L/m, zeroeth-order effective pressure p _e0/10⁵ Pa, zeroeth-order shear stress τ ₀/10⁵ Pa, effective-pressure index b and shear-stress index a. Showing the proportion of positive growth-rate constants binned against each parameter. Abcissa is the parameter value while ordinate is the proportion. Shown for all positive growth rates and for growth rates greater than 10^-5,10^-4 and 10^-3 a^-1 Computations are for U_b = 100 m a^-1; growth-rate constants are proportional to basal velocity in this linearized analysis, so for example for U_b = 1000 m a^-1, the rate constants quoted on the ordinate label should be increased by a factor of ten.

Fig. 5. The coupled flow of ice and till. Same parameter study as for Figure 4. Showing the proportion of negative wave velocities binned against each parameter. Abcissa is the parameter value while ordinate is the proportion. Shown for all negative wave velocties and for growth rates greater than 10^-5, 10^-4 and 10^-3 a^-1. Computations for U_b = 100 m a^-1.

These modes arise from a kinematic-wave mode, which in the uncoupled linear approximation would simply be a travelling wave of sediment. Coupling with the ice turns this into a wave which can grow or shrink depending upon the parameters, and whose small real part value is a conséquence of the fact it is associated with the neutrally stable wave (eigenvalue with zero real part).This also means that the growth rate can vary by several orders of magnitude and provide time constants so long that the instability would never be seen in practise, and the regions of parameter space where glaciologically significant growth rates (i.e. time constant of centuries) can be seen are rather more restricted. This is called the bed mode or bed-wave mode, and the real part of its eigenvalue is denoted by G ^B.

However, it seems that instabilities for reasonable till thicknesses (greater than ten metres) develop so slowly (rate constants ≲ 10^-3 a^-1 that in general, they compare with, or are greater than, expected ice-stream occupation times, and therefore unlikely to have sufficient time to develop. More detailed investigations show that for low ice thicknesses (which imply large slopes because we are keeping the shear stress constant), and low till thicknesses, growth rates are fast; however, this really corresponds to a valley glacier. In this case the basal velocity is 100 m a^-1. If the basal velocity were 1000 m a^-1, growth rates would be ten times greater, but even so they are not large enough in general to be expected to play a significant role in ice-sheet dynamics. This should he contrasted with the growth rates expected for smaller drumlins (Hindmarsh, 1998b) which are sufficiently large to have dynamical consequences. These low growth rates should also be compared with the much faster time constants introduced when there is hydraulic coupling. This is considered in the next section.

The effect of steady response to a 1 m sinusoidal bedrock wave is shown in Figure 6 (r = 1) for the parameters indicated in the figures. For small till thicknesses the response is mainly taken up by the ice, while for larger thicknesses it is taken up by the till. That there could be such a complicated response is deducible from the scale analysis and the magnitudes of the response are consistent with the scale analysis. The response is rather simpler for r = 0, where the interfacial effective pressure is constant in space (not shown in a figure).

Some evolutions are presented in Figure 7. Two initial conditions are considered, one where a 30 m sine-wave perturbation is applied to the ice-sheet surface, the other a 4 m till wave. Figure 7a a shows an evolution with an ice impulse applied to the nearly always stable case r = 1, where inter-facial pressure are statically controlled. After some rapid Nye diffusion, a travelling wave emerges. In general, this wave does not move at the kinematic-wave speed (it will only do so when the ice and till perturbations are uncoupled). Then, decaying bed modes emerge, of very limited amplitude, which in this case move upstream. Upstream-moving waves are also found when coupling with the water system is introduced. Figure 7b shows the response to a till impulse, in a case where interfacial effective pressures are constant. The wavelength is longer and the ice-surface response is very marked, there being an ice-surface wave of amplitude 40 m (i.e. total relief 80 m). This case is unstable, and the till wave grows slowly.

That there is a bigger response at long wavelengths is predictable from scale theory. At high effective pressure, the ice-surface expression of basal-till-thickness variation is rather marked. This is in accordance with the observations by Bindschadler and others (1996) that there is increased roughness in surface topography of Ice Streams D and E associated with areas of increased bed friction.

5.3. Coupling water flow

Now we consider the flow of water through the basal-drain-age–aquifer system. The zeroeth-order configuration is a flow of water down the plane, through the aquifer and inter-facial drainage system with potential gradient –ερ _i g in physical units. Since we are dealing with an infinite plane, no recharge can be modelled.

Firstly, the ice- and till-thickness-evolution equations require extra terms to account for the effect of water pressure on the effective pressure, and become, in dimensionless form,

(42)

(43)

In section 3 it is shown that the water-flow Equations (11) can be written in dimensionless form as Equation (49). Down the inclined plane, the zeroeth-order potential gradient is –1, so we arrive at the dimensionless equations for water-pressure evolution

(44)

The coefficients in this linear equation have a complicated dependence on the aquifer characteristics and the interfacial flow. In Appendix 3 it is also shown how a configuration where there is dependence on only two new parameters can be created (relation in Equation (50)). These parameters are F (a dimensional quantity) and λ, defined in Equation (12). The first is a water-pressure diffusion coefficient (in other words, by varying it, we are seeing how characteristic time-scales of the drainage system interact with the flow of ice and till), while λ represents an index of how much the transmissibility of the drainage system changes with effective pressure. The Γ_ww coefficient is proportional to F, while the coefficients Γ_w, Γ_i and Γ_D are proportional to the product of F and λ.

There is currently one hydrological model which yields the water-pressure diffusivity along the interface, due to Anandakrishnan and Alley (1997b). They suggest that the diffusivity is between 10 and 100 m² s^-1. However, it is doubtful that it is this well constrained, and we consider a range of values necessary to reproduce the cases where the interfacial effective pressure is statically controlled, where it is constant, and intermediate cases.

Fig. 6. The coupled flow of ice and till introducing bedrock typography (i.e. beneath the till). Showing steady profiles for r = 1 (i.e. constant interfacial effective pressure) arising as a consequence of a 1 m sinusiodal variation in bedrock topography at indicated parameter values. Basal axes are position x and the parametric variable D₀, the till thickness. Vertical axes are ice thickness (tl), till thickness (tr), shear stress (bl) and effective pressure (br). Note how as the till thickens the response is taken up by the till thickness rather than the ice thickness.

Fig. 7. The coupled flow of ice and till showing the evolution of ice and till profiles. Basal axes are position x and time t. Lefthand column is ice thickness, righthand column is till thickness. Cases (a) and (b) are for indicated parameters. The initial condition was (a) a 30 m sine wave in the ice, interfacial effective pressures statically determined (r = l) and (b) a 3 m sine wave in the till, interfacial effective pressures constant (r = 0). Case (a) is stable, there is rapid decay of the Nye diffusion mode, while the slower decaying bed mode expresses itself as an upstream-moving wave. Case (b) is unstable and is a case where smaller absolute variations in deforming layer thickness lead to larger absolute variations in the ice thickness. Again, rapid decay of the Nye diffusion mode occurs, while the bed mode grows slowly.

We now expand the water pressure periodically

our final set of mode equations is

(45a)

(45b)

(45c)

(45d)

(45e)

(45f)

The two uncoupled drainage styles can be retrieved as special cases. The case where the effective pressure remains constant corresponds to the righthand side of Equations (45c) and (45f) being zero (e.g. very long wavelength or low transmissibility). The case where the water pressure does not change occurs when the only non-zero coefficients are s _w and Γ_ww Sufficiently large k ² G _ww/S _W and small λ ensure this. Another special case can be retrieved by setting Q _t, Q _n, Q _D and Q _f to zero, implying no motion of the bed. The basal-ice velocity can be regarded as due to sliding, and the hydraulic system some sort of linked-cavity system; in other words, a hard-bed configuration.

As an example, calculations to investigate how far the parameters F and λ had to be varied to obtain the limiting cases were carried out for L = 3000 m, H ₀ = 1000 m and U _b = 100 m a^-1. These showed that in order to retrieve the bed-mode eigenvalues for r = 1 (interfacial effective pressures statically controlled) that the required values of the hydrological parameters were λ = 0 and F = 10 000 m² a^-1, while in order to retrieve the case r = 0 (interfacial effective pressures constant) λ = 1 and F = 0.01 m² a^-1. The value of 50 m² a^-1 suggested by Anandakrishnan and Alley (1997b) thus falls roughly between the two. This would not be the case for other length-scales, and a scale analysis could be usefully carried out.

The dynamics of the system can be analysed by computing the eigenvectors and eigenvalues of the system. This yields three complex conjugate pairs of eigenvectors and eigenvalues. One pair of modes always has a dominant projection onto the bed profile. This is called the bed mode, and the real part of its eigenvalue is denoted by G ^B. The other two mode pairs have significant projections onto the water pressure and the ice thickness. One of these mode pairs always has negative real part of its eigenvalue, and this one is defined as the Nye mode, and the remaining mode pair is called the pressure mode, with the real part of its eigenvalue denoted by G ^b. On occasion, this mode pair has a stronger projection onto the ice thickness than the other mode pair but in these cases the maximum difference is less than 3%. In general, the mode pair with negative real part eigenvalue projects more strongly onto the ice thickness.

5.4. Computations for the linearized ice–till–water system

In this section results are given in dimensional units. Parameter space has been explored for Equation (45) by sampling at approximately 20 000 points for a mobile bed and a fixed bed (i.e. Q _i = 0). The set of sample points S _wtp comprised the direct product of the sets λ € {0.1, 1, 10}, F/m²a^-1 € {10^-4, 10^-2, 1, 10, 10⁴, 10⁶, 10⁸, 10¹⁰}, D ₀/m € {1, 3, 10, 30, 100}, H ₀/m € {100, 300, 1000, 2000}, U _b/m² a^-l € {100, 200, 500, 1000}, L/m € {300, 1000, 3000, 10 000}, p _e0/10⁵ Pa € {0.1, 0.2, 0.5}, τ ₀/10⁵ Pa € {1}, b € {2.03}, a € {1}, where L = 2π/k. Excluded from S _wtp are the points in phase space where L < H ₀, as in this region the shallow-ice approximation is not valid. Again, varying the basal velocity is accomplished by varying the deformation rate factor A _d. Till sliding was not investigated in this study. The stability properties of this coupled system are somewhat more complex then the uncoupled case. There, the Nye diffusion mode generally has a strong projection onto the ice thickness only, and in general the other mode projects strongly onto the till thickness. In this coupled case, there are unstable water–ice modes with extremely last growth-rate constants. These occur at a rather specific range of water diffusivities.

Figure 8 shows the proportion of cases with bed-mode growth-rate constant greater than indicated amounts, binning along each parameter in turn. This is a rough indication of how the parameter affects stability, and how this qualitative demonstration of sensitivity varies with the growth-rate constants. Figure 9 shows a similar diagram for the pressure modes, for the case of a mobile bed. A similar diagram (not included) yielding fairly similar results occurs for the case of a fixed bed. All of these diagrams show the complexity of the dependence of the stability properties of the system upon glaciological parameters. The dependence of these parameters on the hydraulic diffusivity F is shown in Figure 10, which reinforces the idea that the stability properties of the system are very complex. The distribution of fast unstable modes is different from the distribution of all unstable modes, being concentrated in areas of higher diffusivity. At higher diffusivities still, the system stabilises, becoming equivalent to the static controlled uncoupled case (i.e. r = 1) when λ = 0. For the most part, the fixed bed (Q _i = 0) is more stable than the mobile bed. There are significant numbers of unstable modes in the parameter range around F = 50 m² a^-1 argued by Anandakrishnan and Alley (1997b) to be appropriate for ice streams.

Figure 11 shows the steady response to a 1 m sinusoidal variation in bedrock topography, at indicated parameter values. The parametric variable is λ, the exponent representing the degree to which the transmissibility of the system changes with effective pressure. Of particular note is the fact this drainage parameter affects the ice-surface response.

Figure 12 shows the evolution of ice, till, water-pressure and effective-pressures profiles, for indicated parameters; in particular, note that F = 10⁸ m² a^-1. The initial condition was a 0.1 m sine wave in the ice. This case has a very high growth-rate constant for the water-pressure mode ( = 8 a^–), which is reflected in the short time-span over which the evolution is shown. Over this short period the waves in the ice and till move backward, while a typical feature of the pressure oscillation, the fact that it remains fixed in space, is seen. This is known as a “breather” mode. The change in effective pressure is “large” and indeed is about to float the ice.

Fig. 8. The coupled flow of ice, till and water showing the parameter study of bed-mode growth rates on the following parameters: Non-linear transmissibility exponent λ, till thickness D₀/m, ice thickness H₀/m, wavelength L/m, zeroeth-order-effective pressure p_e0/10⁵ Pa, zeroeth-order shear stress and τ₀/10⁵ Pa. The effective-pressure index b = 2.03 and shear-stress index α = 1. Case with evolving (i.e. prognostic) water pressure and bed, showing the proportion of positive bed-mode growth-rate constants binned against each parameter. Abcissa is the parameter value while ordinate is the proportion. Shown for all positive growth rates, and for growth rates greater than 10^-3, 10^-2 and 10^-1 a^-1. Dependence on hydraulic diffusivity F is shown in Figure 10.

Fig. 9. The coupled flow of ice, till and water, the same parameter study as in Figure 8, showing the proportion of positive pressure-mode growth-rate constants binned against each parameter. Abcissa is the parameter value while ordinate is the proportion. Shown for all positive growth rates and for growth rates greater than 10^-3 10^-1 and 10 a^-1. Dependence on hydraulic diffusivity F is shown in Figure 10.

Fig. 10. The coupled flow of ice, till and water, same parameter study as in Figures 8 and 9, showing the dependence, of growth rates on the bed hydraulic diffusivity F/m² a^-1. Cases of bed modes and pressure modes for cases where the bed is mobile (Q_i ≠ 0) and fixed (Q_i = 0) are shown. The legend applies to all graphs. Note flow distribution of fast unstable modes is different from distribution of unstable modes, being concentrated in areas of higher diffusivity.

Fig. 11. The coupled flow of ice, till and water, showing steady profiles arising as a consequence of a l m sinusiodal variation in bedrock topography, at indicated parameter values. Basal axes are position x and the parametric variable λ, the exponent representing the degree to which the transmissibility of the system changes with effective pressure. Vertical axes are ice thickness (tl), till thickness (tr), shear stress (ml), effective pressure (mr), water pressure (bl) and effective pressure + water pressure (br). This drainage parameter affects the ice-surface response.

Fig. 12. The coupled flow of ice, till and water showing the evolution of ice, till, water-pressure and effective-pressures profiles, for indicated parameters. Basal axes are position x and time t, vertical axes are ice thickness (tl), till thickness (tr), water pressure (bl) and effective pressure (br). Initial condition was a 0.1 m sine wave in the ice. This case has a very high growth-rate constant in the water-pressure mode ().

Fig. 13. Study of the coupled flow of ice, till and water showing the evolution of ice, till water-pressure and effective-pressure profiles, for indicated parameters. Basal axes ore position x and time t, vertical axes are ice thickness (tl), till thickness (tr), water pressure (bl) and effective pressure (br). Initial condition was a 50 Pa sine wave in the water pressure. Bed mode and pressure mode are nearly neutrally stable.

Figure 13 shows the evolution of another set of ice, till, water-pressure and effective-pressure profiles. Here the hydraulic-system diffusivity is F = 1 m² a^-1 and the initial condition was a 50 Pa sine wave in the water pressure. For this case, the bed mode and pressure mode are nearly neutrally stable, but note amplification of this nearly stable response to a very small initial perturbation, causing “large” changes in the effective pressure.

6. Conclusions

This paper has tried to address the issue of how ice, till and water flow couple, asking in particular whether the coupling of flow can explain the texture and variability seen in ice-stream flow. At the moment, there are two separate issues, whether it can and whether it does.

Whether the present theory can explain these phenomena is related to the technical areas explored by this paper. This paper has been concerned with the effect and consequences of scale upon effective-pressure distribution. On small scales, interracial effective pressures are strongly affected by static pressure variations, while on the longer scale there is significant coupling between hydraulic and static gradients, some of which arises as a result of the likely effective-pressure dependence of basal-drainage systems.

Variations in ice thickness cause variations in effective pressure; on the small scale, an ice mound causes an increase in effective pressure, while on the larger scale, this increased overburden is associated with an increase in water pressure which keeps the basal-drainage system transmissibility high. Similar effects are found with mounds in basal topography.

An analysis of the kinematic and shock-formation properties of deforming till is made. Kinematic waves can move up and downstream, and shocks facing upstream or downstream moving upstream or downstream can occur. The upstream-motion property of waves and shocks requires interfacial effective pressures to be significantly affected by static pressure variations. Thicknesses can occur where the velocity of a wave is equal to the shock speed for the jump between that thickness and zero thickness, leading to very persistent migratory forms.

Coupling between ice and till flow is examined at long wavelengths (i.e. greater than the ice is thick). Where interfacial pressures are statically controlled, small changes in ice thickness can significantly affect the effective pressure, and the effect is to reduce the ice-surface expression of basal variation in topography or frictional variations. Where interfacial pressures are controlled by drainage, small variations in till thickness can result in significant variations in ice-surface topography, not obviously related to basal topography. Length scales over which effective pressure is controlled statically may often be smaller than the ice-sheet thickness.

Coupled linearized dynamic analyses are performed. The coupled system is more or less stable when interfacial pressures are statically controlled, and shows instability in 50% of the phase space volume explored when interfacial effective pressures are maintained constant by the basal-drainage properties. There is always a fast stable mode, the Nye diffusion mode, and there is a slow mode, sometimes stable and sometimes unstable, which produces growing till and ice-surface relief. This slow mode is associated with a coupled till and ice wave moving downstream or upstream, but the volume of parameter space where there are unstable upstream-moving waves is very small or zero. In areas of parameter space occupied by ice streams from large ice sheets, the time constants are rather greater than what are believed to be characteristic occupation times for ice streams.

Introduction of coupled water drainage though a system with “linked-cavity properties” (i.e. transmissibility and storage decreasing with effective pressure) introduces, in certain regions of parameter space, unstable modes with time constants comparable with the stable Nye diffusion mode. This time constant seems usually to be associated with the pressure mode, depends on the (linearized) diffusivity of the hydraulic system, and can vary from being too slow to be of practical consequence to highly significant. The model produces unstable upstream- as well as downstream-moving waves in ice and till, although pressure oscillations appear to remain fixed in space (i.e. show “breather” solutions). For the most part, bed modes have long time-scales, although there are restricted areas in phase space where the growth rate becomes very fast (greater than 0.1 a^-1).

Clearly these instabilities are not sustainable. It is expected that shock formation will play a role in the quenching of bed-mode instabilities. The implication of this is that wave-like forms can grow and migrate. Other non-linear phenomena are likely to occur. For example, areas of high water pressure may connect to others, so reducing the water pressure and altering the dynamics. Kamb (1991) considered coupling to the heat equation, though Kamb’s results do not suggest that introducing thermal coupling will stabilise the system.

Many of the issues raised in this paper can be addressed by modelling the non-linear systems. Ensuring that such models represent the inherent instability of the system without introducing further numerical instabilities or removing real instabilities represents a major challenge.

Whether the mechanisms described in this paper do explain ice-stream texture and variability depend upon future observations. The coupling of ice, till and water flow produce a system with stability properties which depend strongly upon the location of the system within parameter space. One would strongly suspect that as the system moves through phase space, its stability properties will change, as is found in chaotic systems or systems in self-organised criticality. If this is the case, the question of prediction and model verification become statistical issues, with the need to observe and monitor whole flow fields. Theory can help us determine how much data we need; at the moment, theory tells us we do not have enough. Satellite data should be able in give us evolving topography and velocity fields and better definition of the bed and internal layers (which should help resolve whether variability is induced by theological or basal heterogeneity) are achievable. Observations of variations of deforming-till thickness of sufficient spatial resolution do not seem possible in the foreseeable future, while observations of basal-water pressure are straightforward but expensive.

A useful and achievable first step is the identification of the migration velocity of topographic features in the ice. Are they standing, moving upstream or moving downstream? Are they growing in amplitude or shrinking? Theory should be able to tell us whether the observed behaviour is due to internal or basal heterogeneity.