2.1. Governing equations for the hybrid description of ice flow

To simplify the Stokes equation, we introduce a parameter, ε = [d]/[l], where [d] and [l] are the vertical and horizontal length scales associated with an ice sheet, with typical values provided in Table 1. Since the lateral extent of an ice sheet is much greater than its thickness (e.g. Reference FowlerFowler, 1997), ε = [d]/[l] is taken to be very small, as in lubrication theory.

Table 1. List of notation with values of constants and scales. Typical ice-sheet scales provide a natural choice of depth and horizontal length and velocity scales. Remaining scales are chosen by balancing terms in the governing equations

For equilibrium at the ice-sheet bed, the basal drag, τ _b, resists the local driving stress, as assumed in the SIA. If, however, there is basal sliding, it is necessary to consider additional stress gradients in the model. In light of this we introduce a new dimensionless parameter, λ, which is the ratio of the magnitude of the basal stress to the viscous shear stress scale, providing an estimate of slip,

When the flow is shear-dominated (basal stresses are of similar magnitude to the vertical shear stresses) while for λ ≪ 1, basal stresses are small compared with vertical shear, and the flow is sliding-dominated. λ is never >1 since τ _b < η[u]/[d] at all times, and equal if, and only if, the ice sheet is frozen to the bed and no slip is occurring (i.e. if the static basal stress is locally in balance with the driving stress, as in the SIA).

This parameter is similar to the traction number of Reference Schoof and HindmarshSchoof and Hindmarsh (2010) and the slip ratio of Reference GudmundssonGudmundsson (2003). More specifically, if we define a deformational velocity scale due to vertical shear, [u _d] ∼ A[τ _b] ⁿ [d] (Reference Schoof and HindmarshSchoof and Hindmarsh, 2010), then our parameter, λ, can be written as λ = [u _d]/[u], which is the inverse of the slip ratio used by Reference GudmundssonGudmundsson (2003) and is closely related to the stress ratio used by Reference Schoof and HindmarshSchoof and Hindmarsh (2010). An advantage of our formulation is that 0 < λ ≤ 1, making it easier to consider the whole range of sliding regimes, as is the purpose of using a higher-order formulation.

With ε and λ as the two control parameters, we derive a set of non-dimensional equations for h, u _b and v _b, as outlined in Appendix A (Reference FowlerFowler, 2011). We make the thin-film assumption, meaning that we neglect the horizontal gradients of the vertical shear stresses. We assume that there are no abrupt changes in basal or surface elevation, as is consistent with observations of pure ice streams. Moreover, membrane stresses are accurately approximated by using the sliding velocity to compute them; if the horizontal velocities, (u, v), change rapidly in space, then this must be due to rapid changes in sliding velocity, (u _b, v _b). These assumptions are explained more comprehensively in Appendix A.

The governing equation for conservation of mass is

(1)

where a is the accumulation rate,

(2)

and S is the resistive stress tensor (Reference Van der VeenVan der Veen, 1999; Reference HindmarshHindmarsh, 2012), given by

(3)

Note that in the specific case where , γ can be approximated as

(4)

in Eqn (1), as is evident from an expansion of |γ| ^n–1 retaining terms greater than .

Equation (1) is a depth-integrated equation for conservation of mass. It shows that a change in depth of ice at some point is caused by divergence of the ice flux and the accumulation rate at that point. If λ → 0 then ice flux approximates to h u _b (plug flow) whereas if λ = 1 the flux corresponds to that for free-surface lubrication flow (as in the SIA).

Conservation of momentum is given by

(5)

which describes the force balance at the bed. The three terms correspond to basal stress, driving stress and membrane stress contributions, respectively. For , the membrane stress contribution is small and the balance is largely between the basal and driving stress; for λ ≪ 1 the membrane stress correction plays a more significant role in the stress balance.

2.2. Constitutive sliding law

Sliding laws of the form τ_b = f(u _b) are long-established. However, because of the strong dependence of basal friction on the effective pressure at the bed, sliding laws of the form τ _b = f(u _b, N), where N is the effective pressure at the bed, provide a better description. In the current framework we do not model the basal water pressure; hence we use the first kind, although some dependence on effective pressure is implicitly included in our choice of f.

More specifically, in the first part of this work we consider a triple-valued sliding law. Reference Sayag and TzipermanSayag and Tziperman (2009, Reference Sayag and Tziperman2011) developed the ideas behind such a sliding law and we adopt their implementation, which is obtained analytically from a simple one-dimensional (1-D) force balance, based on the cross-flow structure of an ice stream. Figure 1 illustrates the resulting non-dimensional relationship, given by

