Basal conditions of two Transantarctic Mountains outlet glaciers from observation-constrained diagnostic modelling

Abstract

We present a diagnostic glacier flowline model parameterized and constrained by new velocity data from ice-surface GPS installations and speckle tracking of TerraSAR-X satellite images, newly acquired airborne-radar data, and continental gridded datasets of topography and geothermal heat flux, in order to better understand two outlet glaciers of the East Antarctic ice sheet. Our observational data are employed as primary inputs to a modelling procedure that first calculates the basal thermal regime of each glacier, then iterates the basal sliding coefficient and deformation rate parameter until the fit of simulated to observed surface velocities is optimized. We find that the two glaciers have both frozen and thawed areas at their beds, facilitating partial sliding. Glacier flow arises from a balance between sliding and deformation that fluctuates along the length of each glacier, with the amount of sliding typically varying by up to two orders of magnitude but with deformation rates far more constant. Beardmore Glacier is warmer and faster-flowing than Skelton Glacier, but an up-glacier deepening bed at the grounding line, coupled with ice thicknesses close to flotation, lead us to infer a greater vulnerability of Skelton Glacier to grounding-line recession if affected by ocean-forced thinning and concomitant acceleration.

Introduction

Basal conditions of glaciers and ice sheets are largely unknown, except where subglacial access is possible (e.g. through drilling or tunnelling); even then only limited information is available at a small proportion of the glacier bed. Yet basal properties, such as thermal regime, the presence, absence and volume of meltwater, and the existence and type of sediment layer overlying bedrock, all influence glacier dynamics (Truffer and others, 2000; Smith and others, 2013; Van der Wel and others, 2013) and their response to environmental changes (Pattyn, 1999; Bamber and others, 2007). Consequently, geophysical techniques (e.g. ground-penetrating radar (GPR) and seismic sounding) are commonly employed to image glacier beds, but these too are logistically difficult to apply over large areas. Although reliant on a number of assumptions, numerical modelling offers an alternative approach to the investigation of subglacial properties. Such an approach has been employed previously in Antarctica to estimate the rate of basal meltwater production (Joughin and others, 2003, 2004), or to model the effect of changes in the basal regime in terms of ice-sheet response (Pattyn, 1996). Here we adopt a similar method to explore likely subglacial conditions beneath two outlet glaciers of the East Antarctic ice sheet (EAIS). The EAIS represents the largest body of fresh water on the planet, with a sea-level equivalent volume of 53.3 m (Fretwell and others, 2013), yet our knowledge of the glacier systems that discharge from its margins is extremely limited and their potential susceptibility to changing environmental boundary conditions therefore remains poorly quantified. With ~15% of EAIS discharge draining through Transantarctic Mountains (TAM) outlet glaciers, knowledge of their subglacial conditions, in terms of the rate and extent of basal sliding, is important for assessments of ice-sheet mass-balance variability and potential contribution to sea level under different climate scenarios.

In this paper we use a one-dimensional (1-D) glacier flowline model to simulate grounded portions of Beardmore and Skelton Glaciers under present-day conditions, upstream of their confluences with the Ross Ice Shelf (RIS; Fig. 1). Although surface velocity measurements of Beard-more Glacier and a seismic survey of the floating portion of Skelton Glacier were made in the 1960s (Swithinbank, 1963, 1969; Crary, 1966), no published data for the subglacial conditions of these two glaciers currently exist. To address this data gap, we focus on predicting likely basal ice temperatures and velocity fields that arise from basal sliding and creep deformation. For constraint, our diagnostic model uses new observational datasets of surface velocity and ice thickness, together with published gridded datasets. The velocity data comprise field-based differential GPS measurements, distributed velocity analyses from high-resolution TerraSAR-X (TSX) satellite data, and regional surface velocities from the MEaSUREs compilation (Rignot and others, 2011). Ice-thickness data were obtained from field surveys (25 MHz pulseEKKO system), from the CReSIS airborne-radar survey and from the Bedmap2 compilation (Fretwell and others, 2013). Surface climatologies were constructed from empirical data (see below), and geo-thermal heat flux values were taken from Shapiro and Ritzwoller (2004). Together these time-invariant data allow an iterative approach to be followed, from which an internally consistent flow regime emerges that is characterized by variable rates of creep deformation and of sliding at the glacier bed. These simulations enable inferences to be made about the dominant mechanisms of flow in these outlets, and whether differences exist between the two glaciers. Finally, we set these findings in a wider context by considering their implications for glacier responses to environmental perturbations, such as rising sea level or atmospheric or oceanic warming.

Fig. 1. Location maps of (a) Skelton Glacier and (b) Beardmore Glacier, showing ice-surface velocities from TSX speckle tracking, and from InSAR (interferometric synthetic aperture radar; Rignot and others, 2011). Locations of GPS installations, radar flight lines and model nodes of the simulations described in the text are also shown. Inset shows the continental context of the two glaciers. (RIS: Ross Ice Shelf; ASAID: Antarctic Surface Accumulation and Ice Discharge.)

Our results indicate that relatively thick (>1 km) sections of these glaciers are at the pressure-melting point and therefore have the potential to slide in these areas, if other conditions are conducive. Furthermore, we suggest that small (less than one order of magnitude) variations in ice deformation rates may account for along-flow differences in observed velocities in areas where sliding is absent. Considered in the context of their respective locations, bed geometries and ice thicknesses, we speculate that the inferred basal conditions may make the northern branch of Skelton Glacier vulnerable to thinning of the adjoining RIS, with grounding-line recession likely to accelerate across a deep subglacial basin. This sensitivity does not exist in the main Skelton Glacier outlet, where the grounding line is currently located downstream of a bedrock high. Beardmore Glacier is the fastest-flowing of our study glaciers, yet its location further south than Skelton Glacier, together with the considerable buttressing effect exerted by the RIS, most likely buffers it from changes taking place near the ice-shelf margins.

