1. Introduction
Evidence of periodic deposition (Reference HeinrichHeinrich, 1988) of ice-rafted debris (Heinrich events) in North Atlantic deep-sea sediments suggests that the Laurentide ice sheet may have surged regularly with a time of order 104 a between surges (Reference Andrews and TedescoAndrews and Tedesco, 1992; Reference ClarkClark, 1994). The primary candidate for such surges is the ice stream discharging down Hudson Strait (Reference BondBond and others, 1992), although other ice streams may have been involved (Reference Bond and LottiBond and Lotti, 1995). Reference MacAyealMacAyeal (1993) suggested a mechanism whereby the ice sheet might behave in this way. Briefly, if the ice is thin, then it is cold-based, sluggish and so thickens; this causes the base to melt, and if the resultant water production is sufficient to allow fast sliding to occur, then a surge may result, drawing the ice down to a thin state once more.
Reference Fowler and JohnsonFowler and Johnson (1995) presented a physically based model of ice-sheet motion controlled by basal sliding, which itself depends on basal water production, They showed that the basal water produced by frictional heat at the base could lead to multiple ice-flux/ice-thickness relationships, and hence oscillatory behaviour by what they called hydraulic runaway, by analogy with thermal runaway that has been previously suggested as an instability mechanism for ice sheets (Reference Clarke, Nitsan and PatersonClarke and others, 1977). Fowler and Johnson were only able to suggest the possibility of surging in their model by means of a crude “lumped” version which essentially reduced the model to a zero-dimensional one. Although they later (Reference Fowler and JohnsonFowler and Johnson, 1996) offered some comments on the possible form of spatially extended surges (comments which we vindicate here), a proper numerical solution of the spatially extended model remains to be done. This is the purpose of the present paper, which also allows us to show for the first time how the magnitude and duration of the surge are controlled by the drainage-response time-scale, an observation that is by no means obvious.
At the outset, it is important to enunciate the philosophy of our approach. Our aim is to validate the idea that ice-sheet dynamics, coupled with basal sliding which incorporates a canal-type basal hydrological system, can generate surging behaviour. It is not our aim to produce the most realistic model incorporating as much realism as possible; nor even is it our intention to produce the best numerical approximation. In fact, to make the study as simple as we can, we intentionally strip the model of what we consider to be inessential detail. It is a matter of debate what constitutes inessential detail, but our approach is to strip the engine down and get it working, then add the trimmings later. This is an established modelling procedure.
Having said that, it turns oui that some unexpected trimmings do need to be added, and this dictates the progress of the paper, as follows. In section 2, we recapitulate the earlier model of ice-sheet flow and drainage due to Reference Fowler and JohnsonFowler and Johnson (1995); section 3 describes a numerical method to solve this model, but also reveals that the model has fundamental shortcomings, and that in fact it is not well posed in a mathematical sense. Hence, in section 4 we introduce realistic and essential modifications which allow a successful numerical method, described in section 5, to produce coherent results. As in the case of valley-glacier surges (Fowler, 1987b), the resulting surges are initiated by the passage of activation and deactivation waves, which propagate upstream and downstream as sharp fronts, and in sections 6 and 7 we provide a “shock structure” type of analysis to deduce the character and propagation speed of these fronts. A discussion follows in section 8.
2. The Fowler-Johnson Ice-Sheet Model
The model is described in more detail by Fowler and Johnson (1995, 1996) and is reviewed for completeness here. We consider a two-dimensional ice sheet on a flat base of thickness h(x, t), where x is a horizontal coordinate, and we specifically assume that 0 < x < l, with hx (= ∂h/∂x) = 0 at x = 0 (the divide) and h = hl > 0 at x — l. We have in mind the transition at x = l from ice sheet to ice shelf, and though a more realistic boundary condition ought to be prescribed, this is not essential to our purpose (indeed, we chose hl > 0 simply to avoid possible numerical difficulties associated with a singularity at x =l if h = 0).
The basic equation is that of mass conservation,
where subscripts denote partial derivatives, a is the accumulation rate and u is a depth-averaged horizontal velocity. Again, to be realistic we would allow for both shearing within t he ice and basal sliding, but the former is relatively small, so we begin by assuming that u is essentially a sliding velocity (at least while the base is at the melting point), and we pose a sliding law relating shear stress τ to velocity u of the form
(Reference Budd, Keage and BlundyBudd and others, 1979; Reference BindschadlerBindschadler, 1983; Reference FowlerFowler, 1987a). Reference KambKamb’s (1991) pseudo-plastic rheology would be accommodated by choosing s = 1, r 1. Here r and S are exponents which may be between zero and one, and c is a measurement of the stickiness of the bed. The basal shear stress is simply
where ρ is ice density and g is the gravitational acceleration. The extra variable N is the effective pressure (ice-overburden pressure minus water pressure), introduced into basal sliding laws by Reference LliboutryLliboutry (1968); it must be described by a basal drainage theory. It is Fundamental to this theory of surges that a drainage law of the type predicted by Walder and Fowler (1991} should be applicable. This was based on the physics of drainage over wet, deformable sediments, such as the detrital carbonates of Hudson Strait, and in its simplest form can be written as-
where c * is a constant relating to the deformable-till properties, and Q is the water flux (per channel) in a distributed system of canals which drain the bed. A number of simplifications are built into Equation (2.4), and we return to these later. In particular, Equation (2.4) is an equilibrium theory (steady drainage), and it cannot apply as Q0.
The precise form of Equation (2.4) is inessential; the important point is that N decreases with increasing Q, and this property would be shared by a drainage system consisting of a patchily distributed film (Reference AlleyAlley, 1989).
The final relation determines the water flux Q. For steady drainage, this is given by
where G is geothermal heat flux, τuis the frictional heat generated by the basal ice motion, q is the cooling rate to the cold ice above, ρ w is water density, L is latent heat and w d is the mean inter-canal spacing. Reference Fowler and JohnsonFowler and Johnson (1996) parameterised the cooling rate to the ice by means of a thermal-boundary-layer description:
where
-ΔT is the ice-surface temperature, c p is specific heat and k is thermal conductivity. The first of these is a heat-transfer term associated with rapid ice flow; the second is representative of conductive cooling. Adoption of Equation (2.6) avoids solving the two-dimensional temperature equation, again for simplicity: the form of Equation (2.6) is not of primary essence.
This parameterisation of the basal cooling allows the thermal effect of changes in ice thickness to be instantly communicated to the base. Reference Fowler and JohnsonJohnson (1995) examined the effect of including a delay in this cooling term (for the lumped (zero-dimensional) model), and found little qualitative effect. Indeed, in some runs this term was ignored entirely. We shall come back to a discussion of the parameterisation of basal cooling in sections 6 and 7.
Equations (2.1) (2.7) are those posed by Reference Fowler and JohnsonFowler and Johnson (1996) when the basal ice is at the melting point. They must be supplemented by appropriate relationships when the basal ice is cold. We will consider ibis situation further in subsequent sections.
Scaling the model
The model consists of Equations (2.1) (2.7) for the variables h, u, τ, N, Q, q and, with boundary conditions
and an initial condition for h. The model is then non-dimensionalised as described by Reference Fowler and JohnsonFowler and Johnson (1996) in terms of scales [h], [u], [τ], etc., and when the variables are written in terms of these scales, we find the dimensionless form of the model
where the dimensionless parameters γ, β, λ have typical values of order one. For example, we take r = S = 1 /2, and choose
where [a] is the accumulation-rate scale (a in Equation (2.9) is the dimensionless accumulation rate). The value of [N] corresponds to ice-stream conditions in Ice Stream B, Antarctica (Reference BentleyBentley, 1987), and if we choose c correspondingly, then we find
and
The definitions of the parameters arc extremely cumbersome, but it can be shown that they are related to the critical friction coefficient e and drainage coefficient c * via the combination
and then
In view of the lack of real constraint on σ, we use Equation (2.12) as a plausible guideline.
3. Numerical solution
Elimination of N and τ in Equation (2.9) enables the sliding law to be written in the form
where
and the modulus sign allows the application of Equation (2.9) even if hx > 0. We define
and then the equation for h takes the form
and is of non-linear diffusion type, with D = D(h, hx, Q), and Q is obtained by quadrature of Equation (2.9) 4.
We thus solve Equation (3.4) using an implicit diffusion solver. Specifically, we define a scheme as follows:
Here Δt is the time step, Δx is the space step and r a = Δt/(Δx)2. The current time level is indexed by j, and the current space location is indicated by i. The diffusion coefficients D ± are essentially Di ±½ and are specifically defined (for the relevant time step) by
This choice corresponds to type 2 of Reference Huybrechts and PayneHuybrechts and others (1996) and is thought to be less accurate than a type 1 (mass-conserving). For the present purpose, we believe it suffices, although in the sequel we do require very small space steps to ensure accuracy and stability in the solver for Q.
The relaxation parameter θ is zero for a fully implicit method, and equals 1/2 for a Crank Nicolson-type scheme which is, however, less stable (though more accurate). The spatial discretisation allows second-order accuracy, and in order to provide second-order accuracy in time, we define Di as D j-1/2 i , that is to say
Equation (3.5) (with Equation (3.7)) is not in itself second-order accurate in lime, but we expect that if the step is iterated (by analogy with the improved Euler method for ordinary differential equations), then it will give second-order accuracy, providing Equation (3.7) (and in fact Equation (3.8) below) is used. A problem in the implementation of Equation (3.5) is then in the non-linear diffusion coefficient D. To obtain second-order accuracy, we iterate the diffusion step by a simple predictor-corrector strategy, in which at the kth iterative step, we replace Equation (3.7) by
and D (0) i = D j-1/2 i . Thus the first pass corresponds to assigning Di at the preceding time step, with subsequent iterates Following each iterative step for h (k) i , we update Q by computing Q (k) i from Equation (2.9). With (only) second-order accuracy required, we use an improved Euler discretisation. We have also tried using a fourth-order Runge-Kutta solver. With the improved Q (k) i , we then compute h (k) i via central differences on h (k) i (and a three-point second-order, or five-point fourth-order, estimate for hx at x = 1). From this we have D (k) = D(h (k), h (k) x , Q (k)), and the iteration can proceed. The iteration terminates at a specified value of k, usually 2, and then h j i = h (k) i , Q j i = Q (k) i .
Results
Extensive numerical computations reveal that ice-sheet “surges” do occur, but these betray many signs of numerical unpredictability and bad resolution. With a good deal of effort, and with a particular choice of the parameters (particularly R and S), given by γ = 0.3, β = 1.2, λ = 0.75, and with Δt = 10-4, Δx = 10-3, R = 1.1, S = 0.4, we find that a front propagates backwards during a surge, as shown in Figure 1. Figures 2 and 3 show associated plots of u and Q. We see that very high spikes of u and a sharp jump of Q are associated with the propagation of this front. The jumps are in fact smooth for a sufficiently fine mesh, and they suggest the existence of a travelling wave. However, efforts to analyze the local structure of this travelling wave failed, and in fact a closer inspection reveals that the apparent wave speed V is greater than Δx/Δt, which indicates that the algorithm is not resolving the wave. Reduction to values Δx = 10-4, Δt = 10-6, only served to sharpen the profiles and raise the maximum values of u and Q in the front. The implication is that the problem as formulated does not have smooth solutions, and we might conjecture that (weak) solutions in which h is discontinuous and u has delta-function-like behaviour do exist, though this is a matter of conjecture.
4. Modifications to the model
We believe that the solution type represented in Figure 1 is realistic, but the model is unable to be solved satisfactorily, because of the absence of relaxational processes which allow for a switch in basal conditions. We believe that ice sheets can flow in either of two stable flow regimes (Reference Fowler and JohnsonFowler and Johnson, 1995), which depend on the amount of basal water present. At low Q, slow flow prevails, while for higher Q, fast flow occurs. The occurrence of multiple flow states can lead to self-sustained oscillations, as explained by Fowler (1987b), but in order to prevent discontinuities in the flow, it is necessary to include relaxation processes between the different regimes, Reference FowlerFowler (1989), in a model very similar to this, describing glacier surges, showed that inclusion of such terms allowed smooth solutions involving the rapid propagation of wave fronts in the surge, and we follow this suggestion here. There are other realistic glaciological effects which could be introduced to alleviate the apparent singularity of the Fowler-Johnson model. In particular, inclusion of longitudinal stresses is known (Reference FowlerFowler, 1989) to act as a smoothing term in the transition between fast and slow flow. We doubt that such damping effects will alleviate the problem here, which is fundamentally due to the implicit assumption of an instantaneous adjustment of the hydraulic pressure to equilibrium. In the lumped model of Reference Fowler and JohnsonFowler and Johnson (1995), this adjustment is the cause of the discontinuous changes in ice flux, and Figure 1 simply shows its spatial analogue.
Drainage relaxation
The description of time-dependent drainage is a complicated problem (Reference NyeNye, 1976; Reference Fowler and JohnsonFowler and Ng, 1996), and we choose to represent the response time of the drainage system by modifying Equation (2.9) 5 as follows:
The term in δ is analogous to that introduced by Reference FowlerFowler (1989), and is a simple representation of the relaxation of N towards equilibrium, δ is the ratio of the drainage-response lime (perhaps a month or less) to the ice-sheet-growth time-scale, so that values of δ~10-6 can be expected. The term is introduced so that when Q → 0, the pore pressure cannot go below zero (or N cannot go above the overburden ice pressure). In fact, it is convenient to choose (which is small) so that when Q = 0, the equilibrium value of N is such that the sliding law Equation (2.9) 2 describes realistically the small amount of shearing in the ice, which is described below.
Freezing base
We now seek to extend the model more realistically to accommodate what happens when Q = 0. There are two principal basal situations to consider: till-based and sediment-based. By till-based we mean that a layer of basal deformable till lies between the ice and an essentially impermeable bedrock; while sediment-based means that the (permeable) sediments below the ice extend to a significant depth.
For a sediment-based ice sheet, let s(x, t) denote the depth of the frost line below the ice-sediment interface. If the sediment porosity is ?, then the motion of the frost line is approximately governed by the Stefan condition
and when this is scaled as in section 2 (with. s ~ ice depth [h]), we obtain the dimensionless equation
where the Stefan number St and Peclet number defined by
and and κ w = k/ρ w c p. With previous estimates, and taking L = 3.3 ? 105 J kg-1, κ w = 38 m2 a-1, ΔT = 50 K as the surface sub-cooling below freezing and c p = 2 ? 103 J kg-1 K-1, we have St ~ 3.4, Pe ~ 7.9. With ? ≈ 0.4, we then have ?StPe~10.7, which implies that (for times of O(1)) s is small relative to ice depth, and the conductive temperature profile in the frozen sediment implies the following (scaled) thermal relation for the basal ice temperature T 0:
where q is the dimensionless heat flux, which is now taken to be given by We discuss the theory for a till-based ice sheet further below.
The above description is by no means watertight, but represents the simplest coherent description of basal freezing, along the same lines as the description of the basal cooling q. The fact that s is small is what motivates the assumption in Equation (4.2) that the heat flux at the freezing front is the same as that at the ice-sheet base, and (equivalently) that the temperature profile in the sediments is linear (Equation (4.5)); this is due to the small aspect ratio of the frozen sediments. Since Equation (4.3) then implies that sub-basal freezing and thawing occur on the (slow) convective time-scale, in fact the description here is likely to be reasonable. The main shortcoming might lie in the fact that q changes rapidly during a surge. However, since it is the heal flux rather than the basal temperature which changes, Equation (4.3) apparently does not imply a consequent rapid change in s.
Shear flow in the ice
When the ice is cold at the bed, we may take the sliding velocity to be zero, but there is a small component of ice velocity due to shearing within the ice. Reference FowlerFowler (1992) showed that the shear is concentrated in a shear layer near the base which gives an effective plug-flow velocity in the ice of (with the present definition of the scales)
where is the activation exponent defined by Fowler (1992; denoted there as γ), = QΔT/RT 22, and
where A is the flow exponent. Taking A = 0.2 bar-3 a-1 (Reference PatersonPaterson, 1994), n = 3, [τ] = 0.2 bar, [h] = 1500 m and [u] = 250 m a-1, we have ~ 10-2, so that even if T 0 ≈ 0, the shearing velocity is very small (note also that ≈ 11).
We define
so that Equations (4.5), (4.6) and (4.7) give, approximately (since s h and u/ ½ 1 in Equation (4.6)),
while Equation (4.3) is (since u 1)
The term τu is (inessentially) introduced in Equation (4.11) purely to facilitate the numerical solution (see Equations (4.13)). When Σ > 0, the base is frozen, and u is given by Equation (4.10); so long as λ * τ n+1 is small, then τu is small and the term makes little difference. This is no longer true if λ * τ n+1~1 (i.e. shear stresses of order 1 bar), which may be approached near the front or (briefly) during surge activation. We have no reason to believe that its inclusion has any significant effect, but we have not explored the point.
Equations (4.10) and (4.11) apply when Q = 0. The term in ν is retained in Equation (4.11) because if h reaches a steady stale with Σ > 0, then Equation (4.11) indicates that s reaches an O(1) steady stale over the longer time-scale of O(1/ν). For the oscillatory solutions we find below, the term is not important.
In solving the model, it is convenient to choose so that the equilibrium value of N as computed (-1/3 when Q = 0) causes the sliding law to give the shearing velocity (Equation (4.10)). We find that this is the case if we choose
which we promptly do.
Water-flux adjustment
We may also modify the equation for Q by including a term δ′Qt on the lefthand side of Equation (2.9) 4. In this case, a stable numerical method is to define the time step as Δt = δ′ Δx, so that integration from (i, j) to (i+1, j+1) is along the characteristics, and is implemented with an improved Euler method.
Summary
The model we now solve is the following:
with given by Equation (4.12). This reduces to Equation (2.9) if δ,δ′,λ * and ν are put to zero, i.e. if we revert to the assumption of equilibrium drainage hydraulic pressure and water flux, and lake zero ice shearing and negligible frost penetration.
5. Numerical Method
The model is solved for h as before. As mentioned, Q is solved along characteristics by choosing Δt = δ′ Δx; both N and Q are solved with the improved Euler method, and we also time-step Σ forward using an improved Euler method from (i, j — 1) to (i, j). The treatment of cold temperate transition points requires a good deal of care, and here we describe possibly the simplest way of dealing with this point; more sophisticated methods could certainly be devised.
The basic problem is illustrated in Figure 4. Each of the cases 1Aa, 1Ab, 1B, 2A, 2Ba and 2Bb represents two successive time steps j — 1 and j, with the spatial grid i indicated horizontally. Solving for Σ or Q is straightforward in Equations (4.13) except at points x CT which define cold temperate transition points (that is, Σ > 0 on one side of x CT and Q > 0 on the other). At each point of the ice-sheet base, we need to solve either for Σ (frozen base) or for Q (warm base), and the difficulty arises because we integrate Q from (i, j -1) to (i + 1,j), but Σ from (i, j-1) to (i, j).
At time step j — 1, there are two main cases to consider:
Case 1 corresponds to frozen bed upstream, case 2 to frozen bed downstream.
Each of these has two separate subclasses:
Sub-case A corresponds to thawing conditions (x CT) moving towards the frozen bed), and sub-ease B corresponds to freezing conditions (x CT) moving towards the warm bed).
In Figure 4, the location of x CT at time levels j - 1 and j is indicated by open circles. To one side, Σ > 0 and Σ can be integrated to step j (at the same value of i). On the other, Q > 0 and it can be integrated to step j at the next value of i. Thus, for the various cases, we may sometimes be able to
compute values of both Σ and Q at some gridpoints at step j, or sometimes neither.
For example, in case 1A, with f > 0, then Σ j-1 i decreases; if Σ j i > 0 (case 1Aa), then x CT remains in the interval (x i, x i+1), and the value of Q j i+1 is computed using a weighted value of the step length Δx. Since this corresponds to backwards propagation of x CT (f > 0, so melting is occurring), and since (see section 6) we expect that the speed of x CT is much less than 1/δ′, we can assume that x CT traces an approximately vertical path in Figure 4 (case 1Aa), and also that Qx ≈ f, whence the distance xi+1 — x CT is approximately Q j-1 i+1. In computing Qj i+1 we replace Δx by this approximation (for case 1Aa). Case 1Ab is similar, except that Σ j i+1 < 0, and x CT crosses into (x i-1, xi ). In this case, xi becomes warm, and we assign the value Qj i ≈ 0. If, on the other hand, f < 0 (case 1B), then x CT must propagate forward at a rate larger than 1/δ′, and this is illustrated in Figure 4; in fact, from the definition of X CT (Σ = Q = 0 on x = x CT, we find that Σx~δ′ and
As indicated in Figure 4 for case 1B, one or more nodes change from Q(> 0) to Σ(> 0). A simple algorithm is to compute Qj k (if Q j-1 i-1 > 0), and if Qj k < 0, then for such k, we assign a value Σ j k ? 0. In practice, we choose Σ j k = 0 which should be accurate for these interpolated nodes, since Σ x ~ δ′.
Cases 1A and 1B represent the two important cases that occur in practice. For completeness we describe two further cases. Case 2A occurs when Q j-1 i > 0, Σ j-1 i+1 > 0 and f j-1 i+1 > 0. In this case, x CT represents a shock, insofar as Q and Σ suffer discontinuities at x CT. A conservation argument then implies that
where Q is evaluated on x CT-, and Σ on X CT+, and we represent this algorithmically in a crude way by selecting the computed value of Qj i+1 if Q j-1 i /Σ j-1 i+1 > 1/δ′, and Σ i i+1 if Q j-1 i /Σ j-1 i+1 < 1/δ′. This has the effect of pushing x CT forwards to the next grid space if CT > 1δ′ = Δx/Δt. Again, a more sophisticated method would be more accurate, but this is sufficient here, where this situation does not in fact arise.
Lastly, case 2B is illustrated in the two cases 2Ba and 2Bb. If the computed Qj i+1 > 0 (case 2Ba), then a shock is propagating, and is dealt with as in case 2A. If Qj i+1 < 0 (case 2Bb), then the front is no longer a shock and may propagate backwards. In either of the illustrated examples, the algorithm simply retains Σ j k (= 0) if Qj k < 0 (as for case 1B).
Results
In the following figures we show the course of a surge using a value of δ = δ′ = 10-2, with space step Δx 4 ? 10-5 and time step Δt = 4 ? 10-7. Consistent results are obtained with coarser resolution, though in certain regions there is then some numerical instability. A remnant of this is in fact visible in Figure 6. Although we have supposed realistic values of δ and δ′ to be of order 10-6, we can anticipate from the form of the equations that such small values will lead to an extremely stiff system, and require inaccessibly small lime steps to capture the anticipated rapid transitions in N and Q. This expectation is validated in sections 6 and 7, where we are able to obtain analytic estimates for the fast time-scales involved, and we then use these to extrapolate the numerical results to more realistic conditions. Figures 5 and 6 show the evolution of the central thickness h max = h(0, t) and the outlet velocity u out = u(1, f) as functions of time. We see a pattern of surging at regular intervals, with a period of the order of the convective time-scale. The surge is initiated by an activation wave which propagates backwards from the margin at a rate controlled by the drainage-pressure relaxation time (δ), and which enables the transition from cold-based to temperate-based dynamics. The propagation of this wave is shown in Figures 7-9.
When the activation wave reaches the divide, a deactivation wave is formed which propagates forwards (at a faster speed which is controlled by the water-flux relaxation time(δ′)) causing the base la freeze again. The outlet ice flux at the margin is shut down very rapidly when this wave travels towards the front, and the slow build-up begins again. The propagation of this deactivation wave is shown in Figures 10-13. The propagation speeds of the activation and deactivation waves are controlled by the relaxation parameters δ and δ′. In sections 6 and 7, we analyze the two fronts, and derive expressions for the propagation speeds, and the form of the variables in the fronts, by methods analogous to those used in studying shock structure. An interpretation of the numerical results in the light of these analytic results is then given in section 8.
6. Surge Activation
During a surge, a wave front propagates backwards at a rate controlled by the water-pressure response time δ. Suppose the wave illustrated in Figures 7-9 is located at x = x s(t). We define local shock variables by rescaling as follows:
where V ~ O (1) is to be determined. The rescaled equations to determine the variables within the activation front become, to leading order (assuming ε 1),
providing varepsilon; 1, and where
and we have chosen
which is consistent with, δ′ 1 if, for example, δ = δ′ 1. Note that the condition δ′ underlies the algorithm used in case 1 in section 5. From Equations (6.1) and (6.4), we have
whence the wave speed is of order
If we choose r = s = 1/2, λ * ~ 10-3, then ~ 10-1 δ -2/3, ε ~ 10-1 δ 1/3. Thus for δ ~ 10-2, the predicted shock width is ~0.02 and speed is ~2, while for the more realistic δ ~ 10-6, we have ~ 103, ε ~ 10-3. From Figure 7, the apparent front speed is about 6, and the front width is about 0.4. This suggests that V ≈ 3, so that the wave speed is about 3, and the shock width is about 20ε. The parameter Λ in Equation (6.3) is defined by
e.g. Λ ~ 109 for δ = 10-2, Λ ~ 1012 for δ = 10-6.
We now seek solutions to Equation (6.2) in which Q > 0 at X > 0, and Σ > 0 at X < 0, so that the cold-temperate transition is at X = 0. In X > 0 (Q > 0), we neglect O(Λ-1/3), so that (if h = h 0 and N → N 0 at X = 0)
(since u 1 at X = 0), and therefore
The depth profile in X > 0 is therefore determined by
and the (scaled) shock speed V is then determined by
where h + is the depth ahead of the front (h + < h 0).
A crude estimate for V is obtained by taking h, h + h 0 in Equation (6.11), so that
If we take h 0 ≈ 2.5, N 0 ≈ 0.2 from Figures 7 and 9 (note that N has been rescaled with λ *-r/s in Figure 7, just as in Equation (6.1)), and r = s = 1/2, then this gives V ≈ 3.15, consistent with the propagation rate in Figure 7. Furthermore, Figure 9 suggests an e-folding decay distance of about (using a ruler on the t = 0.86 profile in Figure 9) 0.077. This is predicted to be X = V, or x – x s = Vε . With ε = 0.022 from Equation (6.5), this would give V ≈ 3.6, again consistent with the propagation rate. These results are added confirmation of the validity of the analysis.
This gives the basic wave profile ahead of the cold-temperate transition point, although the approximation involved in deriving it breaks down as X → 0 (as then ΛQ → 0). More specifically, we see that
as X → 0, and we select the region near X = 0 where Q ~ 1/Λ by writing
To leading order, Equations (6.2) become
In terms of x, this region is of thickness
For r = s = 1/2, the exponent of λ * is one, so that ε * ~ δ 1/2 λ *, and with λ * = 10-3, ε * ~ 10-4 for δ = 10-2, ε = 10-6 for δ = 10-6. Even for the higher values of δ used numerically, this region is still practically at the limit of resolution.
Within this region,
whence
where 0 is an integration constant, and thus
and M is determined from Equation (6.14). By choosing N 0 = Q *-1/3, we can see that M will be zero at X = 0, and thus N is smooth. This upper portion of the front must be joined smoothly to the upstream part where Σ > 0. In fact, we see that in this upper portion,
and if r = s = 1/2, this is (λ *2/δ) r/(r+1) . This is small for λ * = 10-3, and even δ = 10-6, and even for δ = 10-6 it is O(1). This suggests that the portion of the front where Σ > 0 is of little significance, and this is borne out by the numerical results.
7. Surge Deactivation
When the activation wave reaches the divide, a deactivation wave propagates rapidly downstream, thus switching the surge off. We can see from Figures 10-13 that this deactivation wave is associated with the propagation of the cold-temperate transition point x CT downstream. We have seen in section 4 (Equation (4.14)) that this front moves with the speed
When f > 0, the x CT point represents a “warming” front, and the wave speed is less than the characteristic speed, while for f < 0 it is a “cooling” front, and moves faster than the characteristic speed. In the latter case, the front can be considered to be a weak shock. In any event, there is no shock structure associated with this front propagation (at least in this model), and from the numerical solutions, and also Equations (4.13), we see that h experiences a jump of O(δ) across the front, and N decreases from O(-1/3) to O(1). Without any smoothing mechanism for Q, the solutions for N (and thus also h and u) exhibit numerical judders near X CT on a small scale, though these are not apparent in the figures. It can be seen from Figure 11 that for δ′ = 10-2, the speed of the front (between t = 0.906 and t = 0.909) is v ≈ 1.7 ? 102 = 1.7/δ′, consistent with this result.
8. Discussion
The model of ice sheets incorporating canal-type drainage described by Reference Fowler and JohnsonFowler and Johnson (1995) exhibits surges on a regular periodic basis, with a period between surges controlled by the accumulation rate. When the ice is thin and cold, then h thickens towards its equilibrium hu = ax. In the present model, if Σ 1, then u ≈ λ * τn , so λ * h(-hhx )n → ax, and the equilibrium profile is given by
where h m is the marginal value of h.
So long as the marginal ice remains cold, then Σ will tend to an equilibrium (~ 1/ν), so long as f = 0 for some Σ > 0. By consideration of Equation (4.13) 5, we can see that the condition for this not to occur is that
Since h (and thus u and thus ) is given at x = 1, this is a messy but explicit criterion which indicates melting will be initiated at the front if γ is large enough. For O(1) values of γ, λ, and β, Relation (8.2) will inevitably be satisfied since λ * is small.
Physically, the condition to maintain a freezing base requires (assuming u is then very small) that the geothermal heat flux γ can be removed by conductive cooling, which in a steady state is just λ/(h + vΣ). If then γ > λ/h, then this is not possible for a positive frost depth Σ, and the criterion for the initiation of melting is simply Relation (8.2), or (with small u) h > λ/γ.
When melting is initialed, an activation wave propagates backwards at a dimensionless speed of order 3 = 3λ *r/r+1 δ -1/r+1. The surge is active during the passage of this wave, which therefore lasts a time t surge ~ 0.3δ 1/r +1 λ *-r/r+1, and the ice velocity (and thus the discharge flux at the margin) are of order 1/. When the melting front reaches the divide, the dome collapses, and a freezing front deactivates the surge on a time-scale 1/δ′. The timing and magnitude of the surge is thus critically dependent on the relaxât ion times for the basal drainage system.
Dimensional results
The basic units used by Reference Fowler and JohnsonFowler and Johnson (1995) are given by Equation (2.11). We identify the connective time-scale [t], the ice-sheet depth scale [h] and the ice velocity scale [u]. The drainage time-scales can be identified in terms of two processes, pressure relaxation and water-flux relaxation. If we compare Equation (4.13) with, for example, Reference NyeNye’s (1976) hydrological drainage model, then we can identify a convective water-flux time-scale t flux = l/ where is a typical water velocity in the drainage system; then
The pressure-relaxation time-scale is less obvious. In terms of Nye’s hydrological model, it is given by the closure time-scale of ice (and would thus be {A[N]n}-1 where A is the flow-law constant and [N] is the effective-pressure scale). For Walder and Reference Walder and FowlerFowler’s (1994) canal drainage system, the relevant closure time-scale t clos might be given by
It can also be shown (Reference NgNg, 1998) for the Nye model that
where Δ? is the driving hydraulic head, and L is latent heat. If we suppose Equations (8.4) and (8.5) are valid, then with values L ~ 3.3 ? 105 J kg-1, [N] ~ 1 bar and η till ~ 1010 Pa s (Reference Alley, Blankenship, Bentley and RooneyAlley and others, 1986), we would have
and with [t] = 8000 a = 2.5 ? 1011s, we would deduce
These are merely indicative. In particular, it is likely that t clos will be longer than 1 d, if account is taken of the time taken for lateral water transport between flow pathways. That is to say, Equation (8.4) gives a response time for pressure adjustment at a canal, whereas what is of principal concern is the response time of the effective pressure in till away from the supposed drainage canals. In what follows, we suppose always t flux/t clos ≈ 0.06, and give estimates for quantities of interest for various values of t clos.
The surge period is of order [t] ~ 8000 a, entirely consistent with Heinrich-event intervals of 5000-10000 a. Other quantities of interest are:
surging ice velocity, ≈ 7.4 [u];
surge duration, ≈ 0.3[t]/;
deactivation wave speed, ≈ [u]/δ′;
peak surge basal water flux, ≈ [Q]/w d.
According to Reference Fowler and JohnsonFowler and Johnson (1996), [Q] is the water-flux scale per channel, [Q] ~ 2 m3 s-1, while w d is the channel spacing, w d ~ 4 km. This last value is badly constrained, but the quantity [Q] is the water flux per unit width, which does not depend on w d (since in fact [Q] ? w d.) is defined by Equation (6.6), and in our numerical results r = 1/2, δ = 10-2, λ * = 10-3), ≈ 2.15, and this together with the results in Figure 6 (velocity readies 16, duration is 0.15), Figure 8 (maximum Q is about 2) and Figure 10 (deactivation time is 0.01) determines approximately the numerical coefficients above. In Table 1 we show values of these quantities for various choices of t clos (assuming t flux/t clos = 0.06), and using [t] = 8000 a, [u] = 250 m a-1 and [h] = 1500 m, and with λ * (given by Equation (4.9) 4) = 10-3, r = s = 1/2.
The violence of the surge is sensitively dependent on the pressure-relaxation time-scale, ranging from benign if t clos = 30 a to catastrophic if t clos = 1 d. Note, however, that the total ice discharge (which is proportional to the product of velocity and duration) is independent of t clos. For an outlet channel of width 200 km, and an exit depth of 500 m (as used here), the net ice discharge is 4 ? 105 km3.
The activation wave is controlled by the closure time-scale. As the cold--temperate transition point x CT moves back, one can think of the rate of propagation being determined by the rate at which the drainage system is set up, and this is controlled by the balance between (for a canal) till creep and sediment removal. On the other hand, the de-activation wave is controlled by the convective water-flux time-scale, and is effectively instantaneous by comparison.
Even for the most vicious surges, the water flux is not too enormous. The most violent jökulhlaups in Iceland have maximum peak discharges of 105 m3 s-1 km-1 (Reference TómassonTómasson, 1996), although only for several hours. The regular Grímsvötn jökulhlaups have discharges of order 102 m3 s-1 km-1 over several days (Reference BjörnssonBjörnsson, 1992). Total discharges are thus less than about 109 m3 km-1. For the ice-sheet surges, however, the total discharges are larger. Independently of the surge duration, Table 1 suggests a total discharge of 3.8 ? 1010 m3 km-1, or (for a 200 km wide outlet ice stream) 7.6 ? 103 km3, about 1% of the ice discharge. For a drawdown area of 106 km2, this represents melting of 7.6 m of ice.
Apart from t clos the other imponderable in Table 1 is the value of r in the sliding law τ = cur Ns . The values r = s = 1/2 correspond to the view of Reference Boulton and HindmarshBoulton and Hindmarsh (1987) and Reference BoultonBoulton (1996) that coarse (sandy) tills have values of r, s of O(1). The view of Reference KambKamb (1991) is that for clay-rich marine sediments, the sliding law is more nearly plastic, so that s ≈ 1, r 1. We can see the effect of choosing r 1 in Table 1. The values of in the three columns are then = 0.25 ? 107, 0.25 ? 103 and 0.25 ? 103. Even for the mildest case t clos = 30 a, = 250, the surge velocity is 460 km a-1, the duration is 10 years and the peak water discharge is 120 m3 s-1 km-1.
Variable periodicity
The Heinrich-event intervals are not constant. This is hardly surprising since the intervals are controlled by accumulation, which must be significantly affected by the climatic alterations induced by these surges.
Synchronised surging
Reference Bond and LottiBond and Lotti (1995) indicate evidence that different ice streams surged simultaneously, which would apparently cause difficulty for the notion that their oscillations are controlled internally and independently. This point of view has been supported by Reference Revel, Sinko and GroussetRevel and others (1996), while Reference Gwiazda, Hemming and BroeckerGwiazda and others (1996) indicate that, at least for Heinrich event 2, the Hudson Strait origin is dominant. The idea that ice-sheet surges can be caused (as opposed to influenced) by external forcings is, however, dynamically naive. In fact, it is well known that even weakly coupled non-linear oscillators will readily become synchronised, and there are many examples where this phenomenon can be recognised, despise the apparent minuteness of the coupling: firefly illuminations, human menstrual cycles, and mechanical clocks come to mind. It would in fact be almost surprising if the surging of different ice streams did not tend to become synchronised, or at least harmonically resonant, as they are coupled through the climatic input. A similar comment applies to the minute 100000 year Milankovitch forcing of the ice ages.
Acknowledgements
We wish to thank the European Science Foundation which, through its initiatives on free boundary problems i Mathematical analysis of free boundary problems) and the European Ice-Sheet Modelling Initiative (EISMINT), provided grants to facilitate the visits of E.S. to Oxford.
