1. Introduction

Statistically characterizing the topographic structure of the basal boundary beneath ice masses is increasingly being recognized as an important imperative in glaciological research (e.g. Reference Hubbard and HubbardHubbard and Hubbard, 1998; Reference Taylor, Siegert, Payne and HubbardTaylor and others, 2004;Reference Rippin, Bamber, Siegert, Vaughan and CorrRippin and others, 2006). Subglacial topography at all scales imparts important clues concerning the distribution of subglacial sediments, subglacial geomorphic processes and ice-stream stability (Reference Hubbard, Siegert and McCarrollHubbard and others, 2000;Reference Bingham and SiegertBingham and Siegert, 2009), and is thought to exert a fundamental control on ice dynamics (Reference WeertmanWeertman, 1957; Reference Siegert, Taylor, Payne and HubbardSiegert and others, 2004, 2005). Nevertheless, to date, few studies have examined subglacial topographic patterns and their influences quantitatively, perhaps due to the paucity, until recently, of appropriate bed profile datasets from contemporary glacierized environments. Modern techniques have alleviated this issue: radio-echo sounding (RES) has enabled the acquisition of myriads of bed profiles at the regional scale (decimetres to kilometres) across ice sheets and ice caps (e.g. Reference Lythe and VaughanLythe and others 2001; Reference GogineniGogineni and others, 2007; Reference SunSun and others, 2009), Reference Hubbard, Siegert and McCarrollHubbard and others (2000) have demonstrated the ease with which centimetre- and metre-scale bed profiles may be obtained across recently deglaciated terrain using microroughness meters and electro-optical distance meters respectively, and terrestrial laser scanning systems have shown great potential in surface characterization, especially with respect to threedimensional (3-D) surface characterization (e.g. Reference Bauer, Paarand, Kaufmann, Phillips, Springman and ArensonBauer and others, 2003). As ‘real’ bed data have thus become easier to obtain, and there is an increasing drive to integrate such data into numerical models of ice dynamics and ice-stream behaviour (Reference SolomonSolomon and others, 2007;Reference PattynPattyn and others, 2008;Reference HindmarshHindmarsh, 2009), the need for characterizing subglacial topography in a statistically meaningful manner has arisen. In this context, recent studies, and the work described in this paper, have made use of the concept of subglacial ‘roughness’.

Generally, roughness is regarded as the irregularity of a surface. In glaciology, the initial attempt at expressing subglacial bed roughness was due to the requirement of estimating how a glacier’s basal sliding is controlled by the undulation of bedrock. Reference WeertmanWeertman (1957) considered an idealized case of an undulating bed as cubic obstacles with length, a, distributed periodically with a separation distance, λ, between each obstacle. Put another way, the two parameters, a and λ, in this model denote the amplitude and spatial frequency of the undulations of the idealized bed. Expanding Weertman’s theory to consider bedrock surfaces with randomly distributed undulations, Reference KambKamb (1970) and Reference NyeNye (1970) introduced Fourier transformations (FTs) of bed elevation profiles into their theories. A FT can be used to transform any surface into a sum of several periodically corrugated surfaces, and can thus be used to express the amplitude and spatial frequency of the range of undulations present. Hence, these and follow-up studies have quantified roughness in terms of two key parameters: the vertical irregularity (amplitude) and the horizontal irregularity (spatial frequency); or in an equivalent way, as slope is the bridge between vertical and horizontal, using amplitude and slope (Reference PatersonPaterson, 1994).

Though characterizing roughness by FT works well in theoretical dynamic studies as described above, it is not convenient for showing the spatial distribution of roughness. A basic tenet of cartography is that showing the spatial distribution of any property ideally requires that it is displayed in terms of spatial variations in the magnitude of a single parameter. Unfortunately, a full Fourier spectrum for roughness is too complicated to meet this requirement. Reference Hubbard, Siegert and McCarrollHubbard and others (2000) and Reference Taylor, Siegert, Payne and HubbardTaylor and others (2004) circumvented this problem by defining a single-parameter roughness index defined as the integral of the spectrum within a specified wavelength interval; this method has since been applied to show the spatial distribution of subglacial roughness across Antarctica (Reference Siegert, Taylor, Payne and HubbardSiegert and others, 2004, 2005; Reference Bingham and SiegertBingham and Siegert, 2007, Reference Bingham and Siegert2009). In this paper, we build on this method by additionally considering the FT of bed- slope profiles as well as that (traditionally considered) of bedrock amplitude, and link it to horizontal irregularities. We believe that calculating the FT of both bed slope and amplitude enhances the robustness of spectral roughness techniques in subglacial landscape characterization and process inference.

