2.1. Velocity data

We produced a map of ice-flow velocity (Fig. 1) for Pine Island and Thwaites Glaciers using the same interferometric synthetic aperture radar (InSAR) and speckle-tracking methods used in earlier studies (e.g. Reference JoughinJoughin, 2002). Unlike lower-accumulation areas in West Antarctica where RA-DARSAT works well (e.g. Reference Joughin, Tulaczyk, Bindschadler and PriceJoughin and others, 2002), this instrument’s 24 day repeat cycle yields poor coherence for velocity mapping in the high-accumulation areas along the Amundsen Coast. As a consequence, we created our map using data acquired by the European Remote-sensing Satellites (ERS-1 and -2), which provide image pairs separated by 1 and 3 days. Because Pine Island Glacier sped up between 1994 and 1996 (Reference Rignot, Vaughan, Schmeltz, Dupont and MacAyealRignot and others, 2002), we used only the 1996 data in the areas where large changes occurred. Hence, our map represents the 1996 speed of the Amundsen Coast glaciers.

Fig. 1.

Flow speed (color) over the artificially shaded surface (grayscale) of a DEM produced from a combination of laser and radar altimetry (see text) for Pine Island (PIG) and Thwaites (TG) Glaciers. Speed is also shown with a 50 m a⁻¹ contour (purple), 100 m a⁻¹ contours up to 900 m a⁻¹ (thin black), and 1000 m a⁻¹ contours (thick black). Elevations are shown with 500 m contours (gray). Heavy black curves show the approximate locations of drainage divides.

In areas where there were crossing orbits available and where we could successfully unwrap the phase, the velocity estimates are largely determined by the phase data with formal errors of about 3–5 m a⁻¹. In areas without crossing orbits or phase data, speckle tracking with short temporal baselines yields larger uncertainty (5–50 m a⁻¹). Tidal motion leads to large errors on the floating ice (∼100 m a⁻¹), but, at least on Pine Island Glacier (Reference Joughin, Rignot, Rosanova, Lucchitta and BohlanderJoughin and others, 2003b), it is unlikely to be a significant source of uncertainty even when the observation interval is short (J. Scott and others, unpublished information). In addition to these errors, the process for removing vertical displacements can yield slope-dependent errors of up to ∼3% (Reference Joughin, Kwok and FahnestockJoughin and others, 1998).

2.2. Elevation data

The modelling described in the following sections requires accurate estimates of surface elevation (Reference Joughin, MacAyeal and TulaczykJoughin and others, 2004a). In this study, we use a digital elevation model (DEM) produced using a combination of data from the ERS-1 and -2 radar and Ice, Cloud and land Elevation Satellite (ICESat) laser altimeters (Reference Bamber, Gomez-Dans and GriggsBamber and others, 2008; Reference Griggs and BamberGriggs and Bamber, 2008), which is displayed as a shaded relief map in Figure 1.

2.3. Bed elevation data

During the 2004/05 austral summer, the British Antarctic Survey and the University of Texas conducted a joint airborne survey of the Pine Island and Thwaites Glacier catchments (Reference VaughanVaughan and others, 2006; Reference HoltHolt and others, 2007). Figure 2 shows the bed elevation map from this survey along with the flow-speed contours from the map shown in Figure 1, illustrating the close correspondence of the fast-flow features with the deep subglacial troughs (Reference VaughanVaughan and others, 2006; Reference HoltHolt and others, 2007).

Fig. 2.

Bed elevation for Pine Island and Thwaites Glaciers (Reference VaughanVaughan and others, 2006; Reference HoltHolt and others, 2007). Glacier extents are indicated by the contours of glacier speed: 50 m a⁻¹ contour (purple), 100 m a⁻¹ contours up to 900 m a⁻¹ (thin black), and 1000 m a⁻¹ contours (thick black).

3.1. Temperature and melt rates

We solved for the temperature within an ice column lying atop a homogeneous bedrock slab with depth equal to one ice thickness using methods described by Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others (2004b). In doing so, we neglect horizontal diffusion so that temperature within the ice is determined by the following simplified statement of thermal energy conservation:

where u, v and w are respectively the x, y and z velocity components, c is the thermal heat capacity, k is the thermal conductivity, ρ is the density, and W is the strain heating within the ice column due to vertical shear. In this equation the subscript i is used to indicate ice-related parameters, and in the following equation the subscript r is used to indicate bedrock-related parameters. In the underlying bedrock we solve

which includes only vertical thermal diffusion. Following Reference MacAyealMacAyeal (1997), we solve Equations (1) and (2) numerically after recasting them in contour-following vertical coordinates. Using initial conditions described in the following paragraphs, we integrate this model forward in time over several thousand years with a fixed geometry (see section 4) to obtain estimates of present-day englacial temperature, which we also use to compute basal melt rates (Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others, 2004b).

