Revealing the former bed of Thwaites Glacier using sea-floor bathymetry: implications for warm-water routing and bed controls on ice flow and buttressing

The geometry of the sea floor immediately beyond Antarctica's marine-terminating glaciers is a fundamental control on warm-water routing, but it also describes former topographic pinning points that have been important for ice-shelf buttressing. Unfortunately, this information is often lacking due to the inaccessibility of these areas for survey, leading to modelled or interpolated bathymetries being used as boundary conditions in numerical modelling simulations. At Thwaites Glacier (TG) this critical data gap was addressed in 2019 during the first cruise of the International Thwaites Glacier Collaboration (ITGC) project. We present more than 2000 km² of new multibeam echo-sounder (MBES) data acquired in exceptional sea-ice conditions immediately offshore TG, and we update existing bathymetric compilations. The cross-sectional areas of sea-floor troughs are under-predicted by up to 40 % or are not resolved at all where MBES data are missing, suggesting that calculations of trough capacity, and thus oceanic heat flux, may be significantly underestimated. Spatial variations in the morphology of topographic highs, known to be former pinning points for the floating ice shelf of TG, indicate differences in bed composition that are supported by landform evidence. We discuss links to ice dynamics for an overriding ice mass including a potential positive feedback mechanism where erosion of soft erodible highs may lead to ice-shelf ungrounding even with little or no ice thinning. Analyses of bed roughnesses and basal drag contributions show that the sea-floor bathymetry in front of TG is an analogue for extant bed areas. Ice flow over the sea-floor troughs and ridges would have been affected by similarly high basal drag to that acting at the grounding zone today. We conclude that more can certainly be gleaned from these 3D bathymetric datasets regarding the likely spatial variability of bed roughness and bed composition types underneath TG. This work also addresses the requirements of recent numerical ice-sheet and ocean modelling studies that have recognised the need for accurate and high-resolution bathymetry to determine warm-water routing to the grounding zone and, ultimately, for predicting glacier retreat behaviour.

1 Introduction

Knowledge of Antarctica's coastal bathymetry is essential when considering ocean circulation and recent dynamic changes at the ice-sheet margin. Sea-floor bathymetry influences ice–ocean interactions in two ways. First, deep (>500 m water depth) bathymetric troughs and channels provide access routes for warm, salty Circumpolar Deep Water (CDW: ∼0.5–1.5 ^∘C, located below ∼300–500 m water depth; Jacobs et al., 1996, 2013) to present-day grounding zones. The inflow of CDW increases basal melting and ice-shelf thinning (Jacobs et al., 1996; Rignot et al., 2013), leading to ice-shelf disintegration, reduced buttressing, the acceleration of the ice shelves and grounded ice upstream, and ultimately grounding-zone retreat (Schoof, 2007; Joughin et al., 2010; Favier et al., 2014; De Rydt and Gudmundsson, 2016). This effect is particularly significant for ice resting on reverse-slope beds (i.e. retrograde, when the bed slopes down towards the interior of the continent) where grounding-zone retreat may initiate marine ice-sheet instability, a positive feedback that could lead to runaway retreat (Weertman, 1974; Schoof, 2007). Secondly, bathymetric highs can slow ice retreat by acting either as pinning points for floating ice or as “sticky spots” at the grounding zone itself, as well as by partially blocking warm-water access to modern grounding zones (e.g. De Rydt et al., 2014). An ice shelf pinned on a bathymetric high is subject to increased buttressing, and a topographic high at the grounding zone similarly contributes to basal drag that restricts ice flow. Both have the potential to act as stabilising influences (Alley et al., 2007; Parizek and Walker, 2010).

Geophysical mapping at marine-terminating ice-sheet margins is often difficult due to more or less persistent floating ice cover in the form of icebergs, ice tongues, ice shelves and sea ice. This is certainly the case at Thwaites Glacier (TG), West Antarctica, which is one of the two dominant fast-flowing glaciers draining into the eastern Amundsen Sea Embayment (ASE); the other being Pine Island Glacier (PIG; Fig. 1). Together, Thwaites and Pine Island glaciers were responsible for >30 % of the annual ice discharge from the West Antarctic Ice Sheet (WAIS) between 2009 and 2017 (compared with 25 % for 1979–1989; Rignot et al., 2019), and TG and adjacent smaller glaciers accounted for ca. 50 %–55 % of the annual net mass loss from the WAIS since 1992 (Shepherd et al., 2019; Smith et al., 2020). Recent observations and mass balance calculations show that TG is experiencing some of the highest rates of flow acceleration (Mouginot et al., 2014), discharge (Rignot et al., 2019), thinning (McMillan et al., 2014; Milillo et al., 2019; Shepherd et al., 2019) and grounding-zone retreat (Rignot et al., 2014; Milillo et al., 2019) across the entire ice sheet. For example, over the past four decades net mass loss from TG is calculated to have increased from 4.6 Gt yr⁻¹ during the period 1979–1989 to 34.9 Gt yr⁻¹ from 2009 to 2017 (Rignot et al., 2019). Its configuration, on a reverse-bed slope with direct connectivity to the deep WAIS interior (Holt et al., 2006), and its wide marine-terminating ice front (>120 km), with only small, unconfined frontal ice shelves, implies that TG is particularly susceptible to retreat via marine ice-sheet instability (Weertman, 1974; Hughes, 1981; Schoof, 2007; Vaughan and Arthern, 2007). Furthermore, a significant retreat in this system could lead to a much wider WAIS retreat, and the future behaviour of this glacier system is now in the spotlight (e.g. Joughin et al., 2014; Scambos et al., 2017).

Figure 1 Regional bathymetry for the Amundsen Sea Embayment and location of Thwaites Glacier (TG), West Antarctica. Bathymetry is from IBCSO (Arndt et al., 2013); arrows show observed (solid) and inferred (dashed) pathways for CDW across the continental shelf towards the grounding zones of Pine Island, Thwaites and Smith glaciers (following Nakayama et al., 2013; Dutrieux et al., 2014) and along the Dotson–Getz trough (Ha et al., 2014). PIT is Pine Island Trough; PITE is Pine Island Trough East; PITW is Pine Island Trough West; TGT is Thwaites Glacier Tongue; EIS is Eastern Ice Shelf; other ice shelves (I.S.) are also labelled. White outlines north of TG are mapped positions of the B-22A iceberg from 2002, 2010 and 2018, from south to north. Note the more blurry look of the bathymetry in front of TG (in IBCSO this bathymetry is based on the Tinto and Bell, 2011, gravity inversion and interpolation), where ship access has been hampered before austral summer 2018/2019 by persistent fast ice and the presence of the B-22A iceberg.

At present, the eastern and central parts of TG are fronted by two protruding floating ice masses, the Eastern Ice Shelf (EIS) and the Thwaites Glacier Tongue (TGT), which extend for 40–50 km beyond the grounding zone, to the west of which a 20 km wide mélange of icebergs and sea-ice exists (Fig. 2). For simplicity, we shall refer to these floating ice masses collectively as the Thwaites Ice Shelf (Fig. 2a; see Heywood et al., 2016). The TGT extends from the fastest-flowing region of TG (Fig. 2a) and has, on multi-decadal timescales, advanced up to 130 km from the grounding zone before the majority of the floating tongue has calved (MacGregor et al., 2012). As a result, marine areas beyond this have been rendered inaccessible by remnants of the TGT as they have drifted north-northwest, including the very large B-10 and B-22A icebergs that remained grounded on the continental shelf for decades after they had calved before 1965 and in 2002, respectively (Fig. 1) (Ferrigno et al., 1993; Rabus et al., 2003; MacGregor et al., 2012). In contrast, the EIS remains pinned on a sea-floor high that restricts its flow (Rignot et al., 2001; Tinto and Bell, 2011; Jordan et al., 2020) and induces shear between the EIS and TGT. Satellite imagery confirms that since 2006 increased crevassing and fracturing has weakened the shear zone between the EIS and TGT (Kim et al., 2015), with the two ice shelves remaining connected until 2010 (MacGregor et al., 2012). Due to its inaccessibility, few marine observations have been made on the inner continental shelf in front of TG (Jacobs et al., 2012). The existing oceanographic data, along with more comprehensive results from Pine Island Bay, confirm the presence of CDW in the deep troughs east and north of the EIS (Dutrieux et al., 2014; Jenkins et al., 2016) and identify these troughs as potential pathways for warm water to the TG grounding zone (Fig. 2a; Milillo et al., 2019).

Figure 2 (a) New MBES grid for the inner Amundsen Sea Embayment. Ice-velocity data from the MEaSUREs V2 dataset (Mouginot et al., 2019); grounding lines for 1992 and 2011 are from Rignot et al. (2011) and that for 2017 is from Milillo et al. (2019); red arrows delineate CDW pathways following Dutrieux et al. (2014) and Milillo et al. (2019). The black dashed line marks the boundaries of the drainage basin of Thwaites Glacier (Vaughan et al., 2001). (b) NBP19-02 data coverage versus other MBES datasets (Table 1). The dark blue coastline illustrates the ice-shelf and ice mélange extent during survey on NBP19-02 and was digitised from Landsat 8 imagery.

Figure 3 (a) High-resolution bathymetry map of the inner Amundsen Sea shelf in front of TG and the adjacent part of Pine Island Bay showing the large-scale sea-floor morphology including bathymetric troughs (Thwaites Trough, T2–T4; main axes highlighted by black dashed lines) and highs (H1–H3) that form a broad NNE–SSW ridge continuing into a ridge further offshore in Pine Island Bay, northeast of the Eastern Ice Shelf (white dashed line). (b) Mapped sea-floor landforms; streamlined features show the former flow direction of an expanded TG (white arrows). (c) Cross-sectional profiles of the bathymetric troughs; locations of profiles in (a).

Regional bathymetric compilations for the ASE shelf use multibeam echo-sounder (MBES) data where available in offshore regions (Nitsche et al., 2007, 2013; Arndt et al., 2013), with gravity inversions and a limited amount of echo-sounding data from autonomous underwater vehicles (AUVs) for sub-ice-shelf cavities (Jenkins et al., 2010; Tinto and Bell, 2011; Millan et al., 2017; Jordan et al., 2020). These datasets have identified glacially modified depressions on the continental shelf that act as conduits for CDW transport towards present-day ocean-terminating glacier margins in the ASE (Fig. 1) (e.g. Nitsche et al., 2007; Walker et al., 2007; Jacobs et al., 2012; Nakayama et al., 2013). Landward of the MBES data coverage that existed prior to this study, gravity inversions had indicated the presence of a NE–SW-trending ridge on which the EIS is pinned (Rignot, 2001; Tinto and Bell, 2011; Millan et al., 2017). Millan et al. (2017) reported that the ridge is interrupted by at least three channels with water depths between 600 and 1000 m, interpreted to be potential CDW pathways towards the grounding zone.

