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 Reference Walder and FowlerWalder 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 (Reference AlleyAlley, 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 (Reference PatersonPaterson, 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 Reference AlleyAlley (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. Reference Boulton, Menzies and RoseBoulton 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] (Reference Boulton, Menzies and RoseBoulton and Hindmarsh, 1987; Reference Walder and FowlerWalder 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. Reference AlleyAlley, 1993, 1996; Reference Walder and FowlerWalder 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 Reference Boulton, Menzies and RoseBoulton 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 (Reference BearBear, 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., Reference BearBear, 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 Reference AlleyAlley (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. Reference AlleyAlley, 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. Reference BearBear, 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.1. Internal deformation flux relationship

The flux contribution arising from internal deformation can be computed from a postulated viscous relationship for till (Reference Boulton, Menzies and RoseBoulton 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 (Reference Cuffey and AlleyCuffey and Alley, 1996) and a viscous one (Reference HindmarshHindmarsh, 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 (Reference LaxLax, 1973; Reference WhithamWhitham, 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.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.

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; Reference MorlandMorland, 1984; Reference FowlerFowler, 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 Reference BuddBudd (1970) and Reference JóhannessonJó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 (Reference Bindschadler, Vornberger, Blankenship, Scambos and JacobelBindschadler 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

Reference Smalley and UnwinSmalley 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 Reference AlleyAlley (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.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 (Reference NyeNye, 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 (Reference HindmarshHindmarsh, 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 Reference Bindschadler, Vornberger, Blankenship, Scambos and JacobelBindschadler 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.