2. Method

We begin by calculating the single-parameter roughness index of a profile of length, (-l/2,l/2). We first remove the mean, (Z(x)) from the bed elevation profile, Z(x):

(1)

where x corresponds to the horizontal location and Z corresponds to the raw elevation data. We then apply a FT to obtain the spectral power density, S(k):

(2)

where is the FT of Z ₀(x). The roughness index, ξ, is obtained by integrating S(k) in a specified wavelength interval, (k ₁,k ₂):

(3)

When k ₁ is zero and k ₂ is infinity, we obtain a ‘total roughness index’, ξ _t. This is the property used to characterize subglacial roughness across Antarctica in several previous studies (e.g. Reference Siegert, Taylor, Payne and HubbardSiegert and others, 2004;Reference Bingham and SiegertBingham and Siegert, 2009). Figure 1 illustrates this method for four rough surfaces. The surface in Figure 1a is clearly less rough than that in Figure 1b, and this is reflected in their •_t values (0.3 and 1.0 respectively). However, while the surfaces shown in Figure 1b-d all yield the same ξ_t values, they are qualitatively quite different from one another: the surface in Figure 1b is steeper than that in Figure 1d, and that in Figure 1c is the steepest. We suggest, therefore, that while the parameter • has demonstrably been useful for describing roughness, in isolation it may not be sufficient. Here we investigate the utility of considering the extra factor, the FT of bed-slope profile.

Fig. 1. Random surfaces of different roughness, generated using a self-correlation function (Reference ThomasThomas, 1999). (a–d) show the raw elevation profiles, and (e–h) their respective spectral power densities. All {ξ, η} are in units of {ξ_b, η _b}.

Returning to our conception of roughness in terms of vertical irregularities (amplitude) and horizontal irregularities (spatial frequency), we can interpret Figure 1b-d as follows: the amplitudes of the undulations are similar, but the frequencies of the undulations differ. Hence, a factor related to horizontal irregularity needs to be introduced. Given that a bed-slope profile expresses the connection between the vertical and horizontal dimensions, it is reasonable to choose some statistical function of slope as the extra factor. A slope profile, sl(x), can be calculated from the elevation profile, Z ₀(x), by sl(x) = ∂Z ₀(x)/∂x, and its spectral power density, S _sl(k), is S _sl(k), is

(4)

Following Equation (4), we obtain the slope index, ξ_sl, by

(5)

Finally, as slope is the ratio of vertical dimension to horizontal dimension, we define the parameter p, corresponding to horizontal irregularity alone and invariant of vertical stretching of the profile, as

(6)

Our two-parameter roughness index, taking into account both the vertical and the horizontal irregularities, is thus defined as {ξ, η}. Thus, while • describes the magnitude of vertical deviations in the bed slope, η quantifies the horizontal frequency of these deviations. This distinction is expressed graphically in Figure 1 (and is illuminated mathematically in section 3). Here each surface in Figure 1b-d, which give the same value of ξ_t, yields different values of p _t. Figure 1a and b, which have the same η_t, have different ξ_t. Figure 1 therefore demonstrates the utility of {ξ_t η_t} for distinguishing between all four surfaces.

As a further test, we apply our new roughness index, {ξ, η}, to a bed elevation profile acquired from East Antarctica by surface RES obtained as part of the Chinese Antarctic National Research Expedition (CHINARE) program, in 2004 (Reference SunSun and others, 2009). The profile is ˜200km long with a spatial resolution of 70m (Fig. 2). For these discrete data, summing is used in place of integration and a discrete Fourier transform (DFT) is used;a fast Fourier transform (FFT) can further be used in place of the DFT to increase efficiency. In order to consider roughness as a localized property, and to maintain consistency with earlier Antarctic roughness studies (Reference Bingham and SiegertBingham and Siegert, 2009), we calculate {ξ, η} in a moving window of length 1024 points (about 70 km). The results of components ξ and η are shown in Figure 2b and c respectively. It is apparent from Figure 2 that the two parameters describe the surface more comprehensively than each in isolation. For example, the section centred at 35 km has small ξ and small η, reflecting low- amplitude, high-frequency roughness; whereas the section centred at 120km has larger ξ and η, a result of higher amplitudes but a smoother horizontal surface.

