The bedrock topography of Starbuck Glacier, Antarctic Peninsula, as determined by radio-echo soundings and flow modeling

Abstract

A glacier-wide ice-thickness distribution and bedrock topography is presented for Starbuck Glacier, Antarctic Peninsula. The results are based on 90 km of ground-based radio-echo sounding lines collected during the 2012/13 field season. Cross-validation with ice-thickness measurements provided by NASA's IceBridge project reveals excellent agreement. Glacier-wide estimates are derived using a model that calculates distributed ice thickness, calibrated with the radio-echo soundings. Additional constraints are obtained from in situ ice flow-speed measurements and the surface topography. The results indicate a reverse-sloped bed extending from a riegel occurring ~ 5 km upstream of the current grounding line. The deepest parts of the glacier are as much as 500 m below sea level. The calculated total volume of 80.7 ± 7.2 km³ corresponds to an average ice thickness of 312 ± 30 m.

1. Introduction

Continent-wide knowledge of the bedrock topography of Antarctica has recently improved, with the publication of the BEDMAP2 compilation (Fretwell and others, 2013). However, large uncertainties remain at local and regional scales. This is especially true for the Antarctic Peninsula, where the complex topography and its small-scale variability hamper a reliable extrapolation of the ice-thickness measurements available to date. The most obvious way to alleviate this problem is to increase the spatial coverage of direct measurements, although past surveys have shown that both data collection and interpretation can be difficult in the region (e.g. Glasser and others, 2011; Scambos and others, 2011; Farinotti and others, 2013).

Interest in the glaciers flowing into the ice shelves fringing the Antarctic Peninsula has grown following the collapse of the Larsen B ice shelf in 2002 (e.g. Rott and others, 2002; Cook and others, 2005; Cook and Vaughan, 2010). In fact, an increase in flow speed by up to 600% was observed following the removal of the buttressing effect of the floating ice (Rignot and others, 2004; Scambos and others, 2004). This highlights the importance of the dynamical interaction between ice shelves and tributary glaciers (e.g. Gudmunds-son, 2013). Much effort has been invested since the collapse to better understand the processes involved in the disintegration of the ice shelves (e.g. Van den Broeke, 2005; Glasser and Scambos, 2008; Rott and others, 2011). However, our knowledge is still incomplete, hampered as it is by the poorly constrained ice-thickness distribution which often precludes a sound set-up of boundary conditions in model experiments, and forces modellers to seek work-around solutions. With the goal of assessing the contribution of the Antarctic Peninsula to sea-level change until 2200 for example, Barrand and others (2013) overcame the problem by directly imposing grounding-line retreat and ice-shelf break-up scenarios. Better knowledge of the ice-thickness distribution would, moreover, not only benefit future projections, but also improve estimates of current mass fluxes of the Antarctic Peninsula (e.g. Osmanoğlu and others, 2013).

An area that is currently receiving particular attention in this context is Scar Inlet, the remaining part of the Larsen B ice shelf (e.g. Shuman and others, 2011; Berthier and others, 2012). Since the area is expected to break apart in the next few years, it provides an ideal platform for improving the understanding of the mechanisms involved in ice-shelf collapse and associated feedbacks. In this paper, measurements from ground-based radio-echo soundings (RES) are presented for Starbuck Glacier (658370 S, 628230 W; Fig. 1a) one of the tributaries of Scar Inlet, and used to derive a complete bedrock topography. Measurements are compared and cross-validated against ice-thickness data provided by NASA's Operation IceBridge. The work is part of an ongoing effort to better characterize the surface and subsurface topography of the region (Cook and others, 2012; Farinotti and others, 2013).

Fig. 1. (a) Location of Starbuck Glacier (red rectangle and arrow), as seen in the Landsat Image Mosaic of Antarctica (LIMA; Bindschadler and others, 2008). (b) Ground tracks for which RES measurements were collected (black) and successfully interpreted (red), as well as IceBridge flight tracks for 14 November 2011 (grey) and sections with actual ice-thickness data (green). The positions of the temporary GPS base station (orange), and the continuous GPS station (blue) are shown. The black dotted line corresponds to the approximate position of the grounding line (Rignot and others, 2011). Spatial references (crosses) are in polar stereographic coordinates (standard parallel 718 S; origin at South Pole; units km).

