Characterizing bed roughness on the Antarctic continental margin

Abstract

Spatial variability in bed topography, characterized as bed roughness, impacts ice-sheet flow and organization and can be used to infer subglacial conditions and processes, yet is difficult to quantify due to sparse observations. Paleo-subglacial beds of formerly expanded glaciers found across the Antarctic continental shelf are well preserved, have relatively limited post-glacial sediment cover and contain glacial landforms that can be resolved at sub-meter vertical scales. We analyze high-resolution bathymetry offshore of Pine Island and Thwaites glaciers in the Amundsen Sea to explore spatial variability of bed roughness where streamlined subglacial landforms allow for the determination of ice-flow direction. We quantify bed roughness using std dev. and Fast Fourier Transform methods, each employed at local (10⁰ km) and regional (10^1–2 km) scales and in along- and across-flow orientations to determine roughness expressions across spatial scales. We find that the magnitude of roughness is impacted by the parameters selected – which are often not sufficiently reported in studies – to quantify roughness. Important spatial patterns can be discerned from high-resolution bathymetry, highlighting both its usefulness in identifying patterns of streaming ice flow and underscores the need for a standardized way of characterizing topographic variability.

Introduction

Approximately half of the Antarctic Ice Sheet (AIS) lies farther than 5 km from any direct subglacial bed measurement (Morlighem and others, 2020). Interpolation techniques have been used to compensate for unresolved bed topography, but these can lead to the misrepresentation of terrain (MacKie and others, 2021). Moreover, existing topography products that rely on radio sounding systems often fail to identify deep subglacial troughs, which are critical in determining ice stream flow direction (Morlighem and others, 2020). The flow of ice streams is modulated by processes at the ice-bed interface (Stokes and others, 2007; Stokes, 2018), but the extent to which subglacial topography promotes or inhibits ice flow remains uncertain (Favier and others, 2014; Robel and others, 2022). This is due to the complex nature of processes at the ice-bed interface, which makes it difficult to accurately model ice-sheet behavior. The parameterization of basal traction in ice-sheet models is largely reliant on satellite-based observations of the ice-sheet surface (Arthern and others, 2015) and remains a considerable source of uncertainty (Ritz and others, 2015). The lack of direct and high-resolution (i.e. sub-kilometer) observations of subglacial topography limits our ability to separate skin drag and form drag components, often combined when defining basal traction (Kyrke-Smith and others, 2018). The skin drag component of basal friction is impacted by basal meltwater and properties of the uppermost layer of deformable sediments (Iverson and Zoet, 2015), which are not resolvable by topography (i.e. elevation) products. The form drag component, however, which describes the resistance to ice flow that originates as ice deforms around bed obstacles (Weertman, 1964), can be represented by bed roughness measurements. Sliding theories suggest that perturbations at the meter scale can generate enough basal drag to limit sliding (Weertman, 1957; Schoof, 2002; Robel and others, 2022). This is supported by the observation that form drag produced by subglacial roughness can produce significant shearing as grounded ice retreats over rugged topography (Hogan and others, 2020). Thus, the inclusion of high-resolution basal topography as a parameter is essential in producing realistic basal motion (Whillans and van der Veen, 1997; Winsborrow and others, 2010; Morlighem and others, 2020; Law and others, 2023). Bed roughness, defined here as ‘the extent to which terrain varies vertically over a given horizontal distance’ (Rippin and others, 2014), is therefore a useful tool in determining the influence that bed topography exerts on ice-flow velocities (Cooper and others, 2019; Law and others, 2023), though the range of scales at which bed roughness can be quantified is dependent on the spatial resolution of the elevation data.

Studies over large areas (>500 km²) of the AIS use bed roughness derived from radio-echo sounding (RES) to investigate the impact of bed topography on basal processes (e.g. Siegert and others, 2004, 2005; Taylor and others, 2004; Rippin and others, 2006, 2011, 2014; Bingham and Siegert, 2007, 2009; Li and others, 2010). RES provides high along-track resolution, but the transect spacing often exceeds 10 km (Siegert and others, 2004; Bingham and others, 2007; Rippin and others, 2014), which is too wide to capture roughness associated with landform assemblages typical of paleo-ice stream beds (wavelengths between 10¹ and 10² m; Falcini and others, 2018). Results from early studies suggested that variations in bed roughness were spatially organized, where rough beds were found in inland regions of slow-moving ice and smoother beds were found downstream in regions of fast-flowing ice streams (Siegert and others, 2004; Taylor and others, 2004; Bingham and Siegert, 2007). This implies a straightforward relationship where roughness is controlled by, or is a reflection of, ice-flow velocity and distance from the grounding line. However, more recent studies have shown that fast flow is not always associated with a smooth bed (Rippin and others, 2011; Schroeder and others, 2014; Falcini and others, 2018).