In solving Equation (1), we used the InSAR measurements for the horizontal velocity components (Fig. 1). Using this velocity field, we classify each node in the model domain as either fast-moving ‘ice-stream’ flow or slow-moving ‘ice-sheet’ flow, with the boundary separating these node types roughly conforming to the 50 m a⁻¹ flow-speed contour. For the fast-moving nodes, we assume there is substantial sliding or bed deformation so that a constant flow velocity can be assumed throughout the ice column. At ice-sheet nodes where the velocity is depth-dependent, we assume no sliding and compute the normalized variation of velocity with depth for internal deformation (Reference PatersonPaterson, 1994, p. 251). This function is then multiplied by the measured surface velocity to estimate the depth-dependent velocities. We compute a simple linear change in vertical velocity with depth for all nodes, with a basal value of zero and a downward surface value equal to the ice-equivalent accumulation rate, which is based on the average of two gridded accumulation maps (Reference Vaughan, Bamber, Giovinetto, Russell and CooperVaughan and others, 1999; Reference Giovinetto and ZwallyGiovinetto and Zwally, 2000).

Since we use the derived basal temperature gradients to estimate basal melt (see below) and, thus, did not have an a priori basal melt estimate, we neglected the downward advection of ice due to basal melt in our temperature calculations. In regions with a frozen bed, this makes no difference. In the majority of our study area where the melt rates are <25 mm a⁻¹, this will bias the magnitudes of our temperature gradients low by ∼15% or less, which is comparable to the effect of errors in the accumulation rate data used in the model. In the fastest-moving areas where the melt rates approach the accumulation rates in magnitude, caution should be exercised in interpreting the temperature gradients. Because nearly all of the heat used to melt ice in these regions is generated by friction from sliding over the bed, the overall effect on our melt rate results is negligible.

We used the same initialization scheme that we employed for our Siple Coast study (Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others, 2004b), which uses an analytical solution for T that accounts for vertical advection and diffusion only (Reference ZotikovZotikov, 1986). Throughout the simulation, the temperature at the ice–bedrock interface is maintained at the pressure-melting point, T _pmp, for the fast-moving ‘ice-stream’ nodes (i.e. we assume an infinite reservoir of basal water so that the bed never freezes). For the slow-moving ‘ice-sheet nodes’, we assume no flow of meltwater from adjacent nodes and no local meltwater storage so that a node freezes to the bed at any time-step when the basal temperature gradient allows upward conduction of heat that can be balanced by the geothermal heat flux. A node that is frozen to the bed switches to a non-frozen state at any time-step where it reaches T _pmp.

As part of the initialization process, the model must determine which nodes are melting and which are freezing. We assume that a wet bed enables the fast motion of the ice streams and their tributaries, so these nodes are always assumed to be at T _pmp. For the ice-sheet nodes, we use the analytical solution for a frozen bed to determine the bed temperature. If the temperature is below the pressure-melting point, the node is flagged as initially frozen. If the temperature is at or above T _pmp, the node initially is considered melted and we initialize temperature by switching to the analytical solution for a thawed bed (Reference ZotikovZotikov, 1986).

The temperature boundary condition at the ice-sheet surface is the mean annual surface temperature (Reference ComisoComiso, 1994, Reference Comiso2000). At the bottom of the bedrock layer, the boundary condition is ∂T/∂z = −G/k _r, where G is the geothermal flux and k _r is the conductivity estimated for bedrock below the zone of any weathered or altered layers. These upper and lower boundary conditions are used when the bed is frozen, to jointly solve Equations (1) and (2). If there is melting at the bed, Equations (1) and (2) are solved individually, with the additional boundary condition imposed at the ice–bedrock interface that temperature equals the pressure-melting point, T _pmp.

We estimate the basal melt rate, m _r, using

where τ _b is the magnitude of the basal shear stress, U _b is basal speed, k _i is the thermal conductivity of ice, Θ_b is the basal temperature gradient, L _i is the latent heat of fusion and ρ _i is the density of ice (Reference PatersonPaterson, 1994). We calculate Θ_b directly from the just described temperature model and use a spatially homogeneous value of G = 70 mW m⁻² in all our estimates, which corresponds to the geothermal heat flux determined from a borehole at Siple Dome (Reference EngelhardtEngelhardt, 2004). Earlier results estimated 80 mW m⁻² near Ridge B/C (Reference Alley and BentleyAlley and Bentley, 1988) and 60 mW m⁻² at the deep borehole drilled near Byrd Station (Reference RoseRose, 1979). The sensitivity of basal melting is such that a change in G by 10 mW m⁻² yields a corresponding change in melt rate of 1 mm a⁻¹, which is small (<10%) relative to the mean melt rates we estimate for this region. Some indirect geophysical estimates of geothermal fluxes in West Antarctica yield significantly higher values of G (median ∼100 W m⁻² with a range of about 50–150 W m⁻²), but their veracity has not been confirmed by direct observational data (Reference Shapiro and RitzwollerShapiro and Ritzwoller, 2004).

The basal shear heating term, τ _b U _b, encompasses all heating at the bed due to the flow of the overlying ice, and is assumed here to be valid for both sliding across an ‘immobile’ plastic-till slip surface (i.e. failure at the yield stress within a very thin layer near the ice/bed interface) and distributed shear deformation within a viscous layer of finite thickness. We obtain the basal shear stress from the inversions described below. While the data shown in Figure 1 provide the surface velocity, the melt rate calculations require basal velocity, U _b. At the ice-stream nodes, there is little internal deformation, so we assume U _b ≈ U _s. Where the bed is strong beneath the ice sheet and some fast-flowing regions, internal deformation can make a significant contribution so that U _s is much greater than U _b. To compensate for this, we estimate the speed due to internal deformation for each melted node (Reference PatersonPaterson, 1994) and subtract it from U _s to estimate U _b. Since we use τ _b to compute the deformation velocity, this correction only has a significant effect in the regions where the bed is strong.