Skelton Glacier (78°°30’ S, 161°° E)

Skelton Glacier drains a catchment area of ~6000 km², and flows from Taylor Dome through the TAM in southern Victoria Land toward the RIS (Fig. 1a). Between its grounding line and the RIS, Skelton Glacier supports a floating ice shelf, 50 km long, that may be almost grounded (akin to an ice plain) in the 12 km immediately downstream of the grounding line (Crary, 1966). A suite of studies undertaken during the International Geophysical Year (1957–58) reported ice thicknesses, flow and meteorological data from Skelton Glacier, although much of the information concerned the floating portion of the glacier in Skelton inlet. At its confluence with the RIS this floating tongue had a measured thickness of 580 ± 10 to 590 ± 20 m (Crary, 1966; Swithin-bank, 1969) and a velocity of 0.28 m d⁻¹ (Cameron and Goldthwait, 1961), yielding an annual discharge of ~0.8 km³ (Wilson and Crary, 1961). Accumulation ranges from 0.14 m ice eq. a⁻¹ in the névé (Crary, 1966) to 0.06 m w.e. a⁻¹ at the RIS confluence (Wilson and Crary, 1961). Those authors inferred from their calculations that the East Antarctic ice sheet did not contribute ice to Skelton Glacier, but more recent velocity and topographic datasets imply that such an ice-sheet connection does in fact exist (Bamber and others, 2009; Rignot and others, 2011). Mean annual temperatures vary near-linearly with altitude from –22°C on the ice tongue (84 m a.s.l.) to –42°C in the névé (2300 m a.s.l.) (Crary, 1966). A seismic survey revealed that the inner part of Skelton inlet is much deeper than its seaward mouth, reaching a maximum depth of 1933 m b.s.l. (below sea level) and a grounding-line ice thickness of >1500 m in the northern branch of the glacier (Crary, 1966). By contrast, a bedrock high underpins the main (southern) outlet at its grounding line, giving rise to relatively high velocities, where thin ice forms an icefall.

Beardmore Glacier (84°°30’ S, 169°° E)

The first traverse of Beardmore Glacier proved it was connected to the plateau of East Antarctica, and that it represented a significant outlet of the ice sheet in the southern TAM (Shackleton, 1911). The glacier drains a catchment of ~90 000 km² (Marsh and others, 2013) and at its grounding line confluence with the RIS has a mean velocity of ~365 m a⁻¹ (Swithinbank, 1963; Marsh and others, 2014). Velocities at the grounding line oscillate according to both a diurnal and a fortnightly tidal cycle (Marsh and others, 2013). Where it becomes afloat, Beardmore Glacier is ~25 km across, and on its southern side overflows a bedrock high that gives rise to an icefall below Mount Kyffin. Large rifts develop where this section of the glacier becomes buoyant, possibly due to changes in its stress balance as it ‘falls over a submerged cliff’ (Collins and Swithinbank, 1968, p. 109). Whether this cliff is a bedrock feature or is composed of soft sediment, such as might be expected of a grounding zone wedge (Alley and others, 2007), is not known, but ice in the rifted zone was measured to be only half as thick as ice on either flank (Collins and Swithinbank, 1968). At the present grounding line, beyond the rifted zone, both GPR surveys and ice-flexure modelling indicate a characteristic ice thickness of ~900 m (Marsh and others, 2014). Accumulation measurements on Beardmore Glacier are sparse, but Mayewski and others (1990) reported 0.035 m a⁻¹ at 2700 m a.s.l. from the Dominion Range ice-core site (85°15^’ S, 166°10^’ E) near the head of the glacier. Mean annual air temperature at this site is –37.3°C, with a measured temperature of –31.3°C at 230 m depth (Mayewski and others, 1990). Near the mouth of Beardmore Glacier on the RIS, data from the ’Elaine’ automatic weather station (83.1° S, 174° E) indicate a mean annual air temperature of –25°C (Reusch and Alley, 2004), and field observations report an average accumulation rate of 0.28 m a⁻¹ (Denton and others, 1989). Recent airborne-radar data for Beardmore Glacier are not currently available, but a 2010 ground-based radar survey (25 MHz pulseEKKO system) confirmed a floating ice thickness at the ice-shelf junction of ~900 m, thickening to ~1000 m 1–2 km upstream of the grounding line (Marsh and others, 2014). These values are not incorporated into the recent Bedmap2 data compilation (Fretwell and others, 2013), but are consistent with it.

Empirical Constraints

In order to adequately simulate the outlet glaciers with the model described below, our method employs constraining velocity data from two sources. Field measurements of glacier surface velocities were obtained from GPS installations. These data then provided constraints for the calibration of a more extensive velocity field derived from speckle tracking of synthetic aperture radar (SAR) data. Both components are described more fully elsewhere (Marsh and others, 2013), hence only a summary is presented below. On both Beardmore and Skelton Glaciers we are confident that our velocity fields are representative of their mean rate of ice flow, partly because our values compare well with previous studies (e.g. Swithinbank, 1963), and partly because our 11 day TSX repeat cycle averages out high-frequency (e.g. diurnal) velocity fluctuations in the grounding zone, where variability is greatest.