Fig. 2. (a) ˜220 km subglacial elevation profile from East Antarctica, acquired by 21st CHINARE (Reference SunSun and others, 2009). (b, c) Corresponding distribution of ξ and η respectively.

3. Statistical Geometry

To better understand the statistical meanings of ξ and η in a simple way, we consider the total roughness case, ξ_t, where (k₁,k₂) = (0, ∞). According to Parseval’s theorem (Reference RudinRudin, 1987), the value of ξ can be deduced:

(7)

Hence the ξ component of ξ _t is equal to half the mean square of the elevation profile. We calculate the distribution of the mean square of the elevation profile (Fig. 2a), and compare it with the spatial distribution of ξ shown in Figure 2b. As anticipated, both distributions are the same. Similarly, ξ_sl equals half the mean square of the slope profile. This explains why ξ reflects only the vertical irregularity, even though ξ originates from the spectral power density which contains both vertical and horizontal information. While ξ_sl denotes how quickly Z ₀(x) changes, component p, the ratio of • to ξ_sl, actually measures how frequently Z(x) varies about its mean. This can be illustrated mathematically by applying Parseval’s theorem to the definition of η. Consider the total roughness case,

(8)

So η is actually a weighted harmonic average value of (2π/ k)² (square of the wavelength λ), and the weight is S(k). Thus a bed elevation profile with roughness {ξ, η} means it has the same roughness index as a sinusoidal surface Though the horizontal irregularity can partly be shown by ξ in the case of the interval of integration being divided manually into short/medium/long wavelength parts, η is more convenient and suitable for expressing the horizontal irregularity continuously. Where two surfaces have the same total roughness component, ξ_t, the surface with low ξ at high frequencies tends to have larger η, while the surface with low ξ at low frequencies tends to have smaller η. This is shown in Figure 1c, d, g and h. Figure 1h shows that the spectral power density lies mainly in the low-frequency range, so the η of Figure 1d is larger than that of Figure 1c.

At this stage, we note some points concerning the calculation of {ξ, η}. Firstly, although in the case of total roughness, {ξ_t η _t }, it would be more rapid and straightforward simply to derive the mean square instead of following the method outlined in section 2, we prefer our method in all other cases as it helps to discern specific wavelengths of interest. Secondly, in cases where uncertainties in the original slope profile are high, an alternative may be to use a moving mean across a small window rather than using the raw elevation profile. Thirdly, the demonstration that ξ_t is equal to the mean square of the elevation profile is valid when the elevation profile is a continuous function, but in practice this calculation will be based on a discrete set of points, and in this case the relation holds only when these points are equally spaced. Fourthly, it may well be worth noting here that this and previous analyses are all in two dimensions, largely reflecting available data as transects, but that ultimately techniques should be developed for 3-D analysis as empirical data collection improves in this direction.

4. Implications for Glaciology

4.1. Ice dynamics

On the scale of metres, bed roughness exerts a major control on an ice mass’s basal motion. Reference WeertmanWeertman (1957) noted that two mechanisms operate in basal sliding: regelation and deformation. Regelation describes the process by which basal ice melts on the upstream (stoss) side of a basal obstacle and then refreezes on its downstream (lee) aspect. The difference in melting point between the up- and downstream sides results from differences in pressure exerted on each side of the obstacle by the overflowing ice, while the size of the obstacle controls the heat flux between the melting and refreezing ice on stoss and lee. Deformation is the increase in the rate of plastic flow around large bumps due to the additional stresses exerted on the adjacent ice.

The regelation-induced basal sliding velocity, u ₁, deformation-induced basal sliding velocity, u ₂, and total basal sliding velocity, u, are (Reference WeertmanWeertman, 1957)