3.2. Coupled ice-stream/ice-shelf model

The equations that approximate large-scale flow over a weak-bedded ice stream or a floating ice shelf are (Reference MacAyealMacAyeal, 1989)

where x and y are the Cartesian coordinates defining the horizontal plane, u and v are the x and y components of velocity, ρ _i is the density of ice, g is the acceleration due to gravity, z _s is the surface elevation, H is the ice thickness, τ _b,x and τ _b,y are the components of the basal shear stress, and υ is the effective viscosity given by

The rate factor, B, and exponent, n, are from Glen’s flow law with an enhancment factor E (Reference PatersonPaterson, 1994), which is equal to 1 in all our experiments. We solve these equations using finite-element methods similar to those described by Reference MacAyealMacAyeal (1989).

The model’s governing stress-balance equations above include both basal drag and deviatoric stresses acting in the horizontal plane as a means to balance the driving stress. These equations are most applicable when vertical shear is small and ice flow is dominated by basal sliding or bed deformation. The model still produces reasonable results, however, in cases where there is significant deformation, but sliding still dominates (Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004).

In the case of a floating ice shelf, the basal shear stress terms are zero. For grounded ice, we use a general power sliding law of the form (Reference PatersonPaterson, 1994)

where the coefficent, β_m , is squared to ensure a non-negative value in the inversion process. Lacking detailed knowledge of effective pressure at the bed, we have subsumed its effect into the coefficient, β_m , which we determine via a model inversion. If m = 1, these equations correspond to a linear-viscous deforming-bed model (Reference MacAyealMac-Ayeal, 1992), which has been used in previous model studies of Pine Island Glacier (Reference Schmeltz, Rignot, Dupont and MacAyealSchmeltz and others, 2002; Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004). Depending on basal conditions, values of m = 2 or 3 are commonly used to model glacier sliding over a hard bed (Reference PatersonPaterson, 1994).

Studies of till rheology, primarily involving actual till samples, suggest a non-linear relationship between basal shear stress and velocity (Reference KambKamb, 1991; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000a; Reference Kamb, Alley and BindschadlerKamb, 2001). Based on these results suggesting a plastic rheology, an alternate parameterization for the basal shear stress is given by (Reference Joughin, MacAyeal and TulaczykJoughin and others, 2004a)

where direction is determined by the velocity, and the speed-independent magnitude of the basal shear stress is determined by α ². In this parameterization, we assume the till is always in plastic failure (i.e. failure once the basal shear stress reaches the yield stress), so the magnitude of the shear stress is independent of speed (Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000a).

Using an appropriate specification of the basal shear stress and ice rheology, the velocity of the coupled ice-stream/ice-shelf system can be solved numerically using Equation (4). Either dynamic or kinematic boundary conditions are needed to solve these equations. In this study, we use kinematic boundary conditions to specify the velocity (e.g. Fig. 1) on regions of the boundary not immediately adjacent to the ocean. For the forward experiments described below, at the floating ice front, we employ a dynamic boundary condition, which is given by (Reference MacAyeal, Rommelaere, Huybrechts, Hulbe, Determann and RitzMacAyeal and others, 1996)

where n_x and n_y are the components of the outward-pointing unit vector normal (in the horizontal plane) to the ice front. While it was simpler to employ the kinematic condition at the floating ice front for the inversion, the dynamic boundary condition is needed in the forward experiment to allow the speed at the front to vary in response to the simulated ungrounding. The kinematic conditions for the rest of the domain boundary are applied in regions well away from the ice stream (e.g. near divides) where the velocity should be unaffected by the ungrounding.

3.3. Inversion for basal shear stress and flow-law parameter

We invert the forward ice-stream model given by Equation (4) to find the basal shear stress that minimizes the misfit between the model-derived and observed velocities. This inversion relies on ‘adjoint-trajectory methods’ (Reference MacAyealMacAyeal, 1992, Reference MacAyeal1993), which are related to optimal control theory. With this method, an ‘adjoint-trajectory’ version of the model is created to develop an automatic and objective means of evaluating a ‘cost function’ (i.e. least-squares misfit between the modelled and observed velocity) as it varies with changes in the undetermined basal shear stress. As with many models, the ice-stream equations are ‘self-adjoint’, so little modification of the forward model code is needed to solve the inverse problem. Similar methods have been used successfully to study the basal shear stress field of several other West Antarctic ice streams (Reference MacAyeal, Bindschadler and ScambosMacAyeal and others, 1995; Reference Joughin, Fahnestock, MacAyeal, Bamber and GogineniJoughin and others, 2001, Reference Joughin, MacAyeal and Tulaczyk2004a) and the ‘Northeast Greenland Ice Stream’ (Reference Joughin, Fahnestock, MacAyeal, Bamber and GogineniJoughin and others, 2001).

We can apply the inversion procedure to either a linear-viscous model given by Equation (6) with m = 1, or a plastic-bed model given by Equation (7). Because force balance is achieved regardless of the assumed bed model, the inversion results alone cannot determine the most appropriate bed model.