The degree to which bed roughness can identify bed lithology and subglacial bedforms remains underexplored, and the quantification of bed roughness at scales where individual landforms can be resolved has been largely underutilized as a tool to infer bed conditions in lieu of using ice-sheet surface inversions (Taylor and others, 2004; Bingham and others, 2017). While the orientation of elevation transects has been previously considered (Rippin and others, 2014; Bingham and others, 2017; Falcini and others, 2018, 2021; Cooper and others, 2019), not many studies have explored the impact that different elevation detrending scales have on bed roughness at small horizontal scales (10^1–2 m) where local topography is resolved. This is especially problematic as roughness analyses are inconsistently calculated, and are varied in roughness scales of interest across different studies (Smith and others, 2014).

The acquisition of bathymetric data over the deglaciated seafloor around Antarctica presents the opportunity to explore former subglacial bed conditions at higher resolutions and with greater spatial coverage than beneath contemporary ice streams. The seafloor of the Amundsen Sea Embayment (ASE) records the Last Glacial Maximum (LGM) and post-LGM glacial history of the formerly merged ice stream sourced from Pine Island Glacier (PIG) and Thwaites Glacier (TG) (Graham and others, 2016; Fig. 1). During the Quaternary period, these glaciers excavated a trough extending over 500 km in length, while the ice sheet is believed to have reached the shelf edge (Graham and others, 2010). Seismic profiles along the trough indicate that the inner shelf, close to the modern grounding line, predominantly consists of crystalline bedrock, whereas the middle and outer shelf exhibit a younger and unlithified sedimentary substrate (Lowe and Anderson, 2002). On the inner shelf, streamlined bedforms are prevalent in areas of thin sediment cover, though some sediment-filled depressions are also observed near the modern ice shelf front (Fig. 1b; Nitsche and others, 2013). Moving toward the middle shelf, the topography becomes more rugged, with shallow sills and a network of subglacial meltwater channels cutting into the bedrock (Figs 1c, d; Nitsche and others, 2013; Kirkham and others, 2019). The paleo-ice stream beds of PIG and TG converge in the middle shelf, where streamlined bedforms are abundant (Fig. 1d; Graham and others, 2016). Further downstream from the convergence, drumlinized bedforms evolve into mega-scale glacial lineations (MSGLs) at the transition of crystalline bedrock to sedimentary substrate (Fig. 1e; Wellner and others, 2001; Lowe and Anderson, 2002).

Figure 1. (a) Study sites in the Eastern Amundsen Sea and Thwaites Glacier marked by the black, numbered boxes. Arrows show the general direction of paleo-ice flow for Pine Island and Thwaites glaciers, which merged at site 3. (b–e) Multibeam bathymetry of sites in the Eastern Amundsen Sea. Grid cell sizes 35–50 m, from Nitsche and others (2013). (b) Ice-shelf proximal site consists of crystalline bedrock (Cr) mixed with deep pockets of unconsolidated sediment and linear bedforms, i.e. streamlined grooves (SG), crag-and-tails (C-T), and drumlinoid features (Dr). (c) Inner shelf site displaying crystalline bedrock, rugged topography and sinuous channels (Ch). Color ramp as for (b). (d) Site where the Pine Island and Thwaites paleo-ice streams merged, resulting in a change in ice-flow direction. Presence of deep basins (Ba) and channels, a flat topographic high (TH) and grooved crystalline bedrock (SG). Color ramp as for (b). (e) Transition between crystalline bedrock and unconsolidated sediment. (f, g) Swath-radar data from Holschuh and others (2020). (f) Upstream site of the Thwaites bed with MSGLs and bedrock protrusions at shallower depths. (g) Downstream site with streamlining and crag-and-tails either side of large exposed bedrock.

Recent gravity-derived bathymetry beneath the Thwaites ice shelf and ice tongue reveals similarly complex topography with relief comparable to the study sites located in the middle-shelf of Pine Island Bay (PIB; Jordan and others, 2020). Hogan and others (2020) show that high-resolution bathymetry is necessary to capture the spatial variability of bed topography on the inner shelf just offshore of PIG and TG. These two glaciers were responsible for >30% of the annual discharge from the West Antarctic Ice Sheet (WAIS) between 2009 and 2017 (Rignot and others, 2019), and assessing the variability of bed topography and roughness offshore of these glaciers can provide an analog for the subglacial environment of contemporary glaciers and ice streams.

Through high-resolution bathymetry offshore of Pine Island and Thwaites glaciers, as well as the inclusion of elevation models derived from swath-radar underneath TG, this study compares bed roughness results between different methods, orientations and detrending scales to determine the influence of each of these parameters on bed roughness results and the implications for ice behavior. We then compare results from high-resolution elevation models to the coarser BedMachine dataset to assess any potential roughness signatures that might be misrepresented when bed topography is not available at a high spatial resolution. The BedMachine dataset uses a mass conservation method and incorporates various data sources to fill data gaps and provide compatibility with numerical models (Morlighem and others, 2017, 2020). Lastly, we incorporate roughness results from bed topography data obtained from geostatistical simulations conducted on Jakobshavn Glacier by MacKie and others (2021). This allows us to compare roughness statistics between direct observations in the ASE and stochastically simulated topography.

Methods