Considering the use of sea-floor bathymetry at higher spatial scales (than regional compilations), the analysis and interpretation of submarine glacial landforms revealed by MBES datasets provides important information on the dynamics and configuration of former glaciers and ice sheets. On the inner ASE shelf, for example, the bed of an expanded PIG has revealed the past flow direction of a large ice stream that extended more than 400 km to the continental shelf break (Evans et al., 2006; Graham et al., 2010; Jakobsson et al., 2012), as well as evidence for extensive erosion by subglacial water during past glaciations (Nitsche et al., 2013; Kirkham et al., 2019). These offshore areas also contain well-preserved information on the form and composition of the former ice-sheet bed that may, by analogy, shed light on basal conditions under the modern ice sheet (e.g. Clark et al., 2003; Ó Cofaigh et al., 2005). The roughness of an ice-sheet bed, or its variation in the vertical over a certain horizontal distance, is a primary control on basal drag, and therefore ice-flow velocity, and can be analysed for past ice-sheet beds from MBES datasets (offshore) or from satellite-derived digital elevation models (onshore) (e.g. Falcini et al., 2018). Even small (metre-scale) obstacles on the bed have been shown, theoretically, to exert critical basal drag on an overlying ice mass (e.g. Nye, 1970; Schoof, 2002). The great value in analysing MBES datasets for this purpose lies in their higher-resolution and even two-dimensional spatial coverage when compared with radar or seismic-reflection data acquired in over-ice studies (Spagnolo et al., 2017). Previous roughness analyses (from contemporary ice-sheet beds) have associated fast-flowing ice with smoother beds (e.g. Siegert et al., 2004; Rippin et al., 2011); however, recent papers acknowledge that the picture is likely to be much more complex than this with observations of fast flow occurring over even hard, rough beds, including at TG (Schroeder et al., 2014), and acknowledging that processes acting at a variety of spatial scales (including erosion and deposition) will affect spatially varying bed conditions and roughnesses (Jordan et al., 2017; Falcini et al., 2018).

Here, we present the first direct observations of sea-floor bathymetry adjacent to Thwaites Ice Shelf acquired as part of the first International Thwaites Glacier Collaboration (ITGC) cruise on RV/IB Nathaniel B. Palmer in January–March 2019 (cruise NBP19-02). In the first part of the paper, we use these data to investigate the character (bed geometry and substrate composition) of topographic highs that were former grounding zones and ice-shelf pinning points as well as to better resolve sea-floor troughs as potential modified CDW pathways to the modern Thwaites grounding zone. In the second part of the paper, we investigate the roughness characteristics of this palaeo-glacier bed as a potential analogue for the current bed of TG, and we relate flow-line roughness to the drag contribution of an overriding ice mass. We compare bed roughnesses and basal drag contributions from bathymetric profiles from the inner ASE shelf to bed profiles from upstream areas of PIG and TG, as well as to a smooth palaeo-ice-stream bed offshore the nearby Getz A Ice Shelf. To reflect these two components of the paper we describe (1) the observational (geophysical) datasets used and interpret the new bathymetric dataset, which is also provided as a publicly available stand-alone grid; and (2) the spectral approach and how we use it to quantitatively examine the roughness and basal drag contributions from former and modern TG beds. We use sea-floor landform evidence to describe both the flow of an expanded TG over the area as well as the spatial variability in bed composition over a series of topographic highs that once acted as both the grounding zone for TG and as pinning points for its ice shelf. We demonstrate that the offshore area just seaward of Thwaites Ice Shelf is an appropriate analogue for the modern grounding zone of TG, both in terms of its bed characteristics and in the effect of its rugged bed topography on ice flow, by calculating drag contribution for different scales of roughness. Finally, we highlight the importance of high-resolution MBES observational data for constraining gravity inversions and regional bathymetry compilations, which are essential boundary conditions for predictive numerical modelling experiments and for accurately calculating the flux of warm water to the grounding zone.

2 Methods I: multibeam echo-sounder (MBES) datasets

During ITGC cruise NBP19-02, the marine areas in front of TG were unusually clear of sea ice and icebergs, providing a unique opportunity for bathymetric data acquisition at the margin of Thwaites Ice Shelf (Larter et al., 2020). MBES data were acquired using a hull-mounted $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ Kongsberg EM122 echo sounder with 288 across-track beams and an operational frequency in the range 11.25–12.75 kHz. Navigation information and vessel motion information, used to correctly locate depth measurements in real time, were taken from the ship's Seapath 330, a combined GPS and motion reference unit. The MBES was configured with “high-density equidistant” beam spacing, meaning that more than one sounding can be produced per beam (up to 432), effectively increasing across-track resolution, and in “dual-ping” mode, which ensures equal across- and along-track sounding spacing. As an example, in 600 m water depth with maximum port and starboard beam angles set at 60^∘ from nadir (typical conditions and settings on NBP19-02), this results in a sounding spacing on the sea floor of ∼5–7 m; in 1200 m water depth (near maximum survey depths during NBP19-02) this sounding spacing would effectively double to 10–14 m. Here, in addition to this new dataset, we have compiled all available MBES data in the area from UK, US, German, Swedish and Korean expeditions (Table 1; Fig. 2) to produce gridded bathymetric data products. Note that the sounding spacing achievable by each MBES system varies considerably depending on the system setup, with older systems generally attaining lower spatial resolution. For example, at 1200 m water depth and maximum beam angles 60^∘ each side of nadir, the Kongsberg EM120 MBES would achieve an across-track sounding spacing of only 22 m and the SeaBeam 2112 MBES only 35 m. Together, these two systems were responsible for acquiring data from five cruises in the area (Table 1).