Fig. 1. A cross section of the triple-valued sliding law, τ _b(u, v = 0), based on figure 4 of Reference Sayag and TzipermanSayag and Tziperman (2009). The dashed line is the corresponding driving stress (τ_d ).

(6)

where α and β are parameters taken to be −0.9 and 50, respectively.

A multivalued sliding law is expected to be unstable on its middle branch where ∂τ _b/∂u _b < 0 (Reference LliboutryLliboutry, 1968, p.55). Consequently, if the stress is gradually increased while on the lower branch to the first maximum in Figure 1, a transition must then occur to the upper branch, and this transient must involve the rapid variation of some other dynamical quantity that is implicitly involved in the sliding law. Commonly, this hidden variable is taken to be the subglacial effective pressure, which is also presumed to be the underlying cause of the multivaluedness (Reference FowlerFowler, 1987b). We therefore introduce a new variable, v, such that

(7)

where δ is a parameter that describes the timescale on which v relaxes to u _b, i.e. a timescale for the evolution of the subglacial drainage system. This timescale introduces a limit to how rapidly velocities can transition between u ₋ and u ₊.

By replacing u _b with v in the part of the triple-valued sliding law that governs the positioning of the unstable branch (the cubic polynomial), we obtain a relationship between basal stress and velocity with some time dependence in it, defining a timescale over which velocity changes in response to changes in basal stress. In particular, taking δ = 10⁻³ gives a timescale of evolution for the system of 0.4 years, given the other scales in Table 1.

The triple-valued sliding law hence becomes

(8)

We propose that sliding laws of the form τ _b = f(u _b) which express a physically unstable relationship between τ _b and u _b should be replaced with those of the form τ _b = f(u _b, v) .

2.3. Model set-up and numerical solution

Figure 2 shows the 3-D model set-up. There is zero inflow into the domain at x = 0, implemented by setting u _b = 0 and ∂s/∂x = 0; the latter imposes zero flow due to vertical shear at x = 0. Free slip is allowed along the boundary (τ ₁₂ = 0). To avoid considering grounding-line behaviour, we treat x = 1 as an open outflow boundary (Reference Papanastasiou, Malamataris and EllwoodPapanastasiou and others, 1992; Reference Sayag and TzipermanSayag and Tziperman, 2009), where the upper surface is free to evolve and ∂u/∂x = 0. The domain is periodic in the y-direction.

Fig. 2. Schematic illustrating the model set-up. The bed is flat and of constant slope. The surface is also flat initially, before a constant perturbation to the accumulation field is applied in the form of a Gaussian blob, a _p(x, y), as illustrated. There is zero inflow at the x = 0 boundary.

The simulations are initialized in a steady state that is uniform in the cross-flow direction. Downstream of the inflow boundary the ice has a constant thickness, h = 1. There is no variation in velocity in the cross-flow direction, i.e. (u, v) = (u(x), 0). Under these conditions, the velocity field can be computed from Eqns (5) and (8), with v = |u_b| in the steady state. A background ‘balance’ accumulation, a _b(x), is required to maintain the steady state. This balance accumulation is the steady-state solution of Eqn (1), given by

(9)

A forcing is required to produce ice-stream behaviour. A perturbation accumulation, a _p(x, y), is therefore also applied, in the form of a truncated Gaussian at the upper boundary (motivated by Reference Sayag and TzipermanSayag and Tziperman, 2009). More specifically, a _p(x, y) is given by

(10)

where A ₀ is the magnitude of the perturbation. The magnitudes used are necessarily large, but can be thought of as representing the sum of accumulation rates from upstream. Use of such a forcing provides a simple way of reproducing ice-stream behaviour, which allows us to study the resulting downstream flow regimes.

We discretize the governing equations and boundary conditions using a finite-difference scheme with semi-implicit time-stepping. We use a staggered, Cartesian, two-dimensional (2-D) grid with n_i × n_j discrete points (Reference Rommelaere and RitzRommelaere and Ritz, 1996; Reference Sayag and TzipermanSayag and Tziperman, 2009). An upwind difference scheme is implemented at the open outflow boundary, together with the condition ∂u/∂x = 0 (Reference Sayag and TzipermanSayag and Tziperman, 2009). The nonlinear system of equations is solved using a Newton–Krylov solver, provided by the Portable Extensible Toolkit for Scientific Computation (Reference Balay, Gropp, McInnes, Smith, Arge, Bruaset and LangtatenBalay and others, 1997, Reference Balay2012a,Reference Balayb). The fully coupled system is inverted at each time-step, and the solution is accepted when the absolute size of the nonlinear system residual is <10⁻⁸. To aid numerical solution, if this condition is not satisfied after 50 iterations we split the equations into two groups (Eqns (1) and (7) in one, and the components of Eqn (5) in the other) and solve them iteratively using a Picard iteration with a subspace correction (Reference Hindmarsh and PayneHindmarsh and Payne, 1996). The code is tested against two analytic solutions, as outlined in Appendix B. The code converges to the steady-state solutions rapidly and the accuracy of the agreement is good.