GPS data

We instrumented the northern arm of Skelton Glacier with eight GPS stations, aligned parallel to the central flowline and positioned at and upstream of the grounding line (Fig. 1a). We were unable to install instruments along the main trunk of Skelton Glacier due to intense crevassing. On Beardmore Glacier, instrument availability restricted us to four GPS stations, which we installed along the centre line of the glacier 20 and 10 km upstream and 4 and 15 km downstream of the grounding line (Fig. 1b). On both glaciers, data were recorded at 15 s intervals. The measurement period on Skelton Glacier was 38 days; on Beardmore Glacier it was 35 days, although data are not continuous at some of the sites. Data were processed using kinematic Precise Point Positioning algorithms in the NASA Jet Propulsion Laboratory’s GIPSY software (Zumberge and others, 1997; King, 2004). We avoided artefacts and minimized error in horizontal position due to vertical ocean tidal effects, by kinematic processing under loose random-walk constraints with additional ambiguity resolution, although this technique results in some loss of precision compared with static occupation methods (King, 2004). Nonetheless, locational accuracy for individual positions from this method is typically 1.5–2.5 cm, which is more than adequate for our requirements. Diurnal velocity fluctuations were observed at both glaciers, with local velocity maxima occurring at or near the grounding zone and typically coincident with the maximum rate of falling tide (Marsh and others, 2013), in agreement with elastic flexure theory, which predicts a longitudinal extension at the glacier surface during low tide (Holdsworth, 1969). At the grounding line of Skelton Glacier, 6 hour mean velocities vary between 0.25 and 0.35 m d⁻¹ during spring tides, and on Beardmore Glacier between 0.5 and 1.5 m d⁻¹. These diurnal velocity variations are removed by linear regression over the entire dataset to obtain representative annual velocities for the two glaciers.

Satellite data

Speckle tracking makes use of the high pixel noise inherent in coherent imaging systems such as SAR. Speckle tracking has an advantage over optical feature tracking in that it does not require large-scale features (e.g. crevasses) and can be used to map displacement vectors over short periods (Joughin, 2002). TSX speckle tracking is used here to provide an additional spatial representation of the velocity pattern on the glaciers between GPS sites. Full details of the acquisition and processing of these data are presented elsewhere (Marsh and others, 2013), and are not repeated here. For Beardmore Glacier, 16 TSX image pairs with 11 day separation acquired during December 2009 provide full coverage of the glacier, while a further three pairs covering the grounding zone from June 2012 were used to create a stacked velocity map and to smooth the effect of short-term changes in this area. For Skelton Glacier, a mosaic of ten TSX image pairs from May 2012 was used (Fig. 1a and b).

CReSIS ice-thickness data

Although recent gridded data compilations of surface and bed elevations provide a good representation of ice thickness in well-studied sectors of Antarctica (Fretwell and others, 2013), many TAM glaciers have sparse data coverage and are thus susceptible to interpolation uncertainties. We collected new airborne-radar depth sounder (RDS) data for Skelton Glacier in order to better characterize its subglacial topography. In the CReSIS RDS system, ice thickness is derived from multichannel waveform data of different pulse durations. An estimated ice index of refraction of ice (√3.15) is used to convert the difference in propagation time between ice-surface and ice-bottom reflections into ice thickness. An assumption of ice homogeneity is made, and no firn correction is applied. Over-flight of the Taylor Dome ice-core location during the Skelton Glacier survey gave an RDS ice thickness of 550 ± 10 m, in very good agreement with that of 554.2 m determined from the ice core (Morse and others, 1999).

Numerical Model

The numerical glacier flowline model we employ is essentially identical to that of Golledge and Levy (2011), but here we use the model in a diagnostic (rather than prognostic, or time-dependent) manner in order to make explicit use of new observational data, as well as recently published gridded data. These data are employed as constraints in an iterative calculation of creep and sliding velocity fields that aim to characterize the present-day behaviour of the two studied glaciers, assuming that each is currently in a state of equilibrium. Below we describe the modifications made to the model to enable diagnostic simulations to be run and, where necessary, associated components of the model that are of particular relevance. For our model of Skelton Glacier we use a horizontal grid spacing of 500 m, but for the much larger Beardmore Glacier we operate at 1000 m resolution. These and all other model constants and variables used below are defined in Tables 1 and 2. We use CReSIS (Skelton Glacier ice thickness), Bedmap2 (Beardmore Glacier ice thickness), GPS-constrained TerraSAR-X (local velocity field) and MEaSUREs (regional velocity field; Rignot and others, 2011) datasets as model inputs, with native resolutions of 30, 1000, 60 and 900 m, respectively.

Table 1. Physical and model constants

Table 2. Model variables

We start with these known values for ice thickness, bed elevation and ice-surface velocities. Surface climatologies are based on the empirical observations described above, using temperature and precipitation lapse rates to account for altitudinal variations. Since our interest lies in establishing likely values of basal ice temperatures, creep and sliding velocities, we hold our prescribed ice geometry constant and iterate only two model parameters: our deformation-rate coefficient, f _d, and our sliding coefficient, f _s. In the first integration we set f _d according to a temperature-dependent relationship (Golledge and Levy, 2011):

(1)