Table 1 Research cruises that acquired MBES data used in this compilation. IEDA MGDS is the Interdisciplinary Earth Data Alliance Marine Geoscience Data System (USA; http://www.marine-geo.org/index.php, last access: 21 January 2020); UK PDC is the United Kingdom Polar Data Centre (UK; https://www.bas.ac.uk/data/uk-pdc/, last access: 21 January 2020); BAS is the British Antarctic Survey; AWI is the Alfred Wegener Institute (Germany); LDEO is the Lamont-Doherty Earth Observatory of Columbia University; KOPRI is the Korea Polar Research Institute.

For NBP19-02, data processing was performed on board using MB-System (Caress and Chayes, 1996; Caress et al., 2020) in order to apply optimal sound velocity profiles (SVPs) to each data file and to remove erroneous soundings. Ray tracing and sea-floor depths were calculated using SVPs generated from conductivity–temperature–depth (CTD) and expendable bathythermograph (XBT) casts made during NBP19-02 (Larter et al., 2020). Most of the other bathymetry datasets were also initially ping edited during each respective cruise; however, minor additional cleaning was performed in MB-System and QPS Fledermaus (Mayer et al., 2000) after the datasets were collated when clear outliers could be easily identified. Ultimately, and to accommodate the different resolutions of the original datasets, the bathymetric sounding data were gridded in MB-System using a Gaussian weighted mean filter algorithm to produce an isometric 50 m digital elevation model (DEM) for the sea floor on the southern ASE shelf. A degree of interpolation was applied to the final grids in data gaps only, filling areas six cell widths away from cells with real soundings, i.e. for a 50 m grid interpolation will fill cells up to 300 m away from a cell with real soundings. A lower-resolution DEM (500 m grid cells) was produced for studies that typically use coarser bathymetric information. Together, the 50 and 500 m DEMs are presented as a stand-alone regional mid-resolution bathymetric dataset (Hogan et al., 2020) In addition, in order to examine the nature of specific sea-floor features (e.g. ice-shelf pinning points), higher-resolution grids (30 m grid cells) were produced where sounding densities allowed (e.g. Fig. 4). Final grids were visualised and analysed in QPS Fledermaus 7.8.6 and ArcGIS 10.6. Sea-floor trough and channel metrics (including widths, depths, symmetry, form ratio, u-/v-shape characterisation) were derived using the methods described in Kirkham et al. (2019); the reader is referred to Fig. 2 of Kirkham et al. (2019) for a graphical depiction of the channel metrics measured.

Figure 4 Detailed maps of the MBES data and their first derivative and slope over the sea-floor highs in front of Thwaites Glacier. Panels (a) and (c) show the H3 high. Panels (b) and (d) show the H2 high and western flank of H1. Red arrows in (a) and (b) point to gullies incised into the seaward flanks of the highs; the white lines in (c) and (d) mark glacial lineations; bl shows the isolated blocks of H3, and black arrows in (a) and (c) denote their possible transport paths, with the black dotted lines in (a) illustrating semicircular indentations; white arrows in (c) point to a channel at the base of one of the blocks; blue arrows in (a) indicate the downslope transport direction of material into sediment fans at the front of H3 highlighted by bulges in the contours (contour levels 1100, 1125 and 1150 m). Black dashed lines in (c) and (d) locate the profiles in Fig. 5a. GZW is grounding-zone wedge.

3 Results I: a new bathymetric compilation for the inner Amundsen Sea Embayment shelf

Our bathymetric compilation includes more than 2000 km² of new MBES data between the EIS and the TGT and west of the TGT that provides near-continuous bathymetric coverage for ∼40 km north of the present-day ice-shelf margin (Figs. 2, 3). Data gaps remain in the western part of the area in front of the ice mélange, east of Crosson Ice Shelf, as icebergs and bergy bits released from the ice mélange persistently covered this region during NBP19-02. In addition, perennial fast ice and huge icebergs, such as the 80 km long and 45 km wide iceberg B-22A, that calved periodically from the TGT and then moved slowly north- and northwest, thereby episodically running aground, have generally prevented survey between TGT and Crosson Ice Shelf (see area with B-22A outlines in Fig. 1).

The sea-floor bathymetry offshore Thwaites Ice Shelf is dominated by an elongate depression oriented NNE–SSW (Thwaites Trough) and a series of topographic highs (H1-H3) along its southern margin (Fig. 3a). The depression is characterised by water depths of 1100–1250 m, which is 200 to 400 m deeper than the sea floor on its flanks; it typically has a relatively flat or gently inclined floor (Fig. 3a, c). Although the depression appears to be continuous for at least 75 km, and connects with areas of deep (>1300 m) sea floor directly north of the EIS and in Pine Island Trough further east (Fig. 2a), its width varies significantly along its length and the flanks are discontinuous in form. As a result, we do not define this feature as a channel, which implies incision by water flow, but rather as a small trough. The topographic highs that make up the southern flank of the trough (H1-H3) decrease in height (above the surrounding sea floor) from NNE to SSW and appear to form a broad (>15 km) but discontinuous ridge, also oriented NNE–SSW. The large, discontinuous ridge is the extension of a bedrock ridge in Pine Island Trough to the east that has the same orientation (Figs. 2, 3a; Nitsche et al., 2013). The most prominent high (H1) occurs immediately north of the EIS with a shallowest recorded water depth of just 82 m and extends at least 40 km in a north–south direction. The northern part of the EIS is pinned, about 45 km downstream of the grounding zone (Rignot et al., 2001; Tinto and Bell, 2011; MacGregor et al., 2012), on the southern part of this high, forming an ice rumple in the EIS (Matsuoka et al., 2015). The shallowest water depths on the two highs H2 and H3 (from NNE to SSW) are 362 and 611 m, respectively. Indeed, a pinning point for TGT on H2 has been identified from remote-sensing data (Rignot et al., 2011), but the TGT must have fully retreated from that point prior to NBP19-02. North of the trough and ridge (north of $\mathrm{74}{}^{\circ }{\mathrm{30}}^{\prime }$ S) is an area of rugged morphology characterised by shallow sea-floor ridges and deep basins (>1400 m); this area merges with similar terrain in Pine Island Trough described by Nitsche et al. (2013) as their “area 2” (Fig. 2). East of the EIS, the eastern flank of H1 is just exposed but the bathymetry is dominated by a deep (1000–1200 m), rugged area of sea floor, bounded on its eastern edge by a bedrock high in Pine Island Bay that continues southeastwards towards the grounding zone (Fig. 2a). The deepest part of this area appears to form a poorly defined bathymetric trough oriented NNE–SSW (T4; Fig. 3a); oceanographic measurements and models confirm that this deep acts as a pathway for CDW towards grounding lines in western Pine Island Bay (Fig. 2a) (Dutrieux et al., 2014; Nakayama et al., 2019).

3.1 Glacial landforms

The large-scale morphology of the sea-floor offshore Thwaites Ice Shelf is overprinted by linear features oriented subparallel to the trough and ridge. Streamlined subglacial landforms occur in areas where bedrock crops out at the sea floor, either on topographic highs or interrupting the smooth trough floor or rugged basin floors (Fig. 3a, b). These features are identified as crag and tails by their tapering form and are 750–4000 m long, 200–2000 m wide and 25–250 m high (Figs. 3, 4, S1a in the Supplement), as well as their morphological similarity to crag and tails from other deglaciated terrains (e.g. Dowdeswell et al., 2016; Maclean et al., 2016; Nitsche et al., 2016). The wider southern ends of the crag and tails are rugged and have significant relief, suggesting bedrock composition, whereas the northern ends are smooth and elongate, indicative of a sedimentary composition. These landforms are produced subglacially as glacier ice flows across bedrock obstacles producing the characteristic morphology through erosion and deposition (Benn and Evans, 2010; Nitsche et al., 2016). The tapering of these features north-northeastwards indicates palaeo-ice flow of TG in this direction towards Pine Island Trough (Fig. 3b); east of the EIS the orientation of these landforms varies from N–S to NW–SE, depicting flow around the H1 high (Fig. 3b). Curved or semicircular moats – crescentic scours (Lowe and Anderson, 2003; Graham et al., 2009; Graham and Hogan, 2016) – occur around the southern (upstream) sides of some crag and tails (Figs. 3b, 4a, S1d) and are suggestive of erosion by meltwater or, alternatively, a till slurry or mobile basal ice, upstream of bedrock obstacles (see Graham and Hogan, 2016). On the H2 and H3 highs, subtle elongate ridges separated by linear grooves also occur (Fig. 3a). These glacial lineations have parallel sides; lack a pronounced and wider end; and are 1000–2500 m long, 100–200 m wide and 2–10 m high; their crest-to-crest spacing is typically 200–500 m (e.g. Fig. S1b). These features, which were produced subglacially, have relatively short elongation ratios ($<\mathrm{10}:\mathrm{1}$) and are oriented parallel to crag and tails also on the tops of the highs but are slightly oblique to crag and tails in the trough (Fig. 3). A distinct but discontinuous scarp is mapped on the top of the H2 high, upstream from its frontal edge (Figs. 3, 4b, d). The scarp has a curved planform and a steep (3–4^∘) northern and gentle (0.5^∘) southern slope; glacial lineations occur on the gentle back-slope of this feature which extends for about 15 km in a SW–NE direction. The asymmetric geometry and lineated back-slope of this landform identify it as a grounding-zone wedge (GZW; Fig. S1c) (Alley et al., 1989; Larter and Vanneste, 1995; Dowdeswell and Fugelli 2012; Jakobsson et al., 2012), i.e. a wedge of sedimentary material that built up at the grounding zone when it was stable for a time on the H2 high. This feature is interrupted, however, by a ∼40 m deep groove or small channel with bedrock ridges visible on either side, suggesting that the sedimentary wedge is not thick enough to fully bury the underlying topography; i.e. it is only tens of metres thick. Discrete 6–10 m deep linear to curvilinear furrows with small berms (4–6 m) were mapped on H3 and north of the GZW front scarp (Fig. 4a, b); these are interpreted as iceberg plough marks.

Sea-floor highs in the area are also variously gouged and streamlined, resulting in a pattern of grooves and bedrock ridges (Fig. S1e). Grooves on the H1–H3 ridges occur on their southern (upstream) parts and exhibit a range of orientations that probably relate to the structure of the underlying bedrock exploited by glacial erosion. The grooves are typically <20 m deep, <200 m wide and <6000 m long. The surfaces of sea-floor highs north of the H1–H3 ridge have a more streamlined appearance resulting from shallow, semi-parallel shallow grooves that are preferentially aligned with the crag and tails (Figs. 2, 3). Streamlined bedrock highs of this type are typical of inner-shelf morphologies around Antarctica (e.g. Wellner et al., 2006; Livingstone et al., 2013) including in the adjacent eastern part of Pine Island Bay (Nitsche et al., 2013). Taken together, the orientations of the streamlined subglacial landforms (crag and tails, glacial lineations, bedrock grooves/ridges) define the former ice-flow directions of an expanded TG (depicted by white arrows in Fig. 3b), confirming that this area is the former bed for grounded glacier ice.

3.2 Trough and channel metrics

Multiple bathymetric troughs and bedrock channels were mapped and analysed in the area beyond Thwaites Ice Shelf as part of this study (Fig. 3a). Troughs and channels in the adjacent part of Pine Island Bay have been described comprehensively by Nitsche et al. (2013) and Kirkham et al. (2019). The larger troughs in our study area, which have also been identified on gravity-derived regional bathymetry maps (Millan et al., 2017; Jordan et al., 2020), are considered important as potential pathways for the transport of CDW towards the grounding zone of TG and warrant full description here. We distinguish these comparatively larger troughs, based on their size, connectivity and variable flank form (as described above), from the channels which are notably smaller in scale (Fig. S2). The latter have continuous, parallel-sided flanks and undulating thalwegs and incise into rugged sea-floor areas interpreted as bedrock (e.g. Fig. S1e). It is widely accepted that channels of the latter type were eroded by pressurised subglacial water flow (Lowe and Anderson, 2003; Nitsche et al., 2013; Kirkham et al., 2019), whereas the troughs likely relate, at least in part, to underlying geological structures such as dykes and tectonic deformations (see Gohl et al., 2013) that have been variously modified by ice, as well as possibly by subglacial water flow.

Cross-sectional analyses of the troughs (n=166) reveal their large scale, with average widths, depths and cross-sectional areas being 2090 m, 90 m and 144 000 m², respectively (Fig. S2). The troughs are typically 10–30 times as wide as they are deep, although we note that there is a significant size difference between the main NE–SW trough (Thwaites Trough) and the remaining troughs (Figs. 3a, S2). By comparison, the bedrock channels (n=822) are on average 520 m wide, 50 m deep with a cross-sectional area of 18 000 m². The channels are generally 5–10 times wider than they are deep. The derived b values, which characterise cross-sectional shape (Pattyn and Van Huele, 1998), suggest that the bedrock channels are between v and u shaped, whereas the larger troughs have no dominant cross-section shape (Fig. S2). In general, the trough floors are flat or inclined in cross profiles (Fig. 3c) and are gently undulating in along-trough profiles.

The main Thwaites Trough is oriented NNE–SSW, which is oblique to the northerly palaeo-ice-flow directions immediately in front of Thwaites Ice Shelf (Fig. 3b). This indicates that, at the time that the subglacial landforms were produced, ice was thick enough not to be fully steered by even major elements of the bed topography. The two southernmost troughs that we have analysed (T2 and T3 in Fig. 3a) are oriented perpendicular to Thwaites Trough (i.e. NNW–SSE), and the troughs north and east of the EIS are generally aligned with palaeo-ice-flow directions (Figs. 2, 3a). Note that the T1 trough is not well covered by our MBES dataset and is not discussed in detail here. The T2 and T3 troughs, whose floors have water depths of 800–900 m, separate the H1–H3 bathymetric highs and are of interest as potential pathways for CDW to the TG grounding zone. Long profiles from the T2–T4 troughs (Fig. S3) identify sill depths along the pathways of the troughs that may act as important pinch points for ocean circulation, in particular if they impede CDW inflow towards the Thwaites grounding zone (cf. De Rydt et al., 2014). T2 has a smooth long profile with a prominent sill at 710 m depth at about $\mathrm{107}{}^{\circ }{\mathrm{3}}^{\prime }$ W, $\mathrm{75}{}^{\circ }{\mathrm{3.6}}^{\prime }$ S, whereas T3 has a rugged profile with three bathymetric sills in its northern (ice distal) part with depths of 750–760 m, as well as several other sills further south (ice proximal) around 780 m water depth (Fig. S3a). The bathymetry around T4, east of the EIS, is generally deeper (>1000 m) than most of T2 and T3, so the main constriction on this trough seems to be where it crosses the NE–SW-trending sea-floor ridge in Pine Island Bay (Fig. 2a). At this location, around $\mathrm{105}{}^{\circ }{\mathrm{24.4}}^{\prime }$ W, $\mathrm{74}{}^{\circ }{\mathrm{35.4}}^{\prime }$ S, there is a sill at 880 m depth (Fig. S3d). Channel widths at these locations are 5000 m for the T3 sills and 2500 m for the T4 sill (although this is not the only interruption in the ridge in Pine Island Bay). Widths were measured at 500 m depth as this is taken to be a reasonable top-CDW depth for the area (based on oceanographic measurements; Bastien Queste, personal communication, 2020). The bathymetry of the highs west of T2 is >500 m depth, meaning that CDW could effectively flood over this topography rather than be constrained to the trough; however, if it was topographically routed (see Nakayama et al., 2019) through T2, then the channel width at the sill is 4700 m (at 640 m water depth).

3.3 Bathymetric highs and ridges

Owing to the importance of sea-floor highs in front of the Thwaites Ice Shelf as barriers to CDW inflow, and as former ice shelf/sheet pinning points, we examine the morphology of the discontinuous NNE–SSW-trending ridge in detail (Fig. 4). The ridge comprises the H1–H3 highs separated by the two troughs described above (T2 and T3; Fig. 3a). The width of the ridge varies significantly, from 6 km in the SW of the study area to at least 40 km over H1, although we acknowledge that data coverage is limited. In some places, the bathymetric highs are strikingly flat topped. These planar features are accentuated in maps of the first derivative of bathymetry, slope, which reveals both low slopes (<2^∘) (Fig. 4c, d) and low roughness over H2, H3 and the western part of H1. The continuation of the ridge further north into Pine Island Trough has a similar surface expression but is generally narrower (Fig. 3a). These areas with low surface slopes are atypical when compared with other bathymetric highs in the area, which have rugged surface morphologies characterised by bedrock grooves and channels (Figs. 2, 4d) (Nitsche et al., 2013; Arndt et al., 2018; Kirkham et al., 2019). Instead, the low slope values are similar to those derived for the base of the troughs in front of Thwaites Ice Shelf (Fig. 4a) and the sediment-filled basins just seaward of the Pine Island Ice Shelf front (Nitsche et al., 2013; Kuhn et al., 2017). At least two distinct levels of flat-topped surfaces occur at 400 and 640 m water depth (Fig. 5). We suggest that this morphology was generated as the highs were overridden and eroded by a formerly expanded TG and Thwaites Ice Shelf (with the necessary ice thickness to reach the depth of the flat-topped surfaces). A value of ∼400–500 m is similar to ice-shelf thicknesses for Thwaites Ice Shelf today (Griggs and Bamber, 2011; Jordan et al., 2020), and a prevalence of flat-topped highs at this depth may, therefore, support recent modification of the sea-floor highs at TG. In contrast, the deeper flat top of H3 (640 m depth) was probably formed at an earlier stage, as was the flat top of the high in Pine Island Trough, as that area is known to have been ice-free (sheet and shelf) for at least the last 10 kyr (e.g. Kirshner et al., 2012; Hillenbrand et al., 2013). The interpretation of erosion or planing off by an ice shelf is supported by the occurrence of glacial lineations on the tops of the highs (Fig. 3b), which are in line with modern ice-velocity vectors (Mouginot et al., 2019) but oblique to the orientation of crag and tails in the troughs, thus indicating a change in flow direction from grounded ice flow to ice-shelf flow over the high. A similar interpretation was made for the lineated surface of a former pinning point of the Pine Island Ice Shelf (PIIS) that has been recently exposed by ice-shelf calving events (Arndt et al., 2018), although that feature was not planed off to form a flat-topped high but rather has a stepped and rugged surface morphology albeit with some gently sloping parts (see their Fig. S3). We note, however, that alternative explanations are possible for this morphology, namely that the flat tops are an inherited feature produced by erosion down to horizontal bedrock strata or that rugged bedrock highs, which are typical of the inner Amundsen Sea shelf (see Nitsche et al., 2013), were mantled by some thickness of glacigenic material that levelled the topography below. The former is relatively easy to discount accepting that the inner shelf of the ASE is composed of crystallite basement with seismic-reflection profiles showing that northward-dipping sedimentary strata only occur on the middle and outer shelf (see Graham et al., 2009; Gohl et al., 2013). In this setting close to the current TG grounding zone, it is perhaps easier to conceive of the latter explanation that rugged bedrock features were mantled by glacigenic material delivered to the area when the grounding zone was located on or near the highs and then flattened by some degree of glacial compaction and/or erosion as it was overtopped by TG and the subsequent Thwaites Ice Shelf. This is consistent with our suggestion for the formation of these flat tops as we cannot tell from our data either what sediment thickness occurs on the highs or how much erosion took place, and we acknowledge that the amount of ice-shelf erosion may have been small, only skimming unconsolidated material from the surface of the highs. The presence of GZWs and glacial lineations on the highs, and sub-bottom profiler data (Fig. S4), confirms that at least some thickness of unconsolidated material occurs on the highs, but seismic-reflection profiles would be required to fully capture the internal structure of these features.

Figure 5 Flat-topped highs in the Amundsen Sea. (a) Cross-sectional profiles over H1–H3 highs at Thwaites Glacier (for location of profiles, see Fig. 4c, d), further offshore in Pine Island Bay (PIB; location of high in Fig. 2a, location of profile in Fig. 5b) and offshore from the Getz A Ice Shelf (location of high in Fig. 1; location of profile in Fig. 5c), showing sea-floor highs planed off at different depth levels. The flat portions of the profiles are marked with grey bars and the depth elevation for that flat top given above, so “H1:400 m” means the flat part of the profile over H1 high is at 400 m water depth. (b) MBES of flat-topped highs part of the discontinuous sea-floor ridge in PIB. (c) MBES of a flat-topped high with glacial lineations (white lines) just in front of the Getz A Ice Shelf (following Nitsche et al., 2016). White arrows show direction of past ice flow based on streamlined subglacial landforms.

No matter the exact formation mechanism, the flat-topped morphology of the highs in our study area is striking and notably rare for sea-floor highs around Antarctica. We note a similar but less pronounced terrain over some highs along the same structural ridge in Pine Island Bay, i.e. east of the EIS (Fig. 5b), and a solitary very flat topped high, with comparable dimensions to those offshore TG, is visible <5 km north of the Getz A Ice Shelf (Fig. 5c). Beyond these rare examples, the best analogy for this morphology probably comes from a set of “iceberg terraces” on terminal moraines at the mouth of a Svalbard fjord, which display remarkably flat topped surfaces at several bathymetric levels. These are interpreted to have formed as tabular, flat-based icebergs overtopped and eroded morainal sediments (Noormets et al., 2016). It should be noted, however, that sediments of this morainal bank complex probably consist of unconsolidated material that has not been overridden (or compacted) by grounded ice, meaning that they are likely more erodible than basement highs in front of Thwaites Ice Shelf. Nevertheless, the flat-topped morphology is suggestive of a sedimentary cap on the pinning points at TG, and the fact that similar features exist in Pine Island Bay and beyond the Getz A Ice Shelf may indicate that pinning points with such sedimentary caps are widespread on the inner shelf in the Amundsen Sea.

New geomorphic information is also revealed by the flanks of the highs (Fig. 4a, b). The northern (ice distal) flanks of the H2 and H3 highs are characterised by subtle downslope-trending gullies that transition into a smooth but inclined sea floor in the troughs at a distinct break of slope. There are also a few, discrete semicircular indentations in the scarp between the surface of the highs and their flanks (Fig. 4a). The gullies have simple non-branching geometries, have small dimensions (widths 150–700 m, depths 5–50 m, lengths <2 km), and typically define broad u shapes in cross section, although some v-shaped forms are present. In addition, their form is consistent with other Antarctic submarine gully systems (e.g. Fig. S1f; Gales et al., 2013; Post et al., 2020). Thus, we interpret the gullies as the result of the downslope mass movement of material from the tops and sides of the H2 and H3 highs via gravitational processes into the small sediment fans at the base of the slope (Fig. 4a). The semicircular indentations may be the headwalls of small slide scars (cf. Noormets et al., 2009; Gales et al., 2013). One 6 km by 2.6 km segment of the H3 high is somewhat detached from other parts of the ridge and appears particularly fragmented on its flanks (Fig. 4a). About 1500 m south of this, on the main H3 high, is a distinct break in slope with the same planform shape as the back of the detached segment that we refer to as a “block” (bl in Fig. 4a, c). Similarly, 2 km west of this block is another somewhat isolated 3 km by 1.8 km block of H3 that is incised by gullies on its northern front and lacks lineations on its surface (Fig. 4a). We consider several interpretations for these features. First, it is possible that they are detached blocks of the H3 high that had been displaced downslope over a short distance (black arrows in Fig. 4a, c), remaining largely intact, but subsequently affected by some gravitational collapse of their flanks. Slide megablocks with similar dimensions (or larger), non-crystalline compositions and degraded flanks are known from, for example, the Hinlopen Slide scar on the northern Barents Sea margin (Vanneste et al., 2006; Hogan et al., 2013). The second possibility is that the blocks are small bedrock highs that have been variously mantled by and surrounded by glacigenic sediment, either deposited subglacially on their surfaces or as it was transported downslope by gravity-driven processes (towards the sediment fans) from the H3 high (see blue arrows in Fig. 4a). If the latter case is true, then the pronounced semicircular indentation on the eastern block (labelled in Fig. 4a) may be an erosional scour mark formed by the recurrent flow of water, as is evidenced by the small channel on the eastern side of the block (white arrows in Fig. 4c).

4 Methods II: bed roughness and basal drag

4.1 Spectral analysis of bed roughness

The idea that the drag experienced by a glacier can be analysed by treating the basal roughness as a superposition of sine waves of different wavelengths (λ) has a long history in glaciology (e.g. Nye, 1970; Kamb, 1970; Hubbard et al., 2000). In these approaches, Fourier methods are used to find the drag contributed by roughness that falls within some range of wavelengths (or band of spatial frequencies). More recently, a similar approach has treated the longer wavelength undulations that can affect flow near the upper surface of the ice (Schoof, 2002). Knowledge about the power spectrum of the basal roughness is an essential input to all of these studies.

The high-resolution bathymetric grid that we have produced (Fig. 2) represents the former beds of expanded TG and PIG. As such, it contains information about bed geometry and roughness that can be used to investigate the behaviour of glacier ice over this terrain, specifically how roughness due to bedrock topography and subglacial landforms at various wavelengths, λ, might generate form drag within the ice (Schoof, 2002). As an example, Fig. 6a shows a bed elevation profile from Pine Island Bay (profile 6 in Fig. S5a) sampled at 25 m intervals along a flow line. To analyse the variance of roughness at each wavelength scale we computed the power spectrum of the bathymetric topography (e.g. Fig. 6b). It is conventional to plot the power spectrum as a function of spatial frequency, defined as $f={\mathit{\lambda }}^{-\mathrm{1}}$, rather than wavelength. This function, P(f), referred to as the periodogram, shows how the variance of roughness is distributed among different frequency intervals. Thus, the variance attributable to roughness within any frequency interval is the integral of the periodogram over that interval. Figure 6b shows the one-sided periodogram P(f_n), evaluated at equally spaced frequencies $f_n=n/a$, where a is the length of a moving window. The periodogram was obtained in MATLAB R2017a using Welch's method (Welch, 1967), with a window of length a=6.4 km and 50 % overlap between consecutive windows along the bed profile. Within each 6.4 km window we removed a linear trend from the bed profile and applied a Hamming window before computing the one-sided power spectrum using the absolute square of the fast Fourier transform, appropriately scaled. We then averaged spectra from the multiple windows along the bed profile to provide results for the flow line as a whole.

Figure 6 (a) Bed elevation (equal to water depth for sea-floor data) versus distance along flow line. (b) Roughness power spectra P_n versus spatial frequency f_n. (c) Scaled basal drag contributions ${\mathit{\beta }}_n/{\mathit{\beta }}^{\ast }$ versus f_n for along-flow profile (6) offshore from Pine Island and Thwaites glaciers (for location see Fig. S5a). The red lines in (b) and (c) are based on an assumption of Brown-noise power spectrum that falls off as the inverse square of spatial frequency. At low spatial frequencies, drag contributions depend on the function F (see Sect. 4.2) with the two limiting cases shown: F₁ (blue), F₂ (black).

Periodograms were computed for a total of 14 palaeo- and modern flow lines for TG and the Dotson–Getz palaeo-ice stream (locations in Figs. S5, S6). For TG and PIG, six palaeo-flow lines were picked manually across the bathymetric grid by tracing lines parallel to subglacial landforms that indicate palaeo-ice-flow directions (e.g. Fig. 3b, profiles 1–6, Fig. S5a). We performed calculations for bed profiles from four additional areas for comparison: (i) modern along-flow bed profiles for TG (profiles 8–9, Fig. S6a), (ii) modern bed profiles for PIG (profiles 10–11, Fig. S6b), (iii) an area of smooth bed topography on the middle continental shelf in the Dotson–Getz palaeo-ice-stream trough (profile 7, Fig. S5b), and (iv) six across-flow profiles from the bathymetric datasets (profiles a–f, Fig. S5a, b). Onshore bed profiles were extracted from the AGASEA (Holt et al., 2006) and Operation IceBridge (OIB; Tinto et al., 2010, updated 2019) airborne radar datasets for TG and PIG beds, respectively. Profiles were selected based on their location, along the central glacier trunk, and their quality in terms of continuity and fewer outliers. The profiles from the Dotson–Getz Trough, offshore from the Getz A Ice Shelf (Fig. S5b), were selected as representative of a sedimentary palaeo-ice-stream bed characterised by mega-scale glacial lineations (MSGL) (Graham et al., 2009; Spagnolo et al., 2014). These were extracted from a MBES dataset fully described by Larter et al. (2009) and Graham et al. (2009).

The power spectrum of natural terrain is often approximated as a power law in frequency (e.g. Jordan et al., 2017). The results of our spectral analysis (Figs. 6 and 7) show that a good approximation can be obtained using an inverse-square power law,

$$\begin{array}{}\text{(1)}& P=A{f}^{-\mathrm{2}},\end{array}$$

where the constant A has units of length. This is the periodogram expected for a Brown-noise or random-walk elevation profile, as produced by taking an uncorrelated random step in the vertical direction for each unit step in the horizontal.

Figure 7 Selection of bed profiles (top), derived power spectra (middle) and basal drag contributions (bottom) for (a, b) along-flow profiles (1) and (3) offshore Thwaites Glacier, (c) onshore along-flow bed profile (8) for Thwaites Glacier, (d) along-flow profile (7) for Dotson–Getz Trough, and (e) across-flow profile “c” offshore Thwaites Glacier. Profile locations are shown in Figs. S5 and S6. The red line in power spectra and drag contribution plots are based on an assumption of Brownian motion (i.e. power decays as inverse square of spatial frequency). At low spatial frequency, drag contributions depend on the function F. Two limiting cases are shown: F₁ (blue) and F₂ (black).

For Brown noise, the parameter A represents the roughness variance per unit length of profile. If we consider a section of profile having length l and restrict our definition of roughness to all wavelengths λ<l, we expect this roughness to have variance obtained by integrating the periodogram P over frequencies $f>l^{-\mathrm{1}}$. For the Brown-noise periodogram (Eq. 1), this integration provides

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}^{\mathrm{2}}=Al.\end{array}$$