Values of constants are given in Table 1. The value of the maximum time-step must always be smaller than the relaxation timescale defined by b so that changes in velocity are not limited by the size of the time-step.

An alternative set-up, together with some basic results, is discussed in Appendix C.

3.1. Reference flow regimes at small λ

In order to assess the suitability of the model for describing ice dynamics in and outside ice streams, we apply the triple-valued sliding law as the basal boundary condition in the set-up described in Section 2.3. The ice rheology is taken to be Newtonian, to allow us to compare results with Reference Sayag and TzipermanSayag and Tziperman (2009).

Figure 3 shows the three distinct flow regimes that occur as a result of varying the magnitude of the perturbation accumulation, A ₀. Other parameter values are given in Table 1.

Fig. 3. Simulation results for three different incoming mass flux values (A ₀ = 15 m a⁻¹, 30 m a⁻¹ and 45 m a⁻¹for (a, d), (b, e) and (c, f), respectively). The plots in the top row are illustrations of the non-dimensional velocity fields, and the plots on the bottom row show time series of the total flux out of the solution domain (blue dots plotted at regular time intervals), with the total incoming flux also plotted (red dashed line).

In Figure 3a, for a small value of A ₀, the solution reaches a steady state in which the velocities all fall on the slow, stable branch of the sliding law and the influx is dispersed across the domain. Figure 3d shows the evolution of the total flux out of the domain, which reaches a steady state after ∼2000 years.

In Figure 3b, for a slightly larger value of A ₀, solutions oscillate between slow-flow and fast-flow behaviour. Initially the velocity field falls entirely on the slow branch of the sliding law. As the accumulation builds up near the inflow boundary, increasing the surface slope, the velocity moves to the fast stable branch of the sliding law and an ice stream forms downstream of the perturbation. The ice stream evolves, drawing down the additional upstream mass accumulation, until there is no longer a driving stress sufficient to maintain the stream. The ice stream narrows before collapsing back to its no-stream state. The cycle then repeats.

Finally, in Figure 3c, further increasing the magnitude of the perturbation accumulation, A ₀, results in sustained flow on the fast, stable branch of the sliding law. A steady-state stream forms, with different regions of the domain on different branches of the sliding law.

These are the same three distinct flow regimes that were reported by Reference Sayag and TzipermanSayag and Tziperman (2009), who used the same set-up, but solved the SSA. The regimes illustrated in figures 5–7 in their work directly correspond to those shown in Figure 3 here. We have neglected to include the depth fields and cross-stream velocity profiles, as the most interesting behaviour is captured in the plots included here. Quantitative differences occur between our results and theirs for a given accumulation perturbation, due to the use of a higher-order formulation in our work, which accounts for vertical variation in horizontal velocities. Particularly in regions of low basal velocity, the added component of velocity at the ice surface due to vertical shear is significant, and results in larger outflow from the domain for a given basal velocity. A further difference arises from the addition of the parameter v into the sliding law (Eqn (8)). It makes our solver convergence more rapid and more robust. Moreover, the jumps in velocity in the solutions are less rapid and the overall behaviour is smoother, as can be seen by direct comparison of Figure 3e and f with figures 7a and 6a, respectively, in Reference Sayag and TzipermanSayag and Tziperman (2009). Setting δ = 0 (i.e. taking v = |u _b|), results in behaviour as seen by Reference Sayag and TzipermanSayag and Tziperman (2009), with rapid jumps in flux when the solution moves across the unstable branch of the sliding law (e.g. in the region 11–12 years in figure 6a of their paper).

3.2. Variations in λ

Having verified (Section 3.1) that our formulation reproduces the three flow regimes obtained by Reference Sayag and TzipermanSayag and Tziperman (2009), we diverge from their study to look at the effect of including vertical variation in horizontal velocity due to gravity-driven shear. We do so by varying the scaling parameter, λ, corresponding to the basal-toshear-stress ratio, and looking at the impact on emergent ice-stream behaviour. A value of λ ≪ 1 corresponds to attributing most of the movement to sliding when considering the stress balance (as in the SSA), while λ = 1 corresponds to assuming that all velocity is due to vertical shear (as in the SIA). Intermediate values correspond to various different proportions of the velocity being attributed to sliding when considering the stress balance.

3.2.1. Shear margin width