incorporating the activation energy necessary for creep, Q =139 kJ mol⁻¹, the universal gas constant, R = 8. 314 J mol⁻¹ K⁻¹, and a thermal parameter, ξ = 5.47 ×10¹⁰ Pa⁻³ a⁻¹ (Hubbard, 2006). Below 263.15 K the value for creep activation energy is reduced to 60 kJ mol⁻¹ (Cuffey and Paterson, 2010) and ξ = 1.14 × 10⁻⁵Pa⁻³3 a⁻¹. A flow enhancement factor, E, is used to account for ice softening by processes not incorporated in the model and is the vertically averaged ice temperature.

The rate of basal sliding, f _s, is initially prescribed a uniform value (f _s0) reduced from published values (Budd and others, 1979; Oerlemans, 1997) to account for the polar (rather than temperate) conditions in our study. This value is subsequently modified by the basal ice temperature, T _b, in a way that decreases sliding speed exponentially below the pressure-melting point, α, in order to allow a numerically smoother transition between sliding and non-sliding areas (Hindmarsh and Le Meur, 2001):

(2)

It is only at the first time step that f _d and f _s are calculated in this way; subsequent solutions arise from the iterative procedure described below. In our model we do not incorporate subglacial water, and so assume that the iceoverburden pressure is equivalent to the full ice thickness, H. The pressure-melting temperature (K) is therefore given by α = 8.7 × 10^-4 H (Hindmarsh and Le Meur, 2001). In the first integration we calculate a steady-state basal temperature field, T _b, from surface temperature, T _s, geothermal heat flux, G, ice thickness, the conductivity of ice, κ _ice, and the pressure-melting temperature (the first three terms in Eqn (3)). Subsequent iterations, however, also include vertical ice velocity, w (from surface accumulation rates), and horizontal velocity solutions to calculate strain and frictional heating, and horizontal and vertical advection:

(3)

In Eqn (3), u _s denotes sliding velocity, is the vertically averaged horizontal velocity of the ice column and τ _b is basal shear stress. Temperature gradients in horizontal and vertical planes are represented, respectively, by the two partial derivatives of the fourth term of Eqn (3). Resulting temperatures are then used to calculate deformation and sliding velocities, u _d and u _s, which combine to give a centre-line velocity, u, for comparison with empirical data:

(4)

where n and p are, respectively, the flow law and sliding exponents. Although a shallow-shelf approximation (SSA) model is not appropriate for mountain glaciers, such as those described here, longitudinal extension and compression in TAM outlets means that gravitational driving stresses, τ _d, alone may not yield a velocity field that accurately replicates observed values (Kavanaugh and others, 2009; Golledge and Levy, 2011). To account for this shortcoming, we implement the scheme of Kamb and Echelmeyer (1986) and Echelmeyer and Kamb (1986), which allows a heuristic adjustment to be made to τ _b, so that they are smoothed according to an ice-thickness-dependent triangular averaging filter:

(5)

where represents longitudinal coupling length, the product of a longitudinal coupling coefficient and ice thickness. The integral spans a distance defined by the longitudinal averaging length, equal to 4 (Kamb and Echelmeyer, 1986). We incorporate a width term in our flux calculations and use a width-dependent shape factor, F, based on tabulated values for a parabolic glacier cross-profile (Paterson, 1994) to account for the effects of horizontal changes in valley width. Since our scheme does not calculate longitudinal stresses per se, but simply allows their effects to be incorporated in a computationally efficient manner, the results are influenced by the choice of parameterization. Thus, to illustrate the effect on our results of changing the coupling values in this scheme, we calculate velocity fields and sliding/deformation-rate parameters for a range of coupling coefficients. In our experiments we use values of 2, 5 and 10, based on those considered appropriate for valley glaciers (1–3) and ice sheets (4–10) (Kamb and Echelmeyer, 1986). Since we are only conducting diagnostic simulations, we are able to compare the effects of modifying the longitudinal coupling coefficient with respect to observed surface velocity gradients, and thus use this parameter in tuning our model to best reproduce the empirical data. Finally, for each of the model configurations we compare the sum of the two velocity components, u, with our observed velocity field, u _in, and obtain the percentage error, u _e:

(6)

Ice at Antarctic temperatures has a high viscosity, so high-magnitude variations in deformation rate are unlikely, whereas large and abrupt changes in sliding velocity are known to occur where basal substrates change or where basal water pressures fluctuate (Stearns and others, 2008; Graham and others, 2009). On this basis we assume that, in areas where basal ice temperatures are close to the pressure-melting point, our initial velocity mismatch may most likely arise from incorrect parameterization of our sliding coefficient, f _s. Thus, for these gridpoints we enter a standard optimization procedure of the type ||x – x*|| ≤ δ in which x – x* defines the misfit from the target, δ. In our implementation, the magnitude of our sliding coefficient is incrementally adjusted so that the resulting sliding velocities evolve toward values in the observed velocity data. This is carried out at each gridpoint according to

(7)

where β represents the velocity misfit as a proportion of u _s:

(8)

After each iteration we calculate an updated sliding velocity from the new sliding coefficient and revise the associated velocity error, u _e (Eqn (6)). Both the rate of basal sliding and the rate of ice deformation are strongly dependent on temperature, and, since a portion of our calculated basal ice temperature comes from glaciological components (strain and frictional heating, Eqn (3)), we next update our basal temperature field according to the new sliding velocities, and recalculate our temperature-dependent deformation rate parameter, f _d. Because of this close interplay between temperature and velocity solutions, the rate of convergence of our iteration is explicitly controlled; in Eqn (7) we use the 10^–2 multiplier to limit incremental sliding parameter adjustments and thereby allow each recalculation of temperature and velocity to evolve smoothly. This iteration proceeds until u _e < tol, where tol represents a prescribed fitting tolerance (10%). At this point, all ‘sliding’ cells are closely fitted to the empirical velocities at each location. Remaining velocity mismatches may be the result of errors in our modelled glacier geometry (thickness, width), or a consequence of any combination of the necessary approximations employed in our simple model. However, it is similarly possible that some, or all, of the remaining velocity mismatch arises from variations in creep deformation rate, as a consequence of differences in debris content of the ice, fabric anisotropy, gas bubble or interstitial water content. To assess the magnitude of creep rate parameter variability necessary to explain all of the remaining mismatch, we rearrange the standard equation for creep velocity (Cuffey and Paterson, 2010) and solve for a new deformation rate value:

