Model investigations of inland migration of fast-flowing outlet glaciers and ice streams

Abstract

Recent observations of increased discharge through fast-flowing outlet glaciers and ice streams motivate questions concerning the inland migration of regions of fast flow, which could increase drawdown of the ice-sheet interior. To investigate one process that could lead to inland migration we conduct experiments with a two-dimensional, full-stress, transient ice-flow model. An initial steady state is perturbed by initiating a jump in sliding speed over a fraction of the model domain. As a result, longitudinal-stress gradients increase frictional melting upstream from the slow-to-fast sliding transition, and a positive feedback between longitudinal-stress gradients, basal meltwater production and basal sliding causes the sliding transition to migrate upstream over time. The distance and speed of migration depend on the magnitude of the perturbation and on the degree of non-linearity assumed in the link between basal stress and basal sliding: larger perturbations and/or higher degrees of non-linearity lead to farther and faster upstream migration. Migration of the sliding transition causes the ice sheet to thin over time and this change in geometry limits the effects of the positive feedback, ultimately serving to impede continued upstream migration.

Introduction

Discharge from the Greenland and Antarctic ice sheets is dominated by flow through large outlet glaciers and ice streams. These outlets may affect ice-sheet mass balance significantly through changes in their length; increased drawdown of the ice-sheet interior should follow from the headward (inland) lengthening of an outlet glacier or ice stream. With respect to ice-sheet mass balance, this raises questions such as: what conditions and processes encourage or limit inland migration of outlet glaciers? what timescales might be associated with inland migration?

The high speeds of outlet glaciers and ice streams are due largely to fast sliding at or near the ice–bed interface. Fast sliding requires lubrication from basal meltwater, which is produced by a combination of geothermal heat and frictional heat. The rate of frictional heating, which is proportional to the product of sliding speed and basal drag, is largest in regions where sliding and basal drag are intermediate (Raymond, 2000). Such regions include the transition from slow to fast sliding, which occurs in ice-stream tributaries and at the heads of outlet glaciers. The speed transition causes longitudinal-stress gradients that modulate the spatial distribution of basal drag, sliding and frictional melting (Price and others, 2002). This opens the possibility for a positive feedback between longitudinal-stress gradients, sliding and frictional melting in the following way: (1) an increase in sliding downstream from the transition leads to longitudinal-stress gradients that require a local increase in basal friction upstream; (2) increased local basal friction causes increased frictional melting; (3) increased meltwater production increases sliding and propagates the causative longitudinal-stress gradient farther upstream. In turn, basal drag and frictional meltwater production increase still farther upstream and the sliding transition migrates upstream over time.

Below, we incorporate this positive feedback mechanism into a flow model that is sufficiently, but not overly, complex. Our goal is to take first steps towards isolating and understanding the impact of this particular feedback on the evolution of a generic ice-sheet drainage system.

Flow Model

Model description

Our flowline model, which is based on the finite-volume method (Patankar, 1980; Versteeg and Malalasekera, 1995), solves the full two-dimensional (2-D) stress-equilibrium equations for ice flow in plane strain. Details of the model are given by Price and others (2007). Here we discuss the governing equations, the general solution method and the parameterizations used to link stresses in the ice with frictional melting and sliding at the bed.

For a viscous fluid in a low Reynolds number flow, conservation of momentum in a 2-D Cartesian reference frame is expressed as

(1)

where x and z are the along-flow and vertical-coordinate directions, respectively, and repeat indices imply summation. The first term on the lefthand side of Equation (1) is the body force: the product of ice density, ρ, and acceleration due to gravity in the i direction, g_i . The second term is the stress divergence where the full stress tensor, σ_{ij ,} is given by the deviatoric stress, τ_{ij ,} minus the pressure, P:

(2)

where δ_ij is the Krönecker delta (or identity matrix). The constitutive relation linking deviatoric-stress and strain rate is given by

(3)

where is the strain-rate tensor,

(4)

and η is the effective viscosity,

(5)

In Equation (4), u_i represents the component of the velocity vector, u in the x direction or w in the z direction. In Equation (5), B(T) is the temperature-dependent rate factor, n is the power-law exponent (taken equal to 3) and is the effective strain rate given by

(6)

Ice is assumed to be incompressible:

(7)

We solve Equation (1) in a boundary-fitted, orthogonal, curvilinear-coordinate system. The transformation between this system and a standard 2-D Cartesian coordinate system is given by Price and others (2007). Model solutions shown and discussed below have been transformed back to a Cartesian coordinate system.

General solution method

Integrating Equation (7) over a single finite volume (one gridcell) gives

(8)

where the subscripts U, D, T and B refer to the upstream, downstream, top and bottom faces of the volume, respectively. The area of the relevant cell face is A: for example, A_U = W_UΔz_U, where W_U is the width of the flowband, here taken as unity, and Δz_U is the height at the upstream volume face. Equation (8) states that the net mass flux into and out of a volume is zero. Using Equations (2–6) in Equation (1), with an initial guess of the pressure field, we solve for estimated velocities u^* and w^*. Inserting these values into Equation (8) gives

(9)

where S is the mass source (or sink) within each volume. In general, the first estimate of the velocity field (based on a guess of the pressure field) does not satisfy continuity: S ≠ 0. We satisfy continuity using an iterative pressure-correction method (Patankar, 1980) in which a non-zero mass source specifies a pressure perturbation that is used to improve the estimated pressure and velocity fields. Through Equation (9), the updated velocity field leads to a further improvement in the estimate for the mass source/sink (i.e. one with magnitude approaching zero) and a further improvement to the estimated pressure perturbation. Simultaneously, the updated velocity field is used to update the estimated effective viscosity through Equations (4) and (5). Iterations continue until the solution has converged.

The converged velocity field is used to estimate the change in shape of the free surface (and thus the change in the domain geometry) at a future time-step. Changes in domain geometry and the redistribution of mass within the x-z plane are accounted for when re-gridding the finite-volume mesh at the start of each time-step.

Basal motion

The frictional-melt rate, , is the product of the sliding speed, u _b, and the basal drag, τ _b:

(10)

where L _f is the latent heat of fusion. The sliding speed is defined as the horizontal velocity at the top of a basal layer of deforming fluid for which the rate of deformation is:

(11)

The rate factor, C, and the exponent, p, are analogous to A and n in Glen’s flow law for ice, and the effective stress, τ _e, is analogous to the effective strain rate in Equation (6). Setting p > 1 specifies a power-law rheology for the basal layer; p = 1 corresponds to a Newtonian-viscous rheology. The inverse-flow law is:

(12)

where D = C^−1/p is an inverse-rate factor, or a ‘stiffness’ parameter, with units Pa a^1/p and is as in Equation (6). For p ≫ 1 the magnitude of D can be thought of physically as representing the yield strength of the basal layer (see Fig. 1).

Fig. 1. Shear strain rate vs shear stress in basal layer. Thick curves are for D = 20 kPa a^1/p and thin curves represent a reduction in that value by 20% (ΔD of 20%) for different rheologies: p = 1 corresponds to a Newtonian-viscous rheology for the basal layer; p > 1 represents a power-law rheology.

An infinite number of values for p are possible. Here, we examine model behavior for p = 1 and p = 3. We choose p = 1 because it offers a physically plausible, conceptually simple link between basal stress and basal motion that will be the least sensitive to perturbations of the stiffness parameter, and because it has historical precedent (e.g. Alley and others, 1987; Raymond, 2000). To explore the effects of a non-linear sliding relation we choose p = 3, which is sufficiently large to introduce non-linear behavior but sufficiently small to allow for numerical stability and efficiency in model computations. Further, we note that numerous studies show support for a spatially variable pattern of alternating strong and weak basal resistance in regions of flow transition, such as ice-stream tributaries (Price and others, 2002; Joughin and others, 2004a; Peters and others, 2006). In this case, basal resistance comes primarily from ice deformation around roughness elements and/or sticky spots, which we expect to enforce p → n (Raymond, 2000), where n = 3.

We assume that : as the melt rate increases, the stress required to produce a given strain rate within the basal layer decreases. This assumption avoids complications associated with the transport and storage of basal water. Support for its use, at least for describing conditions beneath the West Antarctic ice streams, comes from work by Tulaczyk and others (2000). A more general model would specify the strength of the basal layer over time as a function of the total amount of water present at the bed. Such a model would need to include terms for the import, storage and drainage of basal water in addition to its rate of production (below we discuss how including these effects might affect our conclusions).

Here we specify D to be constant at some initial value (discussed below) until a threshold melt rate, , is reached. For increasing melt rate in the range , the stiffness parameter decreases smoothly from D to (D + ΔD) following a half cosine bell curve (Fig. 2b). For , the stiffness parameter is specified constant with value D + ΔD. The ‘sliding transition’ discussed below is the region over which the stiffness of the basal layer changes from D to D + ΔD (i.e. the region for which ).

Fig. 2. Stiffness parameter as a function of frictional melting rate. Schematics of (a) the basal melting-rate, , increase in the along-flow direction; (b) the stiffness-parameter, D, decrease by ΔD for increasing melt rate in the range ; and (c) variation of D in the along-flow direction.

Boundary conditions, initial conditions and assumptions

The model domain is 650 km long and ice thickness decreases from 2000 m at the upstream end to 400 m at the downstream end. The basal topography is assumed to be flat. Length scales were chosen to approximate those of a typical, large-scale ice-sheet drainage, such as an ice stream in West Antarctica. The upstream boundary of the domain is a flow divide at which we specify a zero-flux boundary condition. The surface is specified as stress-free, and a no-slip condition is specified at the bottom of the deforming basal layer. We avoid having to specify the velocity (or flux) at the downstream boundary by instead specifying hydrostatic pressure. Our region of focus is the transition from slow-to-fast sliding that occurs far inland on the ice sheet (>200 × the mean ice thickness upstream from the downstream boundary); the velocity field at the downstream boundary has no significant effect on the results reported here.

All models start from a quasi-steady state in which the accumulation rate, the stiffness parameter, D, the stress exponent, p, and the thickness of the basal layer are held steady and constant, and the rate of elevation change everywhere within the domain is ≤10⁻⁴ m a⁻¹. Because our goal is to investigate how the evolution of frictional melting affects the strength of the basal layer, we specify temperate, isothermal ice. In this case there is no need for calculation of the temperature field and the rate factor, B(T), is held constant and steady at a value of ∼1.7× 10⁵ Pa a^1/3. We assume the pre-existence of basal meltwater that is sufficient to allow sliding everywhere. In this way, the frictional melt rate above the base melt rate defines the evolution of the basal layer strength through changes in the stiffness parameter, . Henceforth, when we discuss the ‘melting rate’ we specifically mean the frictional melting rate, where τ _b in Equation (10) is calculated as the sum of all basal resistances in the along-flow direction (Van der Veen and Whillans, 1989).

We do not account for the effects of drag from valley side-walls (in the case of outlet glaciers) or from slow-moving interstream ridges (in the case of ice streams). Our focus is primarily on the inland and tributary regions of an ice sheet, where we expect the effects of lateral drag to be small (e.g. Price and others, 2002). We acknowledge that, in regions like an ice-stream onset, lateral drag may not be insignificant. Implications of this simplification are discussed further below.

Perturbing equilibrium models

We perturb equilibrium models by introducing an instantaneous reduction of the stiffness parameter (equivalent to an instantaneous jump in sliding speed) over approximately the downstream half of the model domain. Here, ‘instantaneous’ refers to a change initiated over a single model time-step of 1 year. This is sufficiently long to insure that the elastic contribution to the stress-equilibrium equations is negligible,¹ and as in the majority of ice-sheet models we neglect elastic effects. We arbitrarily choose a point along flow near the middle of the model domain and define the melting rate at that point as the threshold melting rate ( in Fig. 2). A specified value of then allows us to define , which determines where along the melting-rate profile the non-zero gradient in the stiffness parameter begins and ends. Starting from the initial steady-state condition, we specify that the stiffness parameter remains unchanged (i.e. D _new = D _ini) at locations where , and at locations where the value of the stiffness parameter is reduced by some fraction of its initial value (D _new = D _ini + ΔD). In between, the stiffness parameter decreases smoothly as a function of , as shown in Figure 2.

We choose D _ini = 20 kPa a^1/p ,which gives surface velocities near the center of the model domain that are representative of those observed near modern-day ice-stream onset regions (≤100 m a⁻¹). In fact, the initial value of D is somewhat arbitrary;² the fractional change, ΔD, relative to its initial value is important for determining model response. Values reported for ΔD below are given as a percent reduction from the initial value (e.g. ΔD = −0.1 × D _ini is reported as ’ΔD of 10%). Below, we present a detailed discussion of the model response to perturbations, ΔD, of 5%, 10% and 15%. The latter value approaches the limit on a perturbation magnitude imposed by issues of numerical convenience and economy (the larger the perturbation, the more difficult it becomes to achieve a converged solution within a practical number of model iterations). Nevertheless, the model response to this range of perturbations is reasonable with respect to recent observations of outlet-glacier acceleration. For example, using p = 3, a ΔD of 15% leads to just under a two-fold increase in sliding downstream from the transition region (similar to that observed in Greenland, as reported by Rignot and Kanagaratnam, 2006).

We examine model response to perturbations for two sets of experiments: one set with the feedback turned ‘off’, and the other with the feedback turned ‘on’. In both cases, the perturbation, ΔD, causes increased meltwater production upstream. However, for experiments with the feedback turned off (‘reference’ models), the stiffness parameter is held fixed after the initial perturbation and as a result the sliding transition does not move over time. The reference models are analogous to model experiments by Payne and others (2004) in which the geometry and stress field in the ice evolved in response to a perturbation in basal resistance but the distribution of basal resistance did not. For experiments with the feedback turned on (‘linked’ models), and evolve in response to the changing stress field and geometry; the sliding transition can move over time. For the linked models, we assume that the patterns of and are in equilibrium with the stress field: during each time-step, we iterate on the value of and the value of while holding the geometry constant until we obtain a consistent solution. Iterations continue until the velocity field no longer changes. Computationally this is not problematic because convergence usually occurs within a few iterations.

Results from experiments with the reference model are used as a control for comparison with those from the linked model. This is convenient because the dominant physical processes are common to both the reference and linked models, but the former avoids the additional complication of a migrating sliding transition.

Model Results

We discuss the results from experiments with the reference and linked models by examining the evolution of the stress fields, the geometry and the basal melting rate after a perturbation has been applied (Figs 3–12). In our discussion of terms in the force balance we follow a commonly used sign convention (e.g. Van der Veen and Whillans, 1989): driving stress is positive in the downslope direction; resistance to the driving stress by friction at the ice–bed interface is taken as positive; a pulling (pushing) force from the downstream direction is taken as a positive (negative) longitudinal-stress gradient.

Fig. 3. Time series of the distribution of longitudinal-stress gradient for the reference model. The along-flow coordinate is given by x, the height above the glacier bed is given by z, ΔD is 10% and p = 3. The vertical dashed line marks the midpoint of the sliding transition. Vertical and horizontal coordinates are scaled by the initial ice thickness, H, at the sliding transition (at x = 0). The longitudinal-stress gradient is scaled by the magnitude of the maximum instantaneous longitudinal-stress gradient (that at t = 0). Note that the greyscale axes span a different range in each panel.

Fig. 4. Stress balance and melting rate for the reference model. Panels in each column cover the same region and represent the same times and values of ΔD and p as in Figure 3. Top row: driving stress (thick solid curve) and basal drag (dashed curve) scaled by the initial driving stress at x = 0. The thin black curve is the steady-state driving stress prior to the perturbation. Middle row: depth-averaged longitudinal-stress gradient. Bottom row: frictional-melting rate scaled by the value of the threshold melting rate. For reference, the two thin curves are the initial (steady-state) melting-rate profile and the instantaneous melting-rate profile after the perturbation.

Fig. 5. Spatial and temporal evolution of the threshold melting rate, , for the reference model (dashed curves) and the linked model (solid curves). Curves track the location of relative to its initial position as a function of time after a perturbation, ΔD, of 10%, with p = 3. For the linked model, the location of is synonymous with the location where the stiffness parameter, D, starts to decrease (see Fig. 2). Zero on the vertical axis coincides with zero on the horizontal axis in Figures 3 and 6. The inset shows details during the first 250 years after the perturbation. Details of the controlling processes during different time intervals (1. longitudinal-stress gradients (LSG); 2. increasing surface slope; and 3. long-term thinning) are discussed in text.

Fig. 6. Time series of the spatial distribution of longitudinal-stress gradient for the linked model. The along-flow coordinate is given by x, the height above the glacier bed is given by z, ΔD is 10% and p = 3. The vertical dashed line marks the midpoint of the sliding transition. Vertical and horizontal coordinates are scaled by the initial ice thickness, H, at the sliding transition (at x = 0). The longitudinal-stress gradient is scaled by the magnitude of the maximum instantaneous longitudinal-stress gradient (that at t = 0). Note that the greyscale axes span a different range in each panel.

Fig. 7. Stress balance and melting rate for the linked model. Panels in each column cover the same region and represent the same times and values of ΔD and p as in Figure 6. Top row: driving stress (thick solid curve) and basal drag (dashed curve) scaled by the initial driving stress at x = 0. The thin black curve is the steady-state driving stress prior to the perturbation. Middle row: depth-averaged longitudinal-stress gradient. Bottom row: frictional-melting rate scaled by the value of the threshold melting rate. For reference, the two thin curves are the initial (steady-state) melting-rate profile and the instantaneous melting-rate profile after the perturbation.

Fig. 8. Spatial and temporal evolution of the threshold melting rate, , for linked model experiments. Curves a–f track the location of relative to its initial position as a function of time after perturbations. a: ΔD of 5%, with p = 1; b: ΔD of 10%, with p = 1; c: ΔD of 15%, with p = 1; d: ΔD of 5%, with p = 3; e: ΔD of 10%, with p = 3; f: ΔD of 15%, with p = 3. The location of is synonymous with the location where the stiffness parameter, D, starts to decrease.

Fig. 9. Shift of melting-rate profile after a perturbation to the melting rate. Solid curves in (a) and (b) represent two different initial melting-rate profiles that experience the same perturbation to the melting rate, , at location x ₀. After the perturbation (dashed curves), the melting rate initially at x ₀ has been displaced upstream by a distance Δx. The horizontal axis encompasses a distance of several tens of ice thicknesses near the location of the melting-rate perturbation.

Fig. 10. Thinning rate as a function of time upstream from the perturbation for linked model experiments. Curves track thinning rates at a distance of 150× the ice thickness upstream from the initial perturbation (near the ice divide). Dashed curves represent the response to ΔD of 5%, dash-dot curves represent the response to ΔD of 10% and solid curves represent the response to ΔD of 15%. Curves for p = 1 and p = 3 are labeled.

Fig. 11. Relative differences in thinning characteristics between linked and reference models. Curves represent the percentage difference in the maximum thinning rate and the timing of maximum thinning rate. The change in the stiffness parameter from its initial value is plotted on the horizontal axis. The percentage increase in the maximum rate of thinning is plotted on the vertical axis (dotted curves), as is the percentage increase in how soon that maximum occurs (dashed curves). Positive numbers indicate that maximum thinning rates are larger and occur sooner in the linked model. Circles and squares represent results for p = 1 and p = 3, respectively.

Fig. 12. Thinning rate as a function of distance upstream from the perturbation for linked model with ΔD of 10% and p = 3 (same conditions as shown in Figs 6 and 7). Curves track thinning rates upstream at times 0, 3, 5, 10 and 20 years after the initial perturbation.

Evolution of the stress field and geometry: reference model

When a perturbation is first applied (t = 0 years; Fig. 3), longitudinal-stress gradients immediately become positive throughout the ice thickness upstream of the sliding transition, and negative throughout the thickness downstream. In response, basal drag is elevated upstream from the transition and depressed downstream relative to the driving stress. The driving stress does not respond immediately (t = 0 years; Fig. 4) because the surface geometry has not yet had time to adjust to the new stress state. For the case shown in Figures 3 and 4, the effects of this initial ‘impulse’ extend ∼20× the ice thickness on either side of the transition. The melting-rate profile everywhere increases relative to the initial condition (t = 0 years; Fig. 4) because downstream from the transition the sliding speed has increased, and both the sliding speed and the basal drag have increased upstream.

After 3–5 years, the longitudinal-stress gradients are more evenly distributed and their magnitudes have decreased (t = 3–5 years; Figs 3 and 4). The surface geometry starts to adjust to the new stress field by getting steeper upstream from the transition in response to relatively high basal drag there. Downstream the surface gets flatter in response to the relatively low basal drag there. As a result, the driving-stress profile evolves to more closely match that of the basal drag. The melting-rate profile reflects these changes as well; relative to the melting-rate profile immediately after the perturbation, it continues to increase upstream of the transition but it decreases downstream (t = 3–5 years; Fig. 4).

The surface geometry continues to adjust to the evolving stress field, and after 10–20 years the longitudinal-stress gradients are very different from the initial pattern (t = 10–20 years; Figs 3 and 4). Downstream from the transition the reduced basal stiffness results in faster sliding. Because the lower portion of the ice column moves relatively faster, less internal deformation is needed to pass the balance flux and, in steady state, the upper portion of the ice column moves relatively slower. The result is that across the transition the pattern of longitudinal-stress gradients along the surface is similar to that at the bed but rotated by 180°; moving downstream it grades from extensional to compressional at the bed but from compressional to extensional at the surface. The reduction in velocity in the upper ice column is of a smaller magnitude than the increase in velocity near the bed. Nevertheless, because it occurs over a large fraction of the ice column, and that fraction is relatively stiff (because of its smaller effective viscosity), the depth-averaged pattern of longitudinal-stress gradients more closely resembles the pattern along the surface than along the bed. The effect of longitudinal-stress gradients is small except within ∼5× the ice thickness on either side of the sliding transition, where positive gradients raise the basal drag slightly. Previous authors have made similar observations with respect to equilibrium patterns of longitudinal-stress gradients: Weertman (1957) noted the 180° rotation pattern at a jump in basal sliding, and Budd (1970) discussed the ‘smoothing’ effect of longitudinal-stress gradients on the basal-drag profile with respect to the driving-stress profile.

Large-scale adjustments to the velocity fields and surface shape are essentially complete after 500 years. While both the driving stress and the basal drag decrease across the transition, the gradient in both terms on either side of the transition is nearly the same as that prior to the perturbation (t = 500 years; Figs 3 and 4). The reduction in driving stress and basal drag across the transition reflects the adjustment of the geometry needed to achieve a new steady state. Upstream, steeper surface slopes are required to pass the steady-state ice flux through a relatively thinner ice column, for which the basal velocity has not changed. Downstream, shallower slopes are required to pass the steady-state ice flux because the basal velocities have increased (owing to the reduction in the basal stiffness parameter). Ice-sheet thinning continues at a decreasing rate for another several thousand years.

Evolution of melting rate: reference model

The dashed curve in Figure 5 shows the location of in the reference model as a function of time after the initial perturbation. It is convenient to divide the trajectory of into three time periods for which the change in the position of is controlled by three different processes: (1) an initial, rapid but short-lived (order of tens of years or less) period of upstream motion; (2) a longer period (up to several hundred years) of more gradual upstream migration; and (3) a long period (several thousand years) during which the position of stabilizes at its most upstream extent and then migrates slowly back downstream.

The initial rapid upstream motion is a response to the sudden change in sliding speed imposed at the transition (interval 1; Fig. 5 inset): tension across the transition increases the basal friction, sliding speed, and hence the melting rate, upstream. This initial, rapid, upstream propagation of is short-lived because the magnitude of the longitudinal-stress gradient decreases with time and so it becomes less effective at raising the melting rate to . The second period of slower but sustained upstream migration is a response to the changing ice-sheet geometry. As thinning propagates upstream from the site of the initial perturbation, the surface slope increases and, in response, the driving stress, basal drag, sliding speed and melting rate continue to increase slowly (interval 2; Fig. 5 inset).

With more time, the geometry continues to adjust; the surface slope flattens, the driving stress, the basal drag and the sliding speed all start to decrease. In this third period, stabilizes its position temporarily before beginning a slow migration back downstream (interval 3; Fig. 5). The slow adjustment is associated with slow thinning of the ice column that is ongoing for several thousand years after the perturbation. By 5000 years after the initial perturbation, the rate of thinning is small and the location of nears a new steady-state position. This position is farther upstream than initially because of changes in the geometry: upstream from the transition, surface slope (and thus driving stress, basal drag and frictional melting) has increased in order to accommodate the steady-state flux through a relatively thinner column of ice (accumulation rate is held constant).

Evolution of the stress field and geometry: linked model

A notable difference between the linked model and the reference model is that positive feedback between longitudinal-stress gradients, basal sliding and frictional meltwater production in the linked model causes the sliding transition to migrate upstream. Although the instantaneous response to the perturbation (t = 0 years; Figs 6 and 7) is the same as for the reference model, differences are significant immediately thereafter: because the feedback allows the basal layer to soften when the melt rate reaches and exceeds , the position of the longitudinal-stress gradient impulse has moved upstream by about 15× the ice thickness after 3 years (t = 3 years; Figs 6 and 7). The driving stress, which is controlled primarily by the surface slope, has started to adjust to the initial perturbation but lags the change in the pattern of basal drag, which has already moved upstream along with the causative longitudinal-stress gradient (t = 3 years; Fig. 7). The magnitude of the longitudinal-stress gradient is larger after 3 years than at the time of the initial perturbation (t = 0 years; Fig. 7). This is because the initial perturbation effectively increases the gradient of the melting-rate profile, which narrows the width of the transition zone and increases the stress gradient associated with ΔD.

After 5 years the sliding transition has migrated upstream by about 20× the ice thickness (t = 5 years; Figs 6 and 7). While upstream migration of the transition zone allows for a sustained perturbation, the longitudinal-stress gradient eventually decreases in magnitude because the geometry of the ice sheet begins to adjust: by 5 years, the driving-stress profile has started to ‘catch up’ with the basal-drag profile (t = 5 years; Fig. 7). After 10–20 years the longitudinal-stress gradient in the transition region has decreased significantly and the geometry has adjusted so that the driving stress and basal drag are nearly equal across the transition region (t = 10–20 years; Fig. 7). Following these adjustments the pattern of longitudinal-stress gradients changes in the same way as those discussed above for the reference case.

As with the reference model, for the next ∼500 years the ice sheet continues to adjust to the new state by slowly thinning. Relative to the reference case, significant differences are: (1) the sliding transition is located ∼20× the ice-thickness upstream from its initial location, (2) the sliding transition is wider than initially, and (3) the associated longitudinal-stress gradient has smaller magnitude (due, in part, to the wider transition). All of these differences are a result of the co-evolution of stresses at the ice–bed interface with the melting-rate profile and with the basal-stiffness parameter.

Evolution of melting rate: linked model

Evolution of the melting rate for the linked model (solid curve, Fig. 5) is similar to that for the reference model, but in the linked model the location of has particular importance with respect to flow dynamics. It is synonymous with the location of the upstream boundary to the sliding transition. Compared to the reference model, there is a longer period during which moves rapidly upstream through longitudinal-stress gradients because the initial perturbation can be sustained for several years as the sliding transition moves upstream. The net effect is that migrates farther upstream in the linked case before slowly moving back downstream. Figure 7 confirms that the initial, rapid upstream motion of is caused by longitudinal-stress gradients and not by the effects of changing ice-sheet geometry: for at least 5 years after the perturbation the change in driving stress lags behind the change in basal drag. During this time, the former is controlled by the changing ice-sheet geometry while the latter is controlled by longitudinal-stress gradients.

Sensitivity to changes in ΔD, p and the melting-rate profile

Figure 8 shows the position of as a function of time after initial perturbations of ΔD varying from 5 to 15%, and with p = 1 and p = 3. Results show that larger initial perturbations cause larger transient longitudinal-stress gradients, which affect the basal friction and sliding speed, and hence the melting rate, over larger distances upstream. A larger initial perturbation comes either from larger jumps in sliding speed across the transition (associated with larger ΔD) or for a basal layer with a higher degree of non-linearity (higher p in Equation (11)).

Upstream migration of the position of also depends on the longitudinal profile of the melting rate; for the same initial perturbation, the smaller the gradient in the melting-rate profile the farther upstream the sliding transition will move. Figure 9 illustrates how the initial gradient influences subsequent upstream migration of the transition. For initial melting-rate profiles (solid curves in Fig. 9a and b), the same perturbation to the melting rate at the sliding transition, , results in new melting-rate profiles (dashed curves). After the perturbation, the position of has moved farther upstream along the new low-gradient profile (Fig. 9b) than for the new high-gradient profile (Fig. 9a); other things being equal, melting-rate profiles with low gradients favor upstream propagation of a perturbation to the melting rate.

The melting-rate profile depends on the downslope profiles of sliding speed and basal drag (Equation (10)). Speed generally increases downslope, so that uniform or decreasing basal drag is needed to maintain a low-gradient profile. Here, that condition will occur for larger values of p. Conversely, a limiting aspect of a melting-rate profile with a small gradient is that a given ΔD will occur over a larger spatial distance (Fig. 2) relative to a melting-rate profile with a large gradient. In this case, the longitudinal-stress gradient associated with a given ΔD will be of smaller magnitude than for a melting-rate profile with a large gradient.

The second period of upstream migration, in which responds to increasing surface slopes, is largely controlled by the magnitude of the initial perturbation. A larger initial perturbation is associated with a larger magnitude of thinning, which has a larger effect on the transient surface slope, driving stress, basal drag, sliding speed and melting rate. Because the associated rate of thinning is larger for a larger perturbation, significant thinning of the ice column takes place relatively sooner and the amount of time spends at its maximum upstream position is relatively shorter: Figure 8 shows that as ΔD and p increase, the maximum of the curve of the vs time becomes more narrow and skewed towards time zero. This and other features of behavior shown in Figure 8 are more easily understood in the context of long-term and far-field responses to the initial perturbation.

Long-term response

After the initial perturbation, the long-term response is ice-sheet thinning; the decrease in the stiffness parameter over a portion of the model domain increases the sliding speed and initially removes mass faster than it is replaced through accumulation. Thinning continues for several thousand years as the geometry and stress fields adjust to the new basal conditions and come to a new equilibrium with the accumulation rate. These adjustments to the geometry of the ice sheet have a strong influence on the long-term trends of melting rates and on the position of the sliding transition. Figure 10 shows the temporal evolution of thinning rates at a distance of 150× the ice thickness upstream from the initial perturbation (near the ice divide) for the same cases as shown in Figure 8. Qualitatively, the curves that track the evolution of thinning rate are similar to those that track the position of . As ΔD and/or p increases, the curves in both Figures 8 and 10 increase in amplitude and are more narrowly peaked and skewed towards time zero. These similarities confirm the tight coupling between the pattern of ice-sheet thinning and the position of the sliding transition.

Response with feedback ‘on’ and feedback ‘off’

Systematic differences between the linked and reference models, evident in Figures 3–7, increase as values of ΔD and/or p increase. For example, the maximum upstream extent of is greater for the linked model than for the reference model, and this difference increases as ΔD and/or p increases. Similar trends but with respect to thinning rates far upstream from the perturbation are summarized in Figure 11, which shows the percentage increase in the maximum thinning rate for the linked model relative to the reference model (dotted curve). Also shown in Figure 11 is the percentage decrease in the arrival time of the maximum thinning rate (plotted as a positive number) for the linked model relative to the reference model (dashed curve). The maximum rate of thinning is larger, and occurs sooner, for the linked models than for the reference models. For p = 1, the largest increase from the reference models to the linked models is ∼8%, and for p = 3 the largest increase is ∼30%. For larger perturbations and higher degrees of non-linearity in the basal sliding relation, larger increases are expected.

The differences arise because, in the linked model, the initial perturbation does not immediately begin to decay but sustains itself through longitudinal-stress gradients for up to 10 years as it propagates upstream. Because the perturbation does not immediately begin to decrease in magnitude, the maximum rate of thinning and the net thinning are larger. As the location of the perturbation moves rapidly upstream for a period of time, thinning at some fixed distance upstream starts sooner than if the site of the perturbation were fixed. Both of these effects are illustrated in Figure 12, which shows thinning rates at different times upstream from the perturbation for the models discussed in Figures 3, 4, 6 and 7.

Discussion

Basal rheology and magnitude of perturbations

Longitudinal-stress gradients that arise immediately in response to a downstream sliding perturbation cause increased basal melting upstream. A positive feedback between longitudinal-stress gradients, basal meltwater production and basal sliding causes the sliding transition to migrate upstream, which causes ice-sheet thinning. Inland migration of the sliding transition and ice-sheet thinning increase for larger ΔD and higher values of p. Here we have restricted our experiments to low-order stress exponents; if subglacial till behaves as a power-law fluid, stress exponents in the range of 5–13 may be more appropriate (Rathbun and others, 2005; Tulaczyk, 2006). Such large values would cause the sliding transition to migrate faster and farther upstream than the results shown here.

While we have applied perturbations to the basal stiffness parameter that are relatively small (ΔD ≤ 15%), we have applied them quickly (over a time-step of 1 year) and over a large portion (∼50%) of the model domain. Real-world perturbations of this type might occur following the rapid movement of ‘pockets’ of subglacial water beneath the ice-stream tributaries in West Antarctica (Gray and others, 2005), or similarly, following the motion of subglacial water associated with the draining and filling of neighboring subglacial lakes, as observed in inland regions of East Antarctica (Wingham and others, 2006). Other ‘pulling’ perturbations might result from, for example, the loss of basal traction provided by an ice plain (e.g. Payne and others, 2004) or the collapse and removal of an ice shelf or ice tongue (Joughin and others, 2004b; Scambos and others, 2004). While such perturbations might be very large locally (i.e. near the grounding line), we expect they would be substantially smaller by the time they have propagated several hundred kilometers inland.

It is not clear how the perturbations applied to our model scale to those experienced by real outlet glaciers and ice streams, but we can attempt some approximate comparisons. For example, Rignot and Kanagaratnam (2006) reported observations in Greenland that indicate a doubling of speed for Jakobshavn Isbræ in the period 1995–2005, and a doubling of speed for Kangerdlugssuaq Gletscher over 5 years (2000–05). In our experiments, a doubling of speed requires a perturbation ΔD ≈ 50% for p = 1, or ΔD ≈ 20% for p = 3, which we apply over much shorter timescales (∼1 year). Observations from Pine Island Glacier, West Antarctica, indicate mean thinning rates along the 200 km trunk of ∼0.75 m a⁻¹ (Shepherd and others, 2001). Mean thinning rates along slower-moving portions of the glacier, including the area upstream of the main trunk, are ∼0.10 m a⁻¹. For comparison, the maximum rate of thinning in our model for ΔD of 15% and p = 3 is ∼0.5 m a⁻¹ (Fig. 12), occurring immediately after and at the location of the perturbation. Thinning rates decrease with time and farther upstream (Fig. 10); maximum thinning rates ∼300 km upstream are 0.06 m a⁻¹,about half the mean value reported by Shepherd and others (2001).

Limits to ongoing inland migration

None of our experiments produced ongoing upstream migration of the sliding transition. Following the initial perturbation, the magnitude of the longitudinal-stress gradient across the transition starts to decrease as soon as the geometry begins to adjust to the new stress state. As the geometry changes, the pattern of longitudinal-stress gradients changes to one less favorable for continued upstream propagation of the sliding transition. To maintain large longitudinal-stress gradients across the sliding transition (necessary to maintain continued upstream migration) the transition region itself must move far enough upstream during each time-step to outpace the effects of topographic diffusion. Here, topographic diffusion, which is fast on ice sheets due to large ice thicknesses and small surface slopes, limits ongoing migration.

Inland migration of the sliding transition is also limited by the upstream gradient of the melting-rate profile. Sliding speed and basal drag, and thus frictional melting, generally decrease upstream; longitudinal-stress gradients arising from a perturbation, ΔD, will then become less effective at raising the melting-rate profile above its ‘background’ level farther upstream. The sliding transition can travel farther upstream in cases where the melting-rate profile has a small gradient. Flow speeds generally increase downstream, so uniform or decreasing basal drag along flow is the one condition that favors a melting-rate profile with a small gradient. Regions on ice sheets where the basal drag is uniform (and small) over many hundreds of kilometers include ice streams and ice plains. In these regions, longitudinal perturbations are transmitted rapidly and over large distances, a conclusion that is supported by several other recent, model-based studies (Payne and others, 2004; Gudmundsson, 2006).

Several of our assumptions could also serve to impede upstream migration of the sliding transition. (1) We specify that the upstream boundary is always a flow divide where the rate of frictional melting is zero. In reality, a wave of thinning propagating upstream would force divide migration and the frictional-melt rate at the former divide location would increase over time as the divide migrated away. The timescale for divide migration is probably longer than the timescale over which the sliding transition can migrate continuously via longitudinal-stress gradients; it is unlikely that our choice of upstream boundary condition is an important factor in limiting upstream migration. (2) We assume that the background melt rate is spatially constant. Conceptually we think of this background melting as coming from an excess geothermal flux over the energy conducted upwards into the ice sheet. If the geothermal flux, surface accumulation rate and temperature are spatially constant, increasing ice thickness moving upstream would act to increase the background melt rate upstream. Relative to the results discussed above, this effect might help to further sustain upstream migration. (3) Experiments confirm that the rate and total upstream migration are larger for p = 3 than for p = 1 and this trend is likely to hold for p > 3. It is therefore possible that continued upstream migration could occur with this physical model for cases where p > 3.

Importance of simplifying assumptions

We do not account for the flow resistance provided by drag against valley side-walls (in the case of an outlet glacier) or against slower-moving interstream ridge ice (in the case of an ice stream). We assume that longitudinal-stress gradients accommodate a reduction in basal drag by increasing it elsewhere so that ultimately the bed still provides all the resistance to flow. In reality, some fraction of the reduced basal drag will be supported at the margins through lateral-stress gradients. This omission means that we overestimate changes in the pattern of basal friction, and hence overestimate increases in the melting rate upstream from the sliding transition. We therefore expect that equivalent experiments conducted with a three-dimensional flow model would give results similar to those shown here, but the magnitude of upstream migration and subsequent changes in geometry would be smaller.

We have assumed that basal resistance is controlled by the basal melt rate. Another possibility is that it depends on the amount of water at the bed, in which case the details of subglacial water transport and storage would be important. For example, for cases where the bed is undersaturated, increases in the melting rate are immediately accommodated by the basal-water system (i.e. the ‘drainage-limited’ state of Raymond, 2000); the timescale for the basal-water system to adjust to an increase in melt rate is shorter than the timescale for the increased melt rate to reduce basal resistance. In this case the feedback discussed above would be impeded. However, if drainage was slow and/or storage capacity was significant, an increase in melt rate might affect basal resistance for some time after the melt rate decreased to or below former levels (i.e. the timescale for the basal-water system to adjust to an increase in melt rate is longer than the timescale for the increased melt rate to reduce basal resistance). In this case, the basal-water system would serve to enhance the feedback discussed above.

Conclusions

A downstream perturbation to basal sliding (a ’pulling’stress) can initialize a positive feedback between stresses in the ice, frictional melting and the basal boundary condition that leads to inland migration of the transition from slow to fast sliding. An initial, short (∼10 years) period of rapid upstream migration is associated with longitudinal-stress gradients and accounts for up to 70% of the total distance that the transition moves. Upstream migration of the sliding transition increases with the size of the perturbation and with the degree of non-linearity assumed in the relation linking basal stress to basal motion. Here, a 15% reduction in the stiffness of the basal layer causes the sliding transition to migrate upstream by about 35× the ice thickness in 250 years.

The feedback between longitudinal-stress gradients, basal sliding and frictional meltwater production has important far-field effects: upstream from the perturbation, maximum thinning rates are larger and occur sooner than for cases where the feedback is not included. These differences increase as the initial perturbation increases and/or as the degree of non-linearity in the relation linking basal stress to basal motion increases. Here, maximum differences are ∼30% for experiments with p = 3. For systems where the stress exponent is even higher (10 or more), much larger differences can be expected.

Our model results and conclusions are conditional on the assumption that basal sliding is a function of the frictional-melt rate. A more realistic assumption might be that basal sliding is related to the total amount of water at the ice–bed interface. Future work should examine how the details of the basal-water system might affect inland migration of outlet glaciers and ice streams.