Basal roughness of the East Antarctic Ice Sheet in relation to flow speed and basal thermal state

Abstract

Basal motion of ice sheets depends in part on the roughness and material properties of the subglacial bed and the occurrence of water. To date, basal motion represents one of the largest uncertainties in ice-flow models. It is that component of the total flow velocity that can change most rapidly and can, therefore, facilitate rapid variations in dynamic behaviour. In this study, we investigate the subglacial properties of the East Antarctic Ice Sheet by statistically analysing the roughness of the bed topography, inferred from radio-echo sounding measurements. We analyse two sets of roughness parameters, one derived in the spatial and the other in the spectral domain, with two roughness parameters each. This enables us to compare the suitability of the four roughness parameters to classify the subglacial landscapes below the ice sheet. We further investigate the relationship of the roughness parameters with observed surface flow velocity and modelled basal temperatures of the ice sheet. We find that one of the roughness parameters, the Hurst exponent derived in the spatial domain, coincides with the thermal condition at the base of the ice sheet for slow flow velocities and varies with flow velocity.

Introduction

Basal motion, which comprises sliding of the ice over the bed and deformation of the bed itself, is a key factor in the control of ice-flow velocities and is thus recognized for its importance in the mass balance of the East Antarctic Ice Sheet (EAIS) (e.g. Shepherd and others, 2006; Ritz and others, 2015; Podolskiy and Walter, 2016). Nevertheless, it is one of the largest uncertainties in ice-flow models and projections of the ice-sheet response to a warming climate (Brondex and others, 2019). Rignot and Mouginot (2012) demonstrated that sliding dominates over the majority of the Greenland Ice Sheet and is the driver for motion already at short distances from ice divides. Sergienko and others (2014) found that sliding occurs over substantially larger areas than previously thought on both ice sheets, exemplifying the importance of this process.

Sliding is a poorly understood process, as the bed is basically inaccessible for in situ observations. The theoretical foundation of today's sliding relations goes back to Weertman (1957) and Lliboutry (1968). Basal velocity is suggested to depend on small-scale roughness of the bed topography (O(~0.05–10 m)) (Weertman, 1957), subglacial hydrology (water pressure, bed separation; e.g. Bindschadler, 1983) and its lithology (Kamb, 1970; Nye, 1970; Winsborrow and others, 2010). Li and others (2010) showed that spectral roughness parameters of the topography can be used to deduce the formation (marine or continental) and erosion history of the respective landscape. This allows us to draw inferences about the bed lithology from measured bed roughness. Furthermore, they propose a sliding relation based on the spectral roughness. Alternatively, spatial roughness, based on self-affine fractal surfaces, has been suggested by Persson (2018) to determine basal friction.

Given that sliding only occurs at a temperate base, Jordan and others (2017) compared the predicted basal thermal state of the Greenland Ice Sheet (as determined by MacGregor and others (2016)) with bed roughness, noting statistical differences between temperate and frozen beds. All influencing factors of basal motion are linked by the roughness parameters, which can be obtained from radio-echo sounding (RES) measurements on regional to continental scales.

Herein we apply two commonly used, but methodologically different derived sets of roughness parameters to an extensive collection of RES data in East Antarctica. The two sets consist of two parameters each. From our analysis in the spectral domain, we obtain two parameters:

1. (1) The vertical roughness ξ, which quantifies the vertical variation of the basal topography; and

2. (2) The parameter η, which is a measure for the horizontal variation of the basal topography, with high values corresponding to a dominance of long wavelengths.

For the same shape of the spectral power density, larger values of ξ indicate larger amplitudes than smaller values for the same value of η for horizontal variation. In contrast, for the same vertical amplitude ξ, the power spectrum with short-scale variations in roughness will contain larger wavenumbers (and thus smaller wavelengths) and will yield smaller values for η than a more smoothly varying surface with missing larger wavenumbers and thus larger η (c. f. Fig. 1 of Li and others, 2010).

From the analysis in the spatial domain, we derive a second set of two parameters:

1. (3) The RMS deviation ν, also known as the Allen deviation, which is a measure of the degree of stationarity of the spatial roughness. For a rough surface with large amplitudes, RMS deviation will yield a large value, whereas a smooth surface will yield a small value. Its interpretation is thus comparable to that of the vertical roughness ξ.

2. (4) The Hurst exponent H, which is a measure of the stationarity of the roughness characteristics as a function of length scale or its degree of fractality, with values between 0 and 1. A Hurst exponent H = 1 indicates a self-similar surface, where the roughness amplitude (e.g. RMS deviation) grows at the same rate in the vertical and horizontal dimension (Shepard and others, 2001).

Tentatively, one can consider that H is a simplified measure for the wavelength distribution at different length scales. For instance, a smaller value of H indicates a comparatively smoother surface than a larger H, as for a small H, the vertical roughness, e.g. the RMS deviation, does not increase as quickly as the spatial length scale considered. Examples can be found in Jordan and others (2017, their Fig. 1). A direct physical example for interpreting roughness values is given by Cooper and others (2019): for an ice stream with a soft, sedimentary bed, a high flow velocity will lead to a smooth bed in the flow direction, thus small roughness amplitude and dominance of long wavelengths, whereas in the across-flow direction, the formation of mega-scale glacial lineations would enhance both the roughness amplitude as well as the dominance of shorter wavelengths.

Methodological details regarding all four parameters are presented in the following section. Afterwards, we discuss the comparison of the large-scale features visible in each set of parameters and relate them to the glaciological and geological settings. In addition to the large-scale overview, we discuss several regions in more detail and investigate spatially anisotropic properties of the roughness parameters, i.e. the dependence of results on the direction of RES profiles with respect to the local direction of surface velocity. To evaluate the potential suitability of each of the four roughness parameters as proxies for basal conditions (i.e. frozen or thawed), we moreover consider the variation of the roughness parameters as functions of basal temperature T and surface velocity v, where we take T from model-based estimates and v from satellite observations. We find that only H indicates a systematic variation with v for low velocities, but the uncertainties are too large to allow conclusive statements.