2. Data

RES data were collected between 4 and 1 1 December 2012, using the British Antarctic Survey Deep-Look Radio-Echo Sounder (DELORES). DELORES is a monopulse radar system operating in the 1–4 MHz range (central frequency 3 MHz) based on designs by Gades (1998) and Welch and Jacobel (2003). The system has been used before in several surveys, mainly focusing on ice streams in West Antarctica (e.g. King, 2009, 2011; King and others, 2009; Woodward and King, 2009). The DELORES transmitter consists of a ±2 kV pulser with a variable firing rate ranging between 1 and 5 kHz. The receiver consists of a chassis computer using a 250MHz digitizing card. Both transmitter and receiver antennas are resistively loaded wire dipoles of 20 m half-length. For the survey, transmitter and receiver were towed by skidoo, separated by a distance of 105 m. Individual radargram traces were produced by stacking 64 radar shots, recorded over a time interval of 15 ms. Survey speed ranged from -10 to 15kmh-¹, resulting in a spatial trace separation of <0.5m. For the days before 8 December, the geographical positions of the individual traces were recorded using differential GPS techniques and Leica 1200 dual-frequency devices. To this end, a base station was installed at 65˚36'25.3"S, 62˚19'19.6"W (orange dot in Fig. 1b) and a second receiver was towed midway between the transmitter and receiver antennas. Horizontal and vertical position accuracy were estimated as <0.2 m and <0.6m, respectively. Due to a receiver firmware failure, positions starting from 9 December had to be recovered from a single-frequency GPS receiver carried on board the towing skidoo. Cross-validation with the data collected prior to 9 December showed that the original positions could be reconstructed with an average horizontal deviation of 2.6 m (75% and 95% of the deviations contained within 3.6 m and 54 m, respectively) and an average vertical deviation of 1.5 m (75% and 95% of the deviations contained within 2.3 m and 9.4 m, respectively). In total, 151 km of radar data were collected (black tracks in Fig. 1b). Bedrock reflections could be successfully detected over 90 km (60% of the total ground covered) along 39 different profiles (red dots in Fig. 1b). The processing of the radar data is described in Section 3.

Additional ice-thickness data, collected on 11 November 2011 by NASA's Operation IceBridge Multichannel Coherent Radar Depth Sounder (MCoRDS; Leuschen and Allen, 2010; Shi and others, 2010), were retrieved from http://nsidc.org/icebridge/portal/ (green profiles in Fig. 1b). These data were used for cross-validation between the two systems and included in the computations.

As described below, the determination of a glacier-wide bedrock topography required the use of additional data besides measured ice thickness. These included a glacier-wide surface topography, and information about ice flow speed and the ice volume flux. A digital elevation model (DEM) of the surface of Starbuck Glacier was extracted from the dataset of Cook and others (2012). Their dataset covers the entire Antarctic Peninsula at 100 m resolution, and was derived from raw ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) global DEM (GDEM) data. Comparisons with elevation data collected by NASA's Ice, Cloud and land Elevation Satellite (ICESat) mission (Cook and others, 2012) revealed a root-mean-square error (RMSE) of ±25 m, which is about two to three times smaller than the RMSE for any of the products available from the SPIRIT (SPOT 5 stereoscopic survey of Polar Ice: Reference Images and Topographies; Korona and others, 2009) project.

Surface ice flow speed was computed from the data collected from a continuous GPS station installed on Starbuck Glacier (blue dot in Fig. 1b). The Leica 1200 station recorded positions at 10s intervals, starting from November 2011 and without interruption. The average flow speed was determined through linear regression of the time series collected until September 2012, and calculated to be 64.88ma-¹ (accuracy <0.01 ma-¹). Seasonal variations were not statistically significant.

3. Methods

3.1. RES data processing

The RES data were processed using ReflexW software. Best results in terms of bedrock-reflection visibility were obtained by applying the following processing sequence:

1. 1. Moving the start time: To ensure consistency in the start time of individual traces, start time was moved to the first detectable break in the radargrams.