The Glen’s flow-law coefficient, B, in Equation (5) can also be inverted for using similar methods (Reference Larour, Rignot, Joughin and AubryLarour and others, 2005). An inversion for both parameters, however, is under-constrained, leading to non-unique results. As a result, for grounded ice we invert only for the basal shear stress, and we determine B using the results from our temperature model and the temperature-dependent relation for B given by Reference PatersonPaterson (1994). This means that the inversion may introduce errors in the basal shear stress distribution that compensate for errors in our assumed values for B, thickness and surface slope in order to achieve the least model–data mismatch. On the floating ice the basal shear stress is fixed at zero, allowing us to invert solely for B (in a strict sense E ^−1/n B) using kinematic boundary conditions on all lateral boundaries. It is important to note that while we are solving for the rheological parameter B, the value determined may include non-rheological effects such as strong rifting and crevassing that contributes to shear-margin weakening (Reference Vieli, Payne, Du and ShepherdVieli and others, 2006).

4.1. Temperature

Using the present-day velocity and thickness distributions and allowing only the temperature to evolve, we ran the temperature model given by Equation (1) for 30 ka to achieve a steady state. In this process, we neglected the influence of temporal changes in surface temperature, accumulation rate (vertical velocity), ice-sheet geometry and ice velocity over this period. The relatively rapid convergence likely can be attributed to our initialization with the vertical advection solution (Reference Joughin, Tulaczyk and EngelhardtJoughin and others, 2003a). Because the latter half of this 30 ka interval falls within the Holocene, it is subject to relatively stable surface temperatures, similar to the present day. Figure 3 shows the final basal temperature gradient from the model, and the temperature gradient for the vertical-advection-only model used to initialize the solution.

Fig. 3.

Modelled basal temperature gradient (a) neglecting horizontal advection and (b) including horizontal advection. Glacier extents are indicated by the contours of glacier speed: 50 m a⁻¹ contour (purple), 100 m a⁻¹ contours up to 900 m a⁻¹ (thin black), and 1000 m a⁻¹ contours (thick black).

Our lack of detailed surface temperature history somewhat limits the accuracy of the results. A more limiting factor may be that the current accumulation rate in this region is poorly known, with errors potentially as large as 100% (Reference Van den Broeke, van de Berg and van MeijgaardVan den Broeke and others, 2006). Finally, the flow speed of Pine Island is actively evolving, and speeds on other glaciers may also have changed over time. Despite similar limitations for the Siple Coast ice streams, our model produced relatively close agreement with measured borehole temperature profiles from that region (Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others, 2004b). No such temperature profiles, however, are available for the Pine Island/Thwaites Glacier region.

Figure 3 illustrates the general pattern of the basal temperature gradient, which determines how rapidly heat is conducted away from the bed. As expected, even with only vertical advection, the temperature gradient is shallower in the deep subglacial troughs. Comparison of Figure 3a and Figure 3b illustrates that horizontal advection has a strong effect on the basal temperature gradient. In areas of convergent flow, where the ice stretches vertically and thickens along flow as it enters the deep main trunks of the glaciers, the magnitudes of the temperature gradients are much smaller than farther upstream. In contrast, as ice is stretched thinner where the troughs shallow near the coast, the temperature gradient steepens substantially relative to the vertical-advection-only results.

4.2. Melt rates

Figure 4 shows our modelled estimates of basal melt, which vary over a wide range. Over much of the area near the grounding lines of Pine Island and Thwaites Glaciers, where the color table saturates (>0.1 m a⁻¹), the melt rate reaches or exceeds 0.4 m a⁻¹. In areas shaded gray, the ice is either frozen to the bed (slower-moving regions outside the main trunk), or the melt rate is negative, indicating basal freeze-on at rates of up to about 5 mm a⁻¹ (faster-moving areas). There are two such areas on Pine Island Glacier’s main trunk, where the weak bed provides insufficient frictional heating to generate melt, allowing freeze-on to occur. In contrast, there is melting beneath the entire fast-moving area of Thwaites Glacier.

Fig. 4.

Basal melt rates for Pine Island and Thwaites Glaciers. Glacier extents are indicated by the contours of glacier speed: 50 m a⁻¹ contour (purple), 100 m a⁻¹ contours up to 900 m a⁻¹ (thin black), and 1000 m a⁻¹ contours (thick black). Although the modelled values may be much higher, the color table saturates at 100 mm a⁻¹ to allow variation in regions of low melt to be distinguished.

The modelled drainage basins of Pine Island and Thwaites Glaciers are nearly equal in size, with areas of 184 000 and 189000 km², respectively. The annual basal melt volume for the Pine Island Glacier catchment is 1.7 km³ a⁻¹, which produces a basin-wide average rate of 9.1 mm a⁻¹. The areal extent of the high-melt region is larger for Thwaites Glacier, yielding nearly twice as much melt at 3.5 km³ a⁻¹ for a basin-wide average melt rate of 18.7 mm a⁻¹.