(9)

in which u _diff = u _in – u. Since we cannot know a priori the relative balance of sliding and creep deformation, we use this iterative approach to calculate the maximum possible sliding velocity and corresponding deformation velocity, as well as a zero-sliding ‘deformation only’ velocity scenario. The difference between these two end members defines the parameter space of uncertainty in our velocity solution. Despite the uncertainties and the simplifications inherent in our approach, our iterated values for f _d and f _s allow a combined velocity field to be calculated that matches observed surface values, and which also yields further information on the possible balance of mechanisms that drive glacier flow.

Results

Glacier geometries and dynamics

Since we constrain our modelled glacier geometries to prescribed values, the principal aim of our iteration procedure is to derive a velocity field that replicates observations of surface flow speed. Our first step is to calculate deformation velocities for each of our domains using a constant (temperature-independent) deformation rate parameter, and preclude sliding. In slow-moving (<20 m a⁻¹) sections of each glacier, velocities can be reproduced relatively closely with a range of values of f _d. However, higher velocities are only accurately simulated in a narrow range of values of f _d, so a close fit along the length of each domain with a single deformation rate parameter value is not possible. A best fit is found with deformation rate coefficients of 4 × 10^-17 for Skelton Glacier and 7 × 10^-16 for Beardmore Glacier, but the considerable scatter in these results (Fig. 2) indicates that spatially variable sliding and/or deformation rate coefficients are necessary to accurately reproduce observed surface velocities. Thus our second approach is to use the optimization scheme to improve the fit of modelled velocities to those used for constraint (Fig. 2). There are outliers, however, in each domain, reflecting an incomplete fitting to constraining data. This arises from our choice of a domain-averaged mismatch tolerance of 10% in the iteration scheme, chosen to avoid over-fitting and the possible introduction of unphysical jumps in along-flow speed. Nonetheless, since these outliers are few, the overall fit of our model to observed data is good and all domain results exhibit correlation coefficients, r ², of 0.99. A spatial view of glacier thickness variation and along-flow modelled velocities of Skelton Glacier is shown in Figure 3, together with velocity data from TSX and in situ GPS measurements, as well as continental gridded velocity data (MEaSUREs; Rignot and others, 2011). High-resolution CReSIS airborne-radar data reveal that the main glacier bed is characterized by an extensive overdeepening in the upper reaches of the glacier (Fig. 3a), where it drains from Taylor Dome into Skelton névé. A second overdeepening occurs 10–25 km upstream of the grounding line, terminated by an abrupt bedrock high ~400 m higher than the deepest part of the basin. More dramatic evidence of subglacial overdeepening is apparent upstream of the grounding line of the northern arm of Skelton Glacier (Fig. 3b). Here the bed depth drops rapidly from near sea level to ~1500 m b.s.l. in a distance of 11 km. This deep subglacial basin was noted by earlier surveys (Crary, 1966), yet the relative paucity of ice-thickness data in this area has resulted in continental interpolations that significantly underestimate ice thickness here (Fretwell and others, 2013). In general, Bedmap2 data fare best in the upper reaches of Skelton Glacier, with increasingly large deviation from CReSIS data in the lower 50 km. TSX and MEaSUREs data both reveal an exponential increase in down-glacier flow speed in the main trunk of the glacier, reaching a maximum of ~268 m a⁻¹ immediately downstream of a bedrock high that gives rise to an icefall. This acceleration due to thinning and steepening of the glacier is consistent with observed and modelled behaviour of other TAM outlets (Kavanaugh and others, 2009; Golledge and Levy, 2011). In the northern arm of Skelton Glacier, peak flow speed is 115 m a⁻¹, much lower than that of the main trunk, and down-glacier acceleration is less uniform, marked by a rapid increase in speed in the upper 10 km of the domain, and a more gradual increase from there to the grounding line.

Fig. 2. Observed versus modelled glacier surface velocities for (a) Skelton Glacier trunk, (b) Skelton Glacier northern arm and (c) Beardmore Glacier. Blue points denote velocity solutions when a constant deformation rate parameter, f _d, is applied to each domain and basal sliding is prevented. Brown points arise from experiments in which both the basal sliding coefficient and the deformation rate parameter are allowed to vary along-flow. Note the different scales in each panel.

Fig. 3. Measured ice surface, bed topography and surface velocity profiles for Skelton Glacier: (a) trunk and (b) northern arm. Also shown in (b) are the locations and magnitudes of surface velocities from GPS installations. Since our TSX velocity data do not extend across the entire domain in (a), we linearly combine these with values from Rignot and others (2011) and use the result as our constraining dataset. Velocities from the iterative modelling procedure are also shown.