2. 2. Interpolation to equidistant traces: Due to varying survey speed, the original data showed heterogeneous spacing between traces. This step ensured a constant spacing, chosen to be 0.5 m.

3. 3. Manual gain: In order to avoid excessive reverberation at large depths in the bandpass-filtering step, a manual gain reduction was applied with a linear decrease from 0 to -20 dB in the time interval [13.0, 14.4] us.

4. 4. Frequency bandpass filter: A first trapezoidal frequency bandpass filter was applied to dewow the data. Lower and upper cut-offs were chosen to be 0.75 and 8.00 MHz, respectively, and the passband was set between 1.50 and 4.00 MHz.

5. 5. Background removal: Persistent noise within individual radargrams was removed using a background-removal filter, i.e. by subtracting the average of all traces.

6. 6. Divergence compensation: In order to compensate for the geometrical divergence losses in signal amplitude with depth, a time-proportional divergence compensation gain was applied. The scaling value was set to 0.01.

7. 7. Resampling of the time increment: In order to reduce computational time in the subsequent migration process, the time increment for the individual traces was resampled from 4 to 10ns.

8. 8. Finite-difference (FD) migration: Starting from a velocity model that prescribes a uniform wave propagation speed of 168 m us–¹ (e.g. Vaughan and others, 2006), a FD migration was performed in order to contract diffractions to a minimum, and recover the correct location of the individual reflectors. The chosen migration scheme was implemented as the ‘steep-dip depth migration solution’ (also known as ‘15° assumption’; Claerbout, 1971), and the central frequency of the RES system (3 MHz) was set as the nominal frequency.

9. 9. Frequency bandpass filter: This second frequency bandpass filter was designed to remove high-frequency noise introduced by the migration step. Filter parameters were kept identical to step 4.

10. 10. Frequency-wavenumber (fk) filter: In many of the profiles, bedrock reflections towards the margins of the glacier were obscured by straight-line events crossing the records with uniform apparent dip. Partial removal of these events, originating from airwave reflections at the side-walls of the glacier, was achieved by applying a fk filter (e.g. Yilmaz and Doherty, 2001, section 6.2). Filter parameters were adjusted according to the apparent velocity of the individual noise events.

11. 11. Trace stacking: In order to increase the signal-to-noise ratio, additional trace stacking was performed as a last step. Three traces were combined in each stack.

Bedrock reflections were then manually picked, and signal travel time converted to ice thickness, assuming a uniform wave velocity of 168mμs-¹. Velocity changes in the firn layer were neglected, since the combination of relatively small accumulation rates (see Section 3.3) and fast firn densification processes, due to relatively high temperatures (e.g. Orr and others, 2008), suggests a relatively shallow and dense firn pack in the surveyed area. In cases where two or more reflector layers were visible in the same radargram, all reflectors potentially deriving from the bedrock were picked in the first instance. The correct layer was then identified by comparing profiles at crossover points (see Fig. 2a for an example). In cases in which the ambiguity could not be overcome, the data segment was discarded.

Fig. 2. (a) Example for two processed DELORES radargrams with a common crossover point, and picked horizons (dashed). (b) Processed DELORES radargram (picked horizon shown in light blue) for the profile coinciding with the IceBridge track. MCoRDS ice thicknesses converted to signal travel time are shown in pink. (c) Location of the profiles shown in (a) and (b).

3.2 Cross-validation with IceBridge MCoRDS data

As shown in Figure 1b one of the ground-based RES profiles coincided with the flight conducted by Operation IceBridge on 11 November 2011. This provided the opportunity to cross-validate the two datasets. For direct comparison, bedrock elevations were computed from the surface and ice-thickness data of both the DELORES and MCoRDS systems. The agreement is excellent (Fig. 2b): median absolute deviation between the two datasets was 17 m (2.8% of the local ice thickness), and the maximal discrepancy was 44 m (7.2%).

Additional comparisons were performed at the five crossover points between the DELORES and IceBridge tracks (Fig. 1b). In this case, median absolute deviation was 27 m (3.7% of the local ice thickness), and maximal deviation was 45 m (11.7%). Although the agreement is slightly less good in this case, the results indicate high accuracy for the MCoRDS data. This is somewhat surprising, since previous airborne surveys (e.g. Holt and others, 2006; Farinotti and others, 2013) highlighted the difficulty in correctly interpreting airborne RES data in regions with prominent topography. The MCoRDS data were included for the further analyses.