Thus, for Brown noise, the variance grows in proportion to the section length considered, and the rms roughness $\mathit{\sigma }=\sqrt{Al}$ grows with the square root of the section length (e.g. Jordan et al. 2017). If we take a longer section of profile, we expect to see larger rms roughness within it. This makes it clear that, at least for Brown noise, the roughness can only be characterised by its rms value if reference is also made to the length scale under consideration. Next, we examine the consequences of roughness for the drag that resists the sliding of a glacier. We pay particular attention to the case when the periodogram follows an inverse-square law (Eq. 1) that is appropriate for a Brown-noise power spectrum.

4.2 Relating bed topography to basal drag

In this section, we will use the power spectra of the high-resolution bathymetry, together with a theoretical expression for form drag (Schoof, 2002, Eq. 67), to assess how roughness at different scales affects the drag that opposes glacier sliding. This allows us to consider the contribution that the observed sea-floor bathymetric roughness at short wavelengths would make to form drag when covered by flowing ice.

The theory of ice flow over an undulating bed (Schoof, 2002) provides an approximate expression for the form drag τ, expressed as a basal shear stress that acts to resist sliding. Here, we use dimensional quantities rather than the non-dimensional quantities (as used by Schoof, 2002). In our notation, Schoof's (2002, Eq. 67) expression for form drag becomes

