Model description

Our flowline model, which is based on the finite-volume method (Reference PatankarPatankar, 1980; Reference Versteeg and MalalasekeraVersteeg 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 Reference Price, Waddington and ConwayPrice 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 Reference Price, Waddington and ConwayPrice 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 (Reference PatankarPatankar, 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. Reference Alley, Blankenship, Bentley and RooneyAlley and others, 1987; Reference RaymondRaymond, 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 (Reference Price, Bindschadler, Hulbe and BlankenshipPrice and others, 2002; Reference Joughin, MacAyeal and TulaczykJoughin and others, 2004a; Reference PetersPeters 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 (Reference RaymondRaymond, 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 Reference Tulaczyk, Kamb and EngelhardtTulaczyk 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 (Reference Van der Veen and WhillansVan 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. Reference Price, Bindschadler, Hulbe and BlankenshipPrice 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 Reference Rignot and KanagaratnamRignot 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 Reference Payne, Vieli, Shepherd, Wingham and RignotPayne 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.

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: Reference WeertmanWeertman (1957) noted the 180° rotation pattern at a jump in basal sliding, and Reference BuddBudd (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.

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 (Reference Rathbun, Marone, Anandakrishnan and AlleyRathbun and others, 2005; Reference TulaczykTulaczyk, 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 (Reference Gray, Joughin, Tulaczyk, Spikes, Bindschadler and JezekGray 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 (Reference Wingham, Siegert, Shepherd and MuirWingham and others, 2006). Other ‘pulling’ perturbations might result from, for example, the loss of basal traction provided by an ice plain (e.g. Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004) or the collapse and removal of an ice shelf or ice tongue (Reference Joughin, Abdalati and FahnestockJoughin and others, 2004b; Reference Scambos, Bohlander, Shuman and SkvarcaScambos 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, Reference Rignot and KanagaratnamRignot 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⁻¹ (Reference Shepherd, Wingham, Mansley and CorrShepherd 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 Reference Shepherd, Wingham, Mansley and CorrShepherd 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 (Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004; Reference GudmundssonGudmundsson, 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 Reference RaymondRaymond, 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.