Introduction
Thwaites Glacier drains part of the Amundsen Sea Embayment of the West Antarctic ice sheet. It is one of the largest and fastest-flowing glaciers in Antarctica and is presently losing mass at a rate that is surpassed only by the adjacent Pine Island Glacier (Reference ShepherdShepherd and others, 2012). Recent changes on Thwaites Glacier include loss of ice-shelf buttressing (Reference MacGregor, Catania, Markowski and AndrewsMacGregor and others, 2012a; Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012), rapid thinning (e.g. Reference Pritchard, Arthern, Vaughan and EdwardsPritchard and others, 2009) and acceleration across its grounding line (e.g. Reference RignotRignot, 2008). Much of Thwaites Glacier is grounded well below sea level and its bed generally slopes downward heading inland (Reference HoltHolt and others, 2006), so it is potentially dynamically unstable if its grounding line retreats significantly (e.g. Reference SchoofSchoof, 2007; Reference Gudmundsson, Krug, Durand, Favier and GagliardiniGudmundsson and others, 2012; Reference ParizekParizek and others, 2013). Understanding the dynamics of Thwaites Glacier and other overdeepened outlet glaciers in the Amundsen Sea Embayment is thus important for improving predictions of the contribution of the West Antarctic ice sheet to future sea-level rise (e.g. Reference Joughin and AlleyJoughin and Alley, 2011).
Although reported acceleration across Thwaites Glacier is currently confined to within ~40 km of the grounding line (Reference RignotRignot, 2008), observations of thinning much farther upstream (Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012) raise questions about the potential vulnerability of its upstream ice-flow configuration to future changes (Reference ParizekParizek and others, 2013). Six distinct tributaries with speeds of <100 m a−1 in the glacier’s drainage basin coalesce into its main trunk, which is up to 200 km wide. The shear margins that delineate the onset regions of these tributaries extend up to 300 km upstream of the grounding line (Fig. 1). Surface speeds within the main trunk are non-uniform; in a central region up to ~40 km from the grounding line, surface speeds exceed 1 km a−1. Slower ice flow (100–500m a−1) outside this central region results in intraglacier zones with large lateral shear strain rates (up to 0.1 a−1; Reference Lang, Rabus and DechLang and others, 2004; Fig. 1).
Fig. 1. (a) 1 km Bedmap2 bed elevation grid across Thwaites and Pine Island Glaciers (Reference FretwellFretwell and others, 2013; World Geodetic System 1984 (WGS84) ellipsoid). Contour interval 100 m. Manually traced shear margins, based on 1996 lateral shear strain rates shown in (b), are outlined by thick solid black lines. The grounding line is derived from the 2008–09 Moderate Resolution Imaging Spectroradiometer (MODIS) Mosaic of Antarctica (personal communication from T.A. Scambos and T. Haran, 2013). (b) Lateral shear strain rate (
; Eqn (A5) in Appendix). Positive (negative) values indicate right-handed (left-handed) shear, i.e. the along-flow surface velocity is decreasing to the right (left) of the flow direction. Inset map shows the location of the study area in Antarctica.
Bed topography and/or spatial variation in resistance to basal sliding largely determine the shear margin positions of many fast-moving glaciers and ice streams (Reference Raymond, Echelmeyer, Whillans, Doake, Alley and BindschadlerRaymond and others, 2001). Resistance to basal sliding depends on several different subglacial properties, including basal hydrology, lithology and roughness (Reference Winsborrow, Clark and StokesWinsborrow and others, 2010). Where abrupt changes in basal lubrication occur, shear margins may be inherently unstable due to thermomechanical feedbacks (Reference RaymondRaymond, 1996; Reference Jacobson and RaymondJacobson and Raymond, 1998). Subglacial geology (e.g. the presence of deformable till) also influences the position of ice-stream margins (Reference Anandakrishnan, Blankenship, Alley and StoffaAnandakrishnan and others, 1998; Reference BellBell and others, 1998; Reference BlankenshipBlankenship and others, 2001; Reference StudingerStudinger and others, 2001). Most of Thwaites Glacier’s tributaries clearly originate within bedrock troughs (Reference HoltHolt and others, 2006; Fig. 1), a pattern similar to that of the ice-stream tributaries in the Ross Sea Embayment (Reference JoughinJoughin and others, 1999). The best example of this topographic control is the main trunk’s western shear margin, which contours around several volcanoes/nunataks (including Mount Takahe and Mount Murphy) (Fig. 1). The role of bed topography is less clear along the eastern shear margin, as observed by Reference HoltHolt and others (2006) (Figs 1 and 2). At its upstream end, along the northern flank of the Byrd Subglacial Basin, the eastern shear margin roughly follows contours of bed elevation, but it is nearly perpendicular to those contours farther downstream. This pattern is not observed for any other of Thwaites Glacier’s tributary shear margins.
Fig. 2. Bedmap2 bed elevation in the vicinity of the eastern shear margin (Reference FretwellFretwell and others, 2013). The format is repeated from Figure 1a. Alternating colored circles and triangles mark the transect– margin intersections, which are used to distinguish these intersections in later figures. Table 1 lists the names of the Airborne Geophysical Survey of the Amundsen Sea Embayment, Antarctica (AGASEA) transects. Transects with low incidence angles (<20°) relative to the shear margin are not included in this study. Also shown is the same shear margin traced using derived from the 2007–08 surface velocity field of Reference Rignot, Mouginot and ScheuchlRignot and others (2011a), which is not substantially different from that derived using 1996 surface velocities.
To better understand the nature of the controls upon Thwaites Glacier’s eastern shear margin, in this study we investigate bed properties there using ice-penetrating radar data. Specifically, we use focused ice-penetrating radar transects to examine the topography and radar reflectivity of the bed. This approach is motivated by the variability in bed conditions across ice-stream shear margins reported by analogous studies of the Ross Sea Embayment. Downstream of their respective tributaries, the shear margins of the Ross ice streams are likely controlled by abrupt transitions in basal resistance to sliding (Reference Raymond, Echelmeyer, Whillans, Doake, Alley and BindschadlerRaymond and others, 2001; Reference Catania, Conway, Gades, Raymond and EngelhardtCatania and others, 2003), although the nature of this putative transition in the basal condition is poorly understood. For example, Reference Peters, Blankenship and MorsePeters and others (2005) found large (>10 dB) and abrupt (<5 km) decreases in bed reflectivity across portions of the downstream shear margins of Whillans and Kamb Ice Streams, which suggest that the bed there transitions abruptly from thawed beneath ice streams to frozen outside the margins. However, Reference Raymond, Catania, Nereson and Van der VeenRaymond and others (2006) found no significant change in bed reflectivity across the same shear margin farther upstream on Whillans Ice Stream.
Table 1. AGASEA transects used in this study
Data
Surface velocity
We use surface velocity fields across Thwaites Glacier to locate the position of the eastern shear margin and to calculate lateral shear strain heating. Several glacier-wide maps of the 1996 surface velocity of Thwaites Glacier and adjacent glaciers have been produced using similar combinations of satellite-based synthetic aperture radar (SAR) interferometry and speckle tracking (Reference Lang, Rabus and DechLang and others, 2004; Reference RignotRignot and others, 2004; Reference JoughinJoughin and others, 2009). Reference Rignot, Mouginot and ScheuchlRignot and others (2011a) recently presented an Antarctic-wide surface velocity map, and the data in the coastal Amundsen Sea Embayment are mostly from 2007 and 2008. Because we focus on the eastern shear margin, which is farther upstream than the largest reported accelerations, and because we wish to maintain internal consistency with the temperature model later used to estimate englacial radar attenuation, here we use the 1996 surface velocity field presented by Reference JoughinJoughin and others (2009).
We identify the position of the eastern shear margin by finding coherent maxima in the absolute value of the lateral shear strain rate field (Fig. 1), which is calculated from the gridded surface velocity field (Appendix). The width of the shear margins is typically <3 km and we estimate the uncertainty in our picks of maxima to be <1 km, based on the typical across-flow width of the peak in . This width is similar to that observed for the Ross Sea ice streams, but measured strain rates there are significantly higher than for Thwaites Glacier (>0.07 a−1 for Whillans Ice Stream versus ~0.01 a−1 for Thwaites Glacier; e.g. Reference Echelmeyer, Harrison, Larsen and MitchellEchelmeyer and others, 1994). This difference could be due to dynamical differences between these glaciers, the finer horizontal resolution of those field measurements compared with the 1 km grid spacing of the remotely sensed surface velocity field we used, or smoothing of the surface velocity field (Reference Joughin, Tulaczyk, Bindschadler and PriceJoughin and others, 2002).
Regional bed elevation
For a regional perspective of bed elevation across Thwaites Glacier (Figs 1a and 2), we use the 1 km resolution Bedmap2 bed elevation grid (Reference FretwellFretwell and others, 2013). In the vicinity of the eastern shear margin, this grid is derived mostly from the 60 MHz ice-penetrating radar transects collected by the Airborne Geophysical Survey of the Amundsen Sea Embayment, Antarctica (AGASEA) (Reference HoltHolt and others, 2006; Reference VaughanVaughan and others, 2006). Those AGASEA bed elevation data are based on bed picks made from unfocused SAR data, i.e. bed echoes that were coherently summed along-track without performing range migration.
Bed properties
The obfuscating effects of surface clutter and volume scattering from crevasses can be substantial across portions of Thwaites Glacier. To reduce these effects as much as possible and to improve along-track resolution, we study changes in radar-detected bed properties across the eastern shear margin using range-migrated SAR-focused airborne radar data collected by AGASEA. Reference Peters, Blankenship and MorsePeters and others (2005, Reference Peters, Blankenship, Carter, Kempf, Young and Holt2007) summarized the radar system’s characteristics and associated processing techniques. The bed was picked in the SAR-focused data using a semi-automatic algorithm that determines the location of the maximum bed echo intensity.
Bed reflectivity (R b) is calculated using the radar equation (e.g. Reference Peters, Blankenship, Carter, Kempf, Young and HoltPeters and others, 2007):
where P
r is received power from the bed (bed echo intensity), P
t is the transmitted power, λ
air is the radar wavelength in air, G
a is the antenna gain, T is the transmission loss at the air–ice interface, L
a is the one-way loss due to dielectric attenuation, Ls
is the total system loss, Gp is the processing gain, h is the height of the aircraft above the ice surface, H is the ice thickness and 
is the real part of the complex relative permittivity of ice. The values of the static system parameters in Eqn (1) are given by Reference Peters, Blankenship, Carter, Kempf, Young and HoltPeters and others (2007) and summarized here in Table 2.
Table 2. Values of static parameters in the radar equation (Eqn (1))
Because of the challenge of calibrating inferred bed reflectivity values across regions of an ice sheet far inland from the grounding line, we instead consider the spreading-and attenuation- ‘corrected’ bed echo intensity 
, which is functionally identical to R
b in Eqn (1) and is also reported in decibels. We assume that 
does not necessarily represent the absolute reflectivity of the bed, but that its spatial variation does accurately capture the spatial variation of R
b. Implicit in this approach is the assumption that the unknown difference between and R
b did not change during the AGASEA survey and that this difference simply represents an unknown system calibration factor that is not presently included in Eqn (1).
For analyses of observed echo intensities using the radar equation (Eqn (1)), the most poorly constrained term is typically the loss due to dielectric attenuation (e.g. Reference Peters, Blankenship and MorsePeters and others, 2005, Reference Peters, Blankenship, Carter, Kempf, Young and Holt2007; Reference MatsuokaMatsuoka, 2011). Radio-frequency dielectric attenuation in polar ice depends primarily on temperature but also on impurity concentrations. To correct observed bed echo intensities for attenuation, we use the radar attenuation model presented by Reference MacGregor, Winebrenner, Conway, Matsuoka, Mayewski and ClowMacGregor and others (2007), which accounts for both temperature and chemistry. Chemistry data are not yet available for the recently drilled deep ice core at the ice divide between the Ross and Amundsen Sea Embayments. The nearest deep ice core for which chemistry data are available is Siple Dome in the Ross Sea Embayment and here we assume that englacial impurity concentrations are uniform across our study area and are equal to the depth-averaged values at Siple Dome (Reference MacGregor, Winebrenner, Conway, Matsuoka, Mayewski and ClowMacGregor and others, 2007). The effect of the spatial variation of impurity concentrations should be small compared with that of englacial temperature (Reference MacGregor, Matsuoka, Waddington, Winebrenner and PattynMacGregor and others, 2012b). Significant uncertainty remains in the parameterization of radar attenuation through ice (Reference MacGregor, Winebrenner, Conway, Matsuoka, Mayewski and ClowMacGregor and others, 2007), but here we ignore this uncertainty in our analysis because temperatures across the eastern shear margin likely vary smoothly and because we do not evaluate absolute values of bed reflectivity in terms of basal conditions.
To calculate temperature-dependent attenuation rates, we use a steady-state temperature model for Thwaites Glacier (Reference JoughinJoughin and others, 2009). This model includes vertical diffusion of heat, horizontal and vertical advection of heat, vertical strain heating, and frictional heating at the bed. However, this temperature model does not include the effect of strain heating due to lateral shear and thus likely underestimates temperatures in the vicinity of the margins that are the focus of this study (Reference RaymondRaymond, 1996; Reference Harrison, Echelmeyer and LarsenHarrison and others, 1998; Reference Jacobson and RaymondJacobson and Raymond, 1998; Reference Raymond, Echelmeyer, Whillans, Doake, Alley and BindschadlerRaymond and others, 2001). To address this limitation, we estimate the temperature increase due to lateral shear strain heating (ΔT
lat) and treat it as a perturbation to the original model. First, we estimate the depth-averaged strain heating rate due to lateral shear 
across Thwaites Glacier as
where 
is the depth-averaged lateral shear stress (Appendix), ρ is the density of pure ice (917 kg m−3) and 
is the temperature-dependent depth-averaged specific heat capacity of pure ice, calculated following Reference Cuffey and PatersonCuffey and Paterson (2010, p.400). Next, we determine the reversed horizontal particle paths p originating at each gridpoint (x,y) using the reversed surface velocity field 
, and then calculate ΔT
lat during the period represented by the reversed particle path as
where dp is the regular horizontal spacing of points along the particle path p(x, y) and (x 0, y 0) is the upstream end of the particle path. This approach ignores diffusion of the additional heat generated by lateral shear (vertical diffusion is included in the original temperature model) and also assumes that ice flow was steady during the period taken for the particle to traverse its reversed path. Surface crevassing is not extensive across the eastern shear margin, so we ignore the effect that cold air pooling in crevasses has upon the near-surface thermal structure of the ice observed in more heavily crevassed margins (Reference Harrison, Echelmeyer and LarsenHarrison and others, 1998).
Based on surface velocities with a horizontal resolution of 1 km, 
is of the order of 10−3 K a−1 within the margins and up to 10−2 K a−1 within the intraglacier shear zones. The upper end of this range is similar to values inferred from borehole thermometry by Reference Harrison, Echelmeyer and LarsenHarrison and others (1998) for a heavily crevassed margin of Whillans Ice Stream, although the horizontal resolution of their measurements was finer (~100–150 m) than that of the velocity grids used here. The total depth-averaged temperature perturbation is <1 K in the eastern shear margin, which is consistent with modeled patterns of heat generation at shear margins (Reference Jacobson and RaymondJacobson and Raymond, 1998). Without the temperature perturbation due to lateral shear strain heating, the mean difference between L
a on the fast and the slow sides of the eastern shear margin is −0.6 ± 1.9 dB. When this perturbation is included, the mean difference is −0.3 ± 2.0 dB. Although this correction is relatively small, we still include it in our attenuation model. Modeled attenuation rates increased from 13 to 18 dB km−1 between the upstream and downstream ends of the eastern shear margin.
Margin Bed Properties
Bed elevation and reflectivity
All 16 radargrams that we examined are shown in Figure 3. An example transect is shown in more detail in Figure 4. The bed echo intensities derived from the radar data are much more variable along-transect compared with the relatively smooth bed topography. This variability is due to a combination of several factors: (1) bed roughness at scales comparable with the radar wavelength (λ
ice = 2.8 m for the 60 MHz data used here) (Reference Grima, Kofman, Herique, Orosei and SeuGrima and others, 2012); (2) clutter from the ice-sheet surface; (3) off-nadir signal; and (4) other noise in the focused data (Reference Peters, Blankenship, Carter, Kempf, Young and HoltPeters and others, 2007). To better distinguish across-margin differences in bed properties, we average values of bed elevation and bed echo intensity (uncorrected and corrected) across 5 km on either side of the margin (
and 
, respectively). We then calculate the best-fit mean difference between the fast- and slow-flowing sides of the margin using a χ
2 minimization that accounts for the standard deviations of values on both sides of the margin. For bed elevation, this relationship is simply
Fig. 3. Radargrams across all transect–margin intersections shown in Figure 2, shown from the upstream (top left) to downstream (bottom right) end of the eastern shear margin. The horizontal and vertical scales are the same for all radargrams. Symbols in the lower right corners of each panel follow Figure 2.
Fig. 4. (a) Example radargram across the eastern margin (transect Y17a; labeled in Fig. 2 and also shown in Fig. 3n). (b) Surface and bed elevation. (c) Surface speed and lateral shear strain rate. (d) Uncorrected and corrected bed echo intensity. Mean values of each property on either side of the margin are shown as horizontal dashed lines.
where Δz
b is the best-fit mean difference between the bed elevation on the fast (inside) and slow (outside) sides of the margin for each transect; identical expressions are used to determine the best-fit values of ΔP
r and 
. We chose a distance of 5 km on either side of the margin over which to average the data because this distance is 1.5–4.0 times the ice thickness across the eastern shear margin and should encompass its zone of thermodynamic influence (e.g. Reference Echelmeyer, Harrison, Larsen and MitchellEchelmeyer and others, 1994; Reference Jacobson and RaymondJacobson and Raymond, 1998) and any transition in bed reflectivity (e.g. Reference Catania, Conway, Gades, Raymond and EngelhardtCatania and others, 2003). For Δz
b, we correct for the deviation of the transect–margin intersection from normal incidence, i.e. the apparent ‘dip’ of the bed. For ΔP
r and , we averaged the data as linear power and converted back to decibels.
On the slower-flowing side of the eastern shear margin, bed elevations are higher for 11 of 16 transect–margin intersections, bed echo intensities are lower for 12 of 16 intersections, and bed reflectivities are lower for 11 of 16 intersections (Fig. 5). Because ice is generally thicker on the faster-flowing side, we might expect lower bed echo intensities on that side due to higher geometric spreading and attenuation losses if bed reflectivity was uniform across the margin (Eqn (1)). Instead, we observe the opposite pattern. After correcting for these losses, bed echo intensity remains higher on the fast-moving sides of the margin (Figs 4d and 5c) and the across-margin difference in the attenuation correction due to lateral shear strain heating is too small to explain this difference.
Fig. 5. (a) Mean bed elevation, (b) bed echo intensity and (c) corrected bed echo intensity for each transect 5 km on each side on the eastern shear margin. Each point represents a different transect–margin intersection, and their labeling follows Figure 2. Vertical (horizontal) error bars represent the standard deviation of each quantity on the slow (fast) side of the margin. Gray dashed lines represent 1 :1 relationships between values on the fast and slow sides of the margin.
The observed across-margin bed echo intensity differences are generally smaller than the predicted difference between frozen and unfrozen beds (10 dB or more; Reference Peters, Blankenship and MorsePeters and others, 2005). These smaller differences could be due to basal conditions that are more complex than simply frozen or thawed, suggesting a more complex subglacial environment that is qualitatively consistent with the predicted style of subglacial hydrology in ice-sheet interiors (e.g. distributed drainage; Reference Carter, Blankenship, Young, Peters, Holt and SiegertCarter and others, 2009). Alternatively, across-margin differences in bed roughness may cause the differences, which we consider in detail below.
Bed roughness
Topography
There are now several methods for quantifying bed roughness from bed topography, mostly using spectral methods to study macroscale (>500 m) bed roughness (e.g. Reference Bingham and SiegertBingham and Siegert, 2009; Reference LiLi and others, 2010; Reference Rippin, Vaughan and CorrRippin and others, 2011). Other methods examine the statistics of bed elevation and the bed echoes themselves (Reference Rippin, Bamber, Siegert, Vaughan and CorrRippin and others, 2006; Reference Oswald and GogineniOswald and Gogineni, 2008). Collectively, these studies have considered a large range of roughness length scales (<1 m to 35 km) and found evidence of power-law relationships between roughness within this range of horizontal scales (e.g. Reference Hubbard, Siegert and McCarrollHubbard and others, 2000), implying a degree of self-affinity. Our focus on across-margin differences in bed properties precludes a straightforward comparison with the spectral roughness estimates from earlier studies because they investigated much larger horizontal scales of both roughness and variability in this larger-scale roughness. Furthermore, if the across-margin differences in bed echo intensity are due to roughness variations, then they will be most sensitive to roughness at scales close to the radar wavelength in ice (Reference Ulaby, Moore and FungUlaby and others, 1982). For a 60MHz radio wave (wavelength in ice of 2.8 m), roughness at this scale cannot easily be measured directly by established spectral methods.
In its simplest form, bed roughness can be calculated as the standard deviation of the linearly detrended along-track bed topography 
within a horizontal range of interest (Reference Shepard, Campbell, Bulmer, Farr, Gaddis and PlautShepard and others, 2001). This form of bed roughness is the root-mean-square (rms) height ζ:
where 
is the mean detrended bed elevation within the horizontal range of interest.
Bed roughness at a specified horizontal length scale Δx is calculated as the rms deviation v of at intervals of Δx (e.g. Reference Shepard, Campbell, Bulmer, Farr, Gaddis and PlautShepard and others, 2001):
The latter approach has also been used to study large-scale bed roughness of ice sheets (Reference YoungYoung and others, 2011; Reference RossRoss and others, 2012). Calculating the rms deviation for a range of horizontal length scales produces a bed roughness ‘spectrum’ analogous to that produced by spectral methods. An obvious lower limit for Δx is the mean horizontal posting of the focused radar data Δx min =16.5 m), which is an upper bound on the horizontal resolution of these data. A suitable upper limit is not as well constrained; Reference Shepard, Campbell, Bulmer, Farr, Gaddis and PlautShepard and others (2001) suggested an upper limit to Δx equivalent to 10% of the length of the profile. For our analysis of 5 km on each side of the margin, this upper limit is 500 m. We calculate v(Δx) across this range (16.5–500 m; Fig. 6) in order to later predict its value at smaller scales, because we are interested in the smaller-scale roughness that may influence bed reflectivity and/or basal drag. The roughness of natural surfaces is often self-affine across a wide range of length scales, i.e. roughness at one length scale often has a power-law relationship to roughness at another length scale (Reference Shepard, Campbell, Bulmer, Farr, Gaddis and PlautShepard and others, 2001). For rms deviation, this relationship is expressed as
where v 0 is the rms deviation at the horizontal unit scale (m), and u is the Hurst exponent, which is a measure of the self-affinity of a surface. The plausible range of u is 0 ≤ u ≤ 1; as u → 1, bump height approaches bump wavelength at all length scales. We calculate the best-fit values of u and v 0 for each side of each transect–margin intersection using linear regression in log space.
Fig. 6. Root-mean-square (rms) deviation of along-transect bed elevation vs horizontal length scale for each transect across the eastern shear margin, with colors following Figure 2. Vertical gray dashed line represents the mean posting between radar traces (Δx min =16.5 m), i.e. the smallest horizontal scale for which rms deviation was calculated directly. Horizontal lines represent the mean values of rms deviation at the mean posting for the portions of the transect inside (solid lines) and outside (dashed lines) the margin.
For the 32 rms deviation profiles shown in Figure 6 (one for each side of 16 transect–margin intersections), the mean Hurst exponent is 0.52 ± 0.15. The mean difference between the measured and modeled values of v(Δx min) is −1 ± 8% of the measured value and 1 ± 22% for v(Δx = 500m). Of these profiles, 78% have data–model differences that are normally distributed, based on an a posteriori 95% confidence χ 2 test, and the mean p value of all data–model differences is 0.06 ± 0.15. This good agreement across the range of horizontal length scales for which we calculated v(Δx) suggests that bed elevation is indeed self-affine across this range.The mean rms deviation at the shortest measured length scale v(Δx min) is 3.7 ± 1.7 m on its fast side and 4.1 ± 1.9 m on its slow side (Fig. 6); these values are not significantly different (p = 0.26). Based on best-fit estimates of v 0 (Eqn (7)), rms deviations at a horizontal scale of 1 m are both predicted to be on average 1 m on either side of the margin, i.e. they are indistinguishable given the ability of our method to resolve bed roughness at this scale. However, any physical interpretation of v 0 assumes that the scale dependencies of the physical processes governing erosion, deposition and ultimately bed roughness at the scale of 1 m are unchanged between 1 m and Δx min
Relationship with bed reflectivity
Because rough surfaces scatter reflected energy over a wider range of angles than smooth surfaces, the amplitude of nadir-reflected bed echoes decreases with increasing bed roughness (Reference Ulaby, Moore and FungUlaby and others, 1982). For unfocused radar-sounding data, most of the recorded energy is received from either nadir or near-nadir (within the first few Fresnel zones). The roughness-induced reflectivity decrease of a reflection from nadir (ρ n) can be modeled using Kirchhoff theory (Reference OgilvyOgilvy, 1991) following Reference Peters, Blankenship and MorsePeters and others (2005) as
where 
and I
0 is the zero-order modified Bessel function. The quantity g can be thought of as the wavelength-normalized rms height 
. Reference Peters, Blankenship and MorsePeters and others (2005) describe this roughness metric (using the notation S) as ‘rms deviation’ (their eqn (8)). However, both their description of this quantity and a review of the origin of our Eqn (8) (Reference BoithiasBoithias, 1987) strongly suggest that rms height across a physically meaningful length scale should be used, not rms deviation at a specific horizontal scale. The subscript D
1 signifies that is calculated across the width of the first Fresnel zone, the radius D
1 of which is calculated as
For the 10 km long portions of the 16 transects that we studied, the mean value of D 1 is 134 ± 7 m.
Figure 7 shows the modeled effect of increasing roughness upon reflectivity based on Eqn (8) for several common ice-penetrating radar frequencies, including that used in this study. This relationship between unfocused radar reflectivity and roughness is consistent with several expectations: (1) reflectivity decreases with increasing roughness; (2) at a fixed bed roughness, the reflectivity decrease is smaller at larger englacial radar wavelengths (i.e. lower frequencies); and (3) where the bed is sufficiently rough, diffuse scattering begins to dominate the returned energy over that which is coherently reflected, so reflectivity does not decrease as rapidly.
Fig. 7. Modeled reflectivity decrease due to roughness for three common ice-penetrating radar frequencies, including the center frequency of the radar data used in this study (60 MHz). Englacial wavelengths (λ ice) for each center frequency are given in the legend. Solid lines represent the modeled reflectivity decrease at nadir for unfocused data (Eqn (8)). The rms height range is prescribed and similar to the observed range across the eastern shear margin. The dashed line is for the focused AGASEA data used in this study (Eqn (9)), assuming an aircraft altitude of 735 m, an ice thickness of 2461 m and a focusing window length of 2 km, which are typical values for the AGASEA data we evaluated. This modeled reflectivity decrease for focused data should not be considered valid for all SAR-focused radar surveys of ice sheets because ice thickness and survey design can vary significantly.
Equation (8) considers only energy reflected from nadir, so it is not strictly applicable to the range-migrated focused data used in this study. This focusing uses matched filters to integrate received energy across a wider aperture (~±120° from nadir) than a purely nadir return, so that energy reflected from a rough bed off-nadir but within the focusing aperture is included in the focused data. In other words, focusing partly corrects for the effect of bed roughness upon observed bed reflectivity. We next seek a relationship between bed reflectivity and roughness for focused radar data that satisfies three physically based assumptions: (1) any appropriate expression should approximate Eqn (8) for nadir-reflected energy within the first Fresnel zone; (2) off-nadir reflected energy is Gaussian-distributed and energy-conserving, following Reference Nayar, Ikeuchi and KanadeNayar and others (1991); and (3) off-nadir energy in the along-track direction that is reflected within the range of angles spanned by the Doppler bandwidth of the focusing matched filter (i.e. its focusing aperture) is completely captured in the range-migrated focused data we used.
Following the above assumptions, this received energy can be modeled as the integral across the focusing aperture for a Gaussian function the width (effectively along-track length) of which is determined by the bed roughness. Increasing the width of the focusing aperture limits the decrease in along-track contribution to bed reflectivity due to increasing bed roughness. However, such mitigation is not possible in the across-track direction for the AGASEA radar system, because the bed is viewed from only one across-track look angle. The across-track contribution to bed reflectivity is therefore simply the integral of returned bed power across the first Fresnel zone (as for both the along-track and across-track directions in Eqn (7)) and is independent of the along-track focusing. The product of these two integrations (along-track and across-track) is the total effect of bed roughness upon its reflectivity. This product can be expressed compactly using the error functions that result from integrating Gaussian functions. For the mean aircraft survey height (735 m), ice thickness (2461 m) and focusing aperture of the range-migrated focused AGASEA radar-sounding data we used, the reduction in focused energy (ρ f) as a function of bed roughness is
where the first error function represents the integration of energy from the along-track focusing aperture and the second error function represents the integration of energy from the across-track angles spanned by the nadir-directed echo from the first Fresnel zone. As expected, Eqn (10) predicts a smaller reflectivity loss at any given roughness compared with Eqn (8) (Fig. 7).
Equation (10) predicts the reflectivity decrease due to roughness compared with a purely specular reflection, not the absolute bed reflectivity of the ice–bed interface. To evaluate these models, we compare Δρ
f (the difference between ρ
f values inside and outside the margin) to the observed differences in 
at the same locations (Fig. 8). These across-margin differences are independent of the absolute bed reflectivity R
b, but should be indicative of relative bed reflectivity differences and hence also corrected bed echo intensity differences 
. For the eastern margin, the correlation coefficient between the observed and modeled reflectivity differences is r = 0.67. A unity correlation between the observed differences and Δρ
f would imply that bed roughness completely controls the across-margin variability in 
. In general, Eqn (10) overestimates the observed reflectivity difference. This implies that, if all the observed reflectivity variability is due to roughness, then the focused radar data are less sensitive to bed roughness than Eqn (10) implies.
Fig. 8. Observed mean corrected bed echo intensity differences between inside and outside the margin vs the bed reflectivity differences predicted by Eqn (8), with labeling following Figure 2. Positive (negative) differences indicate that observed/modeled reflectivity is higher outside (inside) the margin, implying the bed inside (outside) the margin is rougher. Gray dashed line represents a 1 :1 relationship between the observed and predicted reflectivity differences.
Discussion
Nature of across-margin differences
Differences in bed elevation and reflectivity across the eastern shear margin suggest that both bed topography and some other property or combination of properties of the ice– bed interface influence the margin’s position. The relationship with bed topography is not surprising, given the proportionality between ice thickness and driving stress. However, driving stress does not increase consistently across the eastern shear margin. Both sliding and deformation speeds are related nonlinearly to ice thickness. However, calculations using the temperature model applied above suggest that deformation speeds are <10 m a−1 across most of Thwaites Glacier. Thus, the increase in surface velocity across the eastern shear margin (e.g. a doubling in speed along transect Y17a in the vicinity of the eastern shear margin; Fig. 4c) is unlikely to be due to changes in deformation speed. Reference Anandakrishnan, Blankenship, Alley and StoffaAnandakrishnan and others (1998) made a similar argument to support their inference (from seismic data) that changes in subglacial geology affect ice flow across the margin of Whillans Ice Stream.
The broadly negative correlation between bed elevation and reflectivity across the eastern shear margin implies that either the dielectric contrast of the ice–bed interface, its roughness or a combination of these two properties changes there. At 60 MHz, a bed reflectivity change of ~5 dB beneath grounded ice suggests a change between unfrozen till and unfrozen bedrock, while a change of >10 dB suggests a change between a frozen and unfrozen bed, whether till or bedrock (Reference Peters, Blankenship and MorsePeters and others, 2005). We observed bed reflectivity differences both below and above 10 dB, implying that both transitions are present, but there was no clear along-transect pattern (Fig. 5c). On average, these differences are <5 dB, which suggests that the change in the nature of the bed is generally lithological or morphological, rather than thermal. The importance of an unfrozen bed that has a pressurized subglacial hydrological system for inducing significant basal sliding is well established, while the role of bed roughness and smaller changes in lithology is less clear (Reference Cuffey and PatersonCuffey and Paterson, 2010).
A possible change in lithology may be related to the change in bed elevation. Virtually all of Thwaites Glacier is grounded well below sea level (Reference HoltHolt and others, 2006; Fig. 1). During interglacial periods when Thwaites Glacier was thinner than at present, areas outside its margins would have remained grounded longer because they are generally higher. Thus, we might expect thinner layers of potentially deformable marine sediment deposited outside the margins (or even their absence), as has been suggested for the Ross Sea ice streams (Reference JoughinJoughin and others, 1999). However, the small mean difference in bed elevation across the eastern shear margin (43 ± 69 m) suggests that any lithological change across this margin due to this hypothesized spatial variability in depositional environment is unlikely to be significant.
We next consider the possible role of changes in bed roughness in modulating the position of the eastern shear margin, as our analysis suggests that bed roughness changes can explain much of the change in bed reflectivity observed there. Reference Rippin, Vaughan and CorrRippin and others (2011) examined bed roughness at horizontal scales of several hundreds of meters across Pine Island Glacier using spectral methods. They found that ice-flow speed and bed roughness were inversely related on regional scales, and a rougher bed outside at least one tributary. However, they also found smooth beds in some inland areas of slow ice flow. Based on those observations, they argued for an indirect relationship between bed roughness and ice-flow speed. They suggested that smooth beds do not necessarily increase rates of basal sliding, but rather that they lower the basal drag that mostly balances the driving stress. Our results are consistent with their observations and can be used to refine their hypothesis because we focus on bed roughness across large transitions in ice dynamics and we resolve finer scales (down to 16.5 m) of bed roughness that are much closer to the controlling obstacle size that most strongly affects basal drag (<1 m; e.g. Reference Cuffey and PatersonCuffey and Paterson, 2010).
In general, faster-flowing ice overlies smoother beds, but this relationship is apparently nonlinear, as suggested by Reference Rippin, Vaughan and CorrRippin and others (2011) and predicted by theory (e.g. Reference Cuffey and PatersonCuffey and Paterson, 2010, p. 229–238). For the radar transect that crosses the eastern shear margin shown in Figure 4, the surface speed approximately doubles across the margin. Because both sides of the margin are believed to be flowing mostly by basal sliding (based on model estimates of deformation speed), the speed difference probably does not indicate the onset of basal sliding, but rather a change in its magnitude. While bed roughness and reflectivity are overall related across the eastern shear margin, we did not observe significant differences in bed roughness at the shortest resolvable length scales. Using the same radar data as this study, Schroeder and others (in press) also did not observe any clear change in the specularity of bed echoes across the eastern shear margin, further suggesting that bed roughness did not change significantly there. These results emphasize the limitations of our approach to resolving bed roughness and our understanding of its relationship to ice flow, particularly for the soft bed that is likely present along the eastern shear margin.
The magnitude of lateral drag peaks in the vicinity of margins and it reduces the basal drag required to balance the driving stress locally inside the margin. Given their relationship with margin position, the variability in bed roughness we observed is a plausible but not yet definitive mechanism for explaining the expected change in basal drag. Changes in subglacial hydrology or some other bed property (e.g. till rheology) are also plausible and the occasionally large bed reflectivity changes we observed (>10 dB; Fig. 5c) suggest that hydrologic transitions may also be present along short portions of the margin. However, bed roughness provides the most consistent explanation for the across-margin difference in bed reflectivity at the scale we considered.
The basal drag modelling results of Reference JoughinJoughin and others (2009) do not extend beyond Thwaites Glacier’s outermost shear margins, so we cannot discuss our across-margin observations in terms of existing predictions of the spatial variation of basal drag. Within the glacier, Reference JoughinJoughin and others (2009) observed several regions of high basal drag oriented across-flow and separated by regions of very low basal drag. This pattern implies the presence of both hard bedrock and weak till underneath Thwaites Glacier and frames an important question regarding the bed roughness variability we observed: is it a control on, or a consequence of, ice flow? Theoretical studies of the relationship between bed roughness and basal drag are often predicated upon the notion that the bed is hard, but that is unlikely to be the case here. Additional field observations and models are necessary to resolve the relationship between bed roughness and basal drag for soft beds, as also suggested by Reference Rippin, Vaughan and CorrRippin and others (2011). Other West Antarctic glaciers can mobilize basal sediment on decadal timescales (Reference SmithSmith and others, 2007, Reference Smith, Bentley, Bingham and Jordan2012), so it is certainly plausible that Thwaites Glacier has also eroded its bed differentially in a manner that altered basal drag since the last collapse of the West Antarctic ice sheet ~200 ka ago (Reference Pollard and DeContoPollard and DeConto, 2009). We hypothesize that the bed roughness variability we observed is a consequence of ice flow as sediment is remobilized.
Our results also hint at a possible explanation for the perplexing result reported by Reference Raymond, Catania, Nereson and Van der VeenRaymond and others (2006), who observed no change in bed reflectivity across the ‘Snake’ shear margin of Whillans Ice Stream using 2 MHz ice-penetrating radar. Figure 7 shows that 2 MHz radar is not sensitive to bed roughness at vertical scales of <10 m. For Thwaites Glacier, observed rms deviation at a horizontal scale of 85 m (the wavelength of 2 MHz radar in ice) across the eastern shear margin is ~5 m, although it varies between 2 and 20m (Fig. 6). At 2 MHz, an rms deviation of 5m would produce a relatively small reflectivity decrease of ~3 dB in bed reflection power (Fig. 7). If the spatial variation of bed roughness across the Snake is similar to, or lower than, that of Thwaites Glacier’s eastern shear margin, then it is plausible that this roughness pattern could be present without having been detected by 2 MHz radar. Whether this pattern modulates ice flow there or vice versa would be subject to the same uncertainties as described above for Thwaites Glacier’s eastern shear margin.
Implications for margin stability
In part due to their internal dynamics, the Ross Sea ice streams migrate and undergo changes in flow speed (e.g. Reference Hulbe and FahnestockHulbe and Fahnestock, 2007; Reference Catania, Hulbe, Conway, Scambos and RaymondCatania and others, 2012). Comparatively little is known regarding the long-term flow history of Thwaites Glacier, although its modern thinning rates are an order of magnitude faster than those of adjacent glaciers averaged over the last 4.7–14.5 ka (Reference Johnson, Bentley and GohlJohnson and others, 2008). This limited understanding of the long-term flow history of Thwaites Glacier motivates evaluation of its present ice-flow configuration.
Given the across-margin differences we observed and the limited bed controls that they imply, we can make some inferences regarding the stability of Thwaites Glacier’s eastern shear margin. This margin is slightly higher, dimmer (in terms of radar reflectivity) and rougher on its slow side (Figs 5 and 6). This pattern suggests one of two scenarios: (1) the modest across-margin differences in bed properties are sufficient to change the basal drag and hence the driving stress, producing a distinct shear margin; or (2) the edge of the Byrd Subglacial Basin, where the margin’s trace is more clearly normal to the regional bed slope (Fig. 2), is sufficient to initiate the eastern shear margin and maintain a dynamical configuration at odds with the regional bed topography farther downstream, which otherwise suggests a more westerly path. Reference Rippin, Vaughan and CorrRippin and others (2011) found that Byrd Subglacial Basin is relatively smooth at horizontal scales of several hundred meters, but that a rougher bed exists eastward of that basin, including the area crossed by the eastern shear margin. Given the self-affinity of bed roughness (e.g. Fig. 6), their results also seem to suggest a more westerly path for the eastern shear margin. The first scenario (across-margin differences alone produce the margin) would seem to be less stable, because thickness changes across Pine Island and Thwaites Glaciers would likely alter the local force balance more quickly, but we cannot yet determine which is correct.
A qualitative investigation of the internal stratigraphy in the vicinity of the eastern shear margin found the expected dominant pattern, i.e. the internal stratigraphy is smooth and draped over the bed topography (Reference Karlsson, Rippin, Vaughan and CorrKarlsson and others, 2009). Similarly, we observed no clear change in the internal stratigraphy across the eastern shear margin, nor did we observe any features that suggested recent or longer-term migration of its position. However, we did not trace the internal stratigraphy in detail in this study, which could provide more quantitative insight into the strain history of the ice sheet in this region. Firmer conclusions regarding the strain history there must therefore await further research.
The maximum difference between the position of the eastern shear margin derived from the 1996 and 2007–08 surface velocity datasets is ~1 km (Reference JoughinJoughin and others, 2009; Reference Rignot, Mouginot and ScheuchlRignot and others, 2011a), with no clear trend of inward or outward migration (Fig. 2). If we treat that maximum difference as a real migration, then the maximum margin migration rate over this ~12 year period is <0.1 km a−1.This value is within the range of migration rates reported for portions of Whillans Ice Stream (Reference Harrison, Echelmeyer and LarsenHarrison and others, 1998; Reference Stearns, Jezek and Van der VeenStearns and others, 2005), but that ice stream is likely much more vulnerable to margin migration than Thwaites Glacier due to the former’s dynamical configuration (i.e. margin-controlled force balance; Reference Raymond, Echelmeyer, Whillans, Doake, Alley and BindschadlerRaymond and others, 2001). Given that the 1996 and 2007–08 surface velocity fields are gridded at 1.0 and 0.9 km, respectively, the null hypothesis that the eastern shear margin remained stable during this period cannot be invalidated. The small observed differences are likely either gridding artifacts, noise in the surface velocity fields or errors in our margin tracing.
Finally, we note that the downstream end of Thwaites Glacier’s eastern shear margin is <50 km from the onset of an unnamed outlet glacier between Thwaites and Pine Island Glaciers (Fig. 9). This glacier was labeled ‘SW tributary’ by Reference MacGregor, Catania, Markowski and AndrewsMacGregor and others (2012a), although in the strictest sense it is not a tributary of Pine Island Glacier but rather of its ice shelf. Since late 2001, the southern shear margin of Pine Island Glacier’s ice shelf has retreated progressively by ~5 km and rifts have grown along this shear margin (Reference MacGregor, Catania, Markowski and AndrewsMacGregor and others, 2012a). This glacier is also presently thinning (Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012), which is likely due to its recent acceleration (Reference RignotRignot, 2008) and connected to increased sub-ice-shelf melting in the eastern Amundsen Sea Embayment (Reference Jacobs, Jenkins, Giulivi and DutrieuxJacobs and others, 2011). If this ice-shelf rifting and retreat evolves in the coming decades in a manner similar to that of the more dramatic rifting and retreat of this ice shelf’s northern shear margin (a 21 km long rifted zone has developed and disintegrated there over the same period) and other outlet glaciers in this region, then this glacier will eventually lose whatever buttressing is currently provided by Pine Island Glacier’s ice shelf (Reference MacGregor, Catania, Markowski and AndrewsMacGregor and others, 2012a).
Fig. 9. (a) Surface velocity (Reference Rignot, Mouginot and ScheuchlRignot and others, 2011a) and recent coastline history (Reference MacGregor, Catania, Markowski and AndrewsMacGregor and others, 2012a) of Thwaites and Pine Island Glaciers in the vicinity of the downstream end of the eastern shear margin. Background image is the 2008–09 MODIS Mosaic of Antarctica (MOA) and the MOA grounding line is derived from this image (personal communication from T.A. Scambos and T.M. Haran, 2013). The interferometric SAR (InSAR)-derived grounding line is from Reference Rignot, Mouginot and ScheuchlRignot and others (2011b). (b) Expanded view in the vicinity of the southern shear margin of Pine Island Glacier’s ice shelf. The legend gives years of traced coastlines from Reference MacGregor, Catania, Markowski and AndrewsMacGregor and others (2012a), except for the 2012 coastline, which was traced using the same methodology described therein.
An increase in grounding-line ice flux due to buttressing loss can alter ice dynamics far upstream of the grounding line (e.g. Reference Scambos, Bohlander, Shuman and SkvarcaScambos and others, 2004; Reference Scott, Gudmundsson, Smith, Bingham, Pritchard and VaughanScott and others, 2009). Given the weak bed control of the eastern shear margin that we inferred from radar observations, we speculate that the eastern shear margin would migrate outward in response to a substantially increased ice flux, expanding the fast-flow region of Thwaites Glacier. This ice-flux change could originate from the ‘SW tributary’ as described above or the mass loss of Thwaites Glacier proper (its main trunk) could continue and induce such a change. However, as described above, the eastern shear margin has not clearly migrated during a period contemporaneous with these changes.
Conclusions
We found that the dynamic transition between Thwaites and Pine Island Glaciers represented by Thwaites Glacier’s eastern shear margin is coincident with modest radar-observed changes in bed properties there. While the bed topography is generally higher, these differences do not appear to be sufficient to explain the change in ice dynamics. Bed roughness variability can explain the observed bed reflectivity change and is consistent with the pattern of ice flow, but not significantly so. It remains unclear whether this pattern is a result of or a modulator of increased basal sliding, but the likely presence of marine sediment along the length of the margin implies the former.
That we observed little change in the nature of the bed across the eastern shear margin is interesting, especially given its ambiguous relationship with regional bed topography. We did not investigate in detail the bed of the ~100 km wide region between Thwaites and Pine Island Glaciers’ tributaries, although we note that its elevation is generally several hundred meters higher than that at the eastern shear margin and that it is even higher across a dome-shaped rise between two southern tributaries of Pine Island Glacier (Fig. 1). Furthermore, Reference Rippin, Vaughan and CorrRippin and others (2011) found that the bed in this inter-tributary region was rough at larger horizontal scales. Thus, the eastern shear margin may have limited potential to migrate outwards beyond the narrow region we investigated here. There is no evidence of significant margin migration during the most recent two decades, a period during which acceleration and major coastal changes for both Thwaites and Pine Island Glaciers were observed. However, the potential for future increases in ice flux farther downstream indicates that the position of this shear margin should continue to be monitored.
Author Contribution Statement
J.A.M., G.A.C., H.C. and D.D.B. designed this study. D.M.S. developed the relationship between bed roughness and radar reflectivity. I.J. provided the temperature model and 1996 surface velocity data. D.A.Y, S.D.K and D.D.B. collected and processed the radar data. J.A.M. performed the analysis and wrote the manuscript. All authors discussed the results and edited the manuscript.
Acknowledgements
The US National Science Foundation (ANT 0424589, 0739654 and 0739372) supported this work. We thank S.F. Price, T.J. Fudge and T.A. Scambos for valuable discussions, and S.C. Rybarski and S.K. Young for tracing the 2012 coastline of the eastern Amundsen Sea Embayment. We thank the scientific editor (N. Glasser), R. Drews and an anonymous reviewer for constructive comments that improved the manuscript.
Appendix: Lateral Shear Stress
The depth-averaged lateral shear stress is (e.g. Reference Van der VeenVan der Veen, 1999; Reference Price, Bindschadler, Hulbe and BlankenshipPrice and others, 2002)
where 
is the depth-averaged temperature-dependent rate factor, 
is the effective strain rate and n is the flow-law exponent (3). We estimate using the depth-averaged temperature derived from the original temperature model. The effective strain rate is
We define x and y as the along-flow and across-flow directions, respectively; therefore, 
is explicitly zero. Following the incompressibility condition, 
, so the effective strain rate reduces to
We calculate the horizontal strain rates 
and 
from gradients of the surface speed 
projected onto the appropriate surface velocity unit vectors (
for x and 
for y) as
Equation (A5) assumes that there are no longitudinal gradients in the lateral velocity v, i.e. ∂v/∂x = 0.