The value of λ affects the characteristics of an ice stream that forms under this formulation. We consider a simple model of an ice stream where the cross-stream velocity component, v, is zero, and u only varies in the cross-stream direction (Reference Sayag and TzipermanSayag and Tziperman, 2011). Figure 3c shows this to be a fair description for a steady-state ice stream downstream of the onset zone.

Under these simplifications, the only non-trivial, non-dimensional governing equation is

(11)

where

(12)

Assuming that vertical shear is negligible compared with lateral shear in the stream itself, we have τ ₁₃ ≪ (ε/λ)τ ₁₂. This assumption means τ ² can be approximated as , giving

(13)

We now seek the distance, W _m, corresponding to the velocity drop, ΔU, across the shear margin. Assuming a balance between the driving stress and horizontal shear in the margins, we obtain the relation

(14)

Now the value of ΔU is enforced by the triple-valued sliding law (and corresponds to (u ₊ – u ₋) in Fig. 1) . The shear margin must therefore obey

(15)

for constants K ₁ and K ₂.

Taking logarithms of Eqn (15) gives

(16)

suggesting that a logarithmic plot of W _m against λ should give a straight line with gradient −1/(n + 1).

Figure 4 is a plot of margin width against λ from results of simulations run on a 200 × 200 grid. The circles correspond to Newtonian viscosity (n = 1) and the squares to non-Newtonian viscosity (n = 3). Lines with gradients −1/2 and −1/4 are plotted through the Newtonian and non-Newtonian points, respectively. The horizontal gridlines are plotted at intervals of the grid spacing in the simulations. The margin width decreases as λ increases, since membrane stresses play a less dominant role in the stress balance. Without membrane stresses, large jumps in velocity are allowed to occur over a small distance. Moreover, when the margin width is small and spans only two or three gridpoints, it can only take one of a few discrete values, given the choice of resolution. We therefore see that the lines with gradient −1/(n + 1) in Figure 4 are a very good fit to the plotted points, only losing accuracy for numerically unresolved margin widths.

Fig. 4. A plot showing how the shear margin width of an ice stream varies with the scaling parameter, λ. Results from simulations run with n = 1 and 3 are both plotted, where n is Glen’s flow law exponent. Dashed lines with gradient −1/(n + 1) are also plotted.

3.2.2. Effect of λ on flow regime

We also address the question of how the choice of λ affects flow behaviour across all flow regimes with an ensemble of simulations. Values of λ and A ₀ are different in each simulation. Figure 5 maps the resulting flow regimes for different parameter combinations. All three distinct flow regimes occur, each plotted with a different marker; the divides between flow regimes in parameter space are shaded and labelled numerically. We discuss each of these boundaries and explain the physical processes that govern where the boundaries lie in parameter space.

Fig. 5. Phase diagram showing three discrete regions for slow, fast and oscillatory behaviour. Each point plotted represents one simulation that was run with specified values of λ and A ₀.

Boundaries 1 and 2 separate solutions that remain at all times on the slow stable branch of the sliding law (dark circles in Fig. 5) from those that exhibit some oscillatory or streaming behaviour (grey diamonds and white squares in Fig. 5, respectively). The minimum perturbation to the accumulation field large enough to cause the flow to move onto the fast branch of the sliding law occurs at λ ≈ 0.1. For λ < 0.1 the required size of the perturbation accumulation increases as λ decreases further (divide 2), and similarly, for λ > 0.1, the size of a _p must be increased further for there to be any solutions that exhibit fast-flow behaviour (divide 1).

In a bid to explain these observations we consider the necessary conditions for the flow to enter the fast-flow regime. As a result of implementing the triple-valued sliding law, the basal stress must reach a critical value, τ ₊ (Fig. 1) before transitioning onto the unstable branch of the sliding law. Equation (5) describes a balance between basal stress and driving stress, with a correction of due to membrane stresses. For expository simplicity, we write it as

(17)

When λ ∼ 1 the correction for membrane stress effects is , and so to accuracy the balance is between driving and basal stress, as in the SIA. When the driving stress in the onset region reaches the critical value, τ ₊, the transition onto the unstable branch of the sliding law occurs and there is resulting fast-flow behaviour. The correction for membrane stresses, however, becomes more significant for smaller values of λ and so influences how large the driving stress, −h∂s/∂x, needs to become before the basal stress reaches the critical value. In turn, the rate at which the accumulation builds up in the onset zone controls this increase in driving stress. The build-up further depends on the rate at which the incoming ice is drawn downstream. We hence consider the non-dimensional flux out of the domain, which in the Newtonian case, with bed-slope gradient −1 and constant depth field 1 (as is the case for ice streams which form in this suite of simulations), is given by