The majority of the melt (94%) for both catchments is from areas with fast flow (>50 m a⁻¹). Basal shear heating dominates in these regions, so the estimates are relatively insensitive to the geothermal heat flux or the basal temperature gradient. In slow-moving areas, our results show there are large regions with a thawed bed where melt occurs. Figure 3 indicates the melt in these regions is sensitive to the effect of horizontal advection in our estimates of temperature; there would be substantially less interior melting if only vertical advection were considered. The melt rates over much of this area are in the range 1–5 mm a⁻¹, which is comparable to the uncertainty caused by errors in the geothermal heat flux. One such area is at the upper edge of the Thwaites catchment, where we estimate basal melt rates of 2 mm a⁻¹, 30 km from the West Antarctic Ice Sheet (WAIS) Divide drill site (see red star in Fig. 1 for location). Earlier estimates from the other side of the divide also indicate basal melt (Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others, 2004b). While the uncertainty in the geothermal heat flux is large enough to allow for a frozen bed at this location, there is an onset of relatively faster flow in a tributary of Thwaites Glacier about 100 km from the drill site, which suggests a melted bed.

4.3. Basal shear stress

We inverted the ice-stream model for the parameters that determine basal shear stress beneath the fast-moving areas of both glaciers using the methods described above. For the finite-element mesh, we varied the node spacing from ∼1.5 km in fast-moving areas to ∼5 km in slower-moving areas. In the following discussion we refer to the bed as being ‘strong’ where the driving stress is large (>40 kPa) and roughly equal to the basal shear stress. It is important to note that the inversion is not strictly applicable in strong-bedded regions, so the inversion results serve more as a qualitative indicator of bed character in these areas rather than providing quantitatively accurate values (Reference MacAyealMacAyeal, 1993; Reference MacAyeal, Bindschadler and ScambosMacAyeal and others, 1995). We use the term ‘weak’ to describe the regions where the basal shear stress is small (<40 kPa) and likely does not support the full driving stress.

Figure 5 shows basal shear stress for Pine Island Glacier obtained using a viscous-bed model (m = 1). Similar results (not shown) are obtained with a plastic-bed model since force-balance is achieved with either model. Also shown is an image of the same area from the moderate-resolution imaging spectroradiometer (MODIS) Mosaic of Antarctica (MOA), with the contrast stretched to reveal subtle features in the surface topography (Reference Scambos, Haran, Fahnestock, Painter and BohlanderScambos and others, 2007).

Fig. 5.

(a) Inversion of basal shear stress calculated using viscous-bed model for Pine Island Glacier and (b) MOA image. Glacier extents are indicated by the contours of glacier speed: 50 m a⁻¹ contour (purple), 100 m a⁻¹ contours up to 900 m a⁻¹ (thin black), and 1000 m a⁻¹ contours (thick black). The position of the arcuate crevasses discussed in the text is marked ‘AC’.

The inversion reveals a weak area just above the Pine Island grounding line, in the relatively flat area, near flotation, that was earlier identified as an ‘ice plain’ (Reference Corr, Doake, Jenkins and VaughanCorr and others, 2001). The surface slope steepens sharply inland of this region, producing larger driving stress (100–200 kPa), which is resisted locally by the strong bed in this region. Above the steep area, the surface on the main trunk of the glacier and its tributaries flattens out, yielding driving stresses of 10–50 kPa. Over large areas on this trunk, which extends >150 km inland of the strong-bedded region, the bed is weak, offering virtually no resistance (<2 kPa). Instead, the driving stress in this region must be balanced by longitudinal and lateral stress gradients that transfer the stress to isolated ‘sticky spots’, to areas farther upstream and to the bed near the ice-stream margins. This result is consistent with earlier inferences of a weak bed in this region based on geometry (Reference Joughin, Rignot, Rosanova, Lucchitta and BohlanderJoughin and others, 2003b), and analytical (Reference Thomas, Rignot, Kanagaratnam, Krabill and CasassaThomas and others, 2004b) and numerical models (Reference Vieli and PayneVieli and Payne, 2003; Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004). A large band of arcuate crevasses is visible in synthetic aperture radar (SAR) imagery (not shown) (Reference Rabus and LangRabus and Lang, 2003) at the location marked AC in Figure 5, which suggests strong stretching as longitudinal stress gradients compensate for the weak bed. The regions where the shear stress estimates are low correspond to areas in the MOA image with the appearance of more gently undulating surface topography and where flow stripes are visible. This is consistent with the presence of large weak-bedded areas (Reference Gudmundsson, Raymond and BindschadlerGudmundsson and others, 1998).

Figure 6 shows the estimated basal shear stress for Thwaites Glacier. The results indicate a band of strong bed that extends inland for about 80 km from the grounding line. Above this region, there are isolated pockets of weak bed that generally lie between bands of stronger bed. These weak areas likely lie in depressions in the bed topography and may represent areas where till has accumulated as opposed to local highs where it may have been preferentially eroded. These weak areas also agree well with smooth regions visible in the MOA imagery, which suggest a smoother, flatter, ice surface consistent with a weak underlying bed. Relative to Pine Island Glacier, the weak areas seem aligned more across flow than along flow, and there appears to be much more surrounding strong-bedded area, indicating a bed that is substantially stronger on average. This stronger bed likely allows the much steeper surface and higher elevations for the active parts of Thwaites Glacier in comparison with the more drawn-down profile of the upstream regions of Pine Island Glacier.

Fig. 6.

(a) Inversion of basal shear stress calculated using viscous bed model for Thwaites Glacier and (b) MOA image. Glacier extents are indicated by the contours of glacier speed: 50 m a⁻¹ contour (purple), 100 m a⁻¹ contours up to 900 m a⁻¹ (thin black), and 1000 m a⁻¹ contours (thick black).