Although CReSIS data are not currently available for Beardmore Glacier, a 2010 ground-based radar survey across the grounding line gave ice thicknesses consistent with Bedmap2 values (Marsh and others, 2014). Since the greatest errors in the Skelton Glacier bed data appear to occur at or near the grounding line, we are sufficiently encouraged to use Bedmap2 bed topography for our second simulation, accepting that these data may under-represent some of the subglacial variability in bed elevation. We define our domain as extending 281 km up-glacier of the modern grounding line to the margins of the East Antarctic plateau. Ice-surface and bed elevation data reveal a glacier of almost uniform thickness along this length (Fretwell and others, 2013), with the exception of a 300 m bedrock high 167 km from the grounding line, and a shallow (100–200 m deep) depression 10 km upstream of the grounding line (Fig. 4). MEaSUREs and TSX surface velocity data show much higher and more variable flow speeds than evident on Skelton Glacier. A maximum speed of 364. 7 ± 1. 5 m a⁻¹ (Swithinbank, 1963; Marsh and others, 2013) is associated with a steepening and thinning section of the glacier immediately upstream of the grounding line. Further up-glacier, speed decreases approximately exponentially, but is marked by localized accelerations and decelerations that appear to be more closely related to changing surface gradients than to ice-thickness variation (Fig. 4).

Fig. 4. Beardmore Glacier geometry from Bedmap2 data (Fretwell and others, 2013), and velocity interpretations from TerraSAR-X, InSAR (Rignot and others, 2011) and surface GPS installations. Velocities from the iterative modelling procedure are also shown.

Thermal regime and sliding/creep partitioning

Along the main trunk of Skelton Glacier, our model predicts basal ice temperatures close to, or at, the pressure-melting point in the deep upper basin in the névé (78–95 km upstream of the grounding line), in the shallower basin 9–23 km upstream of the grounding line and in the final 3 km adjacent to the grounding line (Fig. 5a, upper panel). In the shorter, northern, arm of Skelton Glacier, the relatively thinner ice characterizing much of the glacier precludes basal melting, except in the lowermost 13 km where ice thickness increases rapidly down-glacier, due to the presence of the very deep basin at the grounding line (Fig. 5b). In both of our Skelton Glacier domains, modelled basal temperatures at the majority of nodes are almost identical to values arising solely from the geothermal heat flux, indicating little or no contribution from internal (glacio-logical) processes. In these areas, strain and frictional heating as well as vertically and horizontally advected heat inputs are minimal. Compared with basal temperatures calculated with a 5 km resolution whole-continent model (Pattyn, 2010), values for the Skelton Glacier centre line compare reasonably well in some sections, but a notable deviation is that our calculation predicts warmer ice in the upper basin. In the northern arm, high-magnitude variations in the Pattyn data contrast with our predicted basal temperature field, that warms gradually towards the grounding line. The consequences of our predicted thermal fields in the two Skelton Glacier domains are shown in the lower panels of Figure 5a and b. Where basal temperatures are close to the melting point, basal sliding becomes possible and modelled glacier velocities no longer arise entirely from creep. Deformation velocities drop where sliding velocities increase, yet the variability in deformation rate required to account for this partitioning is less than an order of magnitude. In both of our Skelton Glacier domains our creep parameter varies little along the length of the glacier, with the exception of much higher (and possibly spurious) values near the ice divides. By contrast, abrupt increases in sliding coefficient of one to two orders of magnitude are necessary to compensate for observed velocity mismatches in areas where sliding is possible. The model is only moderately sensitive to our choice of longitudinal coupling coefficient, with sliding velocities predicted by our range of coupling values varying by much less than an order of magnitude. Importantly, although longer coupling lengths tend to smooth out along-flow velocity gradients, the choice of coupling coefficient does not change the location or extent of basal sliding.

Fig. 5. Surface (T _s) and modelled basal (T _b) ice temperatures, sliding (u _s) and creep (u _d) velocities, and sliding (f _s) and creep (f _d) parameter values arising from the iterative modelling procedure, for (a) the Skelton Glacier trunk and (b) Skelton Glacier northern arm. Grey shading defines the range of possible deformation velocity and creep parameter values that accompany maximum and minimum basal sliding scenarios. Basal ice temperatures from Pattyn (2010) are shown for comparison. Dotted and dashed lines show the relatively low sensitivity of simulated velocities to different longitudinal coupling parameterizations.

In our model of Beardmore Glacier the simulated pattern of flow is quite different to that of Skelton Glacier. The generally thicker, much faster-flowing, Beardmore Glacier has a basal thermal field that is very much warmer than that of Skelton Glacier, with coldest values (–13°C) associated with the thinnest part of the glacier (Fig. 6). Our temperature field contrasts markedly with that of Pattyn (2010), whose values range from –32°C to the pressure-melting point. Due to the much higher predicted ice temperatures in our simulation of Beardmore Glacier, compared with Skelton Glacier, deformation rates are much higher, so the relative contribution made by basal sliding to total (surface) velocities is low (maximum 13%, majority <10%). Our range of longitudinal coupling coefficients yield much greater variation in the predicted deformation rate parameter along the length of Beardmore Glacier than for Skelton Glacier (e.g. with changes of one to two orders of magnitude occurring within a horizontal distance of 10–20 km for a coupling coefficient of 5 (Fig. 6, lower panel)). Since there is no physical reason to suspect such high variability, we infer a greater along-flow smoothing of the overall stress balance of the glacier, most likely due to more effective longitudinal coupling in the thicker, wider and faster-flowing Beardmore Glacier, which we approximate in a parameterized way.