$$\begin{array}{}\text{(3)}& \mathit{\tau }\mathit{\beta }=U.\end{array}$$

In this expression, β is the drag coefficient and U is the speed of ice flow averaged over some horizontal length scale significantly larger than the ice thickness. Here, we choose this length as a=6.4 km, the window length used in our spectral analysis. According to Schoof (2002), the drag coefficient $\mathit{\beta }=\sum _{n=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\beta }}_n$ is the sum of contributions, β_n. Each contribution β_n is caused by roughness that falls within a frequency band of width $\mathrm{1}/a$ centred at frequency, $f_n=n/a$. The spatial wavenumbers corresponding to these frequencies are defined as k_n=2πf_n.

Schoof (2002) provides an expression for the drag contributed by roughness within each frequency band. In our notation, this translates to

$$\begin{array}{}\text{(4)}& {\mathit{\beta }}_n=\mathrm{4}{\mathit{\beta }}^{\ast }{\left(k_n/k^{\ast }\right)}^{\mathrm{3}}F\left(k_n/k^{\ast }\right)\left({\left|{\widehat{h}}_n\right|}^{\mathrm{2}}/H^{\mathrm{2}}\right).\end{array}$$

In this expression, ${\widehat{h}}_n$ is the Fourier component of the bed roughness at spatial frequency, f_n. We estimate ${\left|{\widehat{h}}_n\right|}^{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}{P}_n{a}^{-\mathrm{1}}$ as appropriate for the Fourier series defined by Schoof (2002). The one-sided periodogram P_n can be obtained either directly from the bathymetric observations as described above (Sect. 4.1) or by fitting the inverse-square power law described by Eq. (1) to those periodograms (see Figs. 6 and 7 for examples). The scaling constants are ${\mathit{\beta }}^{\ast }=\mathit{\eta }/H$ and $k^{\ast }=\mathrm{1}/H$, where viscosity η and thickness H are representative values averaged over the length scale a (Schoof, 2002).

We consider two limiting cases for the function F:

$$\begin{array}{}\text{(5)}& {\displaystyle }{F}_{\mathrm{1}}\left(k_n/k^{\ast }\right)={\displaystyle \frac{{\sinh }^{\mathrm{2}}\left(k_n/k^{\ast }\right)-{\left(k_n/k^{\ast }\right)}^{\mathrm{2}}}{\left(k_n/k^{\ast }\right)+\cosh \left(k_n/k^{\ast }\right)\sinh \left(k_n/k^{\ast }\right)}},\text{(6)}& {\displaystyle }{F}_{\mathrm{2}}\left(k_n/k^{\ast }\right)={\displaystyle \frac{\left(k_n/k^{\ast }\right)+\sinh \left(k_n/k^{\ast }\right)\cosh \left(k_n/k^{\ast }\right)}{{\sinh }^{\mathrm{2}}\left(k_n/k^{\ast }\right)}}.\end{array}$$

These are derived by Schoof (2002) as his Eqs. (65) and (66) for small and large bed roughness, respectively. We use dimensional quantities, so $k_n/k^{\ast }{\widehat{h}}_n/H$ and $F\left(k_n/k^{\ast }\right)$ in our notation equate respectively to non-dimensional quantities k_n, $\mathit{\nu }{\widehat{h}}_n$ and f(k_n) in Schoof (2002). For sufficiently small wavelengths $k_n\gg k^{\ast }$, the functions $F_{\mathrm{1}}\left(k_n/k^{\ast }\right)$ and $F_{\mathrm{2}}\left(k_n/k^{\ast }\right)$ both tend to unity. In this case, the contribution to form drag becomes insensitive to ice thickness and to the choice of function used (Schoof, 2002).

If the bed follows a Brown-noise power spectrum, we can use the inverse-square law $P_n=A{f}_n^{-\mathrm{2}}$. Under those circumstances, the drag contribution β_n will grow approximately linearly with frequency at sufficiently high wavenumbers $k_n\gg k^{\ast }$:

$$\begin{array}{}\text{(7)}& {\mathit{\beta }}_n=\mathrm{16}\mathit{\eta }{\mathit{\pi }}^{\mathrm{3}}A{a}^{-\mathrm{1}}{f}_n.\end{array}$$

For the Brown-noise power spectrum, the amplitude of roughness decreases at shorter wavelengths. Nevertheless, Eq. (7) shows that those short-wavelength scales will still be more effective at causing form drag than the longer wavelengths, despite their smaller amplitude. This is because the factor ${\left(k_n/k^{\ast }\right)}^{\mathrm{3}}$ increases faster than the inverse-square law decreases.

As a consequence of the linear increase in the drag contribution with frequency, the total drag would become unbounded, and the sliding would stop, unless the bed of the glacier departs from the Brown-noise assumption and becomes smooth at scales smaller than some wavelength (Nye, 1970). This wavelength, λ_N, provides an upper bound to the spatial frequencies that cause drag, $f_N={\mathit{\lambda }}_N^{-\mathrm{1}}=N/a$. Under this assumption, the total drag can be approximated by truncating the sum:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \mathit{\beta }=\sum _{n=\mathrm{1}}^N{\mathit{\beta }}_n=\mathrm{16}\mathit{\eta }{\mathit{\pi }}^{\mathrm{3}}A{a}^{-\mathrm{2}}\sum _{n=\mathrm{1}}^{N}n=\\ & \mathrm{8}\mathit{\eta }{\mathit{\pi }}^{\mathrm{3}}A{a}^{-\mathrm{2}}N\left(N+\mathrm{1}\right).\end{array}\end{array}$$