4.4. Forward model experiments

As noted earlier, any plausible bed (sliding) model has sufficient degrees of freedom so that it can be inverted to produce a basal shear-stress estimate that achieves a balance of forces consistent with the observed velocity field. To determine whether a power sliding law, plastic rheology or alternative bed model is the most appropriate, observations of how the glacier responds to perturbations in this force balance are required. Here we use the forward ice-stream/ice-shelf model (Equations (4) and (5)) and perturb the modelled grounding line to evaluate the dynamic (diagnostic model) response for different bed models. The results of these experiments are then compared with observed changes. Earlier similar model experiments showed that small changes at the grounding line can produce substantial speed-up far inland of the grounding line, but these studies employed only a linear-viscous bed model (Reference Schmeltz, Rignot, Dupont and MacAyealSchmeltz and others, 2002; Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004; Dupont and Alley, 2005).

In our experiments, we use the bed models given by Equations (6) and (7). With m = 1 in Equation (6), we obtain a linear-viscous model. Other studies often employ sliding laws with values of m = 2 or 3 (Reference PatersonPaterson, 1994). For simplicity, we used a single value of m = 3.

The third model we consider is the plastic-bed model. In areas where the bed consists of weak till, a plastic-till rheology is often considered appropriate (Reference KambKamb, 1991; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000a). Recent theoretical results also suggest that hard-bedded sliding can exhibit similar behaviour, even leading to diminishing resistance with increasing speed (e.g. velocity weakening) in areas with low effective pressure and high sliding rates (Reference SchoofSchoof, 2005). If this conjecture is correct, then of the models we consider, the plastic model should most closely approximate rapid sliding over a hard bed, since it does not increase without bound. We applied the plastic-bed model only on the fast-moving parts of the ice stream. For the slower-moving areas in the plastic-bed experiments, we used Equation (6) with m = 3, which allowed basal resistance to increase in areas outside the glacier to maintain a force balance.

To determine the parameters for the forward model, we simultaneously inverted for the basal shear stress beneath grounded ice and for the flow-law parameter B on floating ice. Using these values, the root-mean-square difference between the InSAR velocity (1996) and our forward reference model (Fig. 7) was 25 m a⁻¹ and 21 m a⁻¹ for the linear-viscous and plastic inversions, respectively. Since we used kinematic boundary conditions for the inversion, differences increased slightly when we used a dynamic boundary condition (Equation (8)) along the front of the ice shelf in the forward model. Since we did not implement an explicit inversion for m = 3, we used the m = 1 shear stress solution to solve Equation (6) for β ₃. As a result, the reference cases for the three bed models shown in Figure 8 differ slightly in speed, but all of them compare well with the InSAR speeds (Reference Joughin, Rignot, Rosanova, Lucchitta and BohlanderJoughin and others, 2003b).

Fig. 7.

Reference model speed (color) for the linear-viscous model. Speed is also shown with 50 m a⁻¹ white contours for values up to 500 m a⁻¹. Inset shows the grounding line at the reference position (white) and after mean retreats of 2.9 and 5.2 km.

Fig. 8.

(a) Modelled speed using different bed models (color) for the reference models (GL0) and grounding line retreats of 2.9 (GL1) and 5.2 (GL2) km. Also shown are InSAR-determined speeds from 1992 and 2000. (b) Speed-up expressed as percentage of the reference model speed (GL0) for grounding line retreat of 2.9 (GL1) and 5.2 (GL2) km.

In a manner similar to earlier experiments (Reference Schmeltz, Rignot, Dupont and MacAyealSchmeltz and others, 2002; Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004), we reduced the basal resistance in areas near the grounding line to zero in order to approximate a grounding-line retreat (Fig. 7 inset). Since it does not actually unground the ice, this produces an area that is steeper and thicker than it would be if it had gradually evolved toward flotation. For Pine Island Glacier, this likely is a reasonable approximation of grounding-line retreat since the ice in this area is near flotation and relatively flat (Reference Corr, Doake, Jenkins and VaughanCorr and others, 2001). We performed experiments with mean grounding line retreats of 2.9 and 5.2 km as shown in Figure 8a. Since the retreat is not uniform across the width of the glacier, we computed the retreat distance as the ungrounded area divided by the glacier width.

Figure 8b shows the speed-ups caused by the specified grounding line retreat expressed as a percentage of the reference-model speed. For the 2.9 km retreat experiment, basal resistance was zeroed over an area of 111 km² from a prior mean basal shear stress of 42 kPa. This produced maximum speed-ups near the grounding line from about 4.5% for the linear-viscous case to just over 16% for the plastic bed. These speed-ups are comparable in magnitude to the actual increase in speed from 1992 to 2000 (Reference Rignot and JacobsRignot and others, 2002). For the 5.2 km modelled retreat, the ungrounded area (203 km²) was just under twice the size with a prior mean basal shear stress of 50 kPa. This produced maximum speed-ups of 11% and 37% for the linear-viscous and plastic-bed models, respectively. In both cases, the m = 3 sliding law produced a response about 50% larger than the linear-viscous case.