3.3. Determination of a glacier-wide bedrock topography

Despite the relatively high density of measurements collected for the main glacier trunk (Fig. 1b) the dataset is not sufficiently dense to determine a glacier-wide bedrock topography based on the RES measurements alone. Thus the ice-thickness estimation approach of Huss and Farinotti (2012) was used. This approach, which develops the methods presented by Farinotti and others (2009), was previously applied to estimate the ice-thickness distribution of every one of the -200 000 glaciers contained in the Randolph Glacier Inventory (Arendt and others, 2012). It is based on principles of ice flow dynamics, and estimates the bedrock topography using information derived at the glacier surface. This is done by disaggregating the glacier into 10 m elevation bands, calculating the mass turnover by estimating the glaciers’ mass-balance distribution, and converting the ice volume flux into ice thickness using Glen's flow law (Glen, 1955). The spatially distributed estimate is then obtained from the ice thickness calculated for individual elevation bands and the information of ‘zero ice thickness at the margins’, using an interpolation scheme that includes local weighting based on surface slope and distance from rock outcrops. The approach explicitly accounts for calving fluxes, and has a set of tunable parameters that can be used for matching direct measurements, if available. In the present case, the method was applied within the ice flow catchment of Starbuck Glacier, defined on the basis of the gradient of the available surface DEM (red dotted outline in Fig. 3), and the parameters of the model were constrained by the ice volume flux, estimated from the available ice flow velocity, and the collected ice-thickness measurements.

Fig. 3. Ice-thickness distribution (color map) and glacier bedrock elevation (contours) for Starbuck Glacier. The red dotted outline represents the ice flow catchment, defined from the gradient of the surface topography. Spatial reference is given analogously to Figure 1b.

In order to estimate the ice volume flux, the ice flow model of Sugiyama and others (2007) was used to compute a flow velocity field for the profile located closest to the continuous GPS station. As input, the model requires the surface and bedrock geometry for the considered cross section (obtained from the available surface DEM and the RES measurements, respectively), and the corresponding average surface slope in the flow direction (derived from the surface DEM, and averaged over a distance of ten times the local ice thickness; Kamb and Echelmeyer, 1986). The parameters of the model, which are the flow-rate factor in Glen's flow law, and a coefficient controlling basal sliding, were adjusted in order to reproduce the measured mean surface velocity. Since the contribution of basal motion to the actual surface velocity is difficult to constrain, two end members were considered for the estimate. The first assumes plug-flow, i.e. that basal velocity accounts for 100% of the observed surface velocity; the second assumes ‘creep only’, i.e. that basal velocity is zero everywhere, and that surface velocity is entirely due to internal ice deformation. According to the model and the two base assumptions, the ice volume flux across the test profile is between 96.4 x 10⁶and 56.7 x 10⁶ m³ ice a-¹ (0.087 and 0.051 Gta-¹, assuming an average density of 900 kg m-³). Assuming the additional ice volume flux yielded by small lateral tributaries downstream of the considered profile (Fig. 1b) is negligible, the estimated volume flux can be considered an estimate of the calving flux of Starbuck Glacier.

So far, no direct measurements exist for the surface mass balance of Starbuck Glacier. However, observations of previous field campaigns suggest that the net surface balance in the region may be close to zero at sea level (personal communication from M. Truffer, 2010, and authors’ own observations). As an approximation, the surface mass balance used as input for the model by Huss and Farinotti (2012) was prescribed by a linear function with altitude. The equilibrium-line altitude was set to sea level, and the accumulation rate at 2000m a.s.l. assumed to be 1.5 mw.e.a-¹ (Turner and others, 2002). This yields a surface mass-balance gradient of 7:5 x 10-⁴mw.e.a-¹ m-¹, and a net surface mass balance of 0.079 Gta-¹, which is in good agreement with the ice flux previously estimated from the surface flow speed.