(9)

(10)

(11)

where a and λ are the average height and length of basal obstacles, C₁ , C ₂ and C are factors related to non-geometric variables such as thermal conductivity and shear stress, and n is Glen’s flow-law exponent (Reference GlenGlen, 1958) whose value is discussed to be ˜3 (e.g. Reference Weertman, Whalley, Jones and GoldWeertman, 1973; Reference Russell-Head and BuddRussell-Head and Budd, 1979).

Based on the statistical meanings of ξ and η discussed in section 3, we can use our basal roughness index as a proxy for the vertical and horizontal length scales that are being used to estimate sliding velocity. Firstly, the interval of integration for {ξ, η} should be in the metre-scale waveband. Then, given that ξ is proportional to the mean square of amplitudes, a can be written as:

(12)

Similarly, λ can be written as

(13)

Hence u ₁, u ₂ and u can be estimated by

(14)

(15)

(16)

This method provides a way of summarily and rapidly evaluating the magnitude and distribution of the aspect of basal sliding directly attributable to basal roughness using simple statistical parameters of a subglacial bed rather than detailed topographic profiles. Other factors such as temperature (hence rate factor), possible presence of basal melt to lubricate the base, lateral stress transmitted from other parts of the ice sheet, and gravitational driving stress relating to surface and basal slopes are also important and contained in c. Because the scattering of electromagnetic waves relates to bed roughness (Reference Tsang, Kong and DingTsang and others, 2000), obtaining simple statistical parameters of a subglacial bed requires fewer RES surveys than for detailed topographic profiles, so for subglacial beds with no detailed topographic profiles and requiring RES surveys, using our method means much fewer RES surveys are needed to evaluate basal sliding. For the areas where detailed subglacial bed profiles have already been obtained, our method requires much less computational time, and input data are needed to generate useful information relevant to basal conditions and motion processes, considered below.

4.2. Geomorphic implications

At the regional scale, Reference Bingham and SiegertBingham and Siegert (2009) proposed that subglacial roughness can be interpreted at least partially in terms of the geomorphic processes operating in subAntarctic and Quaternary glaciated landscapes. Rougher subglacial topography, characterized by larger ξ, was tentatively associated with conditions such as continental settings, mountainous preglacial topography, cold and slow ice flow, and low rates of erosion and deposition. Smoother subglacial topography, described by lower ξ, was associated with marine settings, relatively flat preglacial topography, warm and fast ice flow, and high rates of erosion and deposition. However, many of these factors are interlinked and cannot be distinguished solely by ξ. For example, a continental setting intensively eroded by fast ice flow may have a similar value of ξ to a marine setting with high deposition.

Now, using our two-parameter index, {ξ, η}, it may be possible to characterize and discriminate between subglacial process domains (e.g. erosional as opposed to depositional and bedrock as opposed to marine sediments) more precisely. Figure 3 shows four end-member permutations of the two indices, and their corresponding geo- morphic implications are discussed below.

Fig. 3. Schematic framework for geomorphic interpretation of subglacial landscapes based on {ξ, η}. The subglacial profiles are drawn from Reference Siegert, Taylor and PayneSiegert and others (2005).

Low ξ high η

Over a sufficient length of time, the net effect of subglacial erosion may be assumed to lower peaks and fill valleys, thereby lowering the amplitudes of bed undulations and reducing roughness but leaving the spacing between obstacles largely unchanged. Thus in erosional regions we would expect gradual lowering of • but little change in p. By contrast, large depositional events, most commonly resulting from marine incursions, could, in theory, entirely smother the preglacial landscape, thereby affecting both ξ and p. Making a direct comparison of {ξ, η} between an original landscape and its present ‘deposited’ surface is difficult, because few subglacial landscapes have been imaged below the ice/bed boundary. However, we can understand the difference between an original landscape and its present ‘deposited’ surface by comparing bedrock profiles with overlying englacial internal reflection horizons in RES data. Assuming an area of slow flow, thereby discounting disruption to internal layering caused by ice streaming, we can view the palaeo-landscape (bedrock) and its depositing surface (internal layers just above the bedrock) simultaneously. It is both observed and proven (Reference Hindmarsh, Leysinger Vieli, Raymond and GudmundssonHindmarsh and others, 2006) that the depositing surface tends to parallel the long-wavelength undulations of the underlying bedrock while being insensitive to short-wavelength undulations. In such a situation, deposition demonstrably enlarges η. Therefore, we can deduce that for small-ξ landscapes, those with largeηg are likely dominated by deposition reflecting marine landscapes. This deduction is consistent with previous results (Reference Siegert, Taylor and PayneSiegert and others, 2005): where the bed is below sea level, supporting intensive marine preglacial deposition and postglacial reworking, η is small, especially at high frequencies (implying g becomes large). This subdued topography is often found beneath marine portions of both the East and West Antarctic ice sheets as a consequence of preglacial marine deposition of sedimentary material (Reference Siegert, Taylor and PayneSiegert and others, 2005).