(18)

where U _s is the transversely averaged stream velocity and W _s the width of the ice stream. It is evident that, when , the increased flux due to vertical shear can be significant, whereas when λ is small this effect is insignificant.

For boundary 1 in Figure 5, λ ≳10⁻¹, and hence the membrane stress correction in Eqn (17) is negligible. Equation (18) shows that as λ increases there will be a resultant increase in outflow from the domain for a given basal velocity. For a given basal velocity field, more ice is therefore removed with higher values of λ; hence less ice accumulates in the onset region where the accumulation perturbation is applied. A larger a _p is required at higher λ for the driving stress to reach the critical value to cause transition onto the unstable branch of the sliding law (τ ₊ in Fig. 1). This explains the positive gradient of curve 1 for λ ≳10⁻¹

The curve defining boundary 2 in Figure 5, in contrast, has a negative gradient. For smaller values of λ, the flow velocities remain completely on the slow branch of the sliding law for much larger accumulation rates. At these lower values of λ, λ/3 is small and so the outflow (Eqn (18)) depends only on u _b. It is necessary then to consider the effect of the membrane stress correction in the stress balance (Eqn (17)). While for λ ≳10⁻¹ the dominant balance is between the basal stress and the driving stress, as λ decreases further, the correction in Eqn (17) becomes more significant. The greater accumulation in the onset region causes the driving stress to increase. The increasing velocity in the onset region results in large negative gradients in the longitudinal stress, since the velocities downstream are lower. The correction is therefore negative, meaning that the driving stress needs to reach larger values before the right-hand side of Eqn (17) will reach the critical value, τ ₊, at which the transition to the unstable branch of the sliding law will occur. This effect is only further amplified as λ decreases further, and so larger perturbations to the accumulation field are required to reach the higher driving stresses that result in oscillatory/streaming behaviour. This explains the negative gradient of the curve separating the slow flow behaviour and oscillatory/streaming behaviour for λ ≲10⁻¹. Since lower values of λ correspond to a stress balance where membrane stresses are of increased significance, this result emphasizes the important effect membrane stresses can have in inhibiting the emergence of ice streams.

Finally we consider what governs boundary 3, which separates flow that periodically surges from those solutions that form a steady-state stream. There is no oscillatory flow once this boundary intersects boundary 1, for λ ≳ 0.2. For a steady-state stream to exist, the outflow from the domain must balance the incoming mass accumulation perturbation, ∫ ∫ a _p(x, y) dx dy. Equation (10) gives an expression for a _p that can be integrated:

(19)

giving the total incoming perturbation accumulation.

Equating the outflow (Eqn (18)) and the incoming mass accumulation (Eqn (19)) gives

(20)

Equation (20) gives the necessary condition that must be satisfied for a solution to form an ice stream. If there is enough incoming mass accumulation to maintain a stream of width W _s, the solution will form a steady-state stream; if there is not, the solution will oscillate. Given that the flux in a stream is given by Eqn (18), why then does a solution with a smaller value of λ need more accumulation before a steady-state stream can form? It is due to the fact that a smaller value of λ results in stresses being delocalized and so results in wider shear margins, as discussed in Section 3.2.1. We therefore have a wider strip of the domain flowing at higher velocities and hence greater outflow from the domain at smaller λ. Consequently more incoming mass accumulation is required with smaller values of λ to be able to maintain the resulting wider streams.

4.1. Steady-state ice streams under a double-valued sliding law

To see the effect of using a double-valued sliding law on ice-stream formation, we ran a set of simulations similar to those in Section 3. For small accumulation perturbations these exhibit the same behaviour as with a triple-valued sliding law (as expected, due to the sliding laws overlying each other at low velocities). It is the behaviour on the unstable branch that is of interest; under what parameter combinations can a stable stream-like solution be maintained?

Experimentation shows that under most parameter combinations, a solution to the governing equations resulting in a steady-state stream is not possible. Streams that form increase in velocity rapidly (there are insufficient stresses to hold them back) and since there is limited incoming accumulation a balance is not possible, i.e. Eqn (20) cannot be satisfied, and so the code stops converging unless membrane stresses increase quickly enough to hold back the stream. As such, only under very specific combinations of parameters does a fully coupled nonlinear evolution of the governing equations with the double-valued sliding law (Eqn (21)) reach a steady-state stream solution. Figure 7 gives an example of one solution resulting in a steady-state Newtonian stream.

Fig. 7. An example simulation for which a steady-state stream forms with n = 1, λ = 0.005, τ ^Y = 0.3, A ₀ = 35 m a⁻¹, δ = 0.1 and all other parameter values given in Table 1.