Fig. 6. Surface (T _s) and modelled basal (T _b) ice temperatures, sliding (u _s) and creep (u _d) velocities, and sliding (f _s) and creep (f _d) parameter values arising from the iterative modelling procedure, for Beardmore Glacier. Basal ice temperatures from Pattyn (2010) are shown for comparison. Note the higher variability of the creep parameter than on Skelton Glacier (Fig. 5), and the greater smoothing influence of higher longitudinal coupling lengths. The relatively warm ice of Beardmore Glacier leads to high deformation rates, thus maximum and minimum basal sliding scenarios are almost identical and little variation in deformation velocity and creep parameter occurs. Sliding and deformation rate parameters are, however, more sensitive to the choice of longitudinal coupling coefficient than on Skelton Glacier (Fig. 5).

Acknowledging the simplifications inherent in our scheme, we nonetheless consider a coupling length of ten times the ice thickness, rather than five times, to be more appropriate. Under such a parameterization, down-glacier fluctuation in creep rate parameter values is much reduced, which we consider to be more plausible.

Glacier sensitivity

Our simulations have highlighted key differences in the likely basal properties of our two study glaciers, and marked contrasts in the balance of flow mechanisms that arise as a consequence. Considering our steady-state results in terms of mass conservation, we can make further inferences about possible changes that may arise under perturbed environmental conditions. Thinning at glacier grounding lines currently represents the most significant forcing affecting marine-terminating Antarctic glaciers (Rignot and Jacobs, 2002; Walker and others, 2008). In order to maintain the present flux balance of our two glaciers, a thinning of 5% at their grounding lines would be compensated for by a 5% increase in velocity, and 10% would lead to an acceleration of ~11%. However, all our grounding-line cells are currently close to flotation thickness: the Skelton Glacier centre-line terminal grounded cell is 110 m above the flotation threshold (equivalent to 9% of its present ice thickness), Skelton Glacier north arm appears to be 12 m (1%) below flotation, presumably maintained by flexural rigidity, and Beardmore Glacier is 20 m (2.5%) above flotation. These values are approximate and do not account for the presence of firn, but nonetheless we infer that thinning and acceleration of these glaciers would most likely affect the lightly grounded north arm of Skelton Glacier most significantly, Beardmore Glacier less so, and Skelton Glacier centre line the least.

In terms of perturbations to atmospheric temperatures, we can assess likely impacts in terms of the percentage change in the amount of the glacier bed at (or close to) the pressure-melting point. Assuming no change in ice thickness, we calculate that air-temperature changes at the glacier surface of 0–5 K would lead to quite different responses from each glacier (Fig. 7). As the warmest of our examples, Beardmore Glacier exhibits the greatest potential sensitivity to atmospheric warming, with a steep and almost linear increase in potential basal sliding area with rising air temperature. The centre line of Skelton Glacier shows a similarly linear response, but with a much lower gradient, whereas the north arm of Skelton Glacier illustrates a nonlinear response, characterized by a low gradient at <3 K, and a much steeper slope at >3 K. This most likely reflects the predominantly cold-based configuration of this section of Skelton Glacier, with approximately two-thirds of the bed well below pressure-melting temperatures. The impact of such changes in atmospheric temperatures would, however, be controlled by the rate at which such warming could propagate to the glacier beds. Vertical advection of snow falling under warming conditions would proceed slowly, since accumulation rates are extremely low on both glaciers (0.035–0.28 m a⁻¹), so warming by this mechanism alone would require several millennia to affect the bed of either glacier. Our estimates of glacier sensitivity to atmospheric temperature changes might be considered most relevant over geological timescales, such as glacial/interglacial transitions through the Plio-Pleistocene (5 Ma to present), rather than over coming centuries.

Fig. 7. Percentage increases from present of glacier bed area at pressure-melting point (pmp), for a range of uniform increases in air temperature at the ice surface, assuming static geometry.

Discussion

The experiments described above aimed to characterize basal conditions of two TAM outlet glaciers, and to use these inferences to better understand the controls on their behaviour, both at present and under environmental perturbations. The numerical model takes only one horizontal dimension into account, and is relatively simple in its formulation, making use of many common glaciological assumptions and approximations. Our calculations of basal temperatures depend largely on interpolated geothermal heat-flux fields, which are poorly constrained at a local scale and, along with other input data, may be a potential source of error. Alternatively, errors may arise from the model formulation itself. Like all finite-difference integrations, ours are vulnerable to numerical errors, and our use of a width-dependent shape factor and a longitudinal-averaging scheme may lead to errors in our calculated stress field and the velocities that arise from it. Since our iteration scheme relies on the sequential solution of two sets of equations that continually update sliding and deformation velocities until convergence occurs, we calculate an envelope of possible solutions. These range from a scenario in which no sliding occurs anywhere beneath each glacier, to one where sliding occurs wherever basal temperatures allow. This approach allows a ‘most-likely’ scenario to be estimated, but a further way to gauge the reliability of our results is to compare them with other modelling studies. The three-dimensional simulation of Pattyn (2010) yields basal temperature predictions for the entire Antarctic continent, and therefore provides a wider context for the interpretation of our calculations. Given the greater complexity of the Pattyn (2010) model, it is encouraging that our more simple approach finds agreement in some areas. Where significant mismatch occurs, a number of aspects may be responsible. For example, in the northern arm of Skelton Glacier, the newly acquired CReSIS ice-thickness data indicate a much deeper bed than shown by compilations of continent-wide gridded topography (Fig. 3). Consequently, we would expect warmer basal ice than would be predicted beneath the thinner glacier implied by the gridded, lower-resolution, topography. However, mismatches also occur in our simulations of Beardmore Glacier, where similar datasets were used (Bedmap by Pattyn, 2010, and Bedmap2 in our study). In other areas, such as in the upper basins of Skelton and Beardmore Glaciers, the reason for mismatch may be that our surface accumulation rate (interpreted from observations) may differ from the regional climate model data used by Pattyn (2010). In light of the results described above, it is apparent that the greatest contribution to basal temperatures comes from the prescribed magnitude of geothermal heat flux, so it seems most likely that the differences between our temperature predictions and those of Pattyn (2010) might arise from this source. We use a single value for each of our domains, whereas Pattyn (2010) employs spatially variable values. Assumptions and simplifications made in our model may therefore influence our results to some extent, but conversely, the use of high-resolution empirical data to constrain our models provides a robust foundation for our simulations. Critically, one of our aims is to compare and contrast the two glaciers, and for this purpose we consider our standardized approach valid.

