1. Intoduction

Assessing the probability of rapid grounding-line retreat, and thus of rapid contribution to sea level, is a primary goal of glaciological studies of the West Antarctic ice sheet (WAIS). It is well documented that the grounding line (or zone) has retreated > 1000 km, from the continental shelf edge to its present location, in the time since the Last Glacial Maximum (LGM; Reference Conway, Hall, Denton, Gades and Wad-dingtonConway and others, 1999; Reference Shipp, Anderson and DomackShipp and others, 1999). Of interest now is whether that retreat will continue into the future, and if so, at what rate.

A popular paradigm, and motivation for many West Antarctic studies, is the marine ice-sheet instability hypothesis (Reference MercerMercer, 1978; Reference Thomas, Sanderson and RoseThomas and others, 1979), in which an increase in the rate of ice-shelf discharge causes runaway thinning at, and retreat of, the grounding line. In contrast, laboratory and modeling studies of subglacial till mechanics have been used to suggest that ice streams can thin and flatten to a point beyond which their fast flow will stop (Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000). Observations of the deceleration of the ice plain downstream of Whillans Ice Stream seem to support that prediction (Reference Bindschadler and VornbergerBindschadler and Vorn-berger, 1998; Reference Joughin and TulaczykJoughin and Tulaczyk, 2002). Yet the persistence of Siple Dome through the last glacial cycle suggests that since at least that time ice streams were present, and not much thicker than at the present day (cf., Reference NeresonNereson, 2000). Neither the notion of rapid grounding-line retreat nor sudden shut-down is clearly favored. Moreover, the record would appear to support a perennially thin yet perennially fast-flowing ice sheet.

The most likely source for rapid change in the location of the present-day Ross grounding line is the action of fast-flowing ice streams. The ice streams of West Antarctica progress downstream through a series of distinct flow regimes, from broadly convergent flow in a catchment to extending and accelerating flow down a trunk and across an ‘ice plain’ of very low slope, ending where they go afloat to join the Ross Ice Shelf (RIS) (Reference Alley and WhillansAlley and Whillans, 1991; Reference JoughinJoughin and others, 1999). Ungrounding occurs over the broad expanse of lightly grounded ‘ice plains’ at the confluence of Mercer and Whillans Ice Streams (former Ice Streams A and B), and of Bindschadler and MacAyeal Ice Streams (former Ice Streams D and E; Reference Bindschadler, Stephenson, MacAyeal and ShabtaieBindschadler and others, 1987; Reference Bamber and BentleyBamber and Bentley, 1994). No ice plain exists, at present, downstream of the now-quiescent lower trunk of Kamb Ice Stream (former Ice Stream C).

Much of the ice involved in past variations at the downstream ends of ice streams has moved into the RIS, carrying with it a record of past flow across the grounding line. Flow features preserved in the RIS reveal a dramatic history of discharge variation. Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others (2000) used broad patterns of distortion to streaklines (their ‘flow stripes’) and crevasses in the ice shelf to chronicle some aspects of the last 1000 years of ice-stream discharge and related changes near the Ross grounding line. We use a combination of increased resolution satellite-based imagery and simple ice-flow models to extend earlier studies of the record of change and grounding-line retreat near the present-day Ross grounding zone. Our study consists of three main strands. First, new satellite-based imagery is used to build on the existing record of events at and near the ice-shelf grounding line (Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others, 2000). Second, an ice-shelf flow model is used to test ice-stream discharge scenarios derived from that record. Third, a thermal model is used to illustrate the mechanism behind the discharge cycle apparent in the ice-shelf record. We conclude that this repeated pattern of changing discharge, shown most clearly as a localized slow-down and redirection of flow, is a necessary consequence of rapid ice-stream flow due to the accelerating flow changing the thermal profile in the ice. The range of discharge routes offered by a broad ice plain buffers the discharge against both rapid drawdown and complete shut-down, and may be important in regulating, at least in the case of the Ross Sea embayment, the process of grounding-line retreat.

2. Imagery