Note that there are a number of parameters in the sliding law neither constrained by observations nor directly related to the physics (e.g. the smoothing parameter, δ, the value of τ ₊ at the peak of the sliding law and the magnitude of the accumulation perturbation). The large parameter space means it is difficult to choose the most representative values of all the parameters.

4.2. A simpler model for an ice stream underlain by weak till

To determine the parameter combination needed for a steady-state stream to form under the double-valued sliding law, we consider a simple model of an ice stream and ascertain the balance of terms when the basal stress is not providing full resistance to the stream. The cross-stream velocity component is taken to be zero (v = 0) and u is taken to vary only in the cross-stream direction (i.e. u is a function of y only, as for the shear-margin-width scaling in Section 3.2.1). The domain is taken to be of constant depth, h = 1, and so ∇s = ∇b = (−α, 0). (We use α = 1 for the simulations run in this paper, but we carry out this scaling analysis for the most general case.)

Under these simplifications the only non-trivial governing equation is

(22)

where

(23)

As for the shear-margin-width scaling, vertical shear within the stream is taken to be negligible compared with lateral shear across the margins, meaning that τ ² can be approximated as and

(24)

We also assume a yield stress of the basal sediment beneath the ice stream, τ ^Y. Hence τ ₁₃(b) = τ ^Y in Eqn (22) and it becomes

(25)

Integrating Eqn (25) in y and applying the centre-line condition, ∂u _b/∂y = 0 at y = 0, we get

(26)

Assuming the highest velocity is at the centre of the ice stream, this simplifies to

(27)

and we now have

(28)

This describes a velocity profile across an ice stream under our model formulation with a double-valued sliding law.

Furthermore, defining W _s to be the width of an ice stream, u _b(± W _s/2) = u _∞ to be the basal velocity outside that ice stream (usually close to zero), and taking ΔU = u(0) – u _∞, results in

(29)

as a relationship between the non-dimensional velocity change between the inside and outside of the stream and the non-dimensional width of the ice stream, the basal yield stress of the underlying till, the bed slope and the controlling parameters, λ and ε. This is the well-known result for a parabolic velocity cross-profile of an ice stream (e.g. Reference RaymondRaymond, 1996; Reference Van der VeenVan der Veen, 1999), but also includes the dependence on the controlling model parameters, ε and λ.

For example, given ice with non-Newtonian rheology (n = 3), and taking our typical scalings for an ice sheet, ε = 0.005 and λ = 0.015 (Table 1), we have

(30)

Typical driving stresses for pure ice streams are in the region of 10–20 kPa (Reference RaymondRaymond, 1996; Reference Truffer and EchelmeyerTruffer and Echelmeyer, 2003) with basal stresses of ∼1−10 kPa (Reference Sayag and TzipermanSayag and Tziperman, 2009). Taking |τ ^Y – α| ∼ 10⁻¹, we get the result that ice streams with velocities in the range 10²–10³ m a⁻¹, as in the Siple Coast, will be between ∼30 and 60 km wide, which is in agreement with observations (e.g. Reference PatersonPaterson, 1994).

Furthermore, in Section 4 we experimented with simulations using Newtonian rheology (n = 1), α = 1, τ ^Y = 0.3 and varying values of λ and accumulation magnitude. Equation (29) suggests they therefore should obey

(31)

and we illustrate this behaviour with a plot of ΔU against in Figure 8. The deviation from the scaling is due to numerical inaccuracies, as the width of the ice stream is very sensitive to numerical resolution. The simulations for Figure 8 were run with n_i = n_j = 100.

Fig. 8. Plot of ΔU against , showing the fit of our results to the scaling ( is the line plotted).

It is important to note that the simulations run to create this figure do not all necessarily lie within physical parameter space: the velocities span a much wider range than that found in nature. Moreover, this scaling explains why there is such sensitivity to changes in parameters. Given a certain basal stress underneath the stream, the membrane stresses are only sufficient to hold back the stream if the combination of parameters chosen satisfies Eqn (29).

5.1. Comparison of stress fields under the two sliding laws

We now consider in detail the stress balances in ice streams that occur as a result of contrasting sliding laws. Figure 9 shows cross-stream velocity profiles for streams formed under both triple-valued and double-valued sliding laws, and the corresponding plots of the magnitude of basal stress and lateral shear stress.

Fig. 9. Cross-stream velocity and stress profiles for ice streams formed under a triple-valued sliding law (left column) and a double-valued sliding law (right column). (a) and (d) show the cross-stream velocity profiles. (b) and (e) show the corresponding profiles of basal stress magnitude, and (c) and (f) the corresponding lateral shear stress magnitude, |τ ₁₂|.