Huss and Farinotti (2012) discussed two additional parameters used in estimating the ice-thickness distribution: the flow-rate factor, A, of Glen's flow law and a ‘valley shape factor’, f_s , controlling the general shape of the valley floor (U- or V-shaped), yielded by the spatial interpolation scheme. These parameters, which both only influence the calculated ice-thickness distribution, not the overall mass turnover of the glacier, were calibrated in order to minimize the mismatch between measured and modeled ice thicknesses at all points with RES. Calibration was performed by exploring a parameter space of plausible values ([2.4, 0.02] x10^-24 s^–1 Pa–³ for A and [0.1, 0.9] for f _s), and selecting the parameter combination yielding the minimal average absolute deviation. Best results were found for A = 0.32 x 10-²⁴ s^–1Pa-³ and f _s = 0.65. According to Cuffey and Paterson (2010), the value found for A corresponds to an ice temperature of about -10˚C.

Even with the so-calibrated parameters, the ice-thickness distribution for the measured cross sections cannot be reproduced exactly. This is not surprising considering the simplifications within the model itself. Thus, as a last step, a local a posteriori correction, corresponding to the local difference between measured and modeled ice thickness, was computed for every location with RES measurements. Prior to application, the correction was spatially interpolated using a bilinear interpolation scheme. This ensured the correction also had an influence at locations with no RES measurements, and smoothness of the solution. The accuracy of the calculated bedrock topography, including the effect of the a posteriori correction, is assessed in Section 5.

4. Results

The ice-thickness distribution derived for Starbuck Glacier and the resulting bedrock topography is shown in Figure 3. The total ice volume for the considered surface area of 258 km² is calculated as 80.7 km³, corresponding to an average ice thickness of 312 m. The most prominent feature in the bedrock is the reverse-sloped, relatively wide valley floor that extends upstream from a riegel occurring ~5 km from the grounding line. Ice thicknesses up to ~1020m are found in the deepest parts of the glacier, yielding bedrock elevations as far down as 500 m below sea level. Relatively shallow ice thicknesses are calculated for higher elevations, although the uncertainties are maximal for these regions (see next section). For bedrock elevations >1000m a.s.l., ice thickness rarely exceeds 500 m.

The largest ice thicknesses might be surprising considering the relatively small size of the glacier. However, these values are not dissimilar to those found for nearby Flask Glacier (see Fig. 1a for location), where ice thicknesses of up to 1800 m have been inferred previously (Farinotti and others, 2013). Moreover, it should be noted that the section with the largest ice thicknesses is relatively well covered by RES measurements (Fig. 1b) increasing confidence in the results.

5. Accuracy Estimates

5.1. Accuracy of RES measurements

For RES systems the maximum theoretical resolution, in the sense of the ability to distinguish between two closely spaced reflectors, is a quarter of the pulse wavelength (e.g. Reynolds, 2011). For the central frequency of the employed system, this is ~15m. However, when measuring the distance to a single reflecting interface (e.g. the glacier bedrock), the relevant resolution is given by the so-called range resolution, which is determined by the rise time of the source signal, the bandwidth of the system and the digitization interval, and which can be significantly smaller than the maximum theoretical resolution.

In practical applications, several additional factors limit the accuracy. These include (1) uncertainty in the velocity used for converting signal travel time to depth, (2) the different wave propagation speeds in firn and ice, (3) potential signal distortions introduced through the different filtering steps used during data post-processing, (4) uncertainty in the exact position of the transmitter and receiver during data collection, (5) variation in the geometrical settings of radar antennas when traveling, (6) the assumption that signals in the radargrams originate from reflectors located in the plane defined by the travel line at the glacier surface and a normal vector to that surface, (7) uncertainty in the interpretation of the post-processed radargrams and (8) the accuracy in picking individual reflectors. While points (5)-(8) are likely to result in randomly distributed errors, points (1)-(3) could potentially lead to systematic deviations.

Formally deriving the resulting overall accuracy is not straightforward. However, the presence of a significant bias can be ruled out, in light of the excellent agreement between DELORES and MCoRDS data, whereas the 29 crossover points in the RES profiles for which bedrock reflections could be successfully detected (Fig. 1b) allow a bulk estimate for the stochastic component of the total accuracy. As before, this was done by computing the deviation between the ice thicknesses determined from two different profiles at each crossover point. The median absolute deviation is 11 m, and all deviations are contained within 45 m. Expressed relative to the corresponding local ice thickness, the deviations are 1.4% and 6.1%, respectively. It is remarkable that these numbers are not significantly different from that assessed for the crossover points with the MCoRDS data, which, again, points to the accuracy of the MCoRDS data.