The two glaciers are very different. Skelton Glacier is relatively short, is fed by a local dome of the EAIS, and flows slowly (<100 m a⁻¹) along much of its length. By contrast, Beardmore Glacier is at least twice as long, drains a significant volume of ice from the EAIS, and flows quickly (>100 m a⁻¹) along 80% of its length in our model domain. Applying local climatic conditions but the same initial glaciological parameterization to our two study glaciers has revealed marked differences in their present-day basal regime and in the balance of mechanisms driving flow. A much larger proportion of the bed of Beardmore Glacier is at the pressure-melting point than that of Skelton Glacier, facilitating more extensive basal sliding and higher deformation velocities. The localized ‘slippery spots’ predicted in our model of Skelton Glacier coincide with overdeepenings in the bed, where the greater ice thicknesses favour basal melt. Frictional heating at the bed, and longitudinal coupling within the glacier both appear to be more important in Beardmore Glacier than in Skelton Glacier. Furthermore, Beardmore Glacier appears to be more sensitive to changes in atmospheric temperatures, albeit over millennial time-scales. Thus at a cursory level it seems that Beardmore Glacier might be more responsive to environmental perturbations than Skelton Glacier. However, our simulation predicts that the ice immediately upstream of the Beardmore Glacier grounding line is flowing largely by internal deformation, with minimal basal sliding, due to a reduction in basal temperature arising from thinning ice. This is consistent with inferences from ice-flexure modelling, which found that inclusion of basal sliding was not necessary to explain observed grounding-line velocities (Marsh and others, 2013). This section of frozen bed may, to some extent, buffer the majority of the glacier from oceanic forcings, whilst the higher southern latitude of Beardmore Glacier compared with Skelton Glacier, and consequently lower annual temperatures, may also offset its greater sensitivity to atmospheric warming.

On Skelton Glacier, our model predicts a frozen bed under much of the length of both the main trunk and the northern arm, but a thawed bed, and hence flow by basal sliding, at and immediately upstream of both grounding lines. In the main trunk, sliding occurs within 20 km of the zone of flotation, with the exception of the area underlain by a bedrock high that gives rise to a prominent icefall. The northern arm of Skelton Glacier currently becomes afloat at the seaward edge of a deep basin >1500 m b.s.l. We predict flow by basal sliding in the lower 13 km of grounded ice adjacent to the present grounding line, in an area partly underlain by a reverse-sloping bed. Although the reverse-sloping area is relatively short, the ice thickness here is so large that the horizontal flux through this section is significant. Taking into account conclusions from a previous survey that the ice downstream of here is only lightly grounded (Crary, 1966), as well as our calculation that the current grounding zone is at or even slightly below flotation thickness, we therefore suggest that this section of Skelton Glacier might be vulnerable to changes at its oceanic boundary, specifically, ocean warming or a rise in sea level that might cause sufficient thinning to precipitate rapid grounding-line retreat across the basin, in a manner similar to that observed elsewhere in Antarctica (Jenkins and others, 2010; Park and others, 2013). The long floating tongue of Skelton Glacier, and the narrow fjord in which it is confined, may buttress the glacier, however. Our flowline simulation does not incorporate the lateral drag exerted by valley side-walls on the floating ice tongue, so we are not able to assess the influence of such effects on grounding-line stability, but recognize that they may be significant (Gudmundsson and others, 2012; Gudmundsson, 2013). Overall, however, we consider Skelton Glacier more vulnerable to ocean-forced mass imbalance than Beardmore Glacier, due to its greater proximity to the margin of the RIS and thus the much smaller protection afforded by buttressing effects. Beardmore Glacier is situated far enough south that it is only likely to respond to oceanic changes if a significant portion of the RIS is removed.

Conclusions

We have used field measurements of glacier velocity and ice thickness, together with local climate data, satellite-derived velocity interpretations, continental-scale interpolations of bed topography, ice thickness and geothermal heat flux, to parameterize and constrain a diagnostic model of two outlet glaciers of the East Antarctic ice sheet in the Transantarctic Mountains. Focusing on prediction of their basal thermal regime and the pattern and style of flow that arises as a consequence, we have shown that Beardmore Glacier exhibits more widespread ‘wet-based’ conditions than Skelton Glacier, and is most likely governed by basal sliding and longitudinal coupling to a greater extent than Skelton Glacier. However, after combining our results with inferences from previous studies and interpretations of newly acquired geophysical data, we suggest that although Beard-more Glacier is considerably warmer and faster-flowing than Skelton Glacier, it is the latter whose lower reaches are more likely to be vulnerable to changes in oceanic conditions, due to its closer proximity to the margin of the RIS.