In addition to the increased sensitivity near the grounding line, both the m = 3 and plastic-bed models produced speed-up substantially farther inland than the linear-viscous model (Reference Schmeltz, Rignot, Dupont and MacAyealSchmeltz and others, 2002; Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004). For the linear-viscous bed, force balance was restored largely by increased basal resistance within about 50 km of the grounding line. For m = 3, basal resistance cannot increase within the same distance of the grounding line, so the speed-up extends farther inland. For the plastic bed, force balance can only be restored by lateral and longitudinal stress gradients that distribute the resistance lost near the grounding line to the non-plastic bed areas outside the fast-moving regions. As a consequence, the speed-up extends much farther inland, up to distances of >200 km from the grounding line.

5.1. Melt rates

The basal steady-state temperature gradient estimates show that the inclusion of horizontal advection in the model has a substantial effect. As in the case of the Ross ice streams (Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others, 2004b), the basal temperature gradient becomes more negative in the fast-moving areas where ice thins toward the grounding line. There are several such areas beneath the fast-flowing sections of the Ross ice streams where the weak basal till produces too little frictional heating to avoid basal freezing (Reference Hulbe and FahnestockHulbe and Fahnestock, 2004; Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others, 2004b). For these areas, continued basal lubrication may rely on transport of water generated farther upstream in tributaries, where higher basal shear stresses and shallower basal temperature gradients favour greater melt (e.g. Reference Parizek, Alley and HulbeParizek and others, 2003). In contrast, for Pine Island and Thwaites Glaciers, although the temperature gradient is not conducive to melt, the rapid speeds and high basal shear stress produce strong melt (>0.1 m a⁻¹) near the grounding line. As a result, these high melt rates yield basin-wide average melt rates for Pine Island and Thwaites Glaciers that are roughly 4.5 times larger than for the Ross ice streams.

Limited melt beneath weak-bedded ice streams can provide a thermal feedback whereby an ice stream thins and promotes basal freezing, eventually leading to ice-stream shutdown (Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000b; Reference Bougamont, Tulaczyk and JoughinBougamont and others, 2003a,Reference Bougamont, Tulaczyk and Joughinb). This, along with subsequent reactivation of these ice streams, has led to significant flow variability in the Ross Sea section over the last millennium (Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others, 2000). With the strong melt generated beneath both Thwaites and Pine Island Glaciers, it is unlikely that current or past variability can be attributed to such a feedback (Reference RaymondRaymond, 2000). Furthermore, because melt due to basal shear heating is so large, it is unlikely that temporal fluctuations in geothermal heat flux due to nearby volcanic activity could explain recent changes.

Channels in crystalline bedrock are found on the sea floor in formerly glaciated regions within a few tens of kilometers of the present Pine Island Ice Shelf, which are believed to have been formed by abundant subglacial meltwater (Reference Lowe and AndersonLowe and Anderson, 2003). They estimate relatively modest melt rates of 0.24–0.5 km³ a⁻¹ and conclude additional melt must originate farther upstream. Our basin-wide rate of 1.7 km³ a⁻¹ is considerably higher and almost entirely generated by strong melt (∼0.1 m a⁻¹) over fast-moving regions with a strong bed. These high-melt regions may correspond to regions of crystalline bedrock, similar to that where the channelized topography has been observed. If similar in magnitude to the rates estimated here, local melt in the channelized regions may have been sufficient to produce the observed features.

5.2. Basal Shear Stress

The basal shear stress beneath Pine Island and Thwaites Glaciers is characterized by areas of strong resistance near the grounding line, with patches of weak bed at locations farther inland. These weak areas are more limited for Thwaites Glacier, which may explain its steeper profile. In contrast, the extensive weak area along much of Pine Island Glacier’s main trunk may account for its shallow slope and likely contributes to the relatively low elevation (<1300 m) of the divide between it and Rutford Ice Stream, which also has extensive weak areas (Reference Vaughan, Smith, Nath and Le MeurVaughan and others, 2003; Reference Joughin, Bamber, Scambos, Tulaczyk, Fahnestock and MacAyealJoughin and others, 2006).

Observations from the sea floor in front of Pine Island Glacier reveal areas of exposed crystalline bedrock (Reference Lowe and AndersonLowe and Anderson, 2003), likely similar in character to regions where we calculate high basal shear stress. In other areas near the front of the glacier, piston cores and other geophysical observations indicate areas covered by weak sediment deposits with shear strengths ranging from 1.0 to 6.9 kPa, which is consistent with the magnitudes of shear stresses for the weak areas we find farther inland. Thus, unlike the Ross ice streams, which tend to flow over large expanses of uniformly weak bed (Reference Joughin, MacAyeal and TulaczykJoughin and others, 2004a), our results and ship-borne observations reveal ‘mixed’ bed conditions, alternating between regions of low drag, almost certainly deforming sediments, and those providing greater basal resistance, possibly non-deforming sedimentary rock or even crystalline bedrock of the type identified by Reference Lowe and AndersonLowe and Anderson (2003) on the nearby continental shelf. Over extensive areas near their respective grounding lines, both Pine Island and Thwaites Glaciers are presently grounded on a strong bed, suggesting they have retreated back over weak sediment-laden regions to positions of relative stability.