A total of six study sites were used for analysis where a compilation of multibeam echosounder bathymetry data in the eastern ASE was used to produce the gridded 50 m bathymetric dataset used for analysis of sites 1–4 (Fig. 1a; Nitsche and others, 2013). Elevation data for the bed of TG (sites 5–6) come from swath-radar published by Holschuh and others (2020). The study sites were selected to assess the relationship between topography and the formerly expanded PIG-TG system during and following the LGM (Graham and others, 2010; Nitsche and others, 2013). The diverse set of glacial landforms across the sites allows us to assess and compare roughness values across different relief, bed slopes and geologies. Topographic realizations used for statistical analysis are from simulations by MacKie and others (2021) (Fig. S1).

The site grids were drawn where there was continuous data coverage to ensure that missing data would not impact the roughness results. For each grid, transects oriented parallel and orthogonal to paleo-ice flow direction were inferred from streamlined subglacial landforms, such as grooves and glacial lineations (Figs 1b–g). Based on the width of streamlined features observed in the ASE, the spacing between transects was set to 500 m, which also corresponds to the spatial resolution of BedMachine Antarctica (Morlighem and others, 2020). Elevation values were extracted every 50 m along each transect, to match the horizontal resolution of multibeam bathymetry (Fig. 2). To assess how the configuration of basins and channels may impact roughness measurements, we created a 500 m buffer around the subglacial channels mapped by Kirkham and others (2019) that fell within our study sites in PIB. This buffer was used to compare the roughness values associated with subglacial channels to the roughness of the surrounding area. Once grids for all sites were constructed, a workflow to compare roughness results between different orientations, detrending techniques and scales was implemented (Fig. 3a). Elevation transects were detrended using two methods to remove long-wavelengths components (Taylor and others, 2004): (1) a linear detrend of the entire transect using least-squares regression to assess regional-scale topography and (2) subtracting the mean elevation of a 1.6 km moving window to match the minimum moving window used in the Fast Fourier Transform (FFT) analysis discussed below, which we use to characterize local, kilometer-scale topography. By quantifying roughness at both regional and local scales, we assess the sensitivity of roughness measurements to different detrending methods.

Figure 2. Raw elevation transects and calculated slope transects, both of which have data points at 50 m increments for East Amundsen Sea sites. Left and right columns show transects in the parallel and orthogonal orientations relative to paleo-ice flow direction, respectively. Slope is calculated as the dimensionless ratio of the vertical to horizontal change at every 50 m increment.

Figure 3. (a) Flow chart outlining the steps taken to calculate roughness using the SD and FFT methods described in the Methods section. (b) Example of a single raw elevation profile and corresponding detrended profiles using a local (red) and regional (blue) detrend method. Profile comes from Figure 2a.

Roughness was calculated using a std dev. (SD) method and an FFT method. The SD method provides a metric for the variation of amplitudes in elevation in a straightforward manner, which can be quickly applied to numerous transects with little computational power. SD is commonly used to measure roughness in the Earth Sciences (Smith, 2014), though it is unable to capture the horizontal frequency of undulations. Fourier transformations were introduced in some of the earliest studies on ice sliding over sinusoidal (i.e. idealized) topography, where it was proposed that bed roughness could be described in terms of the power spectrum of the bed elevation (Kamb, 1970). FFT analysis converts bed elevations into a wavelength spectrum to calculate the amplitude and the spatial frequency of undulations present in the bed; the methodology for the FFT calculations used in this analysis follows Li and others (2010). We present a basal roughness index (ξ), which reflects the magnitude of vertical deviations in the bed and is calculated by taking the integral of the spectral power density, S(k), over the moving window (Eqn (1)).

(1)$$\matrix{\xi = \int\!{S( k ) \,{\rm d}k\;\ }} $$

To perform FFT calculations, the convention is to use a minimum of n = 32 data points in each moving window (Taylor and others, 2004); therefore, given the 50 m horizontal resolution of the bathymetry data, a moving window of 1.6 km (50 m × 32) was used to calculate roughness. Both the SD and FFT methods were used to quantify roughness at the local (1.6 km) and regional (20–50 km) scales defined earlier using the two methods of elevation detrending.

To determine the impact of transect orientation on roughness measurements, the directionality of roughness was assessed by comparing parallel- (R _∥) and orthogonal-roughness (R _⟂) values where transects intersect. By implementing the anisotropy ratio (Ω) introduced by Smith and others (2006; Eqn (2)),

(2)$$\matrix{ {\Omega = \displaystyle{{R_\parallel\,{-}\;R_\bot } \over {R_\parallel\,{ + } \;R_\bot }}\;\ } \cr } $$

where R represents the roughness values obtained from the SD and FFT methods, the directionality of roughness values can be compared across sites and methods. Anisotropy ratios approaching 1 suggest R _∥ >> R _⟂, values approaching −1 indicate R _∥ << R _⟂, and values close to  0, suggest an isotropic surface which can represent a smooth or truly random landscape (Falcini and others, 2021). Output results for the analysis are point data, which we use to interpolate and generate raster products containing roughness and anisotropy values.