The subtle surface morphology present on the RIS can be observed to a limited extent in many types of visible-band satellite imagery (e.g. Reference Williams and FerrignoWilliams and Ferrigno, 1988; Reference Casassa, Jezek, Turner and WhillansCasassa and others, 1991; Reference BindschadlerBindschadler, 1993). Fortunately, enhancements made using multiple images can greatly improve our ability to see these features (Reference Scambos, Kvaran and FahnestockScambos and Fahnestock, 1999). Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others (2000) found a complicated array of features in a composite Advanced Very High Resolution Radiometer (AVHRR) image of the ice shelf with a pixel size of 1 km. We take advantage of recent advances in sensor capabilities to produce a multiple-image composite from 250 m pixel size visible-band (band 1) images from the moderate-resolution imaging spectroradiometer (MODIS) sensor on NASA’s Terra platform. The MODIS imagery was obtained from the NASA Goddard Space Flight Center Earth Sciences Distributed Active Archive Center (DAAC) (http://daac.gsfc.nasa.gov/DAAC_DOCS/gdaac_home.html).

At any location, the composite includes up to 10 images spanning a few-week period in November and December 2001, using multiple local solar times from 1000 to 2200h GMT. This range of image acquisition times provides a range of Sun azimuths in the images used. The images were acquired from the DAAC as individual MODIS granules, which were individually remapped into a polar stereo-graphic grid at a pixel size of 250 m, hand-masked for clouds using polygonal areas (the many straight line segments cutting through the composite image are the edges of the polygons used in this masking), and high-pass filtered with a 25km radius to remove the effects of larger-scale slope variations and allow contrast enhancement for smaller-spatial-scale slope variations (subtle topographic ridges, with a short dimension on the order of 1km, are visible). Simple summing of the resultant images produced a composite with a wide (1808) range of solar azimuths that allows curving features and features of any orientation to be well represented (Fig. 1). The range of solar azimuths (illumination directions) requires the viewer to exercise some caution in interpreting the slope directions of the features that are seen in the composite image. A large fraction of the included images are illuminated from the upper left corner of the figure, resulting in slopes that are downward toward the bottom and/or right edge of the figure being dark. A good example of this is the dark line produced by the slope down on to the ice shelf at the grounding line of Kamb Ice Stream. The improved resolution, geolocation and image availability of MODIS relative to the older AVHRR data allow us to show the full extent of features that were only partially visible before.

Fig. 1. Part 1. Composite MODIS image of the RIS and location map with geographic place names. Top right: The black dot in Whillans Ice Stream marks the approximate location of the UpB borehole (about 83828^° S, 138811^° W), and the heavy grey line marks the approximate track of the thermal model. The Kamb grounding line and coast of Siple Dome are collectively called the ‘Siple Coast’, and the Mercer/Whillans grounding zone is the ‘Gould Coast’.

Fig. 1. Part 2. Detail from the composite image showing bowed streaklines that reveal variations in the flow of ice near the grounding line. (a) Tightly folded streaklines indicate recent change on Mercer Ice Stream. (b) Several cycles of flow redirection have occurred in the Mercer/Whillans ice-plain region. (c) Rifts opening at the Kamb Ice Stream grounding line, a stick–slip boundary. (d) Loops downstream of Crary Ice Rise formed on its northern margin and are probably related to the distorted streaklines near the ice-shelf calving front shown in (e) (see also Fig. 3). (f) Loops adjacent to the southern side of Steershead Ice Rise may have been initiated in a former ice plain (see Fig. 2). The labeled geographic locations are SD = Siple Dome, UpB = approximate location of the UpB borehole, NIR = possible new ice rise near the mouth of Mercer Ice Stream, RI = Roosevelt Island, SH = Steershead Ice Rise, and CIR = Crary Ice Rise.

The possible confusion introduced by multiple solar azimuths used in this process is made up for by the greatly enhanced ability to track features as coherent structures for hundreds of kilometers across the ice shelf. Note that the effect of high-pass filtering is to remove the shading produced by slopes across broad features such as Siple Dome and Roosevelt Island (and also Earth curvature); this means that the shading variations in the composite image are dominantly controlled by local slope variations that would be overwhelmed in a simple shaded-relief map of the area.

3. Interpreting Satellite views of Ice-Plain/Ice-Shelf Flow

3.1. Streaklines

In the interpretation of streaklines, misoriented with respect to present-day flow, in terms of past flow events we must first consider how those distortions arise. A model of ice-shelf flow can then be used to constrain the range of time-varying boundary conditions that could have led to the observed distortions. Streaklines may show significant curvature not present in the modern flow field of the ice shelf for two reasons: either the flow direction has varied substantially in the past, or the features have been deformed by shear past a no-slip–slip boundary (Fig. 2). The former can produce broad regions of misalignment between streaklines and modern flowlines while the latter can only occur over limited spatial scales. Shear stresses are naturally large near no-slip–slip boundaries such as the edges of ice streams or ice rises. Reference RaymondRaymond (1996) elegantly demonstrates that the lateral distance over which the resulting shear strain diffuses (the width of the shear boundary layer) depends on the magnitude of the slip resistance at the base of the moving ice. The smaller the slip resistance, the larger the diffusion distance. But even in the case of the nearly zero slip resistance under floating ice, boundary layer width would be at most a few tens of kilometers, with highest rates of shear concentrated near the boundary. Away from sites of no-slip–slip boundaries, rotation rates in the floating ice shelf are small and deformation that would bow streaklines is difficult. For shear to be distributed more widely would require spatial heterogeneity in the basal slipperiness (as appears to be the case on Echelmeyer Ice Stream (Reference RaymondRaymond, 1996)). As we will show using an ice-shelf model in the next subsection, large-radius streakline bows indicate a large redirection angle or a complicated basal boundary condition over a large area, while tight folds indicate significant shear accompanying a smaller-angle flow redirection.

Fig. 2. Sketches of streakline-distortion scenarios. Dark grey indicates grounded, slow-flowing ice, and light grey indicates a lightly grounded ice plain. Heavy black lines indicate lateral margins observed as surface slope breaks, while heavy dotted lines represent past shear margins inferred from crevasses. Light black lines represent streaklines, and the arrows indicate the ice-flow direction. In scenario (a) the regional flow remains nearly constant over time; the ice rise forms, and streaklines are deformed by extremes near the ice-rise margin. In scenario (b) the stream discharges onto an ice plain, and the past flow direction differs from the present direction. The ice rise forms when part of the ice plain becomes dewatered by basal freezing and becomes an obstacle to ice-stream discharge. Streaklines are deformed by the resulting change in flow direction and perhaps also by shear near the ice rise. Distortion would be enhanced over an ice plain with heterogeneous basal traction, a condition that allows lateral shear to be widely distributed.

3.2. Ice-shelf model

We test hypotheses regarding the formation and modification of passive ice-shelf flow features using a map-plane, finite-element model of ice-shelf flow coupled with a Lagrangian tracer routine. The map-plane model generates a velocity field that changes according to changing boundary conditions; this changing velocity field is used to move tracers through the model domain. When connected, sequential positions of the tracers make flow-feature-like traces that can be compared with the streaklines observed in satellite imagery of the ice-shelf surface.

Ice shelves flow by gravity-driven horizontal spreading that transports ice from the grounding line to the calving front. The location of the grounding line is determined by flotation, where the weight of the ice is balanced by buoyancy. In some areas a lightly grounded ice plain, that may be affected by the rise and fall of the tide, lies between the floating shelf and the well-grounded sheet. Resistance to ice-shelf flow is provided by lateral shear at bay walls and ice-rise margins, and by compression upstream of ice rises. It is clear from the observational record that the RIS has experienced repeated changes at or near the grounding line. These may be due to changes in ice-sheet discharge across the grounding line, in the position of the grounding line or an ice plain, or in the geometry of an ice rise.

Changes in ice-stream discharge affect the flowing shelf in two different ways. Changes in boundary velocity propagate nearly instantaneously through the ice shelf while changes in mass flux propagate on the advective time-scale (Reference MacAyeal and LangeMacAyeal and Lange, 1988). Past volume flux is unknown, and we limit ourselves here to the kinematic part of the problem. It may be possible to use the the kinematic history we produce to retrodict ice-stream-mass discharge.

The ice-shelf model derives from the more complete ice-sheet model presented by Reference Hulbe and MacAyealHulbe and MacAyeal (1999). Ice-shelf flow is described by a set of stress-balance equations, simplified by the typical assumption that horizontal flow is depth-independent. The Glen flow law is embodied in an effective viscosity

(1)

in which is an ‘inverse’ rate factor and α is an enhancement parameter, u and v are the two horizontal components of velocity in the x and y directions of a Cartesian coordinate system, and n is the flow-law exponent, taken here to be 3. The rate factor varies in a simple way. Over most of the shelf, ice stiffness is taken to be homogeneous, with a value that replicates the broad-scale flow of the shelf (6.6 × 10⁷ Pa s^1/3 cf. Reference MacAyeal, Bindschadler and ScambosMacAyeal and others, 1995; Reference Hulbe and MacAyealHulbe and MacAyeal, 1999). Variations are imposed near certain boundaries by setting the parameter α to a value other than unity, as discussed in the next paragraph. The stress-balance equations are derived by Reference MacAyeal and ThomasMacAyeal and Thomas (1982) using the dynamic boundary conditions of a stress-free top surface, a zero shear-stress basal boundary, and a depth-integrated sea-water pressure applied along the seaward front of the ice shelf. A no-slip condition is specified at no-flow boundaries such as bay walls, the downstream ends of interstream ridges and the edges of ice rises. The stress boundary conditions may be modified by the application of a non-zero basal shear stress under an ice plain, or by ice softening within shear margins. The former is accomplished by incorporating a basal stress term that is the product of the horizontal velocity and a (scalar) friction parameter into the stress balance (cf. Reference Hulbe and MacAyealHulbe and MacAyeal, 1999). This is equivalent to treating the area as a small basal-friction case of ice-stream flow. The friction parameter used for the ice plain is 10^-7 Pa s, two orders of magnitude smaller than the typical ice-stream values suggested by Reference MacAyeal, Bindschadler and ScambosMacAyeal and others (1995) and Reference Hulbe and MacAyealHulbe and MacAyeal (1999). Shear-margin softening is accomplished by manipulating the effective viscosity near large-shear boundaries, as discussed in the next paragraph. Ross Sea embayment and ice-shelf geometry are taken from the BEDMAP database and its antecedents (BAS, 2001).

The infrastructure of the ice-shelf model has been demonstrated to correctly reproduce the broad-scale features of RIS flow (Reference MacAyeal, Rommelaere, Huybrechts, Hulbe, Determann and RitzMacAyeal and others, 1996). However, several special boundary conditions are required to accurately simulate the pattern of flow around Crary Ice Rise Fig. 3. First, to promote ice flow through the narrows between Crary Ice Rise and the Transantarctic Mountains (TAM), ice softness must be enhanced in shear margins along the coasts of both features. This is reasonable, given the large degree of crevassing along Crary’s margins and the observation that ice from glacier shear margins advects downstream along the TAM. In the model, the ice is softened by arbitrarily adjusting the value of α in Equation (1). A reasonable agreement with observed velocity is achieved with α = 0.2 in narrow (2-3 km) bands along the TAM coast from 170^°W to 180^°W and about 15 km outboard of the southern side of Crary, as observed in satellite imagery of the region and in the Ross Ice Shelf Geophysical and Glacio-logical Survey (RIGGS) velocity measurements (Reference Thomas, MacAyeal, Eilers, Gaylord, Hayes and BentleyThomas and others, 1984; Reference BindschadlerBindschadler, 1993). A shear margin with α = 0.5 is specified along the northern margin of the ice rise (see section 3.5). Second, to reduce what the model would otherwise simulate as a very fast discharge north of Crary Ice Rise, a spatially uniform basal friction is applied at the base of the ice plain as discussed above. As a result, basal shear stress takes up 40-65% of the driving stress over most of the Mercer/Whillans ice plain. The ice-shelf geometry also imposes a technical constraint on the model, as mesh resolution must be fine enough to resolve across-flow stress gradients across the narrows. Although of little consequence to the present study, a narrow shear margin with α = 0.5 also appears to be required around the margins of Roosevelt Island. Other modeling studies, which fail to account for shear margins or the ice plain, or to adequately resolve the narrows between Crary and the TAM, fail to match measured ice speed in this area (e.g. Reference MacAyeal and ThomasMacAyeal and Thomas, 1986). Improvements to the ice-shelf flow pattern could also be made with a more complicated spatial variation in the rate factor, but for our purpose the manipulations yield equivalent results (a discussion of this issue can be found in Reference Scambos, Hulbe, Fahnestock and BohlanderScambos and others, 2000).

Fig. 3. Simulated present-day flow of the RIS (a) without and (b) with shear-margin softening and basal traction in the Whillans ice-plain area. Geographic names are as given in Figure 1. The contour interval is 100 m a^–1. The flow field in part (b) is the ‘present-day’ field for the streakline simulations. Details of the modifications are provided in the text. The important difference for this study is the flow of ice around Crary Ice Rise. Without margin softening, flow between Crary Ice Rise and the TAM is too slow. Present-day discharge speeds for ice stream and TAM glacier outlets into the ice shelf are reported in ma^–1, and the sections of model boundary where the discharges are applied are marked with heavy grey lines. Locations of shear margins, simulated by softening the flow-law rate factor, are drawn in white on part (b). Crary and Steershead Ice Rises are represented by nodes in the model domain so that those areas may be treated as either floating or grounded, as in the Steershead experiment discussed in section 3.3. Roosevelt Island is not included in the model domain, and a no-slip boundary condition is applied around its coastline. This model is used in tracer experiments that attempt to reproduce streaklines observed in the present-day ice shelf.

Thus configured, the model is used to conduct many streakline simulations. Three examples from Figure 1 are discussed here to illustrate our interpretation of the observed ice-shelf features.

3.3. Steershead Ice Rise

Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others (2000) suggested that the widely bowed streaklines south of and abutting Steershead Ice Rise (Fig. 1f) formed by shear deformation after the ice rise grounded. In their proposed scenario, ice-shelf flow over the Steershead region was arrested due to some change in ice thickness at the rise location. Subsequent flow past the obstacle deformed the streaklines by shear (Fig. 2a). The scenario seems unlikely when the shear zone width is compared with the spatial scale of the bowed streaklines (20km or so and nearly 100 km, respectively; see Fig. 1f) but it is possible that over a sufficiently long time the rotation could be accomplished. We test the misorientation-by-deformation hypothesis by reproducing Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others’ (2000) sequence of ice-stream discharge and ice-rise grounding events in the ice-shelf model and advecting tracers through the time-varying flow field. Ice-rise grounding is simulated by imposing a no-slip basal boundary condition at model mesh nodes within the ice-rise area. The changing boundary conditions are summarized in the caption to Figure 4. As may have been anticipated, streaklines are deformed only within several ice thicknesses of the ice rise (the thickness is about 450m in this area). This is a due to the small lateral shear strain rates away from the ice rise and the limited time, about 350 years (Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others, 2000), over which the deformation must occur. Thus, the unusual orientation of the streaklines cannot be the result of passive deformation of the floating ice.

Fig. 4. Model simulations of streakline deformation adjacent to ice rises. Light black lines represent simulated streaklines at the end of the model run. Significant streakline folding is limited to within several ice thicknesses of an ice rise. The simulations run 2000 years, with the following boundary conditions and events: the discharge speeds (m a^–1) of the ice streams are Mercer: 270, Whillans: 800, Kamb: 650, Bindschadler: 650 and MacAyeal: 550; Kamb Ice Stream slows to 10ma^–1 100years ago; Whillans and Kamb Ice Streams narrow at 540 and 340 years ago, respectively; Crary Ice Rise forms 1100 years ago and develops shear margins 200years later; Steershead becomes an ice rumple 540 years ago and grounds firmly 340 years ago; TAM glacier discharge is constant, at observed rates, over the entire simulation. Discharge gates are shaded grey, and small arrows indicate lateral margin jumps at the times noted here. The inset shows streakline folding prior to shear-margin development at Crary Ice Rise. This is the process now underway on Mercer Ice Stream (Fig. 1, part 2). The model domain includes the Whillans ice plain (Fig. 1).

Two past conditions can account for the bowed streaklines south of Steershead. The ice-shelf flow direction in this area may have been significantly different in the past than it is at present and the area may have been a lightly grounded ice plain, probably with a spatially heterogeneous pattern of basal slip. We suggest that the observational record implies both conditions. Both Kamb Ice Stream and the northerly track of Whillans Ice Stream discharged ice into that plain, with the result that flow from Kamb would have been directed more to the north than at present. A more northerly route for Kamb Ice Stream is also indicated by the abandoned shear margin along the downstream, southern side of Siple Dome (Reference Jacobel, Scambos, Nereson and RaymondJacobel and others, 2000). About 350 years ago, Steershead became more firmly grounded than the surrounding ice plain, and ice discharging from Kamb was forced to divert southward, to a direction similar to the present-day flow direction. The result is a misorientation of relict streaklines with respect to the present flow direction. Transverse shear, distributed across the plain south of Steershead as the discharge direction changed, would have promoted streakline distortion, but redirection in a floating ice shelf would not (cf. the Steershead experiment in Fig. 4). We further suggest that, without a supply of new ice into the ice plain, the grounding line must have later retreated through the region. (The evidence is in the form of ice-shelf crevasses, a topic that is revisited in section 3.6.)

3.4. Mercer Ice Stream: active shear deformation

The new MODIS imagery reveals a pattern of streakline deformation not previously observed: relatively tight folds that resemble recumbent folds seen in other rocks (Fig. 1a and Fig. 5). In the narrow outlet to Mercer Ice Stream, the distance over which the streaklines are distorted is comparable to the typical ice-stream shear boundary layer width. Thus, transverse shear deformation resulting from a more modest change in flow direction than implied in the Steershead case can account for the observed distortion. The locus of the redirection (NIR in Fig. 5) is a site where the bed is relatively high, the ice is relatively thin and present-day flow is very slow (cf. the BEDMAP database, BAS, 2001; Reference Joughin and TulaczykJoughin and Tulaczyk, 2002). South of the now slow-flowing ice, the grounding line is embayed. In both the Steershead and Mercer locations, it appears that streaklines in ice formerly flowing over a now-stagnant area were pinned when the ice slowed and that subsequent redirection of ice flow has distorted the streaklines. The difference between the two lies in the angle through which the flow redirects. For Mercer, the angle is small and the streakline folds are tight. For Kamb, the angle is large and the streaklines bow widely.

Fig. 5. Ongoing flow redirection and streakline deformation on Mercer Ice Stream. (a) Composite MODIS image with labeled geographic features: IP, lightly grounded ice plain; NIR, possible new ice rise; RIS, floating Ross Ice Shelf. (b) Traced streaklines. (c) Past flow pattern of Mercer Ice Stream. (d) Present-day flow redirected around an apparent obstruction.

4. A Possible Thermal Mechanism behind Episodes of Flow Redirection and Grounding-Line Retreat

We have argued, using satellite imagery interpreted with the help of an ice-shelf flow model, that evidence of episodic flow redirection and grounding-line retreat is found through-out the area covered by Figure 1a-f. Redirection is actively occurring on Mercer and possibly Whillans Ice Streams. Past episodes are evident in the discharge history of Kamb Ice Stream around Steershead Ice Rise and of Whillans Ice Stream around Crary Ice Rise. In each case, an ice rise is the locus of flow redirection. We suggest that the stagnant ice involved with Mercer’s ongoing redirection event will become an ice rise if grounding-line retreat continues. The simplest explanation we have found for this repeated series of events is based on the downstream thermal evolution of ice near the bed of a fast-flowing ice stream. In this section we illustrate what we mean by the downstream thermal evolution of an ice stream using a simple model of temperature in a ‘column’ of ice advecting down Whillans Ice Stream. Its implications are discussed in section 5.

The ice within an ice stream experiences rapid change as it moves downstream. Upstream, the ice is relatively thick, and thermal conditions favor the production of meltwater. Subglacial sediment prevents adequate drainage from developing, basal water pressure is large and the resulting rapid sliding is constrained largely through drag at the ice-stream margins. As the ice moves downstream, fast longitudinal stretching thins the ice and steepens the basal temperature gradient. In some locations, stream widening also leads to thinning, and a reduction in margin control on the flow. In the thin, downstream reaches of an ice stream, an energy budget that favors sub-freezing basal temperatures is likely (cf. Reference RoseRose, 1979; Reference Parizek, Alley and HulbeParizek and others, 2003).

The downstream trend toward basal freezing in an ice stream due to ice thinning and consequent steepening of the basal temperature gradient can be demonstrated quite simply in a Lagrangian frame. We take as our example Whillans Ice Stream, where a considerable record of flow redirection exists. The thermal evolution of ice between the location where ice temperature is known, the UpB borehole Fig. 1, and the ice-plain grounding line is tracked using a vertical energy-balance model that simulates the downstream progress of an imaginary column of ice.

The vertical temperature profile at any geographic location within the ice sheet depends on boundary conditions at that location, and the flow and boundary condition history of the ice upstream of that point. Our model embodies both by placing the energy-balance equation in a Lagrangian frame. Starting from a known thermal state, boundary conditions on a one-dimensional (vertical) domain are varied such that they simulate advection of the model through the three-dimensional ice sheet. We take as an initial condition the observed depth variation of ice temperature at the UpB borehole (Reference Engelhardt, Kamb and FahnestockEngelhardt and others, 1990), automatically capturing the upstream history of the ice. Boundary conditions on the model are set using modern values, an acceptable simplification over the approximately 300 years required for a column of ice to move from UpB to the downstream edge of the ice plain. Those conditions are: surface accumulation rate from Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others (1999); ice surface temperature based on the datasets of Reference ComisoComiso (1994); and a geothermal flux of 0.072Wm^-2 (based on H. Engelhardt’s Siple Dome measurements) at the base of a bedrock layer.

The temperature evolution of a column of ice advecting downstream can be expressed:

(2)

in which T represents the temperature at any depth, and t, w represents ice motion in the vertical (z) direction, k _i represents the thermal conductivity of ice (2.1WmK^-1), W represents heat generated by shear strain within the ice, ρ represents the ice density (910kgm^-3), and c represents the heat capacity of ice (2009J kg^-1 K^-1). It is assumed, reasonably, that horizontal temperature diffusion is small. Within an ice stream, the vertical shear strain rate is nearly zero and lateral shear strain heating is also small (Reference Hulbe, Wang, Joughin and SiegertHulbe and others, 2003), so the last term in Equation (2) is neglected. Boundary conditions on the vertical velocity are the surface accumulation rate and the basal melting rate (discussed below).

Mass continuity and advection both derive from modern observations. In the case of continuity, this means we assume

(3)

in which is the flow-following horizontal coordinate and û represents the advective speed at any t. The ice thickness, h, is from the BEDMAP digital elevation model (BAS, 2001). The advective speed is calculated using velocity from Joughin’s RADARSAT synthetic aperture radar (SAR) inter-ferometric dataset (Reference JoughinJoughin and others, 1999; personal communication from I. Reference Joughin and TulaczykJoughin, 2002) and from Reference Whillans and van der VeenWhillans and Van der Veen (1993). The assumption of constant thickness and flow speed is acceptable for the short transit time of ice through the boundary condition domain.

The advecting ice is coupled, thermally, to the bed by the heat balance at the interface between the two (following Hulbe and MacAyeal, 1995). Heat flux out of the ice toward the interface is:

(4)

in which k _i represents the thermal conductivity of the ice and δT/δz is evaluated in the two vertical nodes above the interface. Heat flux into the base of the ice is

(5)

where δT/δz is evaluated in the two vertical nodes below the interface, k _r represents the thermal conductivity of the bedrock (3.3 W m K^-1), τ _b represents the basal shear stress, L _f represents the latent heat of fusion of fresh water (3.35 × 10⁵J kg^-1), and b represents the basal melting (freezing) rate. The situation would be modified if sea water infiltrated the basal water system upstream of the grounding line. That could happen as a result of tidal pumping, though the ice-sheet pressure gradient will always tend to flush fresh water downstream. The basal melting rate itself depends on the heat balance at the ice/bed interface:

(6)

If more heat is diffused upward into the ice than is generated at the interface (H _i > H _o), water at the interface must freeze (> 0), and, where the opposite is the case, basal ice melts. In the model, the interdependence of H _i and is accommodated by evaluating the interface heat balance in a sequence of logical steps that make use of temperature gradients at the last location (cf. Reference Hulbe and MacAyealHulbe and MacAyeal, 1999). The validity of this strategy is tied to the validity of our assumption that ice travels quickly through the model domain. It is assumed that where the tendency is toward freezing, water is available.

Frictional heating at the interface between the ice and the bed over which it slides, ûτ _b, has a strong influence on the magnitude of basal melting or freezing. Where basal ice is at the melt temperature, frictional heat goes to melting that ice, but it cannot warm ice already at the melt temperature. Where the tendency is toward basal freezing, heat from friction reduces the freezing rate. In an ice stream, the basal shear stress is some fraction of the gravitational driving stress. Ideally, the computation of basal shear stress would include a thermomechanical model of the subglacial till. Here, we simply partition resistance to the gravitational driving stress between basal shear stress and lateral drag at the (distant) stream margins. Reference Whillans and van der VeenWhillans and Van der Veen (1993) find that the partition is nearly 0:1 in their ‘photo block’ area near the UpB camp. In the model, basal friction is computed using a basal shear stress that varies from 0.1 to 0.5 times the driving stress. In the small-friction (0.1) case, basal freezing sets in almost immediately downstream of UpB. In the case of large basal friction, melting persists nearly to the upstream end of the ice plain, despite a strong tendency toward freezing due to ice-stream thinning. Frictional heat production at the downstream end of the model domain (the ice plain) is small regardless of the partitioning because surface slope, and thus driving stress, is small.

We assume that at each new time-step the interface between the ice and its bed is initially melted. That is, the ice stream is moving over a wet bed with a temperature equal to the local pressure-dependent melt temperature. This is reasonable: water currently exists at the bed of Whillans Ice Stream; ice streams appear to flow fast for millennia, implying a persistent basal lubricant; and without a wet bed the ice stream would cease to slide and our simulation would be at an end. A model designed to predict such slowing would include a depletable basal water reservoir and related till hydrology.

The geothermal flux is used to compute the temperature gradient in a sub-ice bedrock layer with a thermal conductivity typical of sedimentary rocks (3.3WmK^-1). This provides a buffer for the basal ice temperature that tends to mitigate basal freezing. While basal ice is allowed to cool, the assumption of some basal water supply prevents the modeled ice from ‘freezing to its bed’, with the consequence that the accuracy of the geothermal gradient is of minor importance.

The expression for heat balance is solved using a finite-element method. Equation (2) is represented as piecewise linear functions along a one-dimensional, variable-resolution, 26-node domain. Model resolution is greatest near the ice/bed interface, where the temperature gradient is the largest. Node spacing scales to the ice thickness. The trajectory of the ice-column model through the boundary condition space is determined using a combination of observed velocity and, where velocity is missing, surface slope Fig. 1. Once the trajectory is determined, gaps in velocity observations are filled by linear interpolation between known values.

The thermal evolution of ice flowing from the UpB borehole site to the downstream end of the ice plain is depicted in Figure 6. The geographic starting point is selected such that model trajectory lies close to the kinematic center line of the ice stream. The starting point is within a few kilometers of the borehole location. At the UpB location, the ice/bed interface tends to melt, but the high thinning rate, about 50% over the 500km trajectory, dramatically steepens the basal temperature gradient, despite enhanced thermal diffusion that opposes the steepening. Thus, with distance downstream, the ice begins to consume latent heat from basal meltwater (generated at an earlier time or elsewhere in the system). This is simply a consequence of the mismatch between time-scales for the two processes. The point at which basal freezing, and basal water consumption, commences, scales to the proportion of the driving stress taken up by basal drag. But the fundamental point is clear: the rapidly flowing ice stream thins downstream and thus tends toward freezing and the consumption of basal water. In the case of Whillans Ice Stream, basal water consumption may begin, intermittently, as far upstream as the UpB camp and is well established, at a rate of about 0.002ma^-1, by the time the ice stream discharges onto the ice plain. The amount of latent heat the thinning ice must extract from basal water is on the same order as the heat supplied by the geothermal flux (0.072Wm^-2).

Fig. 6. Simulated change in the temperature of a column of ice advected down Whillans Ice Stream, starting from the UpB borehole and ending at the downstream edge of the ice plain. (a) Ice-sheet geometry along a trajectory that includes the model’s flow-following particle path. The UpB borehole, at which the model is initialized, is located by a cross. (b) Surface speed along the profile, from SAR (diamonds) and surface-based measurements (squares). The solid line represents advective speed used in the model calculations. (c) Modeled basal freezing rate for simulations where the basal shear stress is 0.1 times (solid line) and 0.5 times (dashed line) the driving stress. A positive value corresponds to freezing of basal water. The inset panel shows the vertical temperature profile at selected locations along the simulation path for a basal shear/driving stress ratio of 0.1:1. The heavy dashed lines represent the initial (UpB) and final (grounding-line) states of the model. Lighter weight lines between those bounds plot the temperature gradient at intermediate locations; solid lines indicate a tendency toward basal melting, and dashed lines indicate basal freezing. The mean time between calculations is about 4.5 years. For clarity, only every fifth calculation is plotted. These are not steady-state temperature profiles.

Where the fast-flowing ice stream thins the most, latent-heat demand at the ice/bed interface is the greatest. If the basal water system cannot supply enough water to meet the heat demanded by the ice, the local water source will be depleted through freezing, with the probable result that traction between the ice and till will increase. Larger traction could buffer the tendency to freeze by increasing frictional heating. However, even in simulations where basal shear stress takes up half the driving stress, the tendency toward basal freezing dominates in the ice-plain region.

5. Discussion: The Thermal Explanation for Flow Variability and Redirection

Many sets of bowed streaklines and transverse crevasses, revealed by new MODIS images, indicate significant reorganizations of flow and grounding-line retreat along the Siple and Gould Coasts Fig. 1. Some of these features were identified by Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others (2000), who used them to infer a 1000 year history of Ross Ice Stream discharge events. Here, we have used that history as a framework upon which to build model experiments to study the origin of the bowed (deformed) streaklines. Those experiments, and an understanding of the downstream evolution of the basal temperature gradient in a fast-flowing ice stream, lead us to propose a simple thermal mechanism that causes variations in the rate of flow that may determine the pattern of grounding-line migration in the Ross Sea embayment.

As the energy balance at the interface between the ice stream and its bed changes with distance downstream, so too does the stream’s ability to slide fast. The West Antarctic ice streams rest on unconsolidated marine sediments (Reference Studinger, Bell, Blankenship, Finn, Arko and MorseStudinger and others, 2001) that have been shown in the laboratory to deform rapidly when water-saturated (Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others; 2000; Reference Kamb, Alley and BindschadlerKamb, 2001). In areas where basal meltwater is present, the till inhibits its drainage, leading to large basal water pressure (Reference Engelhardt, Kamb and FahnestockEngelhardt and others, 1990). Thus, the saturated till promotes rapid ice motion despite small driving stresses, but the stream must slow when there is insufficient meltwater, either because the loss of water changes till porosity and thus till strength, or simply because the interface between ice and bed is no longer lubricated. Even without a comprehensive model of till thermomecha-nics, it is clear that where basal meltwater is abundant and distributed, ice streams will flow fast, and where basal meltwater is scarce, they will not.

As we have shown, heat flow in the basal ice, determined by the basal temperature gradient, controls the tendency to melt or freeze water at the base of the ice sheet, and thereby must in some way regulate ice-stream flow. Smaller temperature gradients limit the diffusion of heat from the interface between the ice and its bed and promote meltwater production. Larger gradients enhance the diffusion of heat from the interface. Where the heat flux up into the ice exceeds the heat flux to the base of the ice from geothermal and mechanical sources, the basal ice must chill or, if liquid water is available, the deficit must be balanced by latent heat from that source. That is, basal freezing must occur. Where the meltwater source becomes depleted, basal traction will increase and the ice stream will no longer be capable of rapid motion.

As ice flows downstream through a relatively confined channel, the rate of longitudinal stretching increases (Reference Alley and WhillansAlley and Whillans, 1991). Continuity requires a commensurate rate of ice-stream thinning. The more vigorous the ice stream, the faster the rate of downstream extension and thus of downstream thinning (the rate is ~-10^-3a^-1 at the UpB camp on Whillans Ice Stream; Reference Hulbe and WhillansHulbe and Whillans, 1997). In a steady-state ice stream, where most of the motion is due to basal sliding, downstream thinning alone must cause the basal temperature gradient in the advecting ice to steepen. Similarly, if the ice sheet thins or thickens at a fixed geographic location, the basal temperature gradient must change accordingly. Other conditions, such as the surface air temperature and accumulation rate, geothermal heating, dissipation of shear strain energy, and frictional heat generation at the ice/bed interface, add spatial variations to the simple trend.

Ice is relatively thin (the bed is relatively high) at every site around which we infer ice flow to have been redirected Fig. 7. The streakline record indicates that ice once flowed over these now-bypassed locations. Flow-driven thinning of the ice leads to a steepening of the basal temperature gradient. Where the gradient is sufficiently steep, basal freezing will commence. Over time, the thinning, fast-flowing ice may deplete its basal lubricant source, with the result that basal shear stress will increase and the ice stream will decelerate. The decelerating region becomes an obstacle to fast flow, causing discharge redirection. Thus, it may be inevitable that the fastest-flowing regions within the thin, downstream ends of ice streams slow to become the loci of flow redirection.

Fig. 7. Bed elevation and loci of flow redirection. Bed elevation, from the BEDMAP dataset (BAS, 2001), is contoured at a 100 m interval. Sites of flow redirection are shaded dark grey. The image is a subset of the MODIS composite in Figure 1.

An ice-plain kinematic cycle emerges from the combination of our conceptual framework and the RIS record: fast ice-stream flow; local deceleration within the fast-flowing region which creates an obstruction; flow diversion around the obstruction; and eventual grounding-line retreat past the obstruction (producing an isolated ice rise). Some areas that slow may be reactivated by margin migration. Reference RaymondRaymond (1996) has shown that margins deposit heat outboard of the stream, and Reference Echelmeyer and HarrisonEchelmeyer and Harrison (1999) measured active margin migration on Whillans Ice Stream. With either grounding-line retreat or reactivation by margin migration, a complete cycle of flow variability (in both speed and direction) can be established on an ice plain. In a regional sense, the cycle would be seen as repeated episodes of slowdown, flow redirection, and retreat of the grounding line. This view is in contrast to the two other possible scenarios of grounding-line retreat: single episodes of deceleration and stabilization (which seems unlikely in a system that has been retreating for an extended period) or a continuous, gradual retreat of the grounding line as the ice in the interior is drawn down (which cannot explain the complex patterns of flow recorded on the RIS or measured changes in flow speed).

The recent stagnation of Kamb Ice Stream can be understood within this framework. The bathymetry Fig. 7 shows a very shallow sea-floor in the area bounded by the present grounding line, Crary Ice Rise and Steershead Ice Rise. The extensional features in the floating ice in this area suggest a recent ungrounding that occurred across a retreating grounding line. These features also suggest that the grounding line in this area divided ice with large basal traction from floating ice, most notably in the area just to the north of Crary Ice Rise. This raises the possibility that the shut-down of Kamb Ice Stream began on an ice plain that once filled the Kamb embayment. As the Kamb grounding line retreated past Steershead to its present location, the downstream end of the ice stream would have become constricted between the large interstream ridges that have apparently persisted for millennia (Reference NeresonNereson, 2000). Without a broad ice plain, avenues for redirection are limited or none-existent.

Recent changes at the mouth of Whillans Ice Stream, which has slowed over the last three decades (Reference Joughin and TulaczykJoughin and Tulaczyk, 2002; personal communication from P. Vornber-ger, 2002) also fit within our kinematic cycle framework. R. Bindschadler and others (unpublished information), using a combination of velocity and surface elevation measurements, find that the recent changes are spatially heterogeneous, with the discharge azimuth shifting southward, around the southern side of Crary Ice Rise, as ice just upstream of the ice rise slows. In the context of the observational record of repeated flow redirections across the Whillans and Mercer ice plain over the last 1000 years, the recent changes may be interpreted as part of an ongoing discharge cycle, rather than a shut-down of that system.

In each case discussed here, the locus of flow redirection is an area where the bed is relatively shallow Fig. 7. According to our hypothesis, the most important aspect of this correlation is that fast-flowing ice is the thinnest and thus has the steepest basal temperature gradient in these areas. It may also be that the hydraulic potential gradient due to the associated thickness and surface gradients would tend to push basal water away from the area, limiting the latent-heat source. Relatively shallow shoals that extend from the continental shelf edge back to the present-day grounding line (cf. the BEDMAP bathymetry, BAS 2001) could have served as a platform for such events throughout the last 7000 years, allowing the gradual retreat implied by the geologic record (Reference Conway, Hall, Denton, Gades and Wad-dingtonConway and others, 1999). If this is the case, the observed retreat may have occurred at the maximum rate possible, given the thermal limiting of rapid flow we have discussed here. The location of the present grounding line may in large part be determined by the thermal effects of the flow-driven thinning in the lower reach of the ice streams. This mechanism has the potential to mitigate catastrophic retreat.

Discharge that has a history of this episodic variation, such as the record on the RIS shows, may well be an indicator of an ice-stream/shelf system that is in rapid retreat, as ice must thin to the point of freezing-on in order to produce the flow variability that is seen. Thermally driven variability can also be considered in the converse: if the system were stable or advancing, ice would not thin to the point of freeze-on, and redirection and flow variability would not occur. The well-behaved streaklines on the Ronne-Filchner Ice Shelf may indicate that the ice streams flowing into that system do not have actively retreating grounding lines.