Although there are weak areas beneath Thwaites Glacier, these tend to be somewhat isolated and surrounded by nearby areas of relatively strong bed, which likely compensate for the weak regions and may promote some degree of stability. In contrast, there is a >140 km long weak region extending up much of the length of Pine Island Glacier’s main trunk. This region is ‘buttressed’ by the ∼40 km wide band of high basal shear stress located just above the grounding line, much of which is only a few tens to hundreds of meters above flotation. The results in Figure 8 along with those of earlier studies (Reference Schmeltz, Rignot, Dupont and MacAyealSchmeltz and others, 2002; Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004; Reference Thomas, Rignot, Kanagaratnam, Krabill and CasassaThomas and others, 2004b) suggest that the glacier is extremely sensitive to even small changes in loss of this basal resistance. Such results are consistent with observations that this glacier has accelerated by >40% between 1974 and 2006 (Reference Joughin, Rignot, Rosanova, Lucchitta and BohlanderJoughin and others, 2003b; Reference RignotRignot and others, 2008) as areas near the grounding line have thinned by >1 m a⁻¹ (Reference Shepherd, Wingham, Mansley and CorrShepherd and others, 2001), and the grounding line has retreated by several kilometers (Reference RignotRignot, 1998). This retreat seems to be largely confined to the region near the grounding line where we infer a relatively weak bed. Thus, the basal shear-stress distribution shown in Figure 5 suggests the potential for a large instability if the ice continues to thin and retreat over this region, ungrounding ice from the strong-bedded region that appears to be holding back ice from the weak-bedded interior.

5.3. Bed model

Our experiments (Fig. 8) reveal that in addition to the magnitude of the basal shear stress, the form of the bed model plays a large role in the modelled sensitivity of the glacier to changes in basal resistance near the grounding line, affecting both the magnitude and spatial scale of the response. While early studies suggested that a linear-viscous bed was appropriate for flow over weak till (Reference Alley, Blankenship, Bentley and RooneyAlley and others, 1986; Reference Blankenship, Bentley, Rooney and AlleyBlankenship and others, 1986), more recent studies find that a perfectly plastic rheology is a better model (Reference KambKamb, 1991; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000a). Furthermore, a linear-viscous bed model is not appropriate for sliding over bedrock (Reference PatersonPaterson, 1994). Thus, it is likely that the linear-viscous model, which yielded the least sensitivity in our experiments, is an inappropriate model and earlier studies based on this bed model likely underestimate the glacier’s sensitivity to changes in basal resistance near the grounding line (Reference Schmeltz, Rignot, Dupont and MacAyealSchmeltz and others, 2002; Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004).

Although a plastic-bed model may be appropriate for weak till, flow over a hard bedrock typically is represented by sliding laws of the form of Equation (6) with m = 2 or m = 3 (Reference PatersonPaterson, 1994). Our results indicate that a sliding law with m = 3 produces a substantially larger (∼50%) speed-up that extends significantly farther inland than the response for a linear-viscous model to the simulated ungrounding. In addition to this nearly instantaneous response, prognostic ice-stream model studies suggest that the ensuing time-dependent response will propagate inland more rapidly as the degree of non-linearity in the bed model is increased (Reference Price, Conway, Waddington and BindschadlerPrice and others, 2008).

Theoretical work indicates that for a bed with bounded slopes, basal shear stress cannot increase without bound as it does in a power sliding law (Reference IkenIken, 1981; Reference SchoofSchoof, 2005). Instead, these theoretical studies suggest that for sliding at high speed and low effective pressure, resistance initially increases with speed, similar to that for power-law sliding, but reaches a peak after which it decreases with faster sliding. None of the models we used is equivalent to the sliding model proposed by Reference SchoofSchoof (2005), but the plastic model has the characteristic that shear stress does not increase without bound as sliding speed increases. Because shear stress is fixed and does not decline, the plastic model is likely to be less sensitive to geometric changes that alter the force balance than the type of sliding law proposed by Reference SchoofSchoof (2005). Hence, even with the plastic-bed model, our results may underestimate the glaciers’ sensitivity to grounding-line retreat.

The extent of the simulated ungrounding in our model experiments is consistent with the magnitude of observed changes in grounding line position (Reference RignotRignot, 1998). There is enough observational uncertainty to allow plausible adjustment of the modelled grounding line so that any of the three models can match the magnitude of speed-up in the near-grounding-line region. Only the modelled speed-up for the plastic bed is still appreciable at distances of >200 km from the grounding line, which is consistent with the extent over which there is substantial thinning (Reference Shepherd, Wingham and MansleyShepherd and others, 2002). Thus, with a plastic-bed or similar sliding law (Reference SchoofSchoof, 2005), the direct response to a grounding line shift could produce thinning comparable to the observed magnitudes. This model also provides the best apparent agreement with the extent of the observed inland speed-up (Fig. 8b).

When considering the agreement between observed and modelled speed-up, it is important to consider the observations showing speed-up were acquired 8 years apart. Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others (2004) demonstrated how an initial speed-up can steepen surface slopes and diffuse rapidly inland over the course of years or decades to produce inland thinning and speed-up over a much larger area than affected by the initial speed-up. This type of effect undoubtedly has some influence on the InSAR speeds shown in Figure 8. Thus, while the agreement between model and data suggests something approaching a plastic-like bed response, the results are far from definitive. Nevertheless, the agreement is close enough to clearly demonstrate at least the plausibility of this type of behaviour.