The same conditions for the SD method were applied to the BedMachine dataset (Morlighem and others, 2017, 2020) to investigate where roughness values may be under- or over-estimated depending on the spatial resolution of the elevation raster used. Roughness results derived from high-resolution datasets were subtracted from BedMachine results to generate rasters showing where and by how much BedMachine results differ from the ‘true’ roughness of resolvable landforms. We also applied the SD method to one of the 250 topographic realizations from MacKie and others (2021) as well as the corresponding elevation dataset from Greenland BedMachine (Morlighem and others, 2017). The main channel is present in both elevation datasets (Fig. S1), but the BedMachine dataset more closely resembles the average of all topographic realizations from MacKie and others (2021). The outliers discussed in the results and discussion sections are defined as roughness values exceeding ±1.5×IQR/√n, where IQR is the interquartile range.

Results

At site 1, on the inner continental shelf closest to the contemporary Pine Island calving line (Fig. 1b), roughness is relatively consistent at both the local and regional scales (10¹–10² m, 10² m²) for both methods used. The lowest roughness measurements (<5 m) from the SD method are found where multichannel seismic data over the deepest water depths of 950–1050 m reveal a basin infilled with >300 m of unconsolidated sediments (Nitsche and others, 2013) with and without the presence of small-amplitude (<5 m) lineations (Fig. 4a). High roughness values from the SD method (>30 m) are found on the slopes of streamlined landforms such as crag-and-tails and whaleback ridges that taper in the direction of paleo-ice flow, previously identified by Nitsche and others (2013). Other high roughness values (20–30 m) are found within the channel buffer, particularly in the orthogonal orientation, where the mean roughness of the channel buffer is 5 m greater than the mean of the site (Table S1). Roughness values of 45–60 m are found where channels have a cross-sectional area >35 000 m², as mapped by Kirkham and others (2019) and where relief between channel thalweg and surrounding bedrock is >250 m. While the spatial pattern for FFT results is similar to the SD method, areas of extreme relief in the orthogonal orientation create outliers (>7000 m²) that are over two orders of magnitude greater than the median (72 m²; Table S2). Roughness ranges and medians for all transects in the parallel orientation are consistently lower than the orthogonal transects for both spatial scales and roughness calculation methods (Figs 5, S2a, Table S2).

Figure 4. Absolute roughness measurements for parallel transects in sites 1–4 in the Eastern Amundsen Sea showing the difference in spatial distribution between scales and methods. Blue lines are subglacial meltwater channels and black, hatched polygons are relict subglacial lakes (Kirkham and others, 2019). White arrows indicate direction of paleo-ice flow.

Downstream (i.e. seaward) of site 1, exposed crystalline bedrock at site 2 becomes more prevalent and streamlined landforms are less common. The topography is rugged with water depths between 500 and 1500 m and several deeply incised channels (Fig. 1c). The upstream section of site 2 is characterized by parallel crag-and-tails and drumlins (Nitsche and others, 2013), while a deep basin floored by streamlined landforms is located downstream and is flanked by steep slopes (Fig. 1c). The magnitude of median roughness for site 2 is nearly double that of site 1 when using the SD method and nearly three times higher in the parallel orientation. When using the FFT method, the medians for site 2 increased by a factor of 5 and 12 for the orthogonal and parallel orientations, respectively (Fig. 5, Table S2). High outliers (>60 m) from the SD method are concentrated along the length of the walls flanking channels inferred as subglacial in origin and with a cross-sectional area >35 000 m². Conversely, the lowest values (<8 m) are present in areas of shallow topography on top of bedrock highs and in deep areas identified by Kirkham and others (2019) as relict subglacial lakes, where no glacial landforms are observed. The FFT method produced high roughness values, but they were not as widespread as those produced by the SD method. Instead, outliers (>6000 m²) are found in areas where the relief associated with subglacial meltwater channels is >250 m (Figs 4b, S2b).

Figure 5. Distribution of the basal roughness parameter (ξ), employing a 1.6 km moving window across all sites. The boxes represent value points between the first and third quartiles (IQR), and the black horizontal bars indicate the median. Individual outliers are plotted where values exceed±1.5×IQR/√n. (a) Distribution of values employing the SD method. (b) Distribution of values employing the FFT method, only applied to sites 1–4. Note use of logarithmic scale on the y-axis.

Site 3, where the paleo-ice streams of Pine Island and Thwaites glaciers merged during the LGM (Larter and others, 2014), has the greatest relief among the sites studied, with water depths between 375 and 1650 m. Topography here is dominated by deep basins that are up to 300 m below the surrounding seafloor (Nitsche and others, 2013) and a central flat topographic high. The steep slopes at the edges of the basin and the large meltwater channels that connect them generate the highest roughness values across all sites, particularly at the regional scale where values are 2–3 times greater than at the local scale (Table S1, S2). These channels are mostly found in the upstream region of site 3, have a cross-sectional area >35 000 m² and have SD roughness values of 30–80 m. Conversely, the streamlined seafloor of the central topographic high, which is cross-cut by geological structures (Graham and others, 2016), records the lowest roughness values of this site in all methods, orientations and scales. Similar to site 2, outliers in the SD-based roughness are located along the walls of the deep basins, while FFT-based roughness outliers are spatially isolated (Figs 4c, S2c).

Site 4 has gentle relief and marks the transition between exposed crystalline bedrock in the middle-shelf to an unlithified sedimentary substrate in the outer-shelf (Wellner and others, 2001; Lowe and Anderson, 2002). Drumlinized bedforms are absent downstream of the transition, and the landscape in this site is dominated by MSGLs, which are not observed in any of the other sites in the ASE. Site 4 has the lowest median roughness values of all sites (Table S2), and the high roughness values in site 4 coincide with the presence of drumlinized features in the upstream region. These SD values are comparable in magnitude to the ones produced by the low-amplitude streamlined landforms found in sites 1 and 3 at the local scale (10–30 m). The lowest roughness values (<1 m) across all sites are found in the downstream region of site 4, where unconsolidated sediments blanket the sea floor and relief is minimal (Figs 4d, S2d). These roughness values are comparable to the ones seen in the sediment-filled basin present in site 1.

Sites 5 and 6 are located in the subglacial environment of TG and are comprised of MSGLs and crag-and-tails similar to those found in sites 1 and 4. The upstream region of site 5 is dominated by elongated bedforms thought to be the tails of crag-and-tails (Alley and others, 2021), which transition into bedrock protrusions downstream (Fig. 1f). At the local scale, high roughness values (>30 m) are found where these protrusions generate relief and where MSGLs terminate in moats, as described by Holschuh and others (2020). Low roughness values are present in the sedimentary basin where MSGLs are present, particularly in the parallel orientation (Fig. S3a). Site 6 is dominated by large bedrock protrusions or ridges, with lineations present upstream and downstream from said ridges. High roughness values (>30 m) are concentrated where the ridges generate relief (<300 m) in both orientations and in the stoss side of the crag-and-tails present downstream of the topographic high. Low roughness values (<5 m) are found where lineations are present in areas of low elevation (Figs 1g, S3).

Roughness results obtained from the topographic simulations conducted by Mackie and others (2021) exhibit similar patterns to those observed in our study sites: roughness is higher in the orthogonal orientation and at the regional scale. When comparing a single realization to the average of all realizations, the former displayed higher roughness values and had higher mean values than any of our study sites, regardless of scale. Meanwhile, the mean value derived from the average of all topographic realizations, aligns closely with our study sites at both scales. The topographic realizations are dominated by a 70 km-long, 4 km-wide, slightly meandering channel oriented in the direction of ice flow. Roughness measurements from the local scale are insufficient to capture the roughness signature generated by the wide channel. Nevertheless, local-scale roughness measurements still indicate values as high as 180 m, exceeding the maximum roughness observed in our study sites at that scale (Figs 4, 5). In contrast, the regional scale is capable of encompassing the entire channel width, offering a more comprehensive depiction of the roughness associated with the channel. Roughness values at the regional scale can reach up to 450 m, significantly exceeding the maximum value of 295 m observed in our study sites at the same scale. The average of all realizations also exhibits significant roughness associated with the channel, albeit to a lesser extent (375 m). However, it is important to note that this approach introduces artifacts that create spikes in roughness values and generates smoother terrain outside the channel, potentially obscuring the presence of other smaller tributaries (Mackie and others, 2021).

Discussion

Impact of methodology employed

The spatial distribution of high and low roughness measurements remains largely consistent between the SD and FFT methods across all sites, similar to results presented by Falcini and others (2018). Yet, there are important differences when comparing the effectiveness of detrending methods in capturing regional and local roughness (Figs 4, 5, S2). Both the FFT and SD methods yield right-skewed roughness distributions, meaning most roughness values are on the lower end, with a few larger values on the higher side. Notably, in the context of this analysis, we specifically consider high outliers. The FFT method produces more outliers due to the significantly greater magnitude and range of roughness values (Fig. 5b, Table S2), as previously noted by Rippin and others (2014). On average, the percentage of data points considered outliers is higher for the FFT method (11%) than for the SD method (5.8%). The number of outliers in the FFT method remains relatively constant across the two detrending scales considered (0.2% change), while the number of outliers in the SD method increases by an average of 19% when using the regional detrend (i.e. detrending across the whole elevation profile) compared to the local detrend of 1.6 km. The increase in outliers when using the SD method is particularly noticeable at sites 1–3, which have high relief and exposed crystalline bedrock (Lowe and Anderson, 2002). However, at site 4, the most downstream site in the ASE, where gentle relief and streamlining of soft sediments are observed (Figs 2g–h), the number of outliers decreases by 17.5%. Contrary to the expectation that using a regional detrend would increase roughness variability, large depositional environments typically found at the downstream ends of paleo-ice stream beds, may exhibit decreased variability and roughness signatures may be dominated by sediment accumulation and drowning of antecedent topography.

While the spatial distribution of roughness values is similar between the local and regional scales across all sites, the regional scale exhibits greater magnitude (Fig. 5), due to the nature of the detrending method, where the range of detrended elevations is considerably greater (Fig. 3b). The local scale effectively provides roughness characterization for smaller-scale (<1.6 km) features and the regional scale considers larger features (i.e. deep basins and meltwater channels) while still removing long-wavelength trends. The choice of scale for detrending and the moving window used to calculate roughness has a direct impact on roughness results. As such, their interpretation requires careful consideration and should not be directly compared with studies that use different scales (Smith, 2014). When comparing results between the two different scales, we found that the spatial distribution of high roughness values showed minimal variation when the FFT method was used. Specifically, when the regional detrend was applied, the average median roughness increased by 14% with the FFT method, compared to a 64% increase using the SD method. The increase in outliers observed with the SD method when the regional detrend was applied coincides with a more widespread distribution of high roughness values across all sites (Figs 4, S2). The SD method detected a greater spatial coverage of high roughness values than the FFT method.

Ultimately, the FFT method is less susceptible, but still impacted, by the scale used to calculate roughness. Although the SD method is more susceptible to the scale used, it can detect roughness at local scales more effectively than the FFT method. Using the SD method at the local scale, high roughness values are typically observed along the length or width of a specific landform, whereas on a regional scale, high roughness values tend to extend beyond landform boundaries and encompass landform assemblages (Figs 4, S2). The FFT method yields similar spatial distributions of high roughness values at both scales, making the distinction less clear. The size of the moving window used to detrend elevation profiles was ultimately dependent on the spatial resolution of the elevation products available, the method employed and the size of landforms present within the study area. Importantly, the SD method is not limited by the 32 sample points required for FFT analysis and can be implemented over even smaller windows. Although the SD method can quantify roughness using a smaller moving window, we opted to use a 1.6 km window in order to make direct comparisons between the SD and FFT methods. Since bathymetric datasets are available for various sectors of Antarctica at a spatial resolution of 50 m or finer, the 1.6 km moving window used in this study can be used as a local scale to compare SD and FFT methods in other deglaciated regions. However, the SD method may be preferred due to its ease of use and ability to detect roughness patterns at scales smaller than 1.6 km without the extreme variation in results associated with the FFT method. To ensure comparability, it is crucial that roughness studies report the scale used to detrend elevation profiles, and the moving window used to calculate roughness, and keep these consistent across study areas, when possible.

Anisotropy

Measurements across all sites reveal that mean roughness values in the orthogonal orientation are higher than those in the along-flow orientation across all sites, scales and methods (Figs 4, 5). This pattern is expected because landforms constructed in the subglacial environment will be preferentially oriented in the along-flow orientation, with the resulting changes in landform amplitude yielding higher roughness in the orthogonal orientation as demonstrated by earlier roughness studies (Rippin and others, 2014; Bingham and others, 2017; Falcini and others, 2018, 2021; Cooper and others, 2019). On average, the anisotropy ratio was higher when the local scale was used, especially in sites 1 and 4, where small-scale streamlined landforms are abundant (Figs 1b, e) and their topographic variability is better represented in smaller moving windows. The heavily skewed distribution of roughness values for the FFT method described earlier is also observed here, as the mean anisotropy ratio is 77% higher than the mean ratio from the SD method (Fig. 6, Table 1). Since the SD method is unit preserving, and deemed an appropriate measure for directionality analysis (Rippin and others, 2014; Falcini and others, 2021), only the anisotropy ratio from the SD method will be described henceforth.

Figure 6. Anisotropy values calculated at every intersection point between parallel and orthogonal transects at sites 1 (a) and 3 (b) represent directionality of roughness measurements at the local scale. Orthogonal roughness dominates in the purple shades, parallel dominates in the green. White/gray shades indicate isotropic or random surfaces. White arrows indicate direction of paleo-ice flow.

Table 1. Mean anisotropy from bathymetry and BedMachine (BM)

Sites 1 (Fig. 6a) and 4, where drumlinoid features and MSGLs dominate, are the most anisotropic landscapes, with a mean anisotropy ratio of −0.27 for the two sites. The high anisotropy values observed in site 4, downstream of the convergence between PIG and TG, support the idea that the increased flow velocity resulting from the convergence led to landforms with higher elongation ratios (Nitsche and others, 2013), which correspond to more anisotropic landscapes. A similar value for an area dominated by MSGLs is reported by Falcini and others (2018), though it is important to note that the window size used in their analysis differs from ours and, therefore, results between the studies are not directly comparable. Also worth noting is that artifacts present over the MSGLs in the bathymetry of site 4 bring the overall anisotropy ratio closer to zero by introducing roughness noise in the parallel orientation (Fig. 4d).

While the presence of meltwater channels creates a more rugged topography at site 2 (Fig. 1c), the orientation of streamlined landforms observed is fairly consistent and yields an average anisotropy ratio of −0.13. Alternatively, the mean anisotropy ratio for site 3 approaches zero (−0.05) in both scales and methods used (Fig. 6b). The isotropic nature of this site can be attributed to several factors, including the irregular alignment of streamlined landforms and a flat topographic high where the Pine Island and Thwaites paleo-ice streams merged (Fig. 1d). The presence of large sinuous meltwater channels, thought to have formed by pressurized subglacial meltwater (Lowe and Anderson, 2002; Nitsche and others, 2013), typically lead to random or isotropic landscapes. Alternatively, we observe anisotropic patterns where the lack of channels suggests a dry bed, allowing glacial sedimentary processes to take place. In such areas, it is possible to identify landforms, such as glacial lineations, irrespective of the method or scale used. This highlights the usefulness of this approach for determining patterns of streaming ice flow.

A distinct decrease in roughness values is evident following the transition from crystalline bedrock to sedimentary strata. This change in roughness, coupled with the subsequent increase in anisotropy downstream, attributed to the presence of MSGLs, enables the identification of geological variation with the ASE.

Comparison with BedMachine

As widely used elevation products in ice-sheet models, we evaluate the performance of BedMachine Antarctica in the ASE, and underneath TG, as well as BedMachine Greenland at Jakobshavn Glacier. We specifically assess the impact of using a coarser resolution dataset on the accuracy of roughness results compared to high-resolution data. While BedMachine is a downsampled version of the high-resolution bathymetry used in this analysis, its lower spatial resolution limits its ability to account for roughness derived from small-scale landforms (<500 m in Antarctica and <150 m in Greenland), resulting in misrepresentation of roughness values across all sites. The biggest discrepancies in roughness occur at the local scale, particularly in areas of sharp relief, like the stoss- and lee-sides of streamlined landforms, the steep walls flanking meltwater channels, and where multiple landforms are in close proximity (Fig. 7, S5). The misrepresentation of SD roughness is as high as 150 m around deep meltwater channels and is more evident in the inner- and middle-shelf of PIB where relief is greater than in sites 5 and 6 underneath TG. In sites 1–4, the mean roughness differences between elevation datasets are relatively small, with an average of 0.7 m. Differences are higher in sites 5 and 6 at 8.7 and 3.9 m, respectively (Table S2). In contrast, areas of low relief, such as the sediment-filled basin in site 1, the flat topographic high in site 3 and where MSGLs are present in sites 5–6, show good agreement in roughness values between the two elevation datasets used. The average difference of 32.4 m in SD roughness between BedMachine Greenland and the topographic simulation from Mackie and others (2021) is much greater than what we observe in the ASE, with SD roughness being underestimated by as much as 300 m in the walls of the main channel (Fig. S5). Indicating that the interpolation used to generate elevation datasets from radar measurements along limited track lines results in greater roughness discrepancies, compared to downsampled versions of 2D bathymetric sets.

Figure 7. Difference in roughness measurements of parallel transects between high-resolution topography and BedMachine, showing where BedMachine under- and over-estimates roughness (red and green shades, respectively). Roughness for sites 1 and 3 (a, b) is derived from bathymetry; sites 5 and 6 (c, d) is from swath-radar (Holschuh and others, 2020).

The directionality of roughness is not captured by BedMachine, as evidenced by the noisy distribution of anisotropy values (Fig. S6). Mean anisotropy values for the ASE sites are closer to zero (Table 1), indicating that BedMachine cannot accurately capture the anisotropy of these sites and thus misses roughness deviations in different orientations. Furthermore, while site 3 is in fact isotropic, the anisotropy ratio derived from BedMachine in site 4 is close to zero, despite it being the most anisotropic site of all. These discrepancies highlight the limitations of BedMachine in capturing MSGLs in soft sediments and providing insights into basal ice-sheet flow and organization. The similarity in mean roughness values in sites 1–4, despite the coarser resolution of BedMachine, can be explained by the use of a downsampled version of the ASE bathymetry being incorporated into the International Bathymetric Chart of the Southern Ocean (IBCSO), which was used in the creation of the BedMachine dataset. While the average roughness difference between the bathymetry and BedMachine elevation datasets in the ASE is minimal, there are notable discrepancies around key features known to influence ice dynamics. As a result, BedMachine fails to capture the roughness signature of features that are important in determining basal conditions. The elevation datasets for sites 5 and 6 are not included in BedMachine, therefore any roughness measurements derived from BedMachine, or any other continent-wide elevation dataset, in the modern subglacial environment are subject to interpolation, leading to smooth and unrealistic topography at the scales considered (MacKie and others, 2021). We do not draw any geomorphic interpretations from the results derived from BedMachine, but instead, we highlight the limitations that arise when an elevation dataset fails to resolve landforms that are useful for constraining ice dynamics (Greenwood and others, 2021).

Subglacial conditions inferred from bed roughness

Work by Siegert and others (2005) established that subglacial bed roughness is dependent on four factors: (1) the direction of ice flow, (2) ice dynamics, (3) lithology and (4) geological structure. Transect orientation with respect to ice flow proved to be one of the biggest controls in roughness results across all sites tested in the paleo-ice stream bed of PIG and TG. Glacial sedimentary processes (erosion and/or deposition) play a crucial role in modulating topography in the along-flow orientation. These processes contribute to the creation and modification of landforms that exhibit a distinct alignment with the direction of ice flow. These patterns of preferential alignment reflect the cumulative effects of multiple glaciation cycles known to have taken place in the ASE (Graham and others, 2016), and highlight the significant role of glacial sedimentary processes in shaping the landscape.

Low roughness values (<10 m) are observed in regions of observed sediment cover in sites 1 and 4, as well as on top of the topographic high in site 3, where the two paleo-ice streams merged (Fig. 4). This finding highlights that low roughness does not always signify the presence of thick sediment cover and bed lithology cannot be inferred from roughness alone. Consequently, low roughness values can imply two different lithologies. Firstly, it can indicate the presence of unconsolidated sediment, which reduces basal shear stress and facilitates ice flow. Alternatively, low roughness can also indicate exposed crystalline bedrock, where ice dynamics over such a hard substrate would suggest increased basal shear stress and slow ice flow (Bell and others, 1998; Wellner and others, 2001). Despite this, the presence of widespread grooves around the topographic high indicates efficient erosion (Bennet and Glasser, 1996; Nitsche and others, 2013), suggesting the merging of the paleo-ice streams was enough to overcome the basal shear stress and erode the bedrock to create a smooth terrain. In this case, the increase in ice flow velocity is attributed to the merging of the ice streams, rather than the transition from crystalline bedrock to a sedimentary substrate further downstream. The merging of the ice streams, combined with the network of subglacial meltwater channels upstream, would have increased the supply of basal meltwater, reducing skin drag and facilitating ice flow. The association of low roughness values with both fast and slow ice flow suggests that skin drag, particularly influenced by the availability of basal meltwater in this context, exerts a direct influence on ice dynamics. These observations suggest that form drag alone should not be a key determinant in the sliding law, emphasizing the importance of understanding the complex interaction between basal meltwater, sediment properties and ice-flow behavior.

The boundary of geological structures, such as bedrock protrusions in site 1 (Fig. 4a), subglacial meltwater channels in sites 2 and 3 (Figs 4b, c) and the boundary between crystalline bedrock and unconsolidated sediments in site 4 (Fig. 4d) all cause a roughness spike and change in anisotropy in both methods and scales tested (Fig. 6). The complex topography of PIB suggests two possible explanations for these roughness spikes. First, the increase in basal shear stress associated with rugged topography and the presence of bedforms is reflected in the increase of roughness values. Second, roughness spikes of even greater magnitude are predominantly associated with the presence of subglacial meltwater channels. These large channels would have lubricated the bed in certain instances (Nitsche and others, 2013), resulting in enhanced ice flow rather than increased basal shear stress. The formation of these channels is believed to be linked to episodic outburst floods from subglacial lakes during previous glacial periods (Kirkham and others, 2019). It is inferred that the largest and deepest channels likely developed over multiple glaciation cycles and required abundant meltwater. Specifically, the highest roughness values observed in PIB can be attributed to large episodic drainage events and it is likely that these channels increased in size progressively, leading to increases in elevated roughness values.

Roughness analysis from swath radar data shows similarities between the geologies of sites 1–4 in PIB and sites 5–6 underneath TG. The widespread streamlining observed in the bed underneath TG, and their associated SD roughness values (<5 m) are analogous to the thick sediment pockets present in site 1 and the MSGLs present in site 4. The elongated bedforms observed in site 5, described as the tails of crag-and-tails (Alley and others, 2021; Fig. 1f), exhibit the same geomorphic characteristics as the crag-and-tails present in PIB at the transition between crystalline bedrock and unconsolidated sediment (Graham and others, 2016; Fig. 1e). In the parallel orientation, both sets of crag-and-tails have comparable SD roughness values. Upstream, the bedrock knobs exhibit roughness values of 20–25 m, whereas their downstream soft till tails display values of <5 m. Underneath TG, these tails have larger amplitudes (50–100 m) than the MSGLs offshore at site 4 (5–25 m). As a result, orthogonal roughness values for these features are 20–25 m greater than their PIB counterparts. Notably, the MSGLs in PIB have been in an open marine environment for at least about 10 cal. ka BP (Hillenbrand and others, 2013), where post-glacial sedimentation might have reduced their amplitude.

Conclusion

We quantified bed roughness at six different sites: four deglaciated sites offshore of Pine Island and Thwaites glaciers and two glaciated sites underneath contemporary Thwaites Glacier. By measuring roughness in different orientations relative to ice flow, using different methods and elevation datasets, and applying different detrending scales, we assess how various parameters influence roughness results. Transects in the orthogonal orientation consistently yield higher roughness values, the trends of which are obscured when using lower-resolution elevation products. The choice of scale at which roughness is assessed has a significant impact on the resulting roughness values and therefore requires careful consideration. Overall, the SD method provides a robust representation of bed roughness in several ways. The results obtained from the SD method accurately identify spatial patterns of roughness and anisotropy indicative of ice streaming. Additionally, the unit-preserving nature of the SD method allows for more reliable comparisons between different scales and locations, making it a useful tool for assessing bed roughness in deglaciated environments. The limitations of low-resolution topography are more apparent in the sites underneath contemporary Thwaites Glacier, compared to sites in PIB, suggesting that interpretations derived from bed roughness at ice-stream beds may not be entirely reliable and these uncertainties must be considered in any modeling work.