Figure 9a illustrates how the cross-stream velocity profile changes with variation in λ under the triple-valued sliding law. Narrowing of the shear margins at higher values of λ, as discussed in Section 3.2.1, can be seen clearly.

Figure 9b shows how the corresponding distribution of the basal stress changes with variations in λ. As λ decreases, increasing the significance of the membrane stresses in the stress balance, the delocalizing effect of the membrane stresses means the basal stresses do not increase back to the stable value within the width of the ice stream. Only at higher values of λ, when the margins are narrower, do the basal stresses rapidly increase to be in balance with the driving stress. The behaviour of the basal stress within the shear margin itself is complex. The stress field at first increases while moving from the slow-flow region into the stream, corresponding to an increase in τ _b from the intersection point of the stable branch of the triple-valued sliding law with the line τ _b = τ _d to the value of τ ₊ at the peak of the sliding law. It then decreases rapidly as the velocity increases to the inner ice-stream velocity, consistent with a decrease in basal stress along the unstable branch of the sliding law. This provides a physically plausible profile, if we assume that outside the ice stream the bed is frozen, and across the shear margin there is a transition to temperate bed conditions, due to the increased internal and frictional heating as velocities rapidly increase (Reference SchoofSchoof, 2012). The fast-flowing ice in the stream results in large horizontal shear across the margin (Fig. 9c), and the basal friction increases due to the increasing velocity. These two effects will cause increased heating, both internally and at the bed, and so further into the margin the temperature will reach the melting point, resulting in internal and basal melting. The meltwater will lubricate the bed, causing a rapid decrease in basal friction while moving further into the ice stream, and reaching a minimum inside the shear margin.

Finally, Figure 9c shows the corresponding distribution of lateral shear stress through an ice stream with a triple-valued sliding law as the basal boundary condition. There are maxima in the lateral shear in the margins, where the rapid transition from slow flow to streaming flow occurs. At high λ these peaks in longitudinal shear stress magnitude are restricted to narrower bands. As a consequence, the lateral shear stress is of larger magnitude in these bands because the same velocity change occurs over a smaller distance (∂u _b/∂y is larger).

Figure 9d–f correspond to Figure 9a–c, but are results with a double-valued sliding law. The velocity in the ice stream is no longer defined by the sliding law, but instead by the scaling, Eqn (29). The velocity of a steady-state stream increases linearly with λ, given a specified ice-stream width. The choice of λ is therefore critical in governing whether a stable stream can form in the model or not, once the magnitude of the forcing accumulation perturbation has been chosen (which, in turn, governs what width and velocity of ice stream can form, according to Eqn (20)).

Figure 9e shows the basal stress profiles across ice streams formed with the double-valued sliding law. At higher values of λ the peak of the sliding law (given by the value τ _b = τ ₊ in Fig. 6a) is increased and hence there are larger values of basal stress in the margins when compared with the triple-valued sliding law. If the maximum was only as large as for smaller values of λ, the membrane stresses would not be strong enough to hold the entire stream back. In Figure 9b for the triple-valued case, when λ is larger (and the margins therefore thinner), it is the basal stress that acts to stop runaway of the ice stream. In contrast, in the double-valued case there is a constant, very low resistance provided by the bed underneath the stream, and so resistance in the margins must be larger for the ice stream to be in balance.

Figure 9f is a plot of the magnitude of the lateral shear stress throughout the stream. As λ increases, the lateral shear magnitude also increases, due to the larger velocity change occurring over the margin. There is insufficient basal traction to support the stream, so membrane stresses play a more significant role in balancing the driving stress.

In summary, under a triple-valued sliding law, at lower values of λ the lateral shear and the basal drag both play an important role in the force balance of an ice stream; while at higher values the basal stress plays the dominant role in providing the resistive drag. A double-valued sliding law, however, enforces a low basal stress beneath the ice stream and so a restricted range of parameters is allowed, such that membrane stresses are large enough to hold back the ice stream throughout. For the membrane stresses to be large enough at higher values of λ, a larger velocity difference across the ice-stream margins is required (hence larger values of ∂u _b/∂y).

5.2. A suitable value of λ?

Throughout this work we have presented models with various values of λ to explore the controls on stress balances in ice streams. How does the scale, λ, affect the proportion of velocity which is due to sliding both inside and outside an ice stream? Figure 10a is a plot showing the fraction of the ice-stream velocity that is due to sliding at the bed; Figure 10b illustrates the fraction of surface velocity due to sliding at the bed over the entire domain for various values of λ in both the Newtonian and non-Newtonian cases.