5.2. Accuracy of the glacier-wide bedrock topography

The accuracies assessed above apply only to locations at which RES measurements are available. The accuracy for locations in which the bedrock topography was determined using the approach of Huss and Farinotti (2012) is expected to be lower. In general, the accuracy is expected to decrease with increasing distance from RES profiles. In order to provide a spatially distributed estimate of this accuracy, a resampling experiment was performed.

In a first step, only one profile was used for calibrating parameters A and fs in the approach of Huss and Farinotti (2012), and a glacier-wide bedrock topography was calculated according to the procedure described above. The ice thickness calculated at the locations of the RES profiles not used for the calibration was then compared on a point-to-point basis to the measurements, and the resulting deviation expressed as a function of the distance to the next RES point actually used in the calibration. This step was repeated for every one of the 39 RES profiles available.

In a second step, the same procedure was repeated using two RES profiles for parameter calibration, and the remaining 37 profiles for assessing the deviation. As a third step, three profiles were used for calibration and in subsequent steps more profiles were used until 38 profiles were used for calibration and 1 for validation.

For any step, k, in which k out of the n = 39 RES profiles are used for calibration, there exist ( n!)/K!(n - k)! potential combinations for choosing a given sample of profiles. For all 38 steps, this results in a total of ~ 5.5x10¹¹ possible combinations, which is unaffordable in terms of computational cost. In order to overcome this problem, k = 39 random trials were performed at each step, resulting in a total of 39 x 38 = 1482 experiments, from which the accuracy was estimated. The value k = 39 was chosen since (1) it is the maximal number allowing the repetition of the same number of experiments at each step without repeating a given combination and (2) repeating the same number of trials at each step prevents the estimate being dominated by steps with a particular number of profiles.

Figure 4a shows the expected relative ice-thickness deviation resulting from this procedure, expressed as a function of the distance to the next RES point used for calibration. The resulting spatially distributed accuracy estimate is shown in Figure 4b. According to the results, uncertainty in total ice volume and average ice thickness is ±7:2 km³ and ±30 m, respectively (numbers refer to the 95% confidence level).

Fig. 4. (a) Estimated uncertainty in relative ice thickness as a function of the distance to the closest RES point. Median (solid curve) and empirical 9 5% confidence intervals (dotted curves) are shown. (b) Spatially distributed accuracy estimate (95% confidence level), derived from the function shown in (a) and the ice-thickness distribution shown in Figure 3 .

6. Conclusions

A total of 90 km of ground-based RES measurements were collected and successfully analyzed for Starbuck Glacier. Repeated measurements at crossover tracks proved to be valuable for both data interpretation and estimating the accuracy. Cross-validation with data collected by NASA's MCoRDS sensor showed excellent agreement, suggesting a high accuracy for the ice-thickness measurements carried out by Operation IceBridge, even for relatively narrow valley glaciers. The various RES profiles were used to calibrate a method to calculate ice thickness from surface topography (Huss and Farinotti, 2012) in order to derive a glacier-wide ice-thickness distribution and corresponding bedrock topography. Additional constraints were obtained from in situ ice flow-speed measurements, and were used to estimate the ice volume flux near the grounding line. A reverse-sloped glacier bed was detected upstream of a bedrock riegel occurring at a distance of ~5 km of the current grounding line. Ice thicknesses up to ~1020m were found, corresponding to bedrock elevation as low as 500 m below sea level. The total volume of the glacier was estimated to be 80:7 ± 7:2 km³ of ice, corresponding to an average ice thickness of 312 ± 30 m. This contribution is a further step towards an improved knowledge of the subsurface topography of the glaciers of the Antarctic Peninsula, which is a prerequisite for investigating the interactions between outlet glaciers and ice shelves.

The RES data, the surface and bedrock topography, the estimated ice-thickness distribution and the estimated accuracy map are accessible at http://dx.doi.org/p82