Data and methods

In this section, we first describe the process to obtain equidistant elevation models from RES observations of ice thickness. We then present our approach to derive two different sets of characteristics of bed roughness, i.e. based on spectral and spatial analyses.

Bed topography

We use a compilation of RES data from more than 20 years of surveying (1994–2017) by the Alfred Wegener Institute (AWI), complemented by three other surveys. The source signal of the AWI radar system is a burst wavelet with a carrier frequency of 150 MHz. The system can be operated in toggle mode, where it is alternating between shots with 60 and 600 ns burst lengths (Nixdorf and others, 1999). The vertical resolution of each of these bursts in ice is ~5 and 50 m, respectively. The data are recorded with a vertical sampling interval of 13.33 ns and are tenfold stacked horizontally to an average trace distance of 75 m. The bed reflection is usually traced and digitized in the differentiated and band-pass-filtered 600 ns data.

After one iteration of manually picking the bed event in a subset of the differentiated data, an auto-picker routine is applied to obtain the maximum horizontal along-track resolution allowed for by the data. The algorithm detects the maximum signal amplitude (i.e. strongest gradient for differentiated data) within a certain time window around the manual picks and also detects the bed reflection between the manual picks, wherever there are manual picks within a distance of 30 traces. The ice thickness d at the position x along the RES profile is then

(1)$${\rm d}\lpar x \rpar = \displaystyle{{c_{{\rm ice}}t\lpar x \rpar } \over 2}\comma \;$$

with the two-way travel time (TWT) t(x) of the radar signal through the ice. We use a radar-wave speed in ice of c _ice = 168.0 m μs⁻¹. The ice thickness d(x) is subtracted from the Bedmap2 surface elevation (Fretwell and others, 2013) at the respective geographic location to get the bed elevation Z(x) (referring to WGS84 vertical datum). Since only relative elevation differences are important for our method (see Eqn (8) below), we can neglect the higher wave velocity in firn and do not apply a firn correction. The justifying assumption to do so is that the firn thickness does not substantially change within 1 km, which is the maximum point distance we use for our roughness calculations. For the same reason, the uncertainty in c _ice does not affect our results, as c _ice, averaged over the ice thickness, does not change over a horizontal distance of 1 km in the regions of interest. The uncertainty in bed-elevation differences thus depends on the uncertainty in surface elevation, the vertical resolution of the radar system and picking accuracy. To assess the accuracy of the derived bed elevation, we conduct cross-over analyses at locations of crossing RES profiles. The RMS error on the elevation differences, which are equal to the difference in ice thickness Δd at cross-over points, for the AWI data is

(2)$$\sigma = \sqrt {\displaystyle{1 \over N}\mathop {\mathop \sum \nolimits}\limits^N_{i = 1} {\lpar { \left( \Delta {\rm d} \right)_i-\overline {\Delta {\rm d}} } \rpar }^2} = 45\;{\rm m}\comma \;$$

with a mean elevation difference $\overline {\Delta {\rm d}} = 56\,{\rm m}$ from N = 3441 measurements.

The characteristics of all RES data are summarized in Table 1. Additionally to the AWI data, we recalculate the roughness parameters with our method for a number of publicly available RES datasets. The spatial distribution of the different datasets is shown by the coloured lines in Figure 1. The data we incorporate in our evaluation are:

1. (1) The final ice-thickness product of the Antarctica's Gamburtsev Province Project (AGAP)-South data, collected by the Lamont-Doherty Earth Observatory (LDEO) (Bell and others, 2011; Wolovick and others, 2013). Data processing uses c _ice = 169 m μs⁻¹ for depth conversion and the bed was picked with a semi-automatic system using the steepest vertical gradient. The output was low-pass filtered along track at 60 m (Wrona and others, 2018).