Fig. 10. Plots illustrating the fraction of surface velocity due to sliding at the bed under a triple-valued sliding law. In (a) the fraction within the ice stream is plotted with different values of λ for both Newtonian and non-Newtonian rheology. In (b) this value is plotted over the domain for specific cases.

In the ice stream itself, increasing λ results in a smaller proportion of the streamflow due to sliding, since vertical shear is favoured in the stress balance. For , over 99% of the flow in the stream is due to sliding, as has been observed in some Antarctic ice streams (e.g. Reference LüthiLüthi, 2010).

The triple-valued sliding law predefines the basal velocity both within and outside the ice stream. In particular, the sliding velocity outside the stream corresponds to the point where the slow stable branch intersects, τ _b = τ _d = 1 in Figure 1. The high proportion of velocity outside the stream due to sliding for small λ is a consequence of this; the vertical shear adds only a very small proportion of velocity. If, however, there were no slip outside the ice stream (if the ice were frozen to the bed) then there would still be some flow, all due to vertical shearing, outside the ice stream. This is in contrast to the SSA, which would result in no flow.

A clear benefit, therefore, in using a hybrid formulation (such as the one described in this paper) is that it can attribute varying degrees of surface velocity to sliding at the bed over a single domain. The SSA produces vertically uniform, plug-flow solutions, which may be accurate in the ice stream itself but do not accurately represent the transition to slow flow outside the stream. The SSA also cannot be used to model the emergence of ice streams, since it provides an inaccurate description of the ice flow before the stream has developed. Of course, the SIA with added sliding also allows different areas of the bed to have varying ratios of surface velocity due to sliding (as in the high-λ limit in these plots), but it is well known that the SIA rapidly loses its accuracy as sliding increases (e.g. Reference GudmundssonGudmundsson, 2003). The SIA does not include membrane stresses in the stress balance which, as we have shown here, have important stabilizing effects on ice-stream formation.

The model formulation presented in this paper also has the potential to be applied to other regions of an ice sheet, where there are no abrupt changes in basal or surface elevation and the basal stress is non-zero (i.e. λ ≠ 0 to avoid the membrane stress term in Eqn (5) becoming infinite). These are valid assumptions in the region of grounding zones, although not over the grounding line itself, where the transition to a shear-free model occurs.

The main advantage of this formulation, however, is the fact that it can accurately model ice that is sliding rapidly while also accrediting varying degrees of velocity to sliding over the domain; this is most useful in the region of ice streams where there are abrupt changes in the dominant motion.

5.3. Limitations of the model

The model used in this work is idealized so it is possible to identify critical controlling mechanisms associated with ice-stream stress balances, to improve our understanding of ice-stream stability.

An important limitation of the model is that we do not incorporate energy conservation, and so have assumed the ice to be at a constant temperature. It is well known, however, that thermal effects play an important role in ice-stream dynamics, since ice-stream beds have been found to be at melting point (Reference Engelhardt, Humphrey, Kamb and FahnestockEngelhardt and others, 1990; Reference Iken, Echelmeyer, Harrison and FunkIken and others, 1993). Moreover, Reference Clarke, Nitsan and PatersonClarke and others (1977) suggested that streaming can be an inherent feature of ice sheets and glaciers under certain conditions, due to a creep instability from purely thermomechanical feedbacks. Further investigations including ice dynamics coupled with the thermodynamics used the SIA (Reference MacAyealMacAyeal, 1992; Reference PaynePayne, 1995; Reference Payne and DongelmansPayne and Dongelmans, 1997). Reference HindmarshHindmarsh (2006a) then demonstrated that the inclusion of membrane stresses dampens these thermomechanical instabilities at short transverse wavelengths, although a thermoviscous feedback mechanism does still provide a mechanism for stream formation under the SSA and full Stokes. While our work to date has focused on feedbacks relating to the basal boundary condition, we intend to include thermal effects more comprehensively in future work. Including them in the rheology would have significant effects at the margins; we would expect the shear thinning effect to be further increased.

Another important aspect neglected in this work is the stabilizing effect of buttressing from ice shelves (e.g. Reference ThomasThomas, 1979; Reference HindmarshHindmarsh, 2006b; Reference Goldberg, Holland and SchoofGoldberg and others, 2009) and the corresponding destabilizing effect that calving events can have on an ice stream (e.g. Reference HughesHughes, 1992). In this study we have cut off our domain before the grounding line and so include no ice-shelf/ice-sheet interaction. We would not, however, expect these effects to significantly alter our results. The hybrid formulation would still provide an apt description of ice-sheet flow, and while buttressing might dampen instabilities (and similarly calving events might allow an instability to develop), this would only add some irregularity to results and some quantitative changes in the location of divides between flow regimes.