1. Introduction
Thwaites Glacier in the Amundsen Sea sector of West Antarctica has an area of -182 000 km2, and much of it is grounded well below sea level with an inward-sloping bed (Reference HoltHolt and others, 2006). The catchment consists of tributaries that merge into the main trunk (Fig. 1). The glacier is currently losing mass at a rate of - 20 Gt a–1 (Reference RignotRignot and others, 2008); thinning is 4 ma–1 in fast-moving regions near the grounding line, decreasing to 0.1 ma–1 in slower-moving (<100m a-1) regions farther inland (Reference Pritchard, Ligtenberg, Fricker, Vaughan and Van den BroekePritchard and others, 2012). Changes over the past three decades include acceleration (Reference Rignot, Vaughan, Schmeltz, Dupont and DupontRignot and others, 2002), increased thinning (Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009; Reference ShepherdShepherd and others, 2012), loss of ice-shelf buttressing (Reference MacGregor, Catania, Markowski and AndrewsMacGregor and others, 2012) and increased ocean warming (Reference Pritchard, Ligtenberg, Fricker, Vaughan and Van den BroekePritchard and others, 2012). Accumulation in the catchment has changed little over past decades (Reference MedleyMedley and others, 2013).
InSAR velocity map (Reference JoughinJoughin and others, 2009) overlain by radar lines (Reference HoltHolt and others, 2006). Green lines are radar lines collected in the AGASEA (Airborne Geophysical Survey of the Amundsen Sea Embayment, Antarctica) campaign. Purple and orange lines indicate areas where radar-detected features were tracked; orange lines are shown in Figure 2 . The numbers indicate radar lines L1–L8. Yellow circles are shallow core locations from the United States International Trans-Antarctic Scientific Expedition. Red circle is the West Antarctic Ice Sheet Divide ice-core site. A–E are GPS locations. Bold black line is approximate grounding line position. Inset map shows study location.
Over the past two decades, Thwaites Glacier has not accelerated as dramatically as Pine Island Glacier, but has increased its rate of mass loss due to widening of the fast-flow region, particularly on the eastern margin (Reference RignotRignot, 2008). The eastern margin is not well constrained by topography or basal properties and may be susceptible to migration (Reference MacGregorMacGregor and others, 2013). Inferring changes in flow pattern is difficult because observations in the upper regions of Thwaites Glacier span only a few decades; other indicators used for inferring past flow directions, such as flow stripes (Reference Fahnestock, Scambos and BindschadlerFahnestock and others, 2000) or crevassing at shear margins (Reference Shabtaie, Bentley, Bindschadler and MacAyealShabtaie and others, 1988), are not visible on Thwaites Glacier. Radar-detected internal layers can be used to extract information about past ice flow over century timescales; here we adapt methods developed by Reference Ng and NgNg and Conway (2004). The radar-detected layers in fast-flowing ice often show trough-and-crest sequences that are inherited from flow disturbances farther upstream. Although the origin of the layer patterns is often not known, the patterns can be used to define the geometry of past flowbands by tracking sequences between radar profiles perpendicular to the ice flow.
Quantifying the amount of thinning in the interior of Thwaites Glacier helps constrain the rate at which perturbations at the grounding line propagate inland (Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004; Reference Joughin, Smith and HollandJoughin and others, 2010a). For the fast-flowing region approximately 100-200 km inland of the grounding line, satellite radar altimetry indicates thinning of -0.10 ma–1 for 1992-2003 with an uncertainty of -0.1 ma–1 (e.g. Reference ZwallyZwally and others, 2005; Reference HelsenHelsen and others, 2008). Satellite laser altimetry suggests slightly greater thinning of -0.15 ± 0.07 m a–1 (Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009) for 2003-09. Altimetry measurements are sensitive to inter-annual variations in accumulation rate and changes in firn compaction. Flux analysis of flowbands is a potential complementary method for determining inland thinning rates because it is not sensitive to interannual accumulation variability. Shabatie and others (1988) used flux analyses to assess the mass balance of Mercer, Whillans and Kamb ice streams on the Siple Coast. They found a similar spatial pattern of thinning and thickening to that found by flux analysis using satellite measurements of velocity (Reference Joughin, Tulaczyk, Bindschadler and PriceJoughin and others, 2002; Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009).
2. Data Sources and Methods
2.1. Radar-detected ice thickness and internal layers
Airborne geophysical surveys of the Amundsen Sea Embayment by the University of Texas Institute for Geophysics (UTIG) and the British Antarctic Survey provide measurements of ice thickness (Reference HoltHolt and others, 2006; Reference VaughanVaughan and others, 2006; Reference ChenBlankenship and others, 2012) and internal layering. We use radar data collected by UTIG with a 60 MHz, phase-coherent High-Capability Radar Sounder (HiCARS) system (Reference Peters, Blankenship and MorsePeters and others, 2005) that were range-migrated using along-track synthetic aperture radar (SAR) focusing (Reference Peters, Blankenship, Carter, Kempf, Young and HoltPeters and others, 2007). Ice thickness was mapped along radar profiles that were collected on 15.6 km grids over Thwaites Glacier (Fig. 1). Radar-detected internal layers were picked and mapped using a seismic package (Geoframe) with strong crossover control. The conversion of radar travel times to depth assumed a constant wave velocity of 168.374 µs m–1 and did not incorporate a firn layer (Reference HoltHolt and others, 2006).
2.2. Radar-detected flowlines
We mapped past flowlines through the study region by tracking positions of recognizable features (crest and trough sequences) between successive pairs of radar cross-profiles (Fig. 2). Positions of recognizable features were picked independently; that is, when picking features between cross-profiles L2 and L3, the positions of picks between cross-profiles L1 and L2 were not considered. A total of 191 flowlines were tracked (Fig. 3; Reference Conway, Catania and FudgeConway and others, 2010); most were located in the eastern and middle tributaries. Fewer flowlines could be tracked in the western tributaries because the radar profiles often crossed the flow direction at angles greater than 308. Four flowlines in the middle tributary could be tracked over 110 km through eight cross-profiles.
Two radio-echograms with traced internal layers in different colors. Flowlines from matching distinct layer patterns between radio-echograms are shown by vertical red lines. Flux-gate areas of the outer flowlines used in the flux balance calculations are shown in dashed purple.
Radar-detected flowlines shown with black circles connected by thin black lines. Thick black lines mark flowlines tracked for 110 km, with the outside flowlines defining the flowband study area. Yellow letters are GPS station locations. White arrows are GPS velocities. Red arrows are InSAR velocities at GPS locations. Large black numbers are average radar-detected directions minus InSAR directions (8). Small white numbers are GPS directions minus InSAR directions. Small red numbers are GPS directions minus radar-detected directions. Background color is bedrock elevation.
Uncertainty in the flowlines was estimated from the differences between picks of a feature at each radar line. Features were on average identified within 140 m, with a maximum difference of 380 m. This results in a ~3° uncertainty in the direction of a flowline. The radar-detected layers provide information on timescales related to the advection of ice between upstream and downstream cross-profiles (Ng and Conway 2004). The velocities range from 75 to 400 m a-1 in the study area, so the time for flow through the study area ranges from -50 years for a feature identified in one pair of radar profiles (~15.6 km) to -600 years for features tracked through eight radar profiles (-110 km).
2.3. Satellite-derived velocity field
We use a map of 1996 surface velocities derived using interferometric SAR (InSAR) and speckle-tracking methods (Reference JoughinJoughin and others, 2009). In areas of crossing satellite orbits where the phase could be successfully unwrapped, formal uncertainty in the horizontal velocity components is 3-5 ma–1. In addition to the formal uncertainties, an additional - 3 % uncertainty is introduced in the process of removing vertical displacements (Reference JoughinJoughin and others, 2009). A third uncertainty component comes from systematic errors in the SAR baseline estimate used to convert image displacements and phase differences into velocity values. While every effort has been made to minimize these errors, for velocity estimates derived from image pairs with short temporal intervals, these errors can be on the order of 50m a-1, and are correlated over several tens to hundreds of kilometers (Reference Joughin, Smith, Howat, Scambos and ScambosJoughin and others, 2010b). The uncertainties of the InSAR velocities presented in Section 3 are the sum of the formal uncertainty and 3% uncertainty from removing vertical displacements.
An InSAR velocity map is also available from Reference Rignot, Mouginot and MouginotRignot and others (2011). The map uses the same 1996 data as Joughin and others (2009) but also incorporates data from 2007 and 2008 that cover portions of Thwaites Glacier. We primarily use the map of Reference JoughinJoughin and others (2009) because it was developed for detailed investigations of Thwaites Glacier.
We discuss differences between Reference Rignot, Mouginot and MouginotRignot and others' (2011) and Reference JoughinJoughin and others' (2009) velocities in Section 3.3.
2.4. GPS measurements
Dual-frequency GPS units (Trimble 5700, R7, NetRS receivers with Trimble Zephyr geodetic antennas) were installed for - 4 week durations during the 2007/08 Antarctic summer along the two main trunks of Thwaites Glacier and on the surrounding slow-moving, flanking ice. GPS antennas were mounted on rigid metal poles that extended -1-2 m across the ice surface. Our processing strategy generally follows that outlined by Reference KingKing (2004) for glaciological applications. GPS stations on streaming ice were processed kinematically relative to a local base station using differential carrier-phase positioning, with epoch-by-epoch zenith tropospheric delay estimation as implemented in the Track software (Chen, 1998) with very loose constraints on rover site motion. Daily solutions for the local base station on slow-moving ice near Mount Takahe were calculated using differential carrier-phase positioning as implemented in GAMIT/GLOBK. Antenna heights were updated daily assuming a linear trend between measured antenna heights at the beginning and end of the field season. Horizontal velocities were extracted from linear fits for the entire field seasons using iteratively reweighted least squares with bi-square weighting after data were transformed into a polar stereographic projection (central meridian is 08 and latitude of true scale is -71° S). Although formal uncertainties (defined as two standard deviations from the best fit) were very small (~0.001 m a–1), we conservatively estimate horizontal velocity uncertainties at ~0.2 ma–1 based on values of residuals from the best-fit line to account for errors associated with antenna height uncertainty, firn compaction and advection.
2.5. Accumulation rate
Measurements of accumulation in the interior of Thwaites Glacier are sparse. Reference KaspariKaspari and others (2004) report measurements from eight firn/ice cores in central West Antarctica. The one core in the Thwaites Glacier catchment indicated average accumulation of 0.48m a-1 (reported accumulation rates are in ice equivalent) over the period 1922-91. Several studies have used these data to help constrain regional climate models (e.g. Reference Arthern, Winebrenner and VaughanArthern and others, 2006; Reference Monaghan, Bromwich and WangMonaghan and others 2006; Reference Van de Berg, Van den Broeke, Reijmer and Van MeijgaardVan de Berg and others, 2006); estimates of accumulation in the Thwaites catchment range from 0.4 to 0.65 m a–1. Recent work using a combination of radar-detected shallow layers and firn cores found an average accumulation rate of 0.46m a-1 for the entire Thwaites catchment (Reference MedleyMedley and others, 2013). For our study area, the accumulation rate for each segment was calculated from Reference MedleyMedley and others' (2013) accumulation map and had an average of 0.54m a-1. The study area was well covered by radar lines, yielding an uncertainty in the accumulation rate for each segment of 0.05 m a–1.
3. Results
3.1. Ice-flow direction
Ice-flow directions defined by radar-detected features and by InSAR differed by an average of 2.58 in the middle tributary but substantially more in the other tributaries. Figure 3 shows the average direction difference for each tributary. The InSAR directions are rotated eastward (towards the top) by an average of 8.38 in the eastern tributary relative to the radar-defined directions. The disagreement is even larger for the two western tributaries although the number of measurements of radar-defined flowlines is smaller. The InSAR directions are rotated westward by an average of 12.98 in the southwestern tributary and 22.48 in the western tributary.
Differences between the ice-flow direction inferred by tracking internal radar features and the InSAR directions could indicate a change in flow pattern over the past few centuries because the radar-detected directions indicate directions for the past few centuries while the InSAR directions are a modern snapshot. However, the differences might also be due to inaccuracies of one of the methods. To examine potential inaccuracies, we compare the flow directions to those measured at five GPS stations. Stations A–C are located in the eastern tributary, and D and E are in the middle tributary. GPS units were not deployed in the southwestern and western tributaries, where differences are largest. The radar-detected and GPS directions differed by up to 3.48, which is within the combined uncertainty; the differences are not biased in one direction (Table 1). The maximum difference between InSAR and GPS directions is 12.78. The agreement is good in the middle tributary, but consistently rotated eastward in the eastern tributary.
Comparison of ice-flow directions
The agreement between radar-detected and GPS directions in the eastern tributary suggests that the InSAR directions are incorrect. The large differences between the radar-detected and InSAR directions in the western and southwestern tributaries are also likely caused by incorrect InSAR directions. Close inspection shows that radar-detected directions follow the bedrock trough and the direction of increasing ice velocity magnitude. In contrast, InSAR velocity vectors are oriented such that an ice particle starts in the center of a tributary and travels across the margin into slower-moving ice.
3.2. Flux imbalance and uncertainty
Flowbands can be tracked between six to eight radar lines in both the eastern and middle tributaries. This is a distance of 80–110 km and presents the opportunity to determine thinning rates from mass conservation. The large differences between the ice-flow directions determined by InSAR and those determined by radar and by GPS for the eastern tributary prevent a reliable flowband analysis. However, the agreement of all three datasets in the middle tributary suggests the InSAR velocities may be sufficiently accurate to allow the thinning rate to be determined.
Four flowlines were tracked through at least eight radar lines in the middle tributary. We focus on the large flowband defined by the outer two flowlines (Fig. 3) and calculate the flux gate area at the top and bottom of each flowband segment. The flux gate area is the cross-sectional area between two flowlines (Fig. 2). Assuming no ice flow across the sides of a flowband, from continuity the flux through the upstream gate equals the flux out the downstream gate plus any mass changes within the flowband:
where A is the flux gate area, u is the average velocity at the flux gate, dn indicates the downstream gate, up indicates the upstream gate, b is the ice-equivalent accumulation rate, h is the rate of change in ice thickness (positive for thickening), m is the rate of basal melting and SA is the surface area of the flowband. The basal melt rate in the study was calculated by Reference JoughinJoughin and others (2009), and the entire study region is melting at the bed. Most of the study area has a near-zero melt rate, with some radar lines averaging a few centimeters of melt per year. The basal melt rate is always less than the accumulation uncertainty, so we do not consider it in our calculations.
The flowband width does not accumulate uncertainty because the radar-detected features are identified independently at each radar profile. Uncertainty in the flux gate width is then the uncertainty in defining the features at each end, which we estimate to be a total of 600 m in width. The radar lines are not always perpendicular to the ice flow, so we convert the measured width to the width perpendicular to the ice flow by multiplying by the cosine of the difference between the direction of the radar line and the direction of the ice flow. Up to 11% of bed returns for a flux gate cannot be distinguished; the maximum distance between bed returns is 1 km. Depths vary by up to 50 m over 1 km distances, but typically variations are less. We interpolate linearly between measured bed returns and estimate that the lack of returns results in <5 m depth uncertainty. Uncertainty of the radar wave speed in ice largely cancels at the upstream and downstream gates, so we ignore this contribution.
The velocity at the flux gate is determined from the InSAR velocity map. The surface velocity is considered equal to the depth-averaged velocity, i.e. the ice moves entirely by sliding. Using a steady-state temperature profile with a surface temperature of -208°C, an accumulation rate of 0.5 ma–1 and a basal temperature at the melting point, the contribution of internal deformation to the surface velocity is ~5m a-1 in the study area. The depth-averaged velocity is -0.95 of the surface velocity. This results in a possible overestimation of the depth-averaged velocity of <1 ma–1, which is much smaller than the InSAR uncertainty and is not considered here. Such an overestimation would have little impact on the flux balance calculations because each flux gate would be biased in the same direction.
The net flux, average depth, width and velocity at each gate are shown in Table 2, as well as the associated uncertainties. The total flux uncertainty is calculated assuming standard propagation of independent errors. The uncertainty for each flux gate is 5-10% of the total flux and is dominated by the velocity.
Values at flux gates
The net flux imbalance as well as the equivalent change in ice thickness averaged over the surface area is shown in Table 3 for each segment. A positive value indicates more mass is entering the segment, which results in thickening; a negative value indicates thinning. The uncertainty for any segment is often an order of magnitude larger than the calculated thickness change. This arises for short segments (15.6 km) because the flux difference is small and the velocity uncertainties at both flux gates dominate the estimate.
Flowband ice thickness change and uncertainty
The full 110 km flowband allows a more robust estimate of average thinning and is also shown in Table 2. The average thinning rate is 0.47±0.42 ma–1. L1 has a large velocity uncertainty, so calculating the thinning from L2-L8 yields nearly the same thinning rate of 0.49 ma–1 but reduces the uncertainty to 0.34 ma–1.
3.3. Additional velocity uncertainty
Uncertainty in the flux analysis is dominated by the velocity uncertainty. Thus far, we have assumed that the formal and slope-dependent uncertainties sufficiently characterized the velocity uncertainty, as in the Ross Drainage where the InSAR velocities were compared with ground-based measurements Reference Joughin, Tulaczyk, Bindschadler and Price(Joughin and others, 2002). However, Thwaites Glacier is a difficult area for InSAR measurements. High accumulation rates and low radar reflectivity lead to poor coherence and rapid temporal decorrelation between images. As a result, in many parts of the glacier, only a few short-repeat image pairs from the European Remote-sensing Satellite (ERS) are available. These estimates have larger errors and biases than those based on the 23 or 38 day pairs available elsewhere.
The InSAR and GPS velocities are compared in Table 4. In the eastern tributary, where the flow directions do not agree, the InSAR velocities exceed the GPS velocities by 16% at site A and 1 3% at site B, which is more than the formal and slope-dependent uncertainty and indicates that there is additional uncertainty from the short repeat interval. In the middle tributary, the InSAR velocity is 5% greater at site D and slightly outside of the uncertainty; at site E, they agree within the uncertainty. The comparison suggests that uncertainty in the velocity magnitude is reasonable in places where the direction is accurate, and inaccurate directions are an indication that the magnitude may be off as well.
Velocity comparison
At four of the five GPS locations, the InSAR velocity is greater than the GPS velocity. It is unlikely that the difference is caused by the different times at which the velocities were measured; InSAR velocities are from 1996 and GPS velocities are from austral summer 2007/08. Thwaites Glacier is most likely accelerating in response to recent terminus changes, so the GPS velocities would be expected to be greater than the InSAR velocities.
We also compared the GPS velocities with the velocity map of Reference Rignot, Mouginot and MouginotRignot and others (2011). Rignot and others use a combination of the 1996 ERS-1 and -2 satellite passes (from which Reference JoughinJoughin and others' (2009) velocities were calculated) as well as the 2007/08 Phased Array-type L-band SAR (PALSAR) data, which do not cover the entire region. InSAR velocity directions from Reference Rignot, Mouginot and MouginotRignot and others (2011) agree well with those calculated using GPS in the eastern tributary (difference is <2.58; Table 1); however, the directions differ by 9.58 and 5.38 in the middle tributary. This is opposite to those from Reference JoughinJoughin and others (2009), which agree well in the middle tributary but not the eastern tributary. Gaps in Reference Rignot, Mouginot and MouginotRignot and others' (2011) velocity map prevent comparisons with the radar-detected directions in the western and southwestern tributaries.
Reference Rignot, Mouginot and MouginotRignot and others' (2011) velocities are less than the GPS velocities at four of the five GPS locations, and the differences exceed the stated velocity uncertainties at three of the five sites. Velocity differences often exceed the stated uncertainties by more than an order of magnitude. At site D, Reference Rignot, Mouginot and MouginotRignot and others (2011) give a 1 m a – 1 velocity uncertainty, yet the GPS velocity is 28m a-1 (12%) faster; at site E, the GPS velocity is 15m a–1 (14%) faster, although the stated uncertainty is 1 m a–1. Agreement between Rignot and others' (2011) velocities and GPS velocities is better in the eastern tributary; the difference exceeds the uncertainty at only one of the three sites. This suggests that the velocity magnitudes are more accurate where the directions agree, as observed with Reference JoughinJoughin and others' (2009) velocities.
The large differences between the GPS and both InSAR velocities (and even larger differences between the two InSAR maps) suggest the uncertainties that accompany the InSAR data are too small and need to include uncertainty associated with short repeat intervals. Hence the estimated uncertainty of the flux analysis presented above is likely higher than stated, which reduces the confidence in the thinning estimates.
GPS velocities provide the necessary precision but lack spatial coverage. If the velocities at the L2 and L8 gates (from Reference JoughinJoughin and others, 2009) are scaled to the GPS velocities, the inferred thinning is reduced to 0.36 m a–1, in closer agreement with the altimetry methods. The uncertainty is difficult to estimate because while the GPS velocities have an uncertainty of 0.2 m a–1, a constant scaling of all velocities along the flux gate is likely too simplistic. If only the GPS velocity uncertainty is used, the total uncertainty is 0.08m a–1, which is similar to the precision of the laser altimetry techniques; the uncertainty in accumulation then dominates the overall uncertainty estimate. Multiple GPS measurements at the flux gates would be necessary to reduce the velocity uncertainty.
4. Conclusions
Flowlines can be effectively determined using radar-detected features in the interior of Thwaites Glacier where other means of identifying flow directions (e.g. flow stripes and crevasses at shear margins) are not possible. Comparisons of GPS and radar-detected flow directions with InSAR directions indicate that InSAR velocities alone are not sufficient to define flowlines and flowbands. The agreement of the GPS and radar-detected flowlines also reveals that flow directions of Thwaites Glacier tributaries have not changed significantly over the past few hundred years.
A flux analysis in the middle tributary where the InSAR velocity directions agree well with the GPS and radar-detected directions indicates thinning of 0.49 m a–1, but with a large uncertainty of 0.34m a-1. Further, InSAR velocities often differ from GPS velocities by more than InSAR uncertainties, indicating that the flux analysis uncertainty is too low. The large uncertainty of the flux balance analysis prevents inference of thinning estimates with similar precision to satellite altimetry, which indicates thinning of 0.10-0.15 ma–1 with uncertainties of -0.1 ma–1 (Reference ZwallyZwally and others, 2005; Reference HelsenHelsen and others, 2008; Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009). InSAR velocities provide a good overview of the velocity structure of Thwaites Glacier and are sufficient for many glaciological investigations (Reference JoughinJoughin and others, 2009; Reference Rignot, Mouginot and MouginotRignot and others, 2011), but they are not sufficiently accurate for flux analyses in upstream tributaries.
Acknowledgements
We thank Brooke Medley and Joe MacGregor for many discussions about Thwaites Glacier. We also thank the US National Science Foundation for supporting this research with grants ANT-0739372 to H.C. and 0632198 to S.A. We also appreciate the efforts of editors Dorthe Dahl-Jensen and David Braaten.