Low ξ, low η

As an enlarged η is linked to deposition, a subglacial topographic surface with small η implies no significant deposition has occurred, so a low ξ associated with this should be caused by erosion. Such a landscape therefore likely represents a continental setting influenced by intensive erosion. This deduction is consistent with previous results: where subglacial sliding is predicted by models (e.g. Wilkes and Aurora subglacial basins, Antarctica) (Reference Siegert, Taylor and PayneSiegert and others, 2005), ξ is found to be small at all frequencies (implying g is evidently unaffected).

High ξ, low η

Bed roughness with high ξ and low η (i.e. high amplitude, low wavelengths) is characteristic of beds which are alpine in nature yet cold-based so that no subglacial sliding occurs. In Antarctica, such conditions have been observed beneath the ice sheet at Dome A and across the subglacial Gamburtsev Mountain Province, the subglacial landscape recording mountain glaciation during rapid ice build-up ˜34Ma ago, and the subsequent preservation of this landscape as a consequence of cold-based, low-erosion conditions predominating beneath the stable ice dome that has since overlain this region (Reference SunSun and others 2009). Subglacial uplands in Dronning Maud Land and parts of West Antarctica have similar features, their alpine-style topography having remained relatively unmodified as the ice sheet grew and the bed over upland areas became cold- based (Reference Jamieson, Sugden and HultonJamieson and others, 2010).

High ξ, high η

A bed retaining a high ξ value suggests it has suffered neither heavy erosion nor deposition, but an enlarged η implies there may have been some periods during which the subglacial landscape experienced deposition. Consequently, high ξ with high η could be linked to continental settings which have undergone deposition over a short period and have subsequently become protected beneath a cold and slow-flowing basal thermal regime. Such a case is presented by the Dome C region, where the bed was previously described as smooth at high frequencies and rough at the medium and long wavelengths (Reference Siegert, Taylor and PayneSiegert and others, 2005).This previous interpretation for the bed of Dome C is similar to our deduction.

Between all the end-member cases described above, a variety of topographies naturally exist, all of which can be characterized quantitatively by relative values of ξ and η The particular value of quantifying bed roughness in this way is that it comprises an independent quantitative means of estimating former ice-sheet dynamics, independent of numerical ice-sheet modelling reconstructions (cf. Reference Jamieson, Sugden and HultonJamieson and others, 2010), and therefore constitutes a potentially useful tool for validating model-derived reconstructions.

Conclusions

We have extended the current use of FT methods to characterize subglacial roughness by supplementing the ‘traditional’ component, ξ (which effectively characterizes vertical roughness), with a new parameter, η which effectively characterizes horizontal roughness. Both parameters have dimensions of square of length. The method can be applied freely at any scale. Our analysis has shown that while (the traditionally used) •ξ alone forms a useful component for characterizing roughness across subglacial landscapes, the addition of the horizontal component, η, significantly enhances our ability to characterize and distinguish between different types of subglacial environment.

At the metre scale, calculation of {ξ, η} represents an efficient means of characterizing the basal boundary condition for investigations of glacier dynamics, and can be used to guide interpretations of the magnitude and distribution of basal sliding. At the regional-scale, {ξ, η} offers a technique for further distinguishing the factors required for gemorphological interpretation of the subglacial bed, and especially for distinguishing between subglacial landscapes of erosion and deposition.