Introduction
It has long been recognized that rapid surge-like, partial or total gravitational, collapse in all or parts of the Antarctic ice sheet, with a concomitant rapid rise in global sea level, depends on the bed beneath the surging portions passing from a frozen to a thawed condition (e.g.Reference WilsonWilson, 1964; Reference WeertmanWeertman, 1966, Reference Weertman1969; Reference McInnes, Radok, Budd and SmithMcInnes and others, 1985, Reference McInnes, Budd, Smith and Radok1986; Reference Jenssen, Budd, Smith and RadokJenssen and others, 1985; Reference Radok, Brown, Jenssen, Smith and BuddRadok and others, 1986; Reference MacAyealMacAyeal, 1992; Reference Verbitsky and SaltzmanVerbitsky and Saltzman, 1997). Thermomechanical ice-sheet models, beginning with Reference Budd, Jenssen and RadokBudd and others (1970, Reference Budd, Jenssen and Radok1971) and continuing in the past decade (e.g. Reference HuybrechtsHuybrechts, 1990; Reference Hansen and GreveHansen and Greve, 1996; Reference FastookFastook, 1997; Reference Verbitsky and SaltzmanVerbitsky and Saltzman, 1997), aspire to map the thawed portions of the bed. All of these modeling studies show that the thawed fraction of the bed is highly sensitive to the basal geothermal heat flux and to past and present surface temperatures and accumulation rates.
It is unlikely that any of these model inputs will ever be known with the accuracy needed to map the thawed fraction of the bed using even the most sophisticated time-dependent three-dimensional thermomechanical models. Ubiquitous drilling into bedrock under the ice sheet could produce an adequate map of the geothermal flux, and ice cores recovered from the drill holes would yield data on past surface temperatures and accumulation rates, but that would be prohibitively expensive, and ice-core data could be compatible with a variety of surface conditions. Even present-day surface temperatures and accumulation rates are not known with the accuracy or uniform distribution required for these sophisticated models.
Apart from ice-core studies, Antarctic surface temperatures and accumulation rates over time could also be provided as output from climate models. For the output to be reliable, however, the climate models would need to accurately couple the dynamics of the atmosphere, oceans, biosphere, ice sheets and sea ice during the last glaciation cycle. Thus, attaining one modeling objective, mapping the basal thawed fraction under the Antarctic ice sheet, relies on even more problematic modeling enterprises. At least time-dependent three-dimensional ice-sheet models exist. No such climate models exist. The best atmospheric circulation models provide only steady-state simulations.
These difficulties might be sidestepped by an ice-sheet model that computes the thawed fraction of the bed under the Antarctic ice sheet from model input that is independent of the basal geothermal heat flux and is only weakly dependent on surface temperatures and accumulation rates. Such a model, and a map of the thawed fraction computed from the model, is presented here. The model uses ice thickness and ice-surface slope variations along flowlines of the Ant arctic ice sheet to determine ice–bed coupling for sheet flow, as measured by the thawed fraction of the bed. The thawed fraction is only weakly dependent on variations in surface temperatures and accumulation rates along the flowlines. The basal geothermal heat flux contributes to thawing the bed, but need not be known in order to map the thawed fraction using this model.
The Model
The model is simple but not simplistic. It is based on the observation by Reference NyeNye (1952) that ice-surface elevations along an ice-sheet flowline can be calculated if the bed topography and the basal shear stress, equated with the gravitational driving stress, are known. The driving stress τ 0 for sheet flow was obtained from the force balance by BGS (1949):
where ρ I is ice density, gz is gravity acceleration in the vertical z direction, h I is ice thickness and α is ice-surface slope. This is the shallow-ice approximation for small bed slopes (Reference HutterHutter, 1983). Let ice-surface elevation h increase by increment Δh in increment Δx measured along horizontal distance x from the ice-sheet margin, following a flowline, and let h R be the bedrock elevation above or depth below sea level at distance x along the flowline. Then h I = h − h R and α = Δh/Δx, and a numerical solution of Equation (1) gives ice elevation h at Δx step i + 1 if h, and τ 0 are known at step i:
where Δh = hi +1 − hi and i = x/Δx is an integer that denotes the number of equal Δx steps at distance x along a flowline of length L from the ice margin to the ice divide.
Equation (2) is also obtained from a general equation for the surface slope Δh/Δx along a flow band which has variable width, variable surface accumulation or ablation rates and variable bed topography, and along which transitions from sheet flow to stream flow to shelf flow can occur as ice moves from the ice divide to the melting or calving margin of an ice sheet. The general equation is equation (3.25) in Reference HughesHughes (1998), where it is derived on p. 51–55. Ice–bed coupling in the general equation is controlled by two variables, the thawed fraction f of the bed for sheet flow and the ratio P
w/P
I of basal water pressure P
w to ice-overburden pressure P
I for stream flow. Both variables appear at each Δx step along the flow band. For sheet flow, the gravitational driving stress is resisted primarily by basal shear stress τ
0, which decreases when a frozen bed becomes thawed, with f recording the decrease. This assumes that basal water exists primarily as a film, as in the Reference WeertmanWeertman (1957) theory of basal sliding, so that P
w ≪ P
I. Hence, P
w/P
I = 0 for a frozen bed and P
w/P
I ≪ 1 for a thawed bed during sheet flow, so that f determines the surface slope. Equation (3.25) in Reference HughesHughes (1998) reduces to Equation (2) in this case. For stream flow, not considered here, f = 1 and P
w/P
I increases along the ice stream to P
w/P
I = 1 at the grounding line of a freely floating ice shelf. The gravitational driving stress in ice streams is resisted primarily by the longitudinal deviator stress
and the side shear stress τ
s.
Equation (2) requires knowing h at an initial step i = 0, and knowing h R and τ 0 at all steps. A value of h = h R for i = 0 at the ice-sheet margin is not allowed. Grounded margins of the Antarctic ice sheet are typically at calving ice walls or at the grounding lines of floating ice shelves, where h > h R at i = 0 in both cases. Equation (2) also requires knowing τ 0 and h R at every step i along ice flowlines. Values of both h and h R along flowlines can be determined from maps of surface and bed topography contoured at 100 and 250 m intervals, respectively, in the glaciological and geophysical folio of Antarctica (Reference DrewryDrewry, 1983). Equation (2) can therefore be solved numerically to obtain the variation of τ 0 along flowlines.
As illustrated in Figure 1 for a rolling subglacial landscape, transition zones between thawed (wet) and frozen (dry) zones begin as isolated frozen patches that become larger and more numerous until they coalesce, and end as isolated thawed patches that become smaller and less numerous until they vanish. This is a freezing zone when ice moves from a thawed zone to a frozen zone, and a melting zone when ice moves from a frozen zone to a thawed zone. In these transition zones, frozen and thawed regions would tend to correspond to the respective hills and dales in bed topography, because glacial erosion occurs where the bed is thawed. Topographic variations allowing frozen and thawed patches seen in the glacial geology produced by former ice sheets exist on scales ranging from several meters to several kilometers (Reference Kleman, Hättestrand and ClarhällKleman and others, 1999), but they are not resolved on the map of Antarctic subglacial topography (Reference DrewryDrewry, 1983). Variations of thawed fraction f across these transition zones can be calculated nonetheless, where f is defined as the fraction of the bed that is thawed in an area of incremental length Δx = 50 km and mean width
in a flow band.
Fig. 1. Basal thermal zones beneath an ice sheet. Frozen (dry) areas are white, and thawed (wet) areas are black. Ice flow from a frozen zone to a thawed zone crosses a melting zone. Ice flow from a thawed zone to a frozen zone crosses a freezing zone.
Freezing or melting transition zones in this model must not be confused with freezing or melting rates when the bed is thawed in thermomechanical models. Only two values of f exist in these models, f = 0 where the bed is frozen and f = 1 where the bed is thawed. A thawed bed has either a freezing rate or a melting rate, as was calculated in the first thermomechanical flowline model for melting rates (Reference Budd, Jenssen and RadokBudd and others, 1970, Reference Budd, Jenssen and Radok1971) and then for freezing rates (Reference HughesHughes, 1973; Reference SugdenSugden, 1977).
In Equation (2), if variations in bed topography represented by h
R do not fully account for variations in surface elevation represented by h, then variations in τ
0 must account for the remaining h variations when sheet flow predominates. Variations in τ
0 are caused by variations in accumulation rate
, surface temperature T
s and flow-band width w, but the major τ
0 variations are caused by variations in ice–bed coupling as represented by the thawed fraction f of the bed. If f = 0, the bed is frozen and ice motion is determined by the vertically averaged creep velocity u of ice given by:
where τ D = τ 0 for a dry (frozen) bed, n is a visco-plastic parameter for ice and A is the temperature-and-fabric-dependent vertically averaged ice-hardness parameter, as determined from creep experiments (Reference NyeNye, 1953; Reference GlenGlen, 1955). If f = 1, the bed is thawed and ice motion is dominated by the basal velocity u 0 of ice given by:
where τ W = τ 0 for a wet (thawed) bed, m is a visco-plastic parameter that depends on how basal ice slides over bedrock or how sediment or till deforms between basal ice and bedrock, and B depends on bedrock roughness for basal sliding or on the properties of a sediment or till layer. For basal sliding over bedrock, m = (n + l)/2 in the Reference WeertmanWeertman (1957) theory, and for deformation of basal sediment or till, m ≈ 1 for viscous flow (Reference Boulton and HindmarshBoulton and Hindmarsh, 1987), and m ≈ 50 for plastic flow (Reference KambKamb, 1991). Other formulations for u 0 have been proposed for basal sliding (e.g. Reference LliboutryLliboutry, 1968; Reference Budd, Keage and BlundyBudd and others, 1979) and bed deformation (e.g. Reference ClarkeClarke, 1987; Reference RaymondRaymond, 1996; Reference HookeHooke, 1998). Until a general theory that combines these processes has been accepted, Equation (4) remains a purely empirical expression in which the pertinent physics are represented by m and B in an unknown way. Reference HughesHughes (1998, p. 87) multiplies B in Equation (4) by(l − P W/P I)2, where P W is the basal water pressure and P I is the ice-overburden pressure such that P W ≪ P I for sheet flow, and P W increases for stream flow from P W ≪ P I at the heads of ice streams to P W ≈ P I as ice streams become afloat. This gives ice streams their concave surface profile. Our study treats only sheet flow, using the Reference WeertmanWeertman (1957) theory.
The relationship between τ 0 τ D, τ W and f for sheet flow is:
Expressions giving τ
D and τ
D at distance x from the end of a low band of length L for constant accumulation rate
and flow-band width w are (Reference Hughes, Denton and HughesHughes, 1981, p. 243, 244):
When n = 3 and m = (n + l)/2 = 2,
and
. This leads to
for a frozen bed and to
for a thawed bed (Reference Hughes, Denton and HughesHughes, 1981, p. 249). Therefore, ice elevations calculated from Equation (2) are weakly dependent on variations in
along flowlines, compared to the variation of f over the range 0 ≤ f ≤ 1 for sheet flow in Equation (5).
Figure 2 is a plot of Equation (2) for h
R = 0, Δx = 20 km, L = 1000 km,
, n = 3, m = 2, A = 4.0 bar a1/3 for ice at −15°C in the Reference GlenGlen (1955) flow law, and B = 0.02 bar a1/2 m−1/2 for a bed roughness factor of 0.05 in the Reference WeertmanWeertman (1957) sliding law. Figure 2 shows variations of ice-surface elevation along a flowline over a flat horizontal bed that are due solely to variations in ice–bed coupling as ice flows over a succession of thawed (melted), freezing, frozen and melting zones on the bed. A thawed zone is wholly wet, a frozen zone is wholly dry, a freezing zone is a transition from a thawed zone to a frozen zone, and a melting zone is a transition from a frozen zone to a thawed zone, as depicted in Figure 1. As already noted, melting and freezing zones in our model are not to be confused with the melting and freezing rates that thermomechanical ice-sheet models compute when the bed is wholly thawed (wet). As shown in Figure 1, our melting and freezing zones are a mosaic of frozen and thawed patches, based on studies of deglaciated landscapes (e.g. Reference Kleman, Hättestrand and ClarhällKleman and others, 1999). Basal melting and freezing rates are equal in a thawed patch during steady-state equilibrium, with basal melting in the back part of the patch and basal freezing in the front part (see Reference HughesHughes, 1998, fig. 9.16).
Fig. 2. Changes in ice-surface slope due to changes in ice-bed coupling along an ice-sheet flowline for constant flow-band width, surface accumulation rate, bed topography and parameters in flow and sliding “laws” for glacial ice. Top: The flowline profile (solid urve) varies between the profiles for a dry frozen bed (upper dashed curve) and a wet melted bed (lower dashed curve) as the ice flows over a succession of melted, freezing, frozen and melting basal thermal zones. Bottom: The actual basal shear stress τ0 (solid curve) for ice moving over the succession of melted, freezing, frozen and melting basal thermal zones lies between the basal shear stress τD for ice moving over a dry frozen bed (upper dashed curve) and the basal shear stress τW for ice moving over a wet melted bed (lower dashed curve).
When accumulation and ablation rates vary along flow bands of variable width, mass-balance continuity provides the following expressions for τ D and τ W (Reference Fastook, Hansen and TakahashiFastook, 1984), which are obtained by combining the equations for the force balance, the flow and sliding laws and the mass balance, as shown by Reference HughesHughes (1998, p. 146–157):
where variable surface accumulation or ablation rates
along flow bands of variable widths w(x) and lengths L from x = 0 define Γ(x) for the mass balance:
and
Equations (7a) and (7b) reduce to Equations (6a) and (6b) when
and w are constant.
Specifying surface and bed slopes
Calculating τ 0 from Equation (2) requires input values of ice-surface elevation and bed elevation at particular stepped distances Δx along each flowline. Because our data points are from where surface flowlines crossed surface- and bed-topographic contour lines and rarely coincided with these intervals, we had to fit continuous surface and bed profiles to the data and derive new sets of paired surface and bed data from the continuous profiles.
A weighted cubic spline was chosen to represent the surface and the bed between elevation contour lines on the Antarctic folio maps (Reference DrewryDrewry, 1983) used in our study. The cubic spline, a third-degree polynomial program developed by Reference Press, Flannery, Teukolsky and VetterlingPress and others (1986), is smooth and passes through all data points, but it can produce oscillations between data points. We included an exponential weighting function P that suppressed oscillations. A nearly linear connect-the-dot solution occurs when p ≥ 3. We used P = 0.5 to fit bed-topographic contours, and, to prevent dips in the surface profiles that give negative driving stresses, we used P = 1.0 to fit surface-topographic contours. Figure 3 illustrates how different values of P affect the weighted cubic-spline curves connecting known bed elevations along a 500 km segment of flowline B4 in Figure 4. All surface and bed profiles along flowlines were generated with an interpolation interval of 1 km on the spline curves.
Fig. 3. Weighted spline interpolations for a section of glacier bed along a 500 km section of flowline BY4 in Figure 4. Weight factors Pare shown for P = 1.0, 0.5 and 0.25 with an interpolating interval of 1 km. Bed topography is known at the triangles, where subglacial topographic contour lines in folio map 3 of Reference DrewryDrewry (1983) intersect the flowline.
Fig. 4. Flowlines of selected flow bands in East Antarctica and West Antarctica along which the thawed-bed percentage was calculated for comparison with bed topography in Figure 5, using Equations (2), (5), (7a) and (7b). East Antarctic flowlines are for Byrd Glacier (BY2–BY5), Mullock Glacier (MU2 and MU3), David Glacier (DA2–DA4), Oates Coast (OL1–OL5), George V Coast (GV1–GV8) and Wilkes Land (W1–W19). West Antarctic flowlines are for Pine Island Glacier (P1–P7), Thwaites Glacier (T1–T5) and ice streams supplying the Ross Ice Shelf (R1–R13). Domes (D) along ice divides and ice margins are identified. Thin solid lines are ice-drainage divides.
Calculation of basal shear stress using Equation (2) assumes that longitudinal stresses and their gradients are negligible. Reference BuddBudd (1968) maintained that effects of longitudinal stresses are negligible over distances >20 times the ice thickness for sheet flow. Mean ice thickness is about 2.5 km for the East Antarctic ice sheet and about 2 km for the West Antarctic ice sheet. Therefore, a step interval of 50 km was chosen. Calculations are made with respect to the midpoint of each step (Reference Hughes, Denton and HughesHughes, 1981, p. 258–259). Bedrock slopes β along the spline curves are often greater than ice-surface slopes a. To check if bedrock slope β affects our results, we also calculated τ 0 using the formula (Reference PatersonPaterson, 1981; Reference HughesHughes, 1998, p. 144–145):
Our spline-fit program allowed α = Δh/Δx and β = Δh R/Δx to be calculated for Δx steps of 1–50 km. Equation (10) was then solved numerically using both Euler and Runge–Kutta methods to see what maximum Δx gave the same τ 0 for both methods and also provided a negligible contribution from bed slopes, given the reliability imposed by the folio data. Steps of Δx = 50 km satisfied these criteria.
Output generation
At each Δx step along a specified flow-band center line, the bedrock elevation, flow-band width, accumulation rate and thawed-bed percentage are required as input to construct an ice-sheet surface elevational profile by the finite-difference model used to solve Equation (2) when To is given by Equations (5), (7a) and (7b). In our application of Equations (2) and (5), constructed ice-sheet surface elevations are fitted to the actual ice surface as closely as possible by varying the fraction of the bed that is thawed in each Δx step for known accumulation rates and flow-band widths in Equations (7a) and (7b). Initially, the thawed-bed fraction is specified as zero along the length of each flow band.
Flow-band center lines used in this study are the flowlines labeled in Figure 4. Flowline profiles were calculated along the center line of flow bands for which bed topography is mapped well enough to justify our study. Near the ice-sheet margin, flow converges into bedrock troughs occupied by ice streams and diverges across bedrock headlands between ice streams. These flow variations may not always be apparent at the scale of the folio maps. We minimized this problem by beginning flow bands at the 200 m ice-elevation contour. Flow-band widths were obtained from initial estimates of the flow-band shape near the margins, but could be changed in successive reconstructions. Change of flow-band width is appropriate where a good match of the reconstructed surface to the existing surface cannot be obtained by varying the flow parameters, indicating flow that is more strongly convergent or divergent near ice-sheet margins than is indicated on the folio map. As seen on the folio maps of surface and bed topography, radio-echo flight-lines along which this topography was measured are particularly sparse near ice-sheet margins.
Accumulation rates were tabulated from the map of surface mass-balance rates by Reference Giovinetto, Bentley and BullGiovinetto and others (1989). Their map was enlarged to the 1: 6 000 000 scale of flow bands on the folio map. Accumulation rates, in kg m−2 a−1, were converted to m a taking ρ I = 917 kg m−3.
Initial values A = 4.0 bar a1/3 and B = 0.02 bar a1/2 m −1/2 n Equations (7a) and (7b) were adjusted so that driving stresses τ 0 calculated from Equation (2) lie on or between the theoretical τ D and τ W curves, following Figure 2. Then thawed fraction f was calculated from Equation (5) for τ 0 values lying between the τ D and τ W curves. At the same time, the reconstructed ice-sheet surface should match closely the actual ice-sheet surface. Where the two profiles failed to match for reasonable values of A and B, flow-band width may be evaluated and changed, where that is more reasonable, such as near ice-sheet margins. However, large variations in flow-band width cause only minor changes in τ 0, τ D and τ W, so the mismatch with known surface slopes near the ice-sheet margin probably reflects transitions from sheet flow to stream flow. Equations (2) and (5) do not adequately express these transitions. Perhaps P W/P I becomes a significant variable in transitions from sheet flow to stream flow. Even so, these transitions take place in only one or two 50 km Δx steps. When a satisfactory profile was generated, the final dataset was run through the program again, and the screen display and dataset were printed. Final values of A and B for n = 3 and m = 2 are listed in Table 1 for East Antarctic flowlines and in Table 2 for West Antarctic flowlines. These values minimize the population standard deviations between calculated and measured surface elevations along flowlines (DIFF).
Table 1. Best-fit constants in the flow law and the sliding law for East Antarctic flowlines (see Fig 4 for list of abbreviations)
Table 2. Best-fit constants in the flow law and the sliding law for West Antarctic flowlines (see Fig. 4 for list of abbreviations)
Determining ice-hardness and basal-sliding parameters
The flow-band model has been applied to the Antarctic flow bands shown in Figure 4, with f calculated along the central flowlines of these flow bands. Most flow bands shown in Figure 4 dovetail into ice streams. Central flowlines are identified by letters and numbers, as shown. Wide ice-drainage basins, such as the ice-drainage basin for Byrd Glacier, have been divided into several narrower flow bands in order to map f contours more accurately. The f variations are shown in Figures 5–9 for selected flow bands. In these figures, τ 0 was computed from Equation (2) at 50 km steps, shown as triangles, to which a solid τ 0 curve is fitted. The solid τ 0 curve is bracketed by the upper τ D and the lower τ W dashed curves computed from Equations (7a) and (7b), respectively. Variations of f along the central flowline of each flow band are then calculated from Equation (5). In these calculations, n = 3 and m = 2 are used, A and B are obtained by bracketing the observed τ 0 variation between the theoretical τ D and τ W variations, and f is assumed to vary along flowlines much more rapidly than do A and B. In moving from the ice divide to the ice margin along flowlines, a frozen bed exists where τ 0 = τ D, a thawed bed exists where τ 0 = τ W, a melting bed exists where τ 0 moves from τ D toward τ W, and a freezing bed exists where τ 0 moves from τ W toward τ D. A rapid drop of τ 0 near the ice margin occurs when sheet flow becomes stream flow.
Fig. 5. The basal thermal regime along East Antarctic flowlines BT3, MU2 and DA2 toward the Ross Sea embayment. Upperpanel: The actual τ0 variation (solid curve with triangles denoting 50 km steps) lies between the theoretical variations τD for creep above a dry frozen bed (upper dashed curve) and τW for sliding on a wet thawed bed (lower dashed curve), with f being the fraction that τ0 lies between the τD and τW curves at each step in order to satisfy Equation (5). Center panel: Variations in flow-band width (solid curve) and accumulation rate (dashed curve) along the flowline. Lower panel: The spline-fitted bedrock topography (curve through circled dots) that was used with f obtained from the top panel to compute a theoretical ice-surface profile (solid curve) that provides a best fit with the actual spline-fitted ice-surface profile (triangles). The dashed curves are for theoretical isostatic rebound of the surface and the bed. DIFF numbers, in meters, are the population standard deviations between the calculated surface elevations and the measured surface elevations along each flowline. The DIFF number identifies values of flow and sliding parameters A and B, respectively, in Tables 1 and 2.
Fig. 6. The basal thermal regime along East Antarctic flowlines OL2 and OL5 towards Oates Land. See Figure 5 caption for explanation.
Fig. 7. The basal thermal regime along East Antarctic flow-lines GV5 and GV7 toward the George V Coast. See Figure 5 caption for explanation.
Fig. 8. The basal thermal regime along East Antarctic flowlines W7, W12, W16 and W19 in Wilkes Land. See Figure 5 caption for explanation.
Fig. 9. The basal thermal regime along West Antarctic flowlines P4, P5, T3 and T4 toward Pine Island Bay. See Figure 5 caption for explanation.
Basal sliding over bedrock is likely for East Antarctic flowlines because East Antarctica is believed to be a crystalline shield from which sediments have been largely stripped away by glacial erosion over millions of years (Reference Marchant, Swisher, Lux, West and DentonMarchant and others, 1993). In contrast, bedrock in West Antarctica may still be carpeted in many places by marine sediments that are cemented by ice where the bed is frozen but become a deforming substrate where the bed is thawed. Therefore, taking m = 2 in Equation (7b) for ice sliding over bedrock in the Reference WeertmanWeertman (1957) sliding theory may apply generally in East Antarctica, but only when a deforming substrate has been stripped from bedrock in West Antarctica, or when m = 2 describes deformation of the substrate. However, the substrate under Ice Stream B, West Antarctica, deforms at m ≈ 50 (Reference KambKamb, 1991). We assume that a deforming substrate will become concentrated beneath ice streams by converging flow, whereas basal sliding over bedrock is more common where sheet flow prevails, even in West Antarctica.
Flowlines are identified for separate ice-drainage systems in Figure 4. In East Antarctica, BY2–BY5 are for Byrd Glacier, MU2 and MU3 are for Mulock Glacier, DA2–DA4 are for David Glacier, OL1–OL5 are for Oates Land, GV1–GV8 are for the George V Coast and W1–W19 are for Wilkes Land. In West Antarctica, P1–P7 are for Pine Island Glacier, T1–T5 are for Thwaites Glacier and R1–R13 are for ice draining into the Ross Ice Shelf. Letters D locate minor domes, especially West Antarctic ice domes. Most flowlines end in ice streams. The effect of different A and B values on bracketing the observed τ 0 variation along a flowline between the theoretical τ 0 curves obtained from Equation (7a) for creep of ice and from Equation (7b) for basal sliding was determined for many flowlines (see Tables 1 and 2). Although Equations (7a) and (7b) do not permit A and B to vary along flowlines, a sense of the possible variations is obtained by seeing how A and B vary from flowline to flowline in Tables 1 and 2 when DIFF is minimized between observed and calculated surface elevations along flowlines. Since A = 5.0 ± 1.0 bar a1/3 and B = 0.04 ± 0.01 bar a1/2 m−1/2 minimize DIFF for nearly all flowlines, it is tempting to conclude that no drastic changes in A and B occur along flowlines. Drastic changes in A are unlikely because A reflects vertically averaged ice temperatures that vary slowly along flowlines (even though basal ice temperatures may vary irregularly). Drastic changes in B are unlikely because basal sliding for millions of years over bedrock (especially over the Precambrian Shield in East Antarctica) should produce a bed roughness that varies slowly along flowlines.
Reference HindmarshHindmarsh (1990) notes that basal ice temperatures could reasonably change from −30°C to the melting point along an Antarctic ice flowline, which would make A vary by a factor of 100, with A 1/(2n+2) used for comparison with f variations. When A changes by a factor of 100, A 1/8 for n = 3 changes by a factor of two. This is within the 0 ≤ f ≤ 1 range of f variations. However, vertically averaged variations of A along flowlines in our model are minor compared to variations of f because variations of A depend on both ice-temperature and fabric variations that interact to keep A relatively constant. The big variations in both temperature and fabric as a frozen bed thaws, and as a thawed bed freezes, occur near the bed. Simple shear in basal ice is favored when a frozen bed maximizes bed traction. Simple shear causes strain softening that lowers A to offset the higher A produced by cold basal ice. When the bed thaws, bed traction is reduced and the strain-softened ice produced by simple shear is replaced by strain-hardened ice produced by longitudinally extending flow for ice sliding over a wet melting bed and longitudinally compressive flow for ice sliding over a wet freezing bed. The increase of A caused by strain hardening tends to offset the decrease of A caused by the warmer ice over a thawed bed. These off setting effects of ice temperature and ice fabric are discussed in more detail by Reference HughesHughes (1998, p. 119–131). The competing effects of thermal and strain hardening or softening of ice are not fully quantified, so they cannot be incorporated with confidence into ice-sheet models.
Variations of B along flowlines are much more problematic. Whether B is computed from Equation (4) or some other formulation, it produces a τ W curve that falls below the τ D curve. Therefore f becomes a fundamental measure of instability in ice sheets. Many ice-sheet models use the mass-balance continuity equation to compute the fraction of ice velocity that is due to basal sliding, as opposed to creep in ice and bed deformation (e.g. Reference Fastook, Hansen and TakahashiFastook, 1984). The relationship between the thawed fraction and the sliding fraction, disregarding bed deformation, is highly non-linear (Reference HughesHughes, 1998, p. 157–158). When 50% of the bed is thawed, only 10% of the velocity is by sliding, and when 90% of the bed is thawed only 50% of the velocity is by sliding. Therefore, only a few frozen patches are highly effective in retarding sliding, so that runaway sliding occurs only when the last frozen patches become thawed. Since ice-stream surges develop from runaway sliding, and ice streams discharge some 90% of Antarctic ice, calculating f for sheet flow becomes a primary goal for assessing the stability of the Antarctic ice sheet, especially when f = 1 near ice-sheet margins where stream flow develops from sheet flow.
Discussion
Interpretation of basal thawed fraction variations
For steady-state sheet flow, a thawed bed is likely when, in order of importance, (1) a high geothermal heat flux warms basal ice; (2) thick overlying ice depresses the basal pressure-melting temperature and slows conduction of basal geothermal and frictional heat to the surface; (3) warmer surface ice at lower elevations slows conduction of basal heat to the surface; (4) low surface-accumulation rates lower the basal temperature gradient; (5) ice flowing faster down a sloping bed increases basal frictional heat and decreases the basal pressure-melting temperature at bedrock bumps, and vice versa for ice flowing slower up a sloping bed; and (6) rapidly increasing ice velocities due to strongly converging flow at the heads of ice streams slow conduction of basal heat to the surface by advecting cold ice downstream. Therefore, a thawed bed (f =1) can be expected wherever a preponderance of these conditions is present. Conversely, a frozen bed (f = 0) can be expected wherever a preponderance of these conditions is absent. Taking x positive upslope along flowlines, as in Equation (2), basal melting increases when ∂f/∂x increases, and basal freezing increases when ∂f/∂x decreases.
Departures from steady-state flow may be indicated when a preponderance of the six conditions noted above favor a thawed bed but f = 0, or favor a frozen bed but f = 1. This indicates transient flow conditions that cause a positive gradient ∂f/∂x to migrate upslope as ice streams become more active and down-draw more interior ice, or a negative gradient ∂f/∂x to migrate downslope as ice streams become less active and down-draw less interior ice.
These criteria have been applied to the Antarctic flowlines shown in Figure 4 in order to assess the state of equilibrium for various Antarctic ice-drainage systems.
Flow of East Antarctic ice toward the Ross Sea
The Ross Sea embayment was created when the western sector of the West Antarctic ice sheet collapsed during the Holocene (Reference Denton, Bockheim, Wilson and StuiverDenton and others, 1989; Reference Hall and DentonHall and Denton, 1996). Collapse would have steepened the surface slope of East Antarctic ice streams passing through the Transantarctic Mountains, thereby increasing their velocity and down- draw of their ice-drainage systems. The major ice streams affected in our study would have been Byrd, Mulock and David Glaciers. Of these, Byrd Glacier would have been affected most recently because it is the furthest south. It also is the fastest and has the largest ice-drainage system, which is comparable in size to the remaining West Antarctic ice sheet and extends to the ice divide between Dome Argus and Dome Circe (see Fig. 4). Being farther north, Mulock Glacier would have been destabilized somewhat earlier in the Holocene. It has a long narrow ice-drainage system that extends to Dome Circe. David Glacier is furthest north, so it would have been the first to destabilize. Unlike Byrd and Mulock Glaciers, David Glacier is not buttressed by the Ross Ice Shelf, so its instability may have persisted longer. Its ice-drainage system extends to the ice divide between Dome Circe and Talos Dome.
Figure 5 shows τ 0 variations caused by f variations along representative flowlines for the ice-drainage systems of Byrd, Mulock and David Glaciers, flowlines BY3, MU2 and DA2, respectively, in Figure 4. In general, high values of f along these flowlines correlate with conditions of thick ice, low surface temperatures (high surface elevations), low accumulation rates, a down-sloping bed and convergent ice flow (decreasing flow-band widths), and low values of f correlate with the opposite conditions, notably thin ice in the Gamburtsev Subglacial Mountains beneath Dome Argus and above bedrock sills at the heads of ice streams. Departure from these general trends may reflect unknown variations in the basal geothermal heat flux. If so, we can conclude that all three ice-drainage systems are presently close to steady-state equilibrium.
Flow of East Antarctic ice toward Oates Land
Flowlines OL1–OL5 in Figure 4 bring interior sheet flow from Wilkes Subglacial Basin right to the coast, where the bed rises to sea level. The bed is smooth and mostly below sea level. Variations of f along flowlines OL1–OL3 are entirely consistent with nearly steady-state flow, as shown for flowline OL2 in Figure 6. Flowlines OL4 and OL5 enter Cook Ice Shelf. In general, f variations deduced from the τ 0 variations are consistent with steady-state flow. However, the bed falls well below sea level near the ice-shelf grounding line, while f remains relatively low. This is an indication that a large ice stream once entered Cook Ice Shelf, but when it shut down its thawed bed began to freeze. Today, sheet flow extends virtually to the grounding line. These features are shown in Figure 6 for flowline OL5 in Figure 4.
Flow of East Antarctic ice toward George V Coast
Flow toward the George V Coast is shown by flowlines GV1–GV8 in Figure 4. Flowline GV1 enters Cook Ice Shelf, and the τ 0 variation requires a variation of f similar to that for flowline OL5 in Figure 6. However, most of the ice flows into Ninnis and Mertz Glaciers. Flow begins in Vincennes Subglacial Basin beneath Dome Circe, and crosses a series of subglacial highlands and basins before entering subglacial trenches beneath these ice streams. Variations in τ 0 reflecting variations in f are shown in Figure 7 for flowline GV5 entering Ninnis Glacier and flowline GV7 entering Mertz Glacier. In general, the f variations are consistent with nearly steady-state flow for all these flowlines.
Flow of East Antarctic ice in Wilkes Land
Wilkes Land includes most of East Antarctica between the northern Transantarctic Mountains and the Amery Ice Shelf. Flowlines W1–W19 begin at Dome Circe and cross Aurora Subglacial Basin to the coast in Figure 4. Sheet flow extends nearly to the coast in this region, yet Dibble Iceberg Tongue, Porpoise Bay and Dalton Iceberg Tongue are features that could have been associated with major ice streams that are now small or absent. On the map of Antarctica produced by the American Geographical Society in 1970, Dibble Iceberg Tongue extends beyond Dibble Glacier, which is a minor ice stream, and Dalton Iceberg Tongue extends beyond Reynolds Trough, which is a surface depression that may have been down-drawn by a former ice stream. Flowlines W7, W12 and W16 enter Dibble Glacier, Porpoise Bay and Reynolds Trough, respectively, and their longitudinal profiles are shown in Figure 8. Variations of f associated with τ 0 variations are consistent with nearly steady-state flow, except near the coast, where the steady-state correlation between thick ice and high f breaks down. This is consistent with Dibble Iceberg Tongue and possibly with Dalton Iceberg Tongue, which are both about 100 km long, having been produced by recent surges of large ice streams. Now only a small ice stream, Dibble Glacier, exists, and sheet flow with f close to steady state extends across Reynolds Trough. Similarly, a major ice stream, now gone, probably entered Porpoise Bay. Another departure from steady-state conditions seems to exist near the end of flowline W19, as seen in Figure 8, even though this flowline ends in Totten Glacier, a major ice stream now active. These observations are compatible with a recent departure from steady-state flow due to termination of stream flow near the ends of flowlines W7, W12 and W16, making ∂f/∂x negative as basal freezing freezes a thawed bed, and due to a recent initiation of stream flow from near the end of flowline W19, making ∂f/∂x positive as basal melting thaws a frozen bed. Both departures would associate low f with thick ice near the ends of these flowlines, as seen in Figure 8.
Flow of West Antarctic ice into Pine Island Bay
About one-third of the remaining West Antarctic ice sheet flows into Pine Island Bay through the two fastest Antarctic ice streams, Pine Island and Thwaites Glaciers, as seen in Figures 4 and 10. Flowline P4 extends from the lowest saddle (1250 m high) on the West Antarctic ice divide where f ≈ 0, to Pine Island Glacier where f ≈ 1. Flowline P4 in Figure 9 has the concave surface of stream flow, with P W/P I increasing to unity along the length of Pine Island Glacier. Across the saddle, ice flows into Rutford Ice Stream and the Ronne Ice Shelf. It is tempting to conclude that the saddle has been down-drawn by these two ice streams, with surface lowering causing the bed to freeze under the saddle. Flowlines P5–P7 reverse the steady-state correlation between a high f for a largely thawed bed and thick ice, especially over Byrd Subglacial Basin where f is low and ice is thick. Flowline P5 in Figure 9 illustrates these anomalies, which extend almost to the ice divide. It is therefore likely that nearly the entire ice drainage system of Pine Island Glacier is not in steady state. Flowline T3 begins in the Ellsworth Mountains, crosses Bentley Subglacial Trench and Byrd Subglacial Basin, and enters Thwaites Glacier. Thin ice leaving the Ellsworth Mountains flows over a frozen bed (f = 0), the bed becomes about half thawed (f ≈ 0.5) as ice thickens, but then the bed becomes nearly frozen (f ∼ 0) for thick ice over Byrd Subglacial Basin, before thawing (f ≈ 1) as ice enters Thwaites Glacier. This indicates steady-state flow until lowering ice being down-drawn by Thwaites Glacier causes the bed to freeze in Byrd Subglacial Basin. Flowline T4 extends from the highest saddle (1750 m high) on the West Antarctic ice divide (f ≈ 0.8) across Byrd Subglacial Basin (f ≈ 0.5), to enter Thwaites Glacier (f ≈ 1). Across this saddle, ice flows into Ice Stream D, which extends nearly to the ice divide and enters the Ross Ice Shelf. The high f and the almost concave profile along most of flowline T4 are compatible with a departure from steady-state equilibrium that extends virtually to the ice divide, but not substantial as seen along flowlines entering Pine Island Glacier.
Flow of West Antarctic ice into the Ross Ice Shelf
Variations of f were calculated for sheet-flow portions of flowlines R1–R13, which begin at the West Antarctic ice divide and enter the Ross Ice Shelf by way of Ice Streams E (R1 and R2), D (R3–R6), C (R7 and R8) B (R9–R11) and A (R12 and R13) in Figure 4. The concave surface of stream flow (f = 1, with P W/P I increasing downstream) extends nearly to the ice divide along all of these flowlines. Variations of f near the ice divide have no explanation that is consistent with steady-state flow. The most likely explanation is that stream flow has down-drawn so much ice that unstable flow extends clear to the ice divide. Alternatively, active Cenozoic volcanism along the ice divide may cause a highly variable geothermal heat flow, and this is what has destabilized the ice (Reference Blankenship, Bell, Hodge, Brozena, Behrendt and FinnBlankenship and others, 1993).
Linking basal thermal zones to ice-sheet stability
Figures 10, 12 and 13 present in map-plane format the model results that were computed and displayed in longitudinal profiles along flowlines, as represented by Figures 5–9. Figure 10 shows isopleths of basal thawed fraction f, expressed in percentage thawed at 25% intervals. Given the uncertainties associated with not being able to vary parameters A and B along flowlines, and with bed topography being interpolated between radio-echo flight-lines, not to mention the relationship between A and basal temperature and the relationship between B and basal sliding or bed deformation, Figure 10 gives an impression of detailed f variations that is probably not justified. Nonetheless, in East Antarctica there is a broad correlation between a largely frozen bed (low f) and thin ice, especially along ice divides, and a largely thawed bed (high f) and thick ice over subglacial basins. In addition, f tends to be low where ice streams begin at high bedrock sills, notably in the Transantarctic Mountains, and f tends to be high where ice streams begin in deep bedrock troughs. Also, there is a tendency for f to increase on a down-sloping bed and to decrease on an up-sloping bed. This is consistent with largely steady-state flow. Departures from these generalities are seen in ice entering Cook Ice Shelf, Dibble Glacier, Porpoise Bay, Totten Glacier and perhaps Reynolds Trough (which may be an ice surface that was down-drawn by a former ice-stream surge that produced Dalton Iceberg Tongue). In West Antarctica, virtually all flowlines have f variations that are consistent with substantial departures from steady-state flow.
Uncertainties in f are impossible to calculate, since A and B in the respective flow and sliding laws used in our model do not vary along flowlines, use of Equation (1) precludes a rugged bed (Reference HutterHutter, 1983) and the map of bed topography is only marginally accurate on the scale of our calculations (Reference DrewryDrewry, 1983). For these reasons, our confidence in f is limited to first-order variations along ice flowlines. We ignore second-order variations. This is illustrated in Figure 11 for variations of τ 0 along flowline P5. Second-order variations lie between the two thin lines, which give the first-order variations. Using this criterion, the uncertainty in f is about ±10%. This uncertainty can be used to delineate four broad regions of f with some confidence. These regions are for f ≈ 0 where the bed is largely frozen, f ≈ 1 where the bed is largely thawed, f increasing in the flow direction where ice is moving over an increasingly thawed bed (a melting basal thermal zone), and f decreasing in the flow direction where ice is moving over an increasingly frozen bed (a freezing basal thermal zone). These four regions are shown in Figure 12. Figure 12 shows the flowlines in Figure 4, with thin dashed portions showing where f is decreasing in the flow direction, thin solid portions showing where f is increasing in the flow direction, areas marked by parallel horizontal lines showing where f ≈ 1 in the interior and where stream flow replaces sheet flow near ice margins, and areas marked by parallel vertical lines showing where f ≈ 0, usually because ice is thin. This information is used in Figure 11 to denote boundaries between melting zones (M) and freezing zones (F) by thick black lines.
Fig. 11. Estimating uncertainties in thawed fraction f of the bed. Calculated basal shear stresses (triangles) along flowline P5 in Figures 4 and 9 are bracketed by thin lines which allow f to vary by ±10% between the theoretical τ0 variations for wholly frozen and wholly thawed beds (dashed lines).
Fig. 12. A generalization of the results in Figure 10 to show how basal thawed fraction f changes along flowlines in Figure 4. In the direction of flow along flowlines, portions where f decreases are shown as thin dashed lines, and portions where f increases are shown as thin solid lines. Regions where f = 0 are shown as parallel horizontal lines, and regions where f = 1 are shown as parallel vertical lines. Boundaries between basal melting zones (M), where f increases in the flow direction, and basal freezing zones (F), where f decreases in the flow direction, are shown as thick solid lines. These boundaries are related to first-order bed topography by showing regions below sea level as dotted areas.
Figure 12 reveals substantial symmetry in the distribution of basal thermal zones for the strongly divergent primary flow from Dome Circe and the strongly divergent secondary flow from the East Antarctic ice divide at approximately 135° E, 70° S. Dome Circe overlies Vincennes Subglacial Basin, and the bed is nearly wholly thawed (f ≈ 1) beneath thick overlying ice. Primary divergent ice flow then passes over a freezing bed (decreasing f) followed by a melting bed (increasing f) in all directions (west, north and east) on the way to the East Antarctic coast, with the freezing bed crossing subglacial highlands and the melting bed crossing subglacial basins (Wilkes Subglacial Basin to the west, Astrolabe Subglacial Basin to the north and Aurora Subglacial Basin to the east). This is all consistent with steady-state flow. Ice flowing northward from Dome Circe becomes a center for secondary divergent flow from Astrolabe Subglacial Basin (∼135° E, 70.5° S), where thick ice over a melting bed (increasing f) thins as it passes over a freezing bed (decreasing f) in the surrounding subglacial highlands, and thickens again as it passes over a melting bed (increasing f) in the subglacial lowlands along George V Coast and western Wilkes Land. This is also consistent with steady-state flow.
Figure 13 denotes dotted areas within thick black borders for the regions having substantial departures from steady-state flow. This map shows that, for the flowlines modeled in this study, the West Antarctic ice sheet has departed substantially from steady-state flow over a broad region between Pine Island Bay in the Amundsen Sea and the Ross Ice Shelf in the Ross Sea, whereas departures from steady-state flow in the East Antarctic ice sheet are restricted to ice entering Cook Ice Shelf, Dibble Glacier, Porpoise Bay, Totten Glacier and possibly Reynolds Trough. These are probable sites of present, in the case of Totten Glacier, or former ice-stream surges. The former surges that produced Dibble and Dalton Iceberg Tongues must have been quite recent, since iceberg tongues are transient features. It is noteworthy that these sites of unstable East Antarctic flow are for ice leaving Wilkes Subglacial Basin through deep troughs that underlie Cook Ice Shelf, for ice leaving Astrolabe Subglacial Basin through a shallower trough beneath Dibble Glacier and for ice leaving Aurora Subglacial Basin through troughs that underlie Porpoise Bay and Totten Glacier. However, no significant subglacial trough underlies Reynolds Trough on the surface, which therefore may have been down-drawn by a former ice stream that surged to produce Dalton Iceberg Tongue.
Fig. 13. Stability of the Antarctic ice sheet deduced from an assessment of steady-state flow along the flowline in Figure 4. Departure from steady-state flow are greatest in the dotted areas enclosed by the thick black lines. Dashed portions of the lines are uncertain boundaries for unstable flow.
The question arises, could these unstable East Antarctic ice margins expand into their subglacial basins to produce widespread gravitational collapse, as has occurred in West Antarctica during the Holocene and may be continuing in Byrd Subglacial Basin? That such widespread collapse is less likely in East Antarctica for the flowlines modeled in this study is seen in Figures 5–9. These figures show, as dashed lines, the surface and bed topography that would exist for complete glacioisostatic rebound of the bed. In contrast to West Antarctica, the bed would rebound to elevations near or above sea level along nearly all East Antarctic flowlines, including the critical ice margins where the marine ice-sheet instability mechanism investigated by Reference WeertmanWeertman (1974) and Reference ThomasThomas (1977) would be most active. Therefore, gravitational collapse of East Antarctic ice that was triggered by the marine instability mechanism would slow or stop as the bed rebounded toward sea level along these ice-sheet margins. However, other mechanisms for ice-stream surges may be able to cause substantial gravitational collapse in this sector of East Antarctica. The surge that produced Dalton Iceberg Tongue and possibly Reynolds Trough on the surface was apparently unrelated to the marine instability mechanism. On a much larger scale of gravitational collapse, a possible surge of Lambert Glacier has down-drawn >10% of East Antarctic ice to produce the Amery Ice Shelf, even though no broad subglacial basin exists in that sector of East Antarctica. The critical instability for these Antarctic ice streams may be the surface inflection point between sheet flow and stream flow, where the steep surface slope combines with the fast ice velocity to maximize production of basal meltwater, and thereby promote ice–bed uncoupling, that allows ice streams to reach farther into the ice sheet, as Reference HughesHughes (1975) proposed.
Conclusions
A simple model that relates surface slope to ice thickness, bed topography and ice–bed coupling has been presented, and then applied to map basal thermal zones for sheet flow in the Antarctic ice sheet. Four basal thermal zones are deduced: thawed, melting, freezing and frozen. Ice–bed coupling is a measure of basal friction, which we interpret to reflect the thawed fraction of the bed directly. If the correlation is indeed direct, then the model can be used to assess the stability of various sectors of the Antarctic ice sheet. However, a direct correlation applies only if the flow and sliding laws of ice are correct, and if bed deformation can be ignored. In addition, ice-hardness parameter A in the flow law and bed-coupling parameter B in the sliding law, including their exponents, must vary slowly along ice flowlines compared to basal thawed fraction f if our analysis is to hold.
Given the importance of f, what is needed to compute it accurately? As already discussed, values of A and B (including n and m) must be known with confidence. The competing effects of ice creep, basal sliding and bed deformation must also be known. This requires an understanding of subglacial hydrology for both hard and soft beds. A deep-drilling program is therefore necessary. Bed topography and morphology must be mapped in detail by radio-echo sounding and seismic sounding to obtain the extent and thickness of soft beds. The Antarctic folio maps (Reference DrewryDrewry, 1983) used in this study have inaccuracies in radio-echo flight-lines and reflections that allow errors of up to 50% in ice thickness. All of the required improvements must come from improved field data and improved thermomechanical models.
If A and B are not adequately represented by Equations (3) and (4), respectively, then f in Equation (5) is not a measure of the thawed fraction of the bed. It includes the thawed-bed fraction, but it is contaminated by our incomplete knowledge of thermal and strain softening in ice, basal sliding over both hard and soft beds, and deformation of soft beds. Nonetheless, our model calculates f values that are broadly consistent with the frozen (f = 0) and thawed (f = 1) basal regions in East and West Antarctica that are calculated by thermomechanical models, and that are shown to be in states of equilibrium similar to those deduced from our model (e.g. Reference HuybrechtsHuybrechts, 1990; Reference Hansen and GreveHansen and Greve, 1996; Reference FastookFastook, 1997; Reference Verbitsky and SaltzmanVerbitsky and Saltzman, 1997). This gives us confidence that f strongly represents the thawed-bed fraction, and that our model computes f with more than superficial accuracy. It should be emphasized that only two values of f really exist, f = 0 for a frozen bed and f = 1 for a thawed bed. Intermediate values of f are calculated only where there is a mosaic of frozen (f = 0) and thawed (f = 1) patches in basal areas wΔx in our model, within which statistical values 0 ≤ f ≤ 1 can be calculated. The flatter the bed, the less likely these mosaics become. Logically, f can be seen as a measure of basal friction that has two extremes, a high value for a frozen bed and a low value for a thawed bed. Our simple model is to thermomechanical models as a quarter horse is to a thoroughbred horse. It is good for a fast determination of basal conditions, but it lacks the physics (the breeding) to explain why these conditions are as they are.
Acknowledgements
The flow-band ice-sheet model used in this study was programmed by J. L. Fastook. We thank R. Greve and R. C. A. Hindmarsh for detailed reviews that greatly improved our manuscript. This work is a contribution to the Experimental Program to Stimulate Competitive Research (EPSCoR) of the U.S. National Science Foundation.