2. (2) Bed picks of the AGAP-North data (Antarctica's Gamburtsev Province Project – North; Ferraccioli and others, 2011), collected by the British Antarctic Survey (BAS) with the PASIN (Polarimetric Airborne Survey INstrument) radar system and c _ice = 168 m μs⁻¹.

3. (3) Operation IceBridge (OIB) HiCARS 1 L2 Geolocated Ice Thickness product (Blankenship and others, 2011), from data collected as part of the ICECAP (International Collaborative Exploration of the Cryosphere through Airborne Profiling, Young and others, 2011) project. TWT-to-depth conversion was conducted with a wave speed of c _ice = 168.5 m μs⁻¹.

4. (4) Ice thickness from the Old Ice A site (OIA) survey around Dome C in 2016. This survey was part of the US–Australian ICECAP follow-up project ICECAP II (Young and others, 2017). A wave speed of c _ice = 169 m μs⁻¹ was used for depth conversion.

Table 1. Characteristics of the analysed RES data

Range resolution and horizontal sampling distance are given for the processed data. RMS error σ is from cross-over analysis of the ice thickness, individually determined for each dataset.

a After stacking.

Fig. 1. East Antarctica with place names of interest, surface elevation (1000 m interval contour lines in black) and used RES datasets. Different colours show data from different institutes/different datasets, as labelled. Antarctic shapefiles from Depoorter and others (2013).

Our intention is to achieve a consistent processing flow for all datasets, covering a large portion of the EAIS when combined. Some of these data have already been interpreted with respect to basal roughness and radar scattering behaviour (e.g. Young and others, 2011, 2016, 2017; Bellot, 2012; Wright and others, 2012; Creyts and others, 2014). However, we re-interpret them as part of the combined product. Note that the slight differences in applied radar wave speeds and firn corrections are insignificant for this study, as each dataset was self-consistently analysed and only variations in bedrock elevation matter. The systems had slightly different range resolutions and thus uncertainties. We performed a cross-over analysis for each dataset individually, with the resulting RMS values σ listed in Table 1. The RMS cross-over values include a contribution from the system uncertainty as well as from unaccounted physical processes, e.g. a spatial anisotropy of bedrock, leading to systematically different values in one flow direction compared to another (c. f. Cooper and others, 2019). The RMS cross-over values for all systems are of the same order of magnitude (ranging between 45 and 80 m), although they were recorded under different conditions in different regions. We therefore consider the system-specific differences a minor issue and, thus, a joint comparison feasible. The local elevation error within each profile (i.e. set of elevations used to calculate each roughness metric) is assumed to be significantly less than the cross-over error and neglected in the roughness computations. Where subglacial lakes are present, it is not possible to identify the surface of the underlying bedrock from radar data. In those cases, our analysis employs the thickness of the ice column and the elevation of the ice base, respectively. Nevertheless, we use the term bedrock elevation throughout the paper.

Spectral roughness parameters

We follow the methodology of Taylor and others (2004) to calculate the spectral roughness parameters and use a Fast Fourier Transform (FFT) to convert bed elevations into wavelength spectra. As the FFT requires a continuous profile with a fixed spatial sampling interval, we first interpolate the data with an interpolation equidistance of 200 m. We do not interpolate between points with more than 1 km distance, but regard these as data gaps. The calculation of each individual roughness value is based on a predefined number of elevation points, contained in a moving window (MW). For the spectral analysis, the MW comprises 2ⁿ elevation records. Taylor and others (2004) suggested that n should equal 5. We, therefore, use MWs of 32 bed elevation entries which results in MWs of 6.4 km length. Line segments with fewer entries are excluded from the processing.

Before passing an elevation profile Z(x) of an MW to the FFT, it is linearly detrended by subtracting the best fit line (Z _bf) from the raw data,

(3)$$Z_0\lpar x \rpar = Z\lpar x \rpar -Z_{{\rm bf}}\lpar x \rpar .$$

The detrending of the profiles prevents an over-estimation of long-wavelength roughness, when, e.g., only the rise of a hill is included in the window, while the small-scale roughness remains mostly unaffected (Shepard and others, 2001; Siegert and others, 2005).

Li and others (2010) showed that, although calculating roughness based on the FFT of elevation profiles is a useful way of describing roughness, it might not be sufficient as it only accounts for vertical irregularities in a surface while missing the horizontal ones. They proposed a two-parameter roughness, using the profile of the local slope profile ${Z}^{\prime}_{{\rm sl}}\lpar x \rpar$ of the detrended elevation profile as a link between the horizontal and the vertical undulations,

(4)$$Z_{{\rm sl}}\lpar x \rpar = {Z}^{\prime}_0\lpar x \rpar = \displaystyle{{\partial Z_0\lpar x \rpar } \over {\partial x}}.$$

The first roughness parameter ξ is then defined as

(5)$$\xi ^{\lpar l \rpar } = \mathop \int \nolimits_{k_1}^{k_2} S_0\lpar k \rpar {\rm d}k\comma \;$$

with the wavenumbers k, wavenumber limits k ₁ and k ₂ of the spectral band under consideration, and the spectral power density

(6)$$S_0\lpar k \rpar = \displaystyle{1 \over l}\vert \tilde{Z}_0\lpar k \rpar \vert ^2\comma \;$$

where $\tilde{Z}_0\lpar k \rpar$ is the Fourier transform of a detrended elevation profile Z ₀(x) of a given length l of the MW. (In the following, we omit the denominator (l) for the sake of clarity, as the length will be mentioned in the context.)

The second roughness parameter η is defined as the ratio of ξ and the corresponding index ξ _sl, calculated with the spectral power density S _sl of the Fourier transform $\tilde{Z}_{{\rm sl}}$ of the slope profile,

(7)$$\eta = \displaystyle{\xi \over {\xi _{{\rm sl}}}} = \displaystyle{{\mathop \int \nolimits_{k_1}^{k_2} S_0\lpar k \rpar {\rm d}k} \over {\mathop \int \nolimits_{k_1}^{k_2} S_{{\rm sl}}\lpar k \rpar {\rm d}k}}.$$

Whereas the parameter ξ reflects the vertical amplitudes of the profile (i.e. the energy in a certain wavenumber band), η represents its horizontal variation, with high values corresponding to a relative dominance of long wavelengths.

Spatial roughness parameters

We follow the methodology described in Jordan and others (2017), which was adapted from Shepard and others (2001), to calculate the spatial roughness parameters ν(Δx) and H, the Hurst exponent. We derive both parameters statistically from the detrended bed elevation profiles of 10 km MWs with a 1 km along-track sampling. All pairs of points in each window are binned by their distance (lag) Δx. We choose a width of 100 m, so that the smallest bin still contains points for the data with the coarsest horizontal sampling (AWI data with average trace distance of 75 m) and a maximum lag of 1 km. Despite the different horizontal sampling of the datasets, we use the same parameters for all data. The RMS deviation ν(Δx) is calculated for each bin,

(8)$$\nu \lpar {\Delta x} \rpar = \sqrt {\displaystyle{1 \over N}\mathop {\mathop \sum \nolimits}\limits^N_{i = 1} {\lpar {Z_0\lpar x_i \rpar -Z_0\lpar {x_i + \Delta x} \rpar } \rpar }^2} .$$

The elevation differences of all pairs of points (N: number of pairs) in the 10 km window, which have a lag Δx within one bin (e.g. 1–100 m, 101–200 m, etc.), are included in the calculation of ν of this specific lag bin. Results for bins with five or fewer pairs of points are discarded.

The results for ν are then used to calculate the Hurst exponent H for each 10 km window. H is the slope of ν over lag distance in a log–log plot (Shepard and others, 2001). We use all ten lags, i.e. the results for ν from the 50 (1–100 m) to 950 m (901–1000 m) bins, to derive the slope,

(9)$$\nu \lpar {\Delta x} \rpar = \nu \lpar {\Delta x_0} \rpar \left({\displaystyle{{\Delta x} \over {\Delta x_0}}} \right)^H\comma \;$$

with a reference horizontal lag Δx ₀. We only use the H of windows containing more than five measuring points of ν in all ten bins. By investigating random samples of the dataset, we found that the breaking point in H as a function of lag only occurred for lags >10 km. We can therefore safely assume that the variogram is not skewed and that we have a sufficiently large sample size to estimate the variogram exponent H (Shepard and others, 2001).

Given the above definition, it has to be emphasized that ξ, η and H present roughness parameters integrated over a range of wavelengths. In contrast, ν provides just a snapshot for the intercept of the variogram at one wavelength bin, the utility of which will be evaluated in the next section. Nevertheless, this has to be kept in mind when comparing the spatial distributions of all roughness parameters. The range of raw values and basic statistics for each roughness parameter are listed in Table 2. Theoretically, the values of H should be in the interval [0, 1], which is the case for 95% of the values resulting from our analysis. However, the application of our parameter estimation algorithm to our data produces a total of 5% of values outside of this plausible range. These will be discarded in the later analysis and interpretation.

Table 2. Range of roughness parameters

For ν we omitted outliers with values >150 m. Statistical properties for H include 5% of implausible values outside the theoretical range interval [0, 1].

For the sake of comparibility, we normalize all roughness parameters to zero mean and unity std. dev. Histograms of the normalized spectral and spatial sets for all four roughness parameters are shown in Figure 2.

Fig. 2. Normalized value distribution (i.e. zero mean and unity std. dev.) of the four roughness parameters (a) ξ, (b) η, (c) ν and (d) H for the complete dataset. Less than a fraction of 10⁻⁶ of the data is outside the range of five std. dev. Basic statistical values of the distributions are shown in Table 2.

Results

We present the two sets of roughness parameters defined above (spectral: (ξ, η), spatial: (ν, H)) for comparison in Figure 3. We first indicate differences between the two sets before debating possible implications in the ‘Discussion’ section below. For the sake of brevity, we restrict our analysis of spatial anisotropy in relation to flow direction and a possible relation to flow velocity to the spatial roughness values.

Fig. 3. East Antarctica with spatial distribution of the normalized roughness parameters (a) ξ, (b) η, (c) ν for the 250 m bin (Δx = [201 − 300] m) and (d) H. Black downward pointing triangles mark Dome Argus (DA) and the drill sites Dome Concordia (DC), Vostok, South Pole (SP), Dome Fuji (DF) and EPICA Dronning Maud Land (EDML) (c. f. Fig. 1). Coordinates in this and all subsequent figures are polar stereographic.

Roughness parameter distribution

Spectral roughness parameters

The areal distribution of the ξ roughness is shown in Figure 3a. The geographic allocation of derived basal roughness varies widely within the study area and is subject to the availability of recorded radar profiles. Distinct tendencies in basal roughness can more easily be identified in areas with higher densities of tracks, which also allows us to delineate large centres of rough and smooth terrain.

A large area of ~200 × 200 km² with relatively low roughness can be found south-west of Kohnen research station (EPICA Dronning Maud Land (EDML) ice core drill site) (the dominating bright cyan region in Fig. 3a). Additional continuous low roughness can be observed on the coastal ridges of Halvfarryggen and Søråsen, south-east and south-west of Neumayer III station (refer to Fig. 1 for locations).

Extensive, rather homogeneous areas can be observed on the East Antarctic plateau. Especially areas south of the Maudheimvidda mountain range show roughness values in this medium range of normalized ξ roughness values. Similar characteristics over smaller extents can be found in coastal areas bordering Fimbulisen and Riiser-Larsenisen and areas south of Halvfarryggen and Søråsen. Relatively high roughness values are generally more scattered than medium and low roughness values and can predominately be observed around the mountainous areas separating the coastal areas from the East Antarctic plateau. High ξ values are found especially in close proximity to bare rock outcrops distributed over the mountain range. Large consistent regions with high roughness values can be identified south-east of the research station Princess Elisabeth as well as in Maudheimvidda. The rough spot directly neighbouring Princess Elisabeth covers an area of nearly 500 km² and is exposed to high ξ values between 2 and 3 std. dev. above the mean. The region with the highest roughness values is the Gamburtsev Mountains, with maximum ξ values distributed around Dome A.

Figure 3b shows the geographic distribution of η values. It shows more spatially varying values than for ξ. Larger areas of low or high η values cannot be delineated. The largest continuous area of low η values can be found on the ridges of Halvfarryggen and Søråsen and within a 50 km wide coastal stretch in-between both ridges. The large plain south-west of Kohnen research station, previously associated with large homogeneous patches of low ξ values, does not show homogeneous values for η. Only the most southern parts of the 200 × 200 km² stretch show low and rather uniform values for η. In addition, a selective distribution of low η values can be observed in littoral zones along the central and very western part of Dronning Maud Land. Other spread-out homogeneous areas can be identified within Amundsenisen, Wegenerisen, Nansenisen, Gamburtsev Mountains and zones within Maudheimvidda. Areas characterized by a homogeneous and extensive presence of high η values, indicating long-wavelength roughness, cannot be readily observed, except for some selective points, predominately distributed along the Fimbulheimen mountain range.

Spatial roughness parameters

We derive ν for ten different lag bins for the combined dataset. The spatial distribution of the 250 m lag results (lags in the range of 201–300 m) is shown in Figure 3c. The pattern of high and low ν values is similar for all lags, but the absolute values increase with increasing lag. For the sake of brevity, we therefore restrict further discussion to the results for the 250 m lag bin. Regions with different roughness are readily identifiable. The Dronning-Maud uplands, separating the coastal area from the plateau, show comparatively high ν values. The region with the largest vertical roughness amplitudes is the Gamburtsev-Mountain province around Dome A. The region south-west of Kohnen station, characterized by large homogeneous patches of low ξ values, is also reflected as very homogeneous in ν roughness compared to the surrounding regions.

The spatial distribution of the Hurst exponent is shown in Figure 3d. The region south-west of Kohnen again stands out as having very low values of H, at least partly. The Gamburtsev Mountains and the mountain range in Dronning Maud Land show high H values. Overall, the spatial distribution of H values is smoother than its roughness-amplitude value ν. Compared to η, the spectral roughness-wavelength proxy, whose distribution is rather homogeneous, H shows more variation and seems better suited for recognising different subglacial landscape types on the continental scale. We will come back to this comparison below.

Specific regions of interests

We now have a closer look at the roughness parameters in selected regions of interest (ROIs) of specific glaciological characteristics. The first ROI is the Jutulstraumen and especially the onset of faster ice flow of this ice stream. The ice flow in this region is concentrated in a deep trough constrained on both sides by the Dronning Maud Land mountain range. The roughness parameter distributions in this area and contours of the flow velocities are shown in Figure 4. Considerably smaller values of the roughness-amplitude parameters, ξ and ν (Figs 4a and c), can be observed inside the ice stream and catchment areas upstream, compared to the regions outside of the ice stream. The contrast is especially strong for the ν values when going from the west into the ice stream, where the topography changes from the subglacial mountains to the deep trough of Jutulstraumen and all the inland ice draining into Fimbulisen ice shelf is channelized. For the roughness-wavelength parameters, η and H (Figs 4b and d), such a differentiation of regions with slower and faster ice flow is not obvious. If any pattern is visible at all, then it is a slight gradient from shorter to longer wavelengths from the grounding line towards the interior. When looking at the cross-over differences of ν and H in the same region (Fig. 5), the differences are in general slightly larger within the ice stream and the catchment basin in the south-east than outside of the main ice stream, especially for H. This tentatively indicates a dependence on profile direction for both the amplitude and wavelength of the roughness in Jutulstraumen, i.e. a spatial anisotropy of roughness parameters. At a few cross-over points outside the stream and catchment basin, larger directional differences are visible as well, especially again for the H values. However, these differences are less pronounced than those in the region of faster ice flow.

Fig. 4. The four normalized roughness parameters (a) ξ, (b) η, (c) ν and (d) H in the Jutulstraumen region. Contour lines indicate surface flow velocity, inset in (b) indicates surface speed from the MeASUREs dataset for the same region (Greene, 2019), with increasing speed with brightness of greyscale.

Fig. 5. Absolute differences for original values (i.e. not normalized distributions) of (a) ν (for the 250 m lag bin) and (b) H at profile cross-over points in the Jutulstraumen region. Larger differences indicate more pronounced spatial anisotropy of roughness parameters. The flow velocities can be seen in the inset of Figure 4b. In the central part of Jutulstraumen, flow direction is roughly from south (bottom) to north (top).

As a second ROI, we examine the Gamburtsev Subglacial Mountains. This well-preserved mountain range in the centre of East Antarctica is especially interesting for its rough terrain and steep mountain slopes, with large roughness values, which already attracted attention in the overall spatial distributions of roughness parameters presented above (Fig. 3). The distributions of the four roughness parameters in the Gamburtsev Mountain region are shown in Figure 6. As expected, we see high amplitude-roughness values ξ and ν (Figs 6a and c) in this region, with maximum values in a long south-west–north-east stretch, passing slightly west of Dome A. This is the central part of the subglacial mountain range with the highest peaks of ~2500 m bed elevation. North of the Gamburtsev Mountains, at the northernmost part of the AGAP radar survey, roughness-amplitude parameters show overall smaller values. This is the onset area of Lambert Glacier, which drains to the Amery Ice Shelf. The distribution of the roughness-wavelength parameters in the Gamburtsev Mountain region (Figs 6b and d) is more homogeneous and less pronounced than the roughness-amplitude proxies. The normalized Hurst exponent H shows relatively high values, which are close to 1 for the original (unnormalized) values, indicating characteristics close to self-similarity. The locations of the highest mountain peaks and highest roughness amplitudes show up as bands of slightly smaller roughness-wavelength values of η than in the surrounding regions. This indicates a dominance of smaller wavelengths in the roughest parts of the mountains. In the northern part with faster ice flow, roughness-wavelength values also get smaller and spatially less variable, although the change is less pronounced than for the roughness-amplitude values.

Fig. 6. The four normalized roughness parameters (a) ξ, (b) η, (c) ν and (d) H in the Gamburtsev Mountain region. The black triangle in the centre of the radar-survey grid marks Dome A. The triangle at the bottom marks Vostok.

The third ROI investigated in more detail is the large patch of homogeneously small values of ξ, ν and partly H, south-west of Kohnen, already mentioned in the presentation of the large-scale overview (see Fig. 3). As we have seen before, both the two roughness-amplitude and the two roughness-wavelength parameters show comparable distributions, but the spatial parameters seem better suited to identify differences by the more pronounced spatial homogeneous characteristics of H compared to η. The close-up confirms the homogeneously small values for the roughness-amplitude ν (Fig. 7c) in the whole region, which was interpreted as a geologically soft basin (Li and others, 2010).

Fig. 7. The four normalized roughness parameters (a) ξ, (b) η, (c) ν and (d) H in the region south-west of Kohnen station (black triangle). Contour lines indicate flow velocity, inset in (b) indicates surface velocity (Greene, 2019) for the same region (velocity increasing with the brightness of greyscale). The onset of Slessor Glacier is visible to the south-west, towards the lower left corner.

To the north of this region, the roughness amplitudes ξ and ν get larger, indicating the foothills of the Maudheimvidda Mountains. The distribution of H (Fig. 7b) shows medium-to-low values at the margin of Maudheimvidda and in the northern part of the basin. Only in the southern half of the basin towards the onset of the northern Slessor Glacier tributary do H values get very small, indicating the dominance of small-wavelength roughness in this area. It has to be noted that all profiles in this region were rather perpendicular to flow direction, which impedes the investigation of the directional dependence of roughness in this region of ice-stream onset. Thus, there are too few radar lines crossing each other to provide meaningful results in terms of spatial anisotropy of roughness.

Discussion

In general, our results indicate that the two parameters ξ and ν, which can be considered a measure for the vertical amplitude of roughness, yield easily depictable spatial patterns. One could assume that they thus also provide more meaningful results than the two parameters η and H, which rather inform about the wavelength of the roughness. In the following, we discuss the comparison and implications in more detail.

Comparison of spectral and spatial roughness parameters

Previous literature on the relation of space-domain and spectral-domain roughness parameters in the context of glaciology objectives was briefly summarized and discussed by Jordan and others (2017). Although the space-domain variogram and deviogram have an approximate correspondence to the frequency-domain power spectrum, Shepard and others (1995) recommended to estimate the Hurst exponent H for roughness analyses from variograms, as they yield supposedly less noisy results and are less likely to bias slope estimates than the power spectrum method (see also Jordan and others, 2017, p. 1250).

Our histograms of normalized spectral versus spatial roughness values (Figs 2a and b vs Figs 2c and d) show significant differences. The distributions of spatial values are more asymmetric and skewed than those of spectral values. According to Jordan and others (2017) and references therein, spatial parameters are better suited to reveal self-affine scaling behaviour. In the context of applying RES methods on ice sheets, H is also related to radar scattering behaviour, allowing for inferences not only about topography but also small-scale (sub-radar resolution) roughness and condition of the base. Whether this is also applicable for the regions in East Antarctica covered by the RES data available for our study will be investigated next.

Roughness, temperature and flow velocity

For application in ice-flow models, a physically justified parametrization of the subglacial roughness is needed to reliably calculate the basal drag and its influence on flow velocity. To investigate the relationship between roughness and flow velocity, we compare our roughness results with the MEaSUREs ice-flow velocities (Mouginot and others, 2017; Greene, 2019). Those velocities were compiled from different satellite missions and sensors and have very variable precisions. Reference errors range between 3 and 35 m a⁻¹, depending on satellite, orbit and acquisition mode (Mouginot and others, 2017). For this comparison, we bin our data by basal temperature, using modelled temperatures from Van Liefferinge and Pattyn (2013). All temperatures we state below are given relative to the pressure melting point (pmp) in units of °C. The relationship of the spatial set of roughness parameters, ν and H, as a function of flow velocity is shown in Figure 8. The data points of these curves are the mean values of all points in the respective temperature and flow-velocity bin. Bin widths are 2 m a⁻¹. The error bars indicate the respective 95% confidence intervals.

Fig. 8. Variation of (a) ν and (b) H with observed surface flow velocity for four different temperature ranges (T given in °C relative to the pressure melting point); (b) also indicates the value at each major ice-core deep drilling site; (c) and (d) are close-ups of (a) and (b) for low flow velocities, v < 50 m a⁻¹. Basal temperatures are from Van Liefferinge and Pattyn (2013). Surface velocities are from Greene (2019) with a reference precision between 3 and 35 m a⁻¹ (Mouginot and others, 2017). Note the different scales on all axes.

Both roughness parameters show strong variations at high flow velocities (v  >  100 m a⁻¹). At flow velocities ${\rm \lesssim }20\;{\rm m}\;{\rm a}^{{-}1}$, ν shows a slightly decreasing trend for all temperature ranges with increasing flow velocity. This means either that the ice can flow faster over a smoother bed (i.e. smaller roughness amplitude), or that the bed is smoothed by the higher erosion rates exerted by faster ice flow (e.g. Rippin and others, 2014); likely, both processes coexist in a positive feedback. The behaviour of H with increasing flow velocity strongly depends on basal temperature. Figure 8d shows that H decreases with increasing flow velocity at cold-bed locations with temperatures below −5^°C relative to pmp. At locations where the bed is thawed or close to thawing ($T{\rm \gtrsim }-5^\circ {\rm C}$ relative to pmp), in contrast, H rather increases with increasing flow velocity.

Furthermore, these curves show that the mean of H systematically decreases with increasing temperature T at locations with flow velocities below 15 m a⁻¹. Although these findings have to be treated with care, one interpretation could be that for low velocities the roughness distribution is more self-affine (large H) for a very cold bed than beds being more close to the pmp. This could be caused by lower cumulative erosion rates over time in those areas which are still very cold today. The rather prominent distinction of H values as such for different temperature bins as well as their different dependency on velocity potentially could indicate that H might be a useful parameter to further constrain the thermal state of the bed, complementary to other methods. However, the velocity uncertainties of MEaSUREs are considerable, making it difficult to determine the significance of our finding. We will focus and elaborate further on this in the next paragraphs.

Does the Hurst exponent indicate thermal conditions at the base?

At low flow velocities (v < 10 − 15 m a⁻¹), H seems to show a systematic variation with flow velocity, depending on the predominant basal temperature range (Fig. 8d). Specifically, our H − v curves potentially hint at a systematic decrease of H with increasing basal temperature at flow velocities below 15 m a⁻¹. Figure 9 shows the behaviour of H with temperature for low (v < 10 m a⁻¹) and high (v ≥ 100 m a⁻¹) flow velocities. The curve for v  <  10 m a⁻¹ confirms that H decreases with increasing temperature, but only for T > −15^°C. For lower temperatures, the relation shows rather the opposite behaviour, but with large uncertainties for very low temperatures. Speculatively, one could consider this an effect of past erosion being much less active at very cold beds over time than for warmer beds, which might have been at pmp in the past and thus subject to more small-scale erosion and thus a more pronounced dominance of longer wavelengths.

Fig. 9. Variation of H with basal temperature (relative to the pmp) for low (v < 10 m a⁻¹) and large (v ≥ 100 m a⁻¹) flow velocities (MeASUREs from Greene (2019)).

For high flow velocities, H seems to increase with temperature, but variation in H (sudden changes of up to 0.25 per 0.3^°C) is too large for a conclusive statement. Such fluctuations in the dataset could either be caused by a spatial anisotropy in regions of faster flow, which was observed in Greenland by Cooper and others (2019), or the regions with faster flow are possibly too different in their characteristics, e.g. due to preglacial topography and erosion history, to exhibit a common roughness pattern. However, given that our estimates of basal temperatures depend on results from simplified ice-flow models, which do not consider all relevant stress components for Antarctic-wide model domains, it is also possible that the temperature estimates are biased or are subject to implausible variations.

Based on the inferred relations, it seems that at locations with low flow velocities, a small Hurst exponent is likely coincident with a temperate bed, a high Hurst exponent likely with a cold bed. This observation applies to larger clusters of low or high Hurst exponents and not to single measurements, because H shows relatively large fluctuations over more regional domains. An application of this hypothesis to our dataset would predict a temperate bed for the region around Kohnen station, parts of the Aurora Subglacial Basin, the southern of the Oldest Ice Reconnaissance survey around Dome Fuji, the region sparsely covered by crossing profiles close to the South Pole, and to the East of Dome C (compare Fig. 10).

Fig. 10. (a) Mean modelled basal temperature of four runs with different geothermal heat flux input (see text for details) with RES profiles overlain (grey). Note that temperatures are given in °C relative to the pressure-melting point. Outside of the ice-covered regions, surface temperature is plotted. (b) Same temperature means, now only plotted at locations with low flow velocity (v < 10 m a⁻¹) and small Hurst exponent (H < 0.4).

Given the nominally large errors of the satellite-based surface velocities of O(10 m a⁻¹), the systematic variation of H(T) for low flow velocities in each temperature bin is intriguing. We therefore test this relation by investigating a different set of modelled basal temperatures and ice-flow velocities at these locations, which are plotted in Figure 10a. The Parallel Ice Sheet Model (Bueler and Brown, 2009; Winkelmann and others, 2011) is applied at 16 km resolution over the entire Antarctic Ice Sheet. Four different geothermal heat flux (GHF) datasets (Shapiro and Ritzwoller, 2004; Fox Maule and others, 2005; An and others, 2015; Martos and others, 2017) are used to derive the present-day temperature field based on 220 ka-long spin-up simulations with paleoclimatic forcing (see model description of ‘AWI_PISM1Pal’ in Seroussi and others, 2019). The ice-sheet geometry, including bedrock elevation, grounding line and ice-shelf front, can evolve freely during the simulations. Details about the paleoclimatic forcing are given in Sutter and others (2019). The initial ice-sheet geometry for the spin-up is also based on Bedmap2 and is refined in the Recovery Glacier area with additional ice thickness datasets from several recent airborne RES campaigns (Leuschen and others, 2010; Ferraccioli and others, 2018; Humbert and others, 2018). The details of the applied interpolation methods are given in Humbert and others (2018). Additional datasets for the Dome C (Young and others, 2017) and Dome Fuji (Karlsson and others, 2018) regions have not been incorporated into the Bedmap2 dataset for the PISM simulations. We use the mean basal temperature (Fig. 10a) of the four PISM simulations to account for uncertainties in GHF.

Limiting the plotted values of temperature in Figure 10b to those regions with low modelled ice-flow velocities, v < 10 m a⁻¹, and small Hurst exponent, H  <  0.4, confirms that in the majority of cases, thermal conditions at the bed coincide with a temperature at or at least close to the pressure-melting point. The only obvious location where the modelled basal temperatures do not match with a small H is south-west of Dome C (including Little Dome C, the candidate site for the Beyond EPICA drilling project). Despite the cluster of low H values, modelled temperatures indicate a cold bed with 6–7^°C below pmp at this location. On the other hand, there are studies indicating a warm bed and a formerly more dynamic ice sheet with subglacial erosion at Dome C (Siegert and others, 2005; Bingham and Siegert, 2009). Furthermore, there is also some indication of water in the same region covered by our data at present day, e.g. borehole temperatures extrapolated to bed suggest basal melt at Dome C (Lefebvre and others, 2008). Apart from the many lakes in this region revealed by earlier RES surveys (e.g. Oswald and Robin, 1973; Tabacco and others, 2006; Wright and Siegert, 2012), Young and others (2017) detected 54 lakes in the area covered by exactly the same RES survey that we also included in our analysis. Specifically, in the Little Dome C region, a rather sharp transition in basal conditions occurs ~2800 m ice thickness, from cold conditions above to temperate below (Beyond EPICA consortium, 2019).

Further investigations with a simple bivariate analysis of H as a function of the two variables T and v, but ignoring the uncertainties of v, tentatively indicate a bi-linear relationship, as visible in Figure 8d for v < 20  m a⁻¹. However, to obtain robust results on a coincident occurrence of low H and temperate beds, a more extensive database, smaller temperature bins, separation of H values for flight direction parallel and perpendicular to flow and in particular more accurate surface velocities seem necessary, as are currently available with the datasets used here. Careful testing should involve measured basal temperatures and extension to larger regions and more comprehensive datasets (e.g. Morlighem and others, 2020), but is beyond the scope of this study.

Roughness in Greenland and Antarctica

A general goal when using roughness as a proxy for physical properties at the base would be to find a common relationship for the Greenland and Antarctic Ice Sheets. A good portion of the Greenland Ice Sheet also has slow flow velocities below 10 m a⁻¹ and the mean and variance of H values are approximately the same as in our cases. When comparing results from Greenland and Antarctica, however, several aspects have to be accounted for: different parametrizations for calculating roughness variables can introduce artificial differences (i.e. non-homogeneous datasets), e.g. the uncertainties of the Antarctic velocity field are usually larger and might only accidently yield a relationship similar to that presented above, or the actual differences in the glacial or subglacial settings of the investigated regions exist and cause a different scattering behaviour.

For Greenland, Jordan and others (2017) showed that H is significantly higher at thawed-bed locations compared to frozen ones. They propose that higher values of H at thawed-bed locations than at frozen ones in their data could be due to a specific hydrological system, characterized by channels and not having so many deep, large lakes, which leads to more diffuse RES-signal scattering and high H values. The follow-on study by Cooper and others (2019) showed that roughness metrics are rather isotropic for small flow velocities, but that the scaling behaviour of roughness shows strongly different results for RES data acquired parallel and perpendicular to faster flowing regions (v > 150 m a⁻¹). They conclude that the roughness is rather a proxy for geologic properties, i.e. low roughness corresponds to a flat and hard – thus less deformable – bed.

The transfer of such findings from Greenland to Antarctica as a whole is not straightforward, as the specific topographic, geologic or glaciological setting of Greenland and subregions are quite different to those of the bulk of regions covered by our data in East Antarctica. The potential absence of channelized drainage systems in Antarctic regions with temperate basal conditions, and a possible replacement with water being collected in large lakes with specular, abrupt reflection signatures, would diminish H. The drainage system in the Aurora Subglacial Basin, for instance, was also identified as consisting of channels and fjord-like structures by Young and others (2011) and Wright and others (2012). Following this chain of reasoning would imply that any conclusions about basal thermal state remain elusive, as long as the characteristics of the subglacial topography and hydrologic system are not constrained further.

Our tentative relationship for the Hurst exponent H as a function of v and T could therefore possibly indicate a relationship between H and flow velocity v, rather than the thermal state at the base. Its somewhat noisy characteristics in Figure 8 for large flow velocities could moreover be partially caused by spatial anisotropy of roughness. For an investigation of the relationship between H and the thermal state at the base, the influence of flow velocity would then need to be decoupled.

Summarizing the open questions from our discussion, we propose that future studies with larger, more comprehensive and accurate datasets in Antarctica, like the BedMachineAntarctica approach (Morlighem and others, 2020) or the upcoming Bedmap3 compilation, need to investigate the following research questions. How can survey-dependent differences and sampling characteristics be minimized? Do smaller bin widths and maximum lags influence the results? Do roughness parameters exhibit a strong anisotropy for regions of faster flow? Are they related to the basal thermal state or rather to flow velocity? Could any of these aspects eventually lead to improved estimates of the basal properties, especially the thermal state?

Conclusions

Due to the inaccessibility of the bed topography, basal roughness is challenging to analyse, particularly over a large geographic extent, such as the EAIS. It has consequently often been under-represented within glaciological models and studies of mass-balance behaviour. We conducted a detailed analysis of basal roughness over large parts of the EAIS. The applied method led to informative results and identified regional variations of roughness between mountainous areas, fast ice streams, smooth (sedimentary) basins and coastal plains. We partially observe spatial anisotropy of roughness, i.e. depending on profile direction, in areas with sufficient data coverage and exposed to larger variations in altitude and faster ice flow, and more spatially isotropic properties in regions with low ice velocities and more consistent terrain characteristics. Nevertheless, a meaningful differentiation between areas of high and low roughness-wavelength parameters proved to be more difficult. The combination of the two roughness parameters ν and H derived in the spatial domain appeared more straightforward to identify regional differences than the spectral parameters ξ and η. However, the spectral amplitude ξ is still sufficiently useful in itself.

For the investigated areas in East Antarctica, we find a weak relationship between roughness parameters, flow velocity and basal temperatures, i.e. a less rough bed topography in regions of faster ice flow. We noticed that, for regionally slow flow, H shows a substantially different behaviour with flow velocity for different basal temperatures. It is systematically smaller when the basal conditions are temperate. However, given the uncertainties in ice-flow velocities, basal temperature estimates and insufficient separation of anisotropy of roughness values, care has to be taken when interpreting clusters of a homogeneous Hurst exponent in terms of the thermal state at the ice-sheet base. Considering the fact that the EAIS is considerably larger and glaciologically and geologically more diverse than the Greenland Ice Sheet, the observations hitherto made for both ice sheets need further careful joint testing, preferably with standardized approaches, ground truthing with direct measurements of basal thermal conditions for verification and separation of the relationship of flow velocity, flow direction and basal temperature with roughness parameters.