A self-adaptive two-parameter method for characterizing roughness of multi-scale subglacial topography

Abstract

During the last few decades, bed-elevation profiles from radar sounders have been used to quantify bed roughness. Various methods have been employed, such as the ‘two-parameter’ technique that considers vertical and slope irregularities in topography, but they struggle to incorporate roughness at multiple spatial scales leading to a breakdown in their depiction of bed roughness where the relief is most complex. In this article, we describe a new algorithm, analogous to wavelet transformations, to quantify the bed roughness at multiple scales. The ‘Self-Adaptive Two-Parameter’ system calculates the roughness of a bed profile using a frequency-domain method, allowing the extraction of three characteristic factors: (1) slope, (2) skewness and (3) coefficient of variation. The multi-scale roughness is derived by weighted-summing of these frequency-related factors. We use idealized bed elevations to initially validate the algorithm, and then actual bed-elevation data are used to compare the new roughness index with other methods. We show the new technique is an effective tool for quantifying bed roughness from radar data, paving the way for improved continental-wide depictions of bed roughness and incorporation of this information into ice flow models.

1. Introduction

The mechanical flow of ice is a product of the balance between the gravitational driving force that acts to move ice downstream and the resistive forces that counter this motion (Paterson, 1994). This counter force is influenced by the nature of the bed, such as thermal regime, material properties of the bed and roughness of the ice-bed interface. The quantitative description of subglacial topography is therefore important to understanding glacial dynamics (Hubbard and Hubbard, 1998; Taylor and others, 2004; Rippin and others, 2006), which is in meter scale roughness that relates to the basal sliding. Theories suggest that ice dynamics (particularly basal flow) are directly related to bed roughness at a range of spatial scales (Paterson, 1994; Joughin and others, 1998; Schoof, 2002), and this is increasingly regarded as an important research area in glaciology (Cooper and others, 2019). Li and others (2010) showed that bed roughness derived from a frequency domain method applied to ice-penetrating radar data can be used to deduce the formation (marine or continental) and erosion history of a variety of end-case landscapes. Recently, bed roughness has also been used to simulate realistic digital elevation models and to locate subglacial lakes (e.g. Mackie and others, 2020; Goff and others, 2014). Clearly, a method that allows bed roughness to be extracted from radar data accurately and efficiently, with maximum information incorporated within the results, would allow it to become better assimilated into assessments of ice dynamics, glacial history and basal conditions (Eisen and others, 2020).

The main methods of quantifying bed roughness from radar data can be divided into two categories: the electromagnetic scattering properties of the bed-echo waveform; and the statistical properties of the along-track topography (Cooper and others, 2019; Eisen and others, 2020).

Analysis of electromagnetic scattering properties involves (1) waveform analysis (e.g. Cooper and others, 2019), (2) analysis of distribution of peak bed-echo power due to scattering loss (e.g. Grima and others, 2014), and (3) specularity analysis using different synthetic aperture sizes (e.g. Schroeder and others, 2013). Both theoretical predictions and observations demonstrate that specular reflections are suppressed in rougher regions (Jordan and others, 2017).

Methods that use statistical properties can be divided into those utilizing the ‘space domain’ and ‘frequency domain’. Space domain methods include some form of functional parameterization, such as the extraction of root-mean-square (rms) height, rms deviation, rms slope and autocorrelation length (Shepard and others, 2001; MacGregor and others, 2013), from which variograms (rms height versus profile length) and deviograms (rms deviation versus horizontal lag) can be constructed, which have been well applied to glacier beds (MacGregor and others, 2013; Jordan and others, 2017). In these methods, bed-roughness results vary with the lengths of the quantified profile or the spatial lag, often following a power-law trend with these lengths (Shepard and others, 2001). In the frequency domain, Hubbard and others (2000) and Taylor and others (2004) defined a single-parameter roughness index as the integral of the spectrum within a specified interval, which can be used to extract the so-called ‘Hurst exponent’ (Martinez and others, 2013) and the corner wavenumber. This single-parameter roughness index has two flaws, however: (1) it cannot be used directly to quantify slope bed irregularities, which leads to topography with different slope irregularities having similar roughness indices; and (2) it requires a fixed-size moving window (MW), which limits the scale of roughness quantified. Addressing the first issue, Li and others (2010) designed a two-parameter roughness quantization method, which we call ‘the previous two-parameter method’ in this paper, that incorporates both vertical change and slope changes. However, this method is still problematic as it requires a fixed-size MW. The outcome is that large-scale roughness results do not incorporate small-scale terrain undulations, and vice versa. To solve this problem, the ‘wavelet transform’ method was developed. This technique uses different scales of MW, which means a variable-size MW, to divide the roughness signal at each position into components related to scale; the wavelet coefficients are the weight of each component (Boggess and Narcowich, 2015). Our method uses a range of scales like the wavelet transform method rather than a single-scale like the previous fast Fourier transform (FFT) methods.