When λ_N≪a, so that N≫1, this gives the approximation

$$\begin{array}{}\text{(9)}& \mathit{\beta }=\mathrm{8}\mathit{\eta }{\mathit{\pi }}^{\mathrm{3}}A{\mathit{\lambda }}_N^{-\mathrm{2}}.\end{array}$$

Therefore, if the Brown-noise inverse-square law power spectra applies down to the finest wavelength, λ_N, the drag will be determined by that scale, along with the viscosity η, and the coefficient of the power law A that can be recovered from the periodogram.

It is common in sliding theories to define the slip length, $L=\mathit{\eta }/\mathit{\beta }$. The slip length L is an important quantity that allows us to make a distinction between two regimes of ice flow. When the slip length is much larger than the ice thickness H, the drag is too small to induce significant shearing within the ice column, and the ice can be considered to slide over the base as a plug flow having uniform velocity with depth. This is the situation modelled for slippery-based ice streams by MacAyeal (1989). By contrast, when the slip length is much smaller than the ice thickness, the drag is able to induce a substantial amount of shearing through the ice column, so the flow velocity varies significantly with depth.

Using Eq. (9), we obtain the following expression for the slip length under the assumption of a Brown-noise power spectrum, truncated at some frequency, $f_N={\mathit{\lambda }}_N^{-\mathrm{1}}$:

$$\begin{array}{}\text{(10)}& L={\mathit{\lambda }}_N^{\mathrm{2}}/\left(\mathrm{8}{\mathit{\pi }}^{\mathrm{3}}A\right).\end{array}$$

One consequence of this is that if we wish to infer the amount of form drag using Eq. (9), or the slip length using Eq. (10), it is not enough to evaluate the roughness parameter A alone. We must also establish λ_N, the smallest wavelength that is effective at causing drag.

5 Results II: assessing roughness and drag contributions for palaeo- and modern glacier beds

5.1 Bed roughnesses

The power spectra of selected bathymetric profiles are shown in Figs. 6 and 7. Figure 6b shows that the one-sided periodogram P(f_n), computed using the bathymetric profile shown in Fig. 6a (location in Fig. S5a), has no strong peaks at any particular preferred scales of roughness. Instead, the periodogram decreases continuously as spatial frequency increases. This decrease approximately follows the inverse-square power law appropriate for Brown noise, so that the periodogram can be approximated as $P_n=P\left(f_n\right)=A{f}_n^{-\mathrm{2}}$. The red line in Fig. 6b shows this power law, with a value A=0.1 m. The spectra are remarkably consistent across many of the profiles considered (Figs. 7, S7), with the exception of the smoother MSGL area (Figs. 7d, S5b). There, the Brown-noise inverse-square law can still provide a good approximation to the periodogram, but the value of A=0.001 m that is required to provide a good match to the observations is much smaller (Figs. 7, S7p, q).

The power-law approximation to the power spectrum also agrees closely with power spectra of profiles of bed elevation from airborne radar flown over TG (Fig. 7c), so the MBES data provide a good analogue to the subglacial undulations that control the sliding of TG today. Despite improvements in the methodology of high-resolution radar surveys of the active subglacial bed (King et al., 2016; Bingham et al., 2017), the MBES data provide a more detailed view of the shorter spatial scales than airborne or ground-based radar. Comparisons with previous studies of subglacial roughness are not straightforward because individual studies have used different window lengths to investigate roughness (e.g. Jordan et al., 2017; Falcini et al., 2018). However, for the Brown-noise power spectra in Eq. (1), we expect the rms roughness for a section of length l to be $\mathit{\sigma }=\sqrt{Al}$. Most of our spectra are close to the Brown-noise spectra with A =0.1 m (Figs. 6, 7), and for this value, we would expect window lengths from 80 to 1000 m to produce rms roughness estimates in the range 2.8 to 10 m. These values are similar to those reported previously for glaciated terrains (Jordan et al., 2017; Falcini et al., 2018).

5.2 Basal drag contributions

For the example bed profile in Fig. 6a, values of ${\mathit{\beta }}_n/{\mathit{\beta }}^{\ast }$ obtained using the periodogram shown in Fig. 6b and Eq. (4) are plotted against spatial frequency f_n in Fig. 6c. The linear dependence predicted by Eq. (7) for the Brown-noise approximation is shown as a red line in Fig. 6c. In these plots we used a value A=0.1 m and a representative ice thickness H=1 km.

The expression for the basal slip length (Eq. 10) lets us use the bathymetry to make a dynamical distinction between regions of fast sliding with little internal deformation L>H and regions of slow sliding and shearing flow L<H. Using Eq. (10), the ratio of slip length to ice thickness is $L/H={\mathit{\lambda }}_N^{\mathrm{2}}/\left(\mathrm{8}{\mathit{\pi }}^{\mathrm{3}}AH\right)$. For a value of A=0.1 m and a typical ice thickness scale of H=1 km, this suggests that features on scales smaller than λ_N=150 m would provide sufficient drag to induce significant vertical shearing within the ice. Since features on this scale are well resolved by the bathymetric profiles (e.g. Figs. 3, 4) and fall within the range of frequencies where the inverse-square power law applies, we conclude that the form drag produced by the observed subglacial roughness would have produced significant shearing within the flow of the grounded ice as it retreated over the highs and ridges surveyed by the MBES. This suggests that it is important to include the effects of form drag caused by basal roughness over such terrain, and by extension over the extant parts of TG today.

A distinction must be made for the region of MSGLs on the Dotson–Getz palaeo-ice-stream bed (Fig. S5b). Here, the elevation profile is exceptionally smooth. The spectral analysis confirms this (Figs. 7d, S7p, q), and the coefficient A that best fits the observations is some 2 orders of magnitude below the more generally applicable value of A=0.1 m. Repeating the above analysis with A=0.001 m shows that the power law would have to apply down to horizontal scales smaller than λ_N=15 m. Features on this scale are not well resolved in our bed profiles (e.g. the MBES grids have cell sizes of 50 m). This means that, in contrast to the more general case, it remains possible that the MSGL terrain is so smooth that the resulting form drag produced little vertical shearing within ice that flowed over it, making the ice dynamics of this area more akin to the flow described for slippery-based ice streams by MacAyeal (1989). This result is consistent with our understanding of how MSGLs form, i.e. via the self-organisation of deforming sediment at the bed under fast-flowing ice (e.g. Spagnolo et al., 2014). We also repeated the analysis in the direction perpendicular to elongated features (Figs. 7d, S7k–o, q). There is no evidence that ice flowed in this direction, but the theory can nevertheless compute the contributions to form drag that would arise in that hypothetical situation. For most of the across-flow lines (Fig. S7), the power spectra are similar to the along-flow direction, and the drag contributions at each frequency are similar. This suggests that the drag coefficient over much of the surveyed terrain is not especially sensitive to the flow direction. For the MSGL terrain there does appear to be some indication that drag would be higher for ice flow in the direction perpendicular to elongated features.

6 Discussion

6.1 Implications from the new bathymetric data

Our results provide the first observation-based, high-resolution geomorphic characterisation of the coastal bathymetry at TG, a former bed for the glacier. These data allow us to investigate bathymetric controls on ocean circulation towards the modern grounding zone, as well as to identify the locations, water depths and substrate compositions of ice-shelf pinning points and former grounding zones. The dominant bathymetric features, a NNE–SSW-trending trough and landward-flanking discontinuous ridge (Fig. 2a), represent a subtly different morphologic terrain from highly rugged, basin-dominated areas north and east of the EIS (Fig. 3a) or the moderate-relief areas of lineated terrain with fewer bathymetric highs in eastern Pine Island Bay (Fig. 2a) (Nitsche et al., 2013; Arndt et al., 2018; Kirkham et al., 2019). The continuity and orientation of the trough and ridge relate to the structure of basement rocks on the inner shelf. NNE–SSW to ENE–WSW structural lineaments have been identified by previous aeromagnetic surveys (Gohl, 2012; Gohl et al., 2013), and gravity-derived bathymetries all resolve a broad NNE–SSW ridge coincident with H1–H3, as well as deeper troughs on either side of the ridge (Tinto and Bell, 2011; Millan et al., 2017; Jordan et al., 2020). Thermochronological analyses of onshore rock samples also infer a NNE–SSW-trending tectonic rift structure (Spiegel et al., 2016).

We highlight several key differences between our new dataset and the available regional bathymetric compilations (Fig. 8). Note that we do not compare our MBES grid with the newly published gravity inversion of Jordan et al. (2020) as that study utilised the MBES dataset to constrain the inversion. Beyond that study, for the area of new MBES data in front of TG, gravity-derived bathymetry generally underestimates sea-floor depths (average of 119 m; Fig. 8a, b), whereas the IBCSO bathymetry, which is based on real sea-floor soundings but relies on gravity-inversion elevations and interpolation in this area (Arndt et al., 2013), generally overestimates sea-floor depths (average of 65 m; Fig. 8c). All of the regional datasets we examined fail to capture the higher-frequency topographic variability revealed by the new MBES data (e.g. Fig. 8d). Although sea-floor highs are sometimes >100 m shallower than the regional products predict, this effect is most notable for the troughs, which are in reality 100 to 550 m deeper than gravity-derived bathymetries and 50 to 250 m deeper than in the IBCSO dataset (Fig. 8d). When we consider cross sections of three troughs that are potential pathways for CDW to the grounding zone (T2, T3, T4; locations marked by asterisks in Fig. 8a), the depth errors are up to 250, 500 and 400 m from west to east. Using a conservative top-CDW depth of 500 m for the TG area (based on hydrographic data acquired during NBP19-02; Bastien Queste, personal communication, 2020; see Nakayama et al., 2013), we calculate the cross-sectional area that CDW occupies in these troughs from our MBES data and from the grid of Millan et al. (2017). We find that the gravity-derived bathymetry underestimates the cross-sectional areas by 7 %–38 % for two of the three troughs and that trough T3 between H1 and H2 (Fig. 3a) is not resolved at all on the Millan et al. (2017) grid (Table S1; Fig. S8). Taking this one step further, we perform a simple calculation of the oceanic heat flux through T2 for the two cross-sectional areas (Millan et al., 2017; MBES) and utilise oceanographic observations from the ASE for ocean temperatures and flow velocities (see Supplement for methods). The total heat flux through the trough cross section defined by the gravity inversion is ∼0.5 TW, and for the MBES cross section it is 1.1–1.3 TW (Table S2). This equates to an underestimation of the heat flux through T2 based on the gravity-derived bathymetry of 55 %–65 %, or more than a doubling in the heat flux through the trough using the deeper bathymetry provided by the MBES grid. To fully quantify the significance of this for the inflow of CDW to the Thwaites ice-shelf cavity and grounding zone requires the use of an ocean circulation model with the MBES as its bathymetry and that is, ideally, calibrated by CTD data within the troughs. Conversely, the identification (and implementation in models) of critical sill depths along the trough pathways could limit CDW inflow along certain routes. These findings have implications for the numerical modelling of warm-water access to the grounding zone, oceanic heat fluxes, the resultant ice-shelf melting rates, and, ultimately, projected mass losses from TG and the WAIS. However, our first-pass calculations underline the importance of high-resolution observational datasets like MBES for capturing high-amplitude bathymetric variations at short to medium wavelengths (i.e. λ<10³ m), particularly in areas close to ice-shelf cavities and the grounding zone.

Figure 8 Difference maps between regional bathymetric datasets and the MBES grid. (a) Millan et al. (2017) minus MBES grid. Yellow asterisks mark the locations of channels discussed in Sect. 6.1. (b) BedMachine Antarctica (Morlighem et al., 2019) minus MBES grid. (c) IBCSO (Arndt et al., 2013) minus MBES grid. (d) Profile data for two profiles for each of the regional datasets compared with profiles from the MBES grid; profiles are located in (b). Note that the Millan et al. (2017) and the BedMachine Antarctica bathymetries are very similar and thus return near-identical bed profiles in (d); the BedMachine Antarctica grid is included for completeness as the most recent regional dataset to cover the area, and its authors highlight the need for more coastal bathymetric datasets (Morlighem et al., 2019).

A new gravity-derived bathymetry model for the Thwaites, Crosson and Dotson–Getz area, constrained by the NBP19-02 MBES data, produced recently by Jordan et al. (2020) has improved resolution compared to previous gravity-derived models as a result of using a strapdown instrument with closer flight line spacing. Despite the improved resolution of their new model, Jordan et al. (2020) concluded that still higher-resolution observations are necessary in areas where knowledge of the bed at scales of less than a few kilometres is required. The need for high-resolution bathymetry has been underscored by recent predictive modelling studies of Antarctic outlet glaciers, which conclude that the shape of the ice-shelf cavity and knowledge of small, kilometre-scale pinning points are both key to improving predictions of ice-sheet retreat and sea-level change (Berger et al., 2016; Favier et al., 2016). Similarly, the latest high-resolution ocean models demonstrate that warm deep water reaches the grounding zone of TG through topographically constrained pathways, again highlighting the critical need for high-resolution bathymetry in making accurate predictions (Nakayama et al., 2019).

6.2 Implications from sea-floor morphology

The geometry and detailed morphology of the H1–H3 ridge also provide insight on ice-shelf pinning points. Historical grounding-zone positions, as mapped from remotely sensed ice-shelf tidal response, confirm that the Thwaites Ice Shelf is still pinned on high H1 and was pinned on high H2 as recently as 1992 and 2011 (Rignot et al., 2011). By 2011, the area of grounding on H2 had reduced to <0.5 km² (Fig. 3a), but the recent configuration and persistence of the TGT suggests that at least some parts of it remain in (ephemeral) contact with the sea floor. Thus, the exposed H1 and H2 sea-floor highs, and by analogy H3, can be studied as current or recent pinning points for the Thwaites Ice Shelf. Glacial lineations and (or) grounding-zone wedges (GZWs) on the surface of H2 and H3, as well as rare iceberg plough marks (Fig. 4b), confirm that these pinning points are mantled by some amount of unconsolidated sediment that can be ploughed or moulded by ice. Sub-bottom profiles over the H2 and H3 highs support this as they show either an incredibly smooth sea-floor response, strongly indicative of unconsolidated sediment cover, or up to 10–15 m of unconsolidated sedimentary units (Fig. S4); furthermore, coring of the top of H3 recovered several metres of sediment (Larter et al., 2020). Although the upper section of the sediment on H3 is glacimarine, deposited after grounded ice had retreated or lifted-off from the high, the presence of a GZW there and on H2 (Fig. 4b, d) suggests that at least the uppermost part of this high may have been constructed via sedimentation at the grounding zone (see Alley et al., 1989, 2007). The potential effects of this grounding-zone sedimentation are twofold: when the TG grounding zone had retreated onto these highs during the Holocene (i.e. sometime before 10.3 ka cal BP; Hillenbrand et al., 2013), GZW formation could have temporarily slowed further retreat (Alley et al., 2007). Second, continued pinning of an ice shelf on the high and GZW, when most of the grounding line had eventually retreated further landward, would have buttressed the grounded upstream section of TG. The new MBES dataset we present here and the sea-floor landforms it reveals, supported by core recovery and sub-bottom profiles, indicate that more sediment is present in this area than is typical of other Amundsen Sea inner shelf environments that experienced rapid ice-sheet retreat, including the adjacent Pine Island Bay and the Dotson–Getz palaeo-ice-stream trough (e.g. Larter et al., 2007, 2009; Graham et al., 2009; Nitsche et al., 2013, 2016). This is likely because the grounding zone was positioned for a relatively long period of time in this area immediately offshore TG and probably because areas so close to the grounding zones of most other large glacier systems have not yet become accessible for shipborne survey. Furthermore, TG has a much larger drainage basin than the Dotson–Getz palaeo-ice-stream trough and therefore the potential to erode and deliver a greater flux of basal sediment to its grounding zone. Our only constraints on grounding-zone retreat through this area (during the Holocene) are from the core on the H1 high, which shows grounded ice withdrawal from the northern part of that high by 10.3 ka cal BP (Hillenbrand et al., 2013) and grounding zones mapped from satellite-era datasets (Rignot et al., 2011). Thus, the TG grounding zone was most probably located between H1 and the current grounding zone, potentially on the sea-floor ridges identified here, for thousands of years delivering a significant volume of sediment to the area. This retreat history is in line with what we know about deglaciation more generally in the Amundsen Sea, where rapid grounding-zone retreat occurred from 15 to 10 ka to reach near modern limits (Hillenbrand et al., 2013; Larter et al., 2014; Smith et al., 2014); however, more marine dates and terrestrial thinning histories will certainly provide additional clarity and chronological constraints for TG.

The submarine landforms observed on and around the sea-floor highs raise the question of the composition of these features. Landforms on the flanks of the pinning points (gullies, slide scars, isolated blocks; Fig. 4) may indicate that these highs consist, at least in part, of an erodible (soft) material with a probable (hard) bedrock core. In marine settings, slide scars and gullies incise large, pronounced sedimentary scarps like the shelf edge (e.g. Noormets et al., 2009; Gales et al., 2013) or the headwalls of major submarine slides (e.g. Laberg and Vorren, 2000; Vanneste et al., 2006) but do not characterise hard bedrock (crystalline) settings. Further evidence comes from the new observation of flat-topped (compacted or planed-off) morphology of the H2 and H3 highs confirming that the upper part of these (down to the level of flattening) consists of a lithology that is apparently erodible by the motion of an ice shelf (based on the orientation of lineations on the highs). Other examples of flat glacial erosion surfaces planed off by ice shelves or flat-based tabular icebergs from the Arctic all document erosion into sedimentary substrates (e.g. Vogt et al., 1994; Jakobsson et al., 2010; Noormets et al., 2016). Small (<100 m high, <5 km wide), flat-topped mounds in the Ross Sea are also thought to consist of unconsolidated volcanogenic deposits rather than crystalline bedrock, and, interestingly, GZWs have built up on them, indicating that these features slowed grounding-zone retreat in that area (Lawver et al., 2012; Greenwood et al., 2018). Although we suggest flattening of the highs by the action of Thwaites Ice Shelf, we cannot say, from our data, how much erosion may have occurred. It may be that surface sediments were simply skimmed from the tops of the highs and transported towards their seaward flanks, which, in conjunction with instabilities relating to ice-shelf grounding (or ungrounding) on the highs, could have promoted slope failures on the fronts and sides of these features (cf. Bellwald et al., 2019). Despite this caveat, all of the landform evidence presented here, supported by cores and acoustic sub-bottom profiles, suggests that the tops, fronts and sides of the H2 and H3 highs are mantled by some thickness of sediment, probably over a bedrock core. Seismic-reflection profiles would be needed to determine the internal structure of these features and sediment thicknesses. In contrast, the morphology of H1 is consistent with it having a crystalline composition. This very shallow, rugged feature is cross-cut by bedrock grooves and channels typical of hard rock exposures on the inner Antarctic shelf (Lowe and Anderson, 2002; Livingstone et al., 2013), and, although bathymetric coverage over this high is incomplete, it has few planed-off sections and no glacial lineations have been identified on its surface yet (Figs. 3b; S4b). Therefore, it is clear that there is a spatial variability in pinning point morphology and composition at TG, as well as across the wider Amundsen Sea area (Figs. 2, 5). More broadly, we also note the relative scarcity of bedrock channels or other landforms related to subglacial meltwater flow in the TG MBES dataset, with the crescentic scours (H3 only; Fig. 4a) being the exception. As an example, Kirkham et al. (2019) mapped more than 1000 subglacial channels in Pine Island Bay, whereas we map only 175 forms here, albeit over a smaller area. It is not clear whether evidence of previous meltwater routing is buried by sediment in the deep troughs or has been destroyed by ice flow over the highs. Physical-property and geochemical analyses on cores from the area, acquired as part of ITGC, should shed light on the frequency and magnitude of meltwater release during the retreat of grounded ice over the sea-floor highs. Schroeder et al. (2013) identified a transition from a distributed channel network with ponded water behind ridges at the modern grounding zone to a system of concentrated channels downstream. It is possible that a similar configuration for the basal hydrological system occurred in this area as ice retreated over the offshore highs and that evidence is preserved in the marine sedimentary record.

The apparent shaping and fragmentation of the H2 and H3 highs (Fig. 4a, b) highlights a potential feedback mechanism between bed properties (composition and topography) and glacier retreat dynamics. If the substrate of a pinning point is soft enough to be moulded by the flow of an ice shelf and to be susceptible to slope failures, it might be eroded over time. Erosion of material from the surface of a pinning point, as it is planed off, in conjunction with retrogressive failures at its seaward flank and possibly larger slide or slump events may act to reduce its surface height, as well as its surface area, by “eating away” at the frontal/side flanks until it cannot serve as a pinning point for the ice shelf (and glacier ice upstream) any longer. As long as the ice shelf continues to move over the high, this process of unpinning would be exacerbated by any increase in flow velocities (leading to increased erosion) and/or by ice-shelf thinning (leading to ungrounding), due to either flow acceleration or sub-ice-shelf melting. The result, in a setting with soft erodible pinning points, is the potential for increased ice-flow velocities to accelerate pinning point destruction which, in conjunction with simultaneous ice-shelf thinning in response to sub-ice-shelf melting, could promote ungrounding earlier than would occur on a corresponding hard, less erodible pinning point. Needless to say, in order for this feedback to occur, an ice shelf would have to continue to flow quickly over the sea-floor high(s) and not form an ice rise, under which erosion rates are considered to be low (e.g. Matsuoka et al., 2015). At TG we note that, at least for the duration of the observational record (∼55 years), the fast-flowing part of the glacier, which feeds the TGT (Fig. 2), has continued to move over the H2 and H3 highs, periodically extending several tens of kilometres before calving (Ferrigno et al., 1993; Rabus et al., 2003; MacGregor et al., 2012), whereas the ice rumple at the end of the EIS restricts flow over H1, with most ice flow being diverted around the rumple (Rignot et al., 2001). It is perhaps interesting to also consider that upwards of 50 m of relative sea-level fall due to glacio-isostatic uplift is thought to have occurred on the inner Amundsen Sea shelf during the Holocene (Whitehouse et al., 2012) and that the uplift of any pinning points would naturally counter ungrounding. Therefore, although we can only speculate on the exact mechanisms affecting rates of unpinning, we suggest that the composition of sea-floor pinning points may be an important factor in their ongoing ability to buttress large Antarctic ice shelves.

6.3 Implications from the new bed roughness data

One major objective of our research is to assess the deglaciated submarine terrain offshore from TG as an analogue for the modern bed to gain new insights on TG bed characteristics. The consistency of derived power spectra and drag contributions for bed profiles from the inner ASE shelf and for upstream areas of Pine Island and Thwaites glaciers (Figs. 7, S7) indicates that the roughness properties of the offshore and onshore areas are comparable across all resolvable frequencies. Furthermore, observations confirm that recent grounding-zone retreat affecting TG has occurred over a series of bedrock ridges with the loss of pinning points and formation of new cavities (Tinto and Bell, 2011; Milillo et al., 2019; Jordan et al., 2020). Further upstream, about 100 km from the recent grounding zone, analyses of radar specularity suggest that the modern TG bed is characterised by high roughness attributed to bedrock cropping out at the glacier base (Schroeder et al., 2014). Our data reveal that the bedrock ridges and intervening troughs underlying the modern grounding zone (Holt et al., 2006; Morlighem et al., 2019), with length scales of up to tens of kilometres and amplitudes of up to several hundreds of metres, constitute a morphological terrain similar to the coastal bathymetry (Figs. 2, 3). Further, we demonstrate that this rugged terrain would exert the same strong influence on basal drag for an overriding ice mass (assuming no cavitation) (Figs. 7, S7). This is consistent with results from inverse methods that determine basal drag for the modern TG bed (Joughin et al., 2009; Arthern et al., 2015). We also note that crag-and-tail landforms (which form subglacially) extend down to the floors of the deep troughs (e.g. Figs. 3, 4b). This confirms that, at least at the time when these features formed, ice of an expanded TG was grounded in the troughs as well as on the highs and probably experienced high basal shear similar to ice at the present-day grounding zone. The orientation of the crag and tails also confirms that ice flow was not directed along troughs but rather overrode the existing topography; this finding is consistent with cosmogenic exposure data from Bear Peninsula (Fig. 1) showing that the ice-sheet surface rose above the top of this terrain during the last glacial period (Johnson et al., 2017).

For shorter wavelengths of bed topography, we can consider the form of the individual sea-floor highs over length scales of several kilometres. We interpret the morphological characteristics of these features as being consistent with the correlation of morphology with bed type known from onshore glacial–geological studies of crag-and-tail type features (e.g. Benn and Evans, 2010), which more recently has been described from on-ice seismic-reflection profiles both for TG (Muto et al., 2019a, b; Holschuh et al., 2020) and beneath the Rutford Ice Stream (Fig. 4 in King et al., 2016). Still, we recognise that high-resolution seismic-reflection data over our bathymetric highs would be required to confirm this. Specifically, the correlation is between hard beds on the stoss sides of topographic highs, associated with crag-and-tail landforms, and soft or sedimentary beds on the lee sides of these features. This pattern is clearly replicated over the H2 and H3 highs, which have rugged upstream ends with crag-and-tail landforms, glacial lineations over their tops and sedimentary tails on their downstream ends (Fig. 4). The correlation of bed types with sea-floor highs (and ridges) holds true for several other areas of the inner shelf around West Antarctica, where streamlining of bedrock highs has often produced landforms with sedimentary tails on the lee sides of bedrock obstacles (e.g. Larter et al., 2009, Graham et al., 2009; Livingstone et al., 2013; Nitsche et al., 2013, 2016), although this is not always the case. Thus, the variability in bed types on topographic highs in offshore regions may provide useful constraints on bed type variability beneath the modern glacier.

Regarding the spectral analysis of roughness and basal drag contributions presented here (Figs. 6, 7, S7), we acknowledge that these only provide an order of magnitude assessment of the contribution to basal drag from the different wavelength scales resolved by the bathymetric DEM (Fig. 2). Analysis beyond the simple 2D-flow-line theory used here (see Sect. 4.2) would be needed to account for 3D-flow effects, as well as for the non-linear dependence of ice viscosity on stress (Glen, 1955). Here, we have not specified the physical mechanism controlling λ_N, the shortest roughness wavelength that influences basal drag. Candidate mechanisms that might limit the influence of roughness at small spatial scales include cavitation (Fowler, 1986), fracture and plucking of crystalline or sedimentary rocks, the formation of a weak internal shear zone (Liu et al., 2020), bulldozing of unconsolidated sediment, or regelation flow around small obstacles (Weertman, 1957). More sophisticated theories accounting for the potential of ice to form cavities in the lee of obstacles could be deployed similarly, but the drag contribution would then depend also on water pressure (Fowler, 1986; Schoof, 2005). Process models of subglacial hydrology, phase change, fracture and sediment transport could all be incorporated in to a more elaborate analysis using MBES datasets as input.

It is clear from our results that the increased spatial resolution of the MBES data is critical for capturing the high-frequency bathymetric variability on the inner continental shelf seaward of TG, which is necessary to understand warm-water incursions into sub-ice-shelf cavities (Figs. 8, S8; Nakayama et al., 2019). The strong correlation of our observations with interpretations of the present bed conditions of TG and, therefore, the robustness of this deglaciated terrain as an analogue for the modern bed further demonstrates that more information can be gleaned from this type of marine dataset (i.e. near-continuous bathymetry with spatial resolution better than 0.05 km). For example, the 3D nature of MBES (with approximately equal resolution in all directions) means that bathymetric variability could be examined in any direction, not only along survey lines, as has been the case until recently with all onshore radar or seismic-reflection profiles of extant bed topography, and over a variety of spatial scales. These analyses add to our understanding of across-flow contributions to basal drag or hydraulic potential (e.g. Muto et al., 2019a) and allow us to consider the spatial variability of bed types (e.g. sedimentary vs. hard beds), particularly where sea-floor sediments are also imaged by seismic profiles and/or cored for ground truthing. Similarly, the application of theories of subglacial processes as discussed above to high-resolution bathymetric datasets will increase our understanding of ice flow over high-frequency bed roughness, particularly if combined with ultra-high-resolution (sub-metre-resolution) bathymetries from AUV surveys (e.g. Davies et al., 2017; Dowdeswell et al., 2020). Indeed, AUV surveys and (or) a dense grid of seismic soundings (only obtainable from non-crevassed ice shelves) are the only way to determine bed geometry in ice-shelf cavities. New techniques such as swath radar that can image the present glacier bed in 3D (Paden et al., 2010; Jezek et al., 2011), albeit in narrow swaths, have already been employed on TG (Holschuh et al., 2020) and could be used in conjunction with offshore bathymetric data to build a better-informed, more complete and more uniform resolution picture of basal conditions under TG and at its grounding zone.

7 Conclusions

New 3D bathymetric data from just offshore Thwaites Glacier reveal that the coastal bathymetry is dominated by a ∼65 km long, ∼1200 m deep trough and a discontinuous ridge with water depths of 650 to <100 m. Spatial variations in the morphology of the ridge segments/highs and available acoustic sub-bottom profiler data suggest differences in substrate composition along the ridge, with the two southernmost highs having a significant erodible component at least in their upper parts, which are sedimentary in composition. The geometry (flat tops) and landform evidence (glacial lineations, gullies, sediment fans) indicate that the bathymetric highs were planed off and variously eroded by the action of Thwaites Ice Shelf as it flowed over them, presumably reducing the height of these former pinning points over time. A feedback mechanism during unpinning may have occurred, whereby as the ice shelf started to lose contact with the high and frontal buttressing weakened, the resultant increase in flow velocities exacerbated erosion of the high and facilitated further unpinning of Thwaites Glacier.

We present three lines of evidence that this coastal bathymetry provides a good analogue for the modern grounding zone of Thwaites Glacier. First, on length scales of several tens of kilometres, the ridge and trough morphology is consistent with the bed topography of the grounding-zone area based on available DEMs and over-ice geophysical data. Second, our spectral decomposition of roughness and basal drag over this rugged, deglaciated terrain is consistent with similar spectral decompositions, and inversions of basal drag, for profiles from the modern grounding zone area and for areas of the Thwaites bed, where bedrock crops out subglacially (e.g. Schroeder et al., 2014). In contrast, smooth beds, characterised by thick sedimentary substrates and linear glacial landforms, produce distinctly different power spectra and drag contributions. Third, the distribution of landforms and substrate types (unconsolidated sediment vs. bedrock) over the ridge indicates that it displays the same correlation of bed type with topography that has been described for upstream bed areas and inferred for the grounding zone (Muto et al., 2019a, b; Holschuh et al., 2020). As such, further analyses of this deglaciated terrain may provide realistic constraints on across-flow roughness and bed type distribution and should inform geophysical observations of the modern TG bed that will be acquired as part of ITGC.

As discussed above, observational datasets like MBES are required seaward of Antarctic ice shelves in order to capture the high-frequency variability that characterises the bathymetry of nearshore areas. Bathymetry derived from gravity inversions cannot adequately reproduce the kilometre- to sub-kilometre-scale features that are important for accurately calculating inflows of warm ocean water in troughs and for defining the topographic highs that may act as pinning points for ice shelves and as barriers to warm-water incursions.