Self-affine roughness (Jordan and others, 2017) is represented by the ‘Hurst exponent’ as a multi-scale roughness parameter that can accommodate a span of scales to a certain extent. In this method, problems occur in relation to ‘breakpoints’, however, where a power law appears to ‘break over’ and obey a different power law. This requires the investigator to subjectively choose where the beginning and endpoints of each trend occur when a breakpoint is evident (Shepard and others, 2001). Besides that, the ‘Hurst exponent’ results show rapid local fluctuations on maps (Jordan and others, 2017; Eisen and others, 2020), which make it difficult to quantify the large span of terrains.

In this article, we describe a new algorithm to quantify bed roughness, improving the previous two-parameter method (Li and others, 2010) by accounting for different scales of roughness through a self-adaptive approach. In addition, the ‘scale’ we use in this paper is the horizontal length of the bed elevation profile which is quantified as roughness. In other words, the ‘scale’ in our method means the horizontal length of the FFT MW and the ‘scale’ in rms height method means the horizontal length of the quantified profile.

2 Methods

2.1 Outline and flowchart of the method

To solve the scale problem of the previous two-parameter method, we propose a self-adaptive two-parameter roughness quantization technique in four steps, as described in Figure 1.

Fig. 1. Flowchart of the self-adaptive two-parameter roughness quantization method. The ‘scale’ we use in here refers to the horizontal length of the bed that is quantified as roughness. MW is moving window.

2.2 Step 1: selecting the range of scale

Our two-parameter roughness method, and the extraction of the characteristic factors, requires a bed elevation profile with a fixed sampling interval. It is therefore necessary to interpolate and resample bed data to ensure the method can work. The roughness, and the characteristic factors, of each position on the bed is calculated within the variable-size MW at multiple scales. The total number of along-track sample points in the MW is M = 2^N, where N is the window-length exponent. Low-frequency fluctuations are observed on a large scale, and high-frequency fluctuations are observed in the small scale. We therefore calculate the roughness and the characteristic factors at multiple scales, which is analogous to the multi-scale decomposition of the wavelet transform. The scale range we need is equivalent to 10² − 10⁴ m as it includes a variety of glacial bedforms that appear within this interval (e.g. drumlins, U-shaped valleys, crag and tails, etc.) (Ó Cofaigh and others, 2002; Stokes and Clark, 2003).

2.3 Step 2: calculating the two-parameter roughness of different scales

The first parameter ξ_N(x _i) of the two-parameter roughness index, whose unit is m² (Li and others, 2010), describes the magnitude of vertical deviations in the bed slope, and is calculated by:

(1)$$Z_0\lpar x_i\rpar = Z\lpar x_i\rpar - \bar Z$$

(2)$$S_N\lpar k\rpar = {1\over M} {\left\vert {\rm FFT} \lpar Z_0\lpar x_i\rpar \rpar \right\vert ^{2}\over M/2}$$

(3)$$\xi_N\lpar x_i\rpar = \int_{k_1}^{k_2} S_N\lpar k\rpar \, {\rm d}k$$

where x corresponds to the location in the bedrock; x _i is the position of the point in the variable-size MW (1 ≤ i ≤ M); Z(x _i) denotes the bed elevation; $\bar Z$ is the mean of Z(x _i) in the variable-size MW; Z ₀(x _i) is the mean-detrended profile; S _N(k) is the spectral power density of Z ₀(x _i); k is the angular frequency; and k ₁ is zero and k ₂ is infinity.

The second parameter η _N(x _i), whose unit is m² (Li and others, 2010), quantifies the slope frequency of these deviations, and is obtained using the ratio of the first parameter defined in equation (1) to the slope index. The slope index is calculated using the FFT to convert the slope of the bed elevation profile into the average power. We obtain the slope index $\eta _{N_{{\rm sl}}}\lpar x_i\rpar$, by:

(4)$$S_{N_{{\rm sl}}}\lpar k\rpar = {1\over M} {\left\vert {FFT} {\big(\partial \left(Z_0\lpar x_i\rpar \right)\big)} \over {\partial x_i} \right\vert ^{2}\over M/2}$$

(5)$$\eta_{N_{{\rm sl}}}\lpar x_i\rpar = \int_{k1}^{k2} S_{N_{{\rm sl}}}\lpar k\rpar \, {\rm d}k$$

where $S_{N_{{\rm sl}}}\lpar k\rpar$ is the spectral power density of the slope profile ${\partial \lpar Z_0\lpar x_i\rpar \rpar \over \partial x_i}$. Finally, we obtain η _N(x _i) as follows,

(6)$$\eta_N\lpar x_i\rpar = {1\over M} {\xi_N\lpar x_i\rpar \over \eta_{N_{{\rm sl}}}\lpar x_i\rpar }.$$

2.4 Step 3: calculating the weights

Three characteristic factors are used to estimate the shape of the bed elevation profile: (1) the slopes, (2) the skewness and (3) the coefficient of variation. They are calculated simultaneously with the variable-size MW.

The slopes are usually used to indicate how inclined a curve is, which include the leading and trailing edge steepness. The leading-edge steepness ${\rm SL}_N^{\rm l}\lpar x_i\rpar$ of the bed elevation Z(x _i) within the variable-size MW is estimated by:

(7)$${\rm SL}_N^{\rm l}\lpar x_i\rpar = {Z^{\rm l}\lpar x_i\rpar \over x_i}$$

where Z ^l(x _i) is the leading edge of the bed elevation. Similarly, the trailing edge steepness, denoted as ${\rm SL}_N^{\rm t}\lpar x_i\rpar$, is calculated by:

(8)$${\rm SL}_N^{\rm t}\lpar x_i\rpar = {Z^{\rm t}\lpar x_i\rpar \over x_i}$$

where Z ^t(x _i) is the trailing edge. The visualized examples of ${\rm SL}_N^{\rm l}\lpar x_i\rpar$ and ${\rm SL}_N^{\rm t}\lpar x_i\rpar$ are shown in Figure 2. The ${\rm SL}_N^{\rm l}\lpar x_i\rpar$ and ${\rm SL}_N^{\rm t}\lpar x_i\rpar$ quantify the leading and trailing edge steepness, respectively.

Fig. 2. The leading and trailing edge steepness ${\rm SL}_N^{\rm l}\lpar x_i\rpar$ and ${\rm SL}_N^{\rm t}\lpar x_i\rpar$, the skewness SK_N(x _i) and the coefficient of variation CV_N(x _i) for N = 10 of four parts of bed elevation profile of TT’ in Figure 4. The red line represents ${\rm SL}_N^{\rm l}\lpar x_i\rpar$ and the green line represents ${\rm SL}_N^{\rm t}\lpar x_i\rpar$. (a) From 136 to 156 km of TT’. (b) From 278 to 298 km of TT’. (c) From 26 to 46 km of TT’. (d) From 145 to 165 km of TT’.

In general, skewness, which quantifies the degree of asymmetry of a curve (i.e. the larger the value, the greater the asymmetry) (Moore and Kirkland, 2007), is a measure of the direction and degree of deviation of the bed elevation profile from its center. Firstly, we calculate the expectation value, E _N(x _i), and the variance, D _N(x _i), of the bed elevation profile as follows:

(9)$$E_N\lpar x_i\rpar = \sum_{i = 1}^{M} \left[{Z\lpar x_i\rpar \over \sum_{i = 1}^{M} Z\lpar x_i\rpar } x_i \right]$$

(10)$$D_N\lpar x_i\rpar = \sum_{i = 1}^{M} \left[x_i^{2} {Z\lpar x_i\rpar \over \sum_{i = 1}^{M} Z\lpar x_i\rpar } \right]- E_N\lpar x_0\rpar ^{2}$$

where x ₀ denotes the center of the variable-size MW. Next, we calculate the skewness by:

(11)$${\rm SK}_N\lpar x_i\rpar = \left\vert \sum_{i = 1}^{M} \left[\left({x_i - E_N\lpar x_0\rpar \over \sqrt{D_N\lpar x_0\rpar }} \right)^{3} {Z\lpar x_i\rpar \over \sum_{i = 1}^{M} Z\lpar x_i\rpar } \right]\right\vert.$$

The coefficient of variation is a measure of the relative dispersion (Yadolah, 2008), and is calculated by:

(12)$${\rm CV}_N\lpar x_i\rpar = {\sqrt{D_N\lpar x_0\rpar }\over \left\vert E_N\lpar x_0\rpar \right\vert }.$$

The visualized examples of SK_N(x _i) and CV_N(x _i) are shown in Figure 2. As we can see in Figure 2, we chose a mountain, two slopes and a valley as examples. The mountain is symmetrical, leading to low values of both SK_N(x _i) and CV_N(x _i). The slopes are asymmetrical and so SK_N(x _i) and CV_N(x _i) will both be high or high SK_N(x _i) and low CV_N(x _i). The valley can lead to low SK_N(x _i) and high CV_N(x _i).

By evaluating the value of ${\rm SL}_N^{\rm l}\lpar x_i\rpar$, ${\rm SL}_N^{\rm t}\lpar x_i\rpar$, SK_N(x _i) and CV_N(x _i) from Eqn (7–12), we can quantify the shape of terrain (Khan and others, 2011; Ilisei and others, 2018), from which we obtain the weight, W _N(x _i), of each scale. More detail about how to obtain the weights and how the weights change with the three characteristic factors can be found in Appendix A and Appendix B, respectively.

2.5 Step 4: calculating the self-adaptive two-parameter roughness index

Finally, the self-adaptive first and second parameter roughness index, ξ_MS(x _i) and η _MS(x _i), are calculated by weighting and summing the roughness of the different scales in Eqn (1) and (3):

(13)$$\xi_{{\rm MS}}\lpar x_i\rpar = \sum_{N = N_{{\rm min}}}^{N_{{\rm max}}} \left[W_N\lpar x_i\rpar \xi_N\lpar x_i\rpar \right]$$

(14)$$\eta_{{\rm MS}}\lpar x_i\rpar = \sum_{N = N_{{\rm min}}}^{N_{{\rm max}}} \left[W_N\lpar x_i\rpar \eta_N\lpar x_i\rpar \right]$$

where N _min and N _max are the range of N. ξ_MS(x _i) and η _MS(x _i) are then the self-adaptive two-parameter roughness index.

3. Experiments

3.1 Simulated bed elevations

Since the sampling interval of bed elevations that we use later is ~19 m (Cui and others, 2018), and the variable-size MW should be in the range of 10² − 10⁴ m as mentioned above, we use 5 ≤ N ≤ 10 and M = 2⁵ − 2¹⁰(608 − 19456 m) to process the real bed elevation dataset. Hence, the simulated bed elevation that we used here has 1024 sample points.

To prove the effectiveness of the algorithm in quantifying roughness, we compare the quantization results of the simulated beds’ roughness with a variety of irregularities, described in Figure 3. The simulated bed profiles are generated using a self-correlation function (Thomas, 1999; Li and others, 2010). Two controllable inputs, which are standard deviation and autocorrelation length, are used to control the vertical and slope irregularities of the bed, respectively. As we increase the standard deviation but keep the autocorrelation length unchanged, the bed elevation profile transforms from that in Figures 3a to b. Accordingly, our self-adaptive first parameter roughness index ξ_MS(x _i) increases and our self-adaptive second parameter roughness index η _MS(x _i) is unchanged. On the other hand, if the standard deviation is unchanged but autocorrelation length is increased, the bed elevation profile changes from that in Figures 3b to c. In this case, ξ_MS(x _i) is unchanged and η _MS(x _i) reduces. Similarly, when the standard deviation is unchanged and autocorrelation length is reduced, which means the bed changes from that in Figures 3b to d, ξ_MS(x _i) is unchanged and η _MS(x _i) increases. Theoretically, higher values of the first roughness parameter relate to higher vertical roughnesses; whereas lower values of the second parameter roughness correspond with higher slope roughnesses (Li and others, 2010). This simple analysis shows that the method is working as it should, and can now be applied to more complex, real-world examples. The comparison between the results of our self-adaptive two-parameter method, the previous two-parameter method and the Hurst exponent method in a simulated bed elevation profile can be found in Appendix C.

Fig. 3. Random simulated bed elevation profiles with four different roughness. (a) Our self-adaptive first and second parameter roughness index, whose unit is m², ξ_MS(x _i) = 0.34, η _MS(x _i) = 1906.1. (b) ξ_MS(x _i) = 1.35, η _MS(x _i) = 1906.1. (c) ξ_MS(x _i) = 1.35, η _MS(x _i) = 455.0. (d) ξ_MS(x _i) = 1.35, η _MS(x _i) = 7347.3.

3.2 Real bed elevations

The radar profile used in this study was acquired in CHINARE 33 by the Snow Eagle 601, which is a modified DC-3 aircraft fully equipped with a suite of geophysical instruments. The apparatus of relevance to this paper is the airborne ice sounding radar, functionally similar with the High Capability Airborne Radar System (HiCARS), which uses a double flat-plate dipole antenna with a peak power of 8 kW and a center frequency of 60 MHz. The radar can penetrate Antarctic ice to a depth of more than 3 km, and can measure internal reflection horizons with a high detection resolution (Cui and others, 2018).

The radar data were processed using along-track coherent stacking, pulse compression and along-track incoherent stacking. The data resolution in the along-track direction is ~19 m after the process of interpolation and resampling, and the depth resolution is ~10 m in air and ~5.6 m in ice (Cui and others, 2018). The bed-ice interfaces were traced and digitized by a semi-automatic program, using the manually picked lower and upper boundaries of the bed interface as a constraint on automatic bed interface picking (Cui and others, 2020). The program was also used to process the former ICECAP data from the Aurora and Wilkes subglacial basins (Blankenship and others, 2016, 2017).

In this paper, we have analyzed the radar profile ‘LSE_GCX0e_Y75a’, which is referred to as TT’ and has a length of 323.9 km. The radar profile is shown in Figure 4b, and its bed elevation and ice-surface are shown in Figure 4c. The roughness results from our algorithm are shown in Figures 4d and e. The radar profile is located between ~ 77°E and ~ 87°E in an approximate west-east direction, covering a region with low ice-flow speed, ranging from 0 to 20 m~a⁻¹ (Fig. 4a). The bed elevations are high between 0 and 60 km, and elsewhere vary closer to the sea level with a relatively low amplitude of several hundred meters. An obvious mountain occurs between ~35 and 55 km, and two relatively smooth areas appear between ~190 and ~230 km and between ~280 km and the end of the profile.

Fig. 4. The profiles of LSE_GCX0e_Y75a (TT’). (a) Ice velocity map and location of TT’. Ice velocity map is from the MEaSUREs InSAR-Based Antarctica Ice Velocity Map, Version 2 (Rignot and others, 2017). Grounding lines are shown in green. (b) Ice-sounding radar profile. (c) Ice-surface elevation and bed elevation. Our self-adaptive first and second parameter roughness index (d) ξ_MS(x _i) and (e) η _MS(x _i).

3.3 Comparison of results

3.3.1 The Hurst exponent

The statistical methods used to calculate the Hurst exponent (H), and thus to quantify the self-affine scaling behavior, are well established in the literature (Malinverno, 1990; Shepard and others, 2001; Kulatilake and others, 1998; Orosei and others, 2003). A feature of terrain with a higher H is that it tends to appear relatively rough at larger length scales (low frequencies) and smooth at smaller length scales (high frequencies). A feature of lower H terrain is that it tends to have similar roughness across different length scales (Jordan and others, 2017). We used the rms height, R(L), to extract the Hurst exponent:

(15)$$R\lpar L\rpar = \sqrt{{1\over P-1} \sum_{i = 1}^{L} \left[Z\lpar x_i\rpar -\bar{Z} \right]^{2}}$$

where P is the number of sample points within the profile window; and L is the length of the profile window. We then extract the Hurst exponent H by the rms height as follows:

(16)$$R\lpar L\rpar = R\lpar L_0\rpar \left({L\over L_0} \right)^{H}$$

where L ₀ is a reference profile length. As L increases, R(L) also increases, as they have a linear relationship. Thus, the Hurst exponent H is the slope of L and R(L). The Hurst exponent is a multi-scale standardized result that can be compared with the rms height, or slopes of different surfaces, as measured by different tools at similar scales (Shepard and others, 2001).

3.3.2 Roughness comparisons

Profile TT’ crosses the Kronshtadtskiy Glacier. The first 50 km is characterized by high ξ_MS(x _i) and high η _MS(x _i), while the other sections are mainly characterized by low ξ_MS(x _i) and low η _MS(x _i) roughness (Figs 4d and e). We now focus on the mountain at about 40 km (Fig. 5). Comparing the roughness results of the profile at this place, we can see that ξ_MS(x _i) and ξ_N(x _i) each have peaks. But ξ_N(x _i) shows very low values from 39 to 42 km, and ξ_MS(x _i) shows high values representing the span of the mountain. η _MS(x _i) has several very low peaks, indicative of high roughness, in contrast to η _N(x _i). The trend of Hurst exponent, η _MS(x _i) and η _N(x _i) is roughly the same, as described in Figure 5c.

Fig. 5. A part of roughness of profile TT’ (from 32 to 52 km). The marked points locate in the center of the profile. (a) The bed elevation profile. (b) Our self-adaptive first parameter roughness index ξ_MS(x _i) and Li's first parameter roughness index ξ_N(x _i). (c) Our self-adaptive second parameter roughness index η _MS(x _i), Li's second parameter roughness index η _N(x _i) and the Hurst exponent H. The ξ_N(x _i) and η_N(x _i) are calculated when N = 5.

Li et al. (2010)'s first and second parameter roughness index ξ_N(x _i) and η _N(x _i) are sensitive to the MW size. To prove this, we calculate ξ_N(x _i) = 55.1 and η _N(x _i) = 32.9 at N = 5 in the central position of the profile, as shown in Figure 5. We then calculate ξ_N(x _i) and η _N(x _i) in same position while increasing N, revealing differences in both ξ_N(x _i) and η _N(x _i) (i.e. ξ_N(x _i) = 30.1 and η _N(x _i) = 22.9 at N = 6; ξ_N(x _i) = 66.0 and η _N(x _i) = 86.3 at N = 7; ξ_N(x _i) = 431.9 and η _N(x _i) = 123.4 at N = 8; ξ_N(x _i) = 800.4 and η _N(x _i) = 195.2 at N = 9; and ξ_N(x _i) = 611.4 and η _N(x _i) = 174.0 at N = 10). In contrast, our self-adaptive first and second parameter roughness index ξ_MS(x _i) = 332.5 and η _MS(x _i) = 105.8 are stable across the scale range, as shown in Figures 5b and c. The unit of ξ_N(x _i), η _N(x _i), ξ_MS(x _i) and η _MS(x _i) is m² (Li and others, 2010).

Through the last 76 km of TT’, we can see several mountain peaks and valleys (Fig. 6). By comparing the η _MS(x _i) and H of the profile in this place, we can see that η _MS(x _i) has several very low peaks, indicative of high roughness, but the Hurst exponent has too many peaks to resolve (Fig. 6c).

Fig. 6. A part of roughness of profile TT’ (from 209 to 285 km). (a) The bed elevation profile. (b) The roughness index ξ_MS(x _i). (c) The roughness index η _MS(x _i) and the Hurst exponent H.

4. Conclusions

We describe a new method to determine subglacial topographic roughness, called the ‘self-adaptive two-parameter roughness quantization’ method. The algorithm is based on the two-parameter roughness quantization method and is analogous to a wavelet transform. According to the weights calculated from three characteristic factors, the new algorithm can determine the scale needed to quantify the roughness at any position. This eliminates the problem of a fixed-size MW of measurements being unable to cope with roughness at multiple scales.

We use simulated bed elevations to test the new algorithm, demonstrating hypothetically the effectiveness of the new algorithm in quantifying roughness. In cases where real bed elevations are used, we show that the algorithm can achieve a better characterization of roughness than the previous two-parameter method. Comparisons between the results of our self-adaptive two-parameter method, the previous two-parameter method and the Hurst exponent method show the advantage of our approach in capturing different scales of terrain span and demonstrate the effectiveness of the new algorithm in quantifying bed roughness.

This method provides an improved tool for quantifying the roughness of the subglacial environment from a multi-scale perspective, especially in the analysis of complex subglacial topography. The new method offers a means by which multi-scale roughness can be quantified at a continental scale. Such a product would have significant glaciological applications, especially in distinguishing between erosional and depositional landscapes. The method could be applied to higher along-track resolution data if it is available in the future for evaluating the magnitude and distribution of basal sliding at the meter scale. The method works particularly well over subglacial landscapes with horizontal lengths greater than the length of a traditional fixed-size MW (e.g. mountains and valleys), resolving the roughness in all places. It also works well in landscapes that have notable roughness at a scale less than a fixed-size MW, as it is able to incorporate the fine-scale roughness within the wider terrain.

The code/routines of the self-adaptive two-parameter roughness quantization method will be freely available under requests from any researchers. We hope the new method can be extensively used by the community to characterize subglacial roughness in future.