The periodic topography of ice stream beds: Insights from the Fourier spectra of mega‐scale glacial lineations

Ice stream bed topography contains key evidence for the ways ice streams interact with, and are potentially controlled by, their beds. Here we present the first application of two‐dimensional Fourier analysis to 22 marine and terrestrial topographies from 5 regions in Antarctica and Canada, with and without mega‐scale glacial lineations (MSGLs). We find that the topography of MSGL‐rich ice stream sedimentary beds is characterized by multiple, periodic wavelengths between 300 and 1200 m and amplitudes from decimeters to a few meters. This periodic topography is consistent with the idea that instability is a key element to the formation of MSGL bedforms. Dominant wavelengths vary among locations and, on one paleo ice stream bed, increase along the direction of ice flow by 1.7 ± 0.52% km−1. We suggest that these changes are likely to reflect pattern evolution via downstream wavelength coarsening, even under potentially steady ice stream geometry and flow conditions. The amplitude of MSGLs is smaller than that of other fluvial and glacial topographies but within the same order of magnitude. However, MSGLs are a striking component of ice stream beds because the topographic amplitude of features not aligned with ice flow is reduced by an order of magnitude relative to those oriented with the flow direction. This study represents the first attempt to automatically derive the spectral signatures of MSGLs. It highlights the plausibility of identifying these landform assemblages using automated techniques and provides a benchmark for numerical models of ice stream flow and subglacial landscape evolution.


Introduction
Mega-scale glacial lineations (MSGLs) are elongated sedimentary ridges of low relief and represent the landform signature of fast-flowing ice streams [e.g., Clark, 1994;Canals et al., 2000;Stokes and Clark, 2001;Ó Cofaigh et al., 2002;King et al., 2009]. Ice streams are the primary drainage routes for much of the ice flowing to the ocean from the Antarctic and Greenland Ice Sheets, as they were for the Pleistocene ice sheets [Bamber et al., 2000;Bennett, 2003;Rignot and Kanagaratnam, 2006;Rignot et al., 2011]. Ice stream dynamics influence ice sheet mass balance and play a fundamental role in regulating sea level changes [Anderson et al., 2002;Stokes et al., 2016]. Ice stream flow is known to be influenced by both external and internal forcing mechanisms [e.g., Bindschadler et al., 2003;Jamieson et al., 2014;Fisher et al., 2015] but is ultimately controlled by the interaction between fast-flowing ice and the bed, where MSGLs are formed [e.g., Engelhardt et al., 1990;Joughin et al., 2004;Piotrowski et al., 2004]. The genesis and evolution of these landforms can therefore reveal key clues about the processes occurring at the ice-bed interface and on the mechanism(s) of ice stream flow [e.g., Jamieson et al., 2016;Spagnolo et al., 2016]. Although the beds of modern ice streams provide important information [e.g., Smith et al., 2007;King et al., 2009], unimpeded access to paleo ice stream beds in both terrestrial and marine settings can shed light on ice stream behaviors and processes that took place over longer time scales than modern observations permit [e.g., Beget, 1986;Hicock et al., 1989;Dowdeswell et al., 2004;Ó Cofaigh et al., 2005;Stokes and Clark, 2001;Margold et al., 2015a]. Despite the importance of this topic and the multiple hypotheses proposed, there is presently no consensus on the formation of MSGLs [Lemke, 1958;Bluemle et al., 1993;Clark, 1993;Clark et al., 2003;Shaw et al., 2008;Fowler, 2010;Stokes et al., 2013; SPAGNOLO ET AL.
PERIODIC TOPOGRAPHY OF ICE STREAM BEDS 1 Spagnolo et al., 2014]. Through further investigation of these landforms, existing theories can be tested and new theories inspired, thus providing novel insight into the processes acting at the ice-bed interface. One aspect of MSGL formation that has received relatively little quantitative attention is the spatial organization of ice stream bed topography.
MSGLs are typically defined by their morphology, characterized by parallel ridges and troughs that are kilometers long, hundreds of meters wide, and only a few meters high [Clark, 1993;Spagnolo et al., 2014]. Despite their similarity to, and spatial association with, other streamlined subglacial bedforms, especially drumlins [Ely et al., 2016], a distinctive aspect of MSGLs is that they can occupy entire ice stream beds without any discernible topographic break between individual landforms. This, and their generally very low (few meters) amplitude, can make it difficult to identify and delimit individual MSGLs [Spagnolo et al., 2014]. Whereas the lack of a topographical break between adjacent MSGLs is evident in ice stream bed topography (Figure 1), the extent to which MSGLs are spatially organized (i.e., regular) remains unclear. Spatial organization can have important implications for hypotheses of their formation [Barchyn et al., 2016], with some predicting that MSGLs evolve via the amplification of initial perturbations that leads to a topography characterized by dominant wavelengths [Fowler, 2010;Fowler and Chapwanya, 2014].
A combination of manual and semiautomated geographic information system techniques has recently highlighted that MSGLs in different settings are characterized by similar ridge-to-ridge (across-flow) spacing [Spagnolo et al., 2014]. However, that analysis relied on (inevitably subjective) manual mapping and semiautomated measurements and primarily focused on average values at the scale of a whole ice stream bed. A more robust and objective analysis of the subglacial topography has not yet been achieved, limiting our ability to quantitatively test hypotheses regarding MSGL formation and, more generally, to characterize the interface between ice streams and their beds. At regional scales, studies of ice sheet basal topographies in Antarctica have indicated that ice stream beds have low roughness values that are consistent with a depositional environment and the presence of streamlined bedforms [Bingham and Siegert, 2009;Li et al., 2010;Rippin et al., 2014]. However, most previous studies were based on one-dimensional analysis of individual radio echo sounding profiles of coarse spatial resolution and only examined wavelengths larger than typical MSGLs. Characterizing the potential regularity of MSGLs, along with their shape and scale, could therefore contribute to the interpretation of ice stream basal roughness and its parameterization in ice sheet models [Schoof, 2002].
In this paper, we apply two-dimensional (2-D) discrete Fourier transforms to describe the orientation-related roughness of MSGLs and associated topography. We apply the technique to marine and terrestrial locations and to both modern and paleo ice stream beds in Antarctica and Canada, to quantitatively address the following questions: 1. Is the topographic signal of MSGLs in the frequency domain distinct from other topographies? 2. To what extent are MSGLs periodic, and do their dominant wavelengths vary from setting to setting or along the length of a single ice stream bed?
The answers to these questions bear directly on theories for the formation and evolution of MSGLs, and the quantitative, spectral analysis presented here allows differentiation among sites that has not been possible using existing methods. Furthermore, characteristic frequency domain signals can serve as input data or output constraints on numerical models of ice stream flow and landscape evolution. Quantification of MSGL spatial regularity and related amplitudes informs our understanding of the low basal roughness beneath ice streams and could prove useful in methods for the automated identification of MSGLs. These results and their implication for ice stream processes are presented and discussed below.
Because topographic amplitude is typically correlated with wavelength [Sayles and Thomas, 1978] and we seek to identify anomalously large-amplitude features, we calculate the median topographic amplitude from all orientations. Amplitudes at each orientation are then normalized by the median amplitude to find the dominant orientations and wavelengths; i.e., those characterized by relatively high-amplitude periodic variations. Our approach enables us to identify and quantitatively characterize the topographic signals with the most regular wave-like expressions.

Calculation of Fourier Transforms
Following customary procedures for spectral analysis, each sample DTM was detrended by removing a best fit plane from the topography before applying a 2-D raised cosine (Hann) window as a weighting function to the entire 100 km 2 DTM [Perron et al., 2008]. Detrending prevents spurious long-period spectral peaks that would arise through the assumption that finite topographic samples are periodic (an assumption intrinsic to the Fourier transform). Windowing the data maintains sharp spectral resolution (the ability to resolve features with different wavelengths) and minimizes spectral leakage (the smearing of high-amplitude peaks into adjacent frequencies that arise from using discretely sampled data to characterize continuous topography). Detrending with a plane, rather than a higher order polynomial, allows us to analyze all topographic wavelengths, without making subjective judgments regarding the removal of low-frequency topographic variations. The result of these steps is a map of elevation deviations from 0 m (via detrending) that smoothly tapers at the edges to 0 m (via windowing; Figure 2c). Two-dimensional fast Fourier transforms were then performed on the detrended and 2-D Hann-windowed DTMs prior to application of a scalar correction factor to the Fourier transform amplitudes to account for the reduction of amplitudes associated with windowing the data. Finally, resulting spectra were shifted to center the zero frequency, and, because the frequency space is 180°rotationally symmetric for the realvalued DTMs, we discarded their bottom halves. Figure 3 illustrates an example of the result of these steps for DTM PI_MSGL_1, showing the mean amplitude (half-height) of features with specified easting and northing spatial frequencies. A different reading of these results is to consider them in polar, rather than Cartesian coordinates. In this view, the amplitude with f x easting and f y northing frequencies is equivalent to the amplitude of a periodic feature with wavelength and periodic topography oriented at As defined, θ is the orientation along the crests of corrugations such as MSGLs, measured clockwise from north. This definition allows us to describe the strike of the MSGLs and the ice-flow direction, in contrast with the traditional approach of defining orientation perpendicular to periodic features. Furthermore, while the wavelengths identified through 1-D Fourier transforms would depend on the angle between MSGL crests and extracted topographic profiles, the wavelength and amplitude of our 2-D analysis are robust to the  Figure 3, orientations with the largest normalized amplitudes are colored and labeled by their trends in the legend (orientations relative to north), whereas smaller normalized amplitudes are dotted gray. This representation shows that the MSGLs with a À12.6°W of N trend reach amplitudes up to 42 times higher than the median amplitude at a wavelength of 690 m.
Journal of Geophysical Research: Earth Surface 10.1002/2016JF004154 orientation of the reference frame. In 1-D, the calculated wavelengths of an MSGL field will lengthen as the angle between MSGL orientations and the orientation of the 1-D profile departs from 90°. However, in 2-D, as the two orthogonal axes (initially northing and easting) are rotated over an MSGL field, the periodic frequency of one axis will increase as that of the second axis decreases, such that λ and their amplitudes are independent of the reference frame orientation. Thus, in 2-D, periodic features are identified only along their along their specific orientations ( Figure S1 in the supporting information).

Quantification of Dominant Wavelengths and Orientations
We examined how topography varies for the full range of possible θ (i.e., orientations from À90°to 90°) by extracting 1-dimensional slices through the 2-D spectra, equivalent to calculating a 1-dimensional topographic spectrum with specific orientations. Our 2-D spectra (e.g., Figure 2d) exhibit consistent "speckle" across the entire frequency space. This speckle represents noise that can be reduced by averaging over both orientations and wavelengths. We extracted amplitude profiles logarithmically sampled in wavelength along every orientation at 0.2°intervals. We then calculated 5°running means across orientations and smoothed the amplitudes across wavelengths with a locally weighted scatterplot smoothing filter implemented in Matlab. We found that a smoothing span of 1/5 of a decade in λ (1/5 of an order of magnitude; e.g., Figure 3a) offered a qualitatively good balance between noise suppression and wavelength resolution. The topographic amplitude as a function of wavelength for orientations spaced every 5°is shown in Figure 4a, centered on the orientation with peak amplitude.
Topography is naturally characterized by large amplitudes at long wavelengths and smaller amplitudes at shorter wavelengths [Sayles and Thomas, 1978]. For each of our DTMs, we identified the background relationship between wavelength (λ) and amplitude by calculating the median amplitude as a function of λ across all orientations, θ. We refer to these median amplitudes as the background spectra. Consistently oriented topographic features that exhibit periodicity (i.e., a high degree of spatial regularity) have amplitudes that rise above these median, background, topographic amplitudes. To specifically examine these periodic features, we normalized the amplitudes for each θ by the background spectra ( Figure 3b). Because the normalized amplitudes are highly sensitive to the background spectra, the median values that serve as the basis for normalization were smoothed with a span twice that used for each θ in the non-normalized plots (e.g., Figure 3a). This approach allows us to rigorously and automatically identify the orientations and wavelengths of periodic, spatially regular topography.
Strongly periodic features within the DTMs have large median-normalized amplitudes ( Figure 3b). To emphasize those orientations with the greatest normalized amplitudes, we identified all orientations with normalized amplitude reaching at least 50% of the largest normalized amplitude among all orientations and attaining normalized amplitudes >3 (these threshold choices are arbitrary and only for visualization). These orientations are colored consistently in Figures 2d and 3. All other orientations are dotted gray. Several examples detailing the results of our approach with simple, synthetic data are included in Figures S1-S6.

Resolution of the Analysis
Our ability to resolve the spectral characteristics of topography is constrained by the limitations of the input DTMs. As the input data are posted at (grid) horizontal resolutions ranging from 2 to 30 m (Table 1), the highfrequency limit beyond which we cannot identify unaliased topographic variations, the Nyquist frequency, varies from 0.25 to 0.017 m À1 , equivalent to λ = 4 to 60 m. At the other extreme, for long wavelengths, we approach a limit beyond which only one or a few λ fit within the extent of the 10 km DTMs. For this reason, we limit our discussion to those features with λ < 3 km; our interpretable range is therefore 60 m < λ < 3 km.
The vertical resolution of our multibeam bathymetry and LiDAR datasets is, in some cases, as coarse as 3 m (Table 1). However, the spatial averaging implicit within the spectral analysis allows us to resolve periodic features with much greater precision. Therefore, we use empirical assessments to demonstrate the vertical resolution of our analyses. Many of the amplitudes for orientations not associated with MSGLs cluster tightly around each other and the median amplitude and follow consistent slopes in the spectral log-log plots for λ < 500 m (e.g., Figure 3a). Similar, consistent slopes are commonly reported for a wide range of surfaces [Brown and Scholz, 1985;Sayles and Thomas, 1978]. We interpret these consistent slopes as indications that the spatially aggregated amplitudes of the DTMs are accurately characterized even to scales <0.01 m (less than the specified vertical resolution of an individual elevation measurement). The nature of the bed data collected from the extant Rutford Ice Stream limits our ability to fully analyze their spectra, as above. These bed-elevation data were collected along radio echo sounding profiles at 7.5 m intervals across the crests of the MSGLs, allowing for spectral analysis across crests of features with a wavelength as fine as 15 m. However, profile spacing along the crests of the MSGLs was 500 m. Although we used a cubic interpolation to resample these data to a 10 m grid, the original data resolution prevents us from characterizing periodicity with wavelength λ < 1 km for orientations not directly aligned with the MSGLs (normal to the radar profiles). In the absence of complete 180°orientation sampling, we are unable to normalize the spectra and identify maxima in λ in the same manner described above. Thus, the Rutford Ice Stream spectra are presented only in their unnormalized forms (Figure 5c).

MSGL Periodic Wavelengths and Orientations
The discrimination of MSGLs is enabled by comparison of topographic amplitudes in one direction with those of all other directions (Figure 3b). In the case of PI_MSGL_1, the most obvious MSGLs have a crest to crest spacing of 700 m and trend slightly west of due north (Figure 1b). Our spectral analysis bears this out, indicating that the dominant wavelength is 690 m and is oriented 12.6°W of N ( Figure 3b). However, our method also reveals subtlety in the DTM, which is not immediately apparent. Overprinting the MSGLs with 690 m spacing, we also find strong periodicity at 310 and 420 m wavelengths. The peaks and troughs in these normalized wavelengths are well defined, with obvious separation among the multiple peaks. A similar pattern of peaks and troughs in the normalized spectra is also found in adjacent orientations, with above-average amplitudes bracketing the dominant MSGL orientation by ±5°.
In a different MSGL setting (MT_MSGL_2), our analysis identifies a dominant wavelength of 310 m oriented at 131.4°E of N (Figure 4a). More subtle peaks in the spectra are found at longer and shorter wavelengths. In contrast, the LAR_fluvial dataset, a terrestrial landscape with no MSGLs but abundant fluvial erosion features, has no clear maxima or minima in periodic orientation or wavelength. Despite several-hundred-meter characteristic widths for the channels apparent in the DTM (Figure 1c), these features are not periodic and all median-normalized amplitudes are less than 3 (Figure 4b).
When we examine the full set of MSGLs, we find that multiple wavelength peaks within a single topographic dataset are common. These peaks in MSGL amplitude reach up to 42 times higher than the median amplitudes from all orientations (background amplitudes) (Table 1 and Figure 3) and are generally sharp, with widths at half-height on the order of the width of our low-pass filter (Figure 4c). Clear gaps in the normalized amplitudes exist between the peaks. Neither the peaks nor the gaps in amplitude are found at the same wavelengths among the different MSGL settings (Figures 4c and S7). The dominant wavelengths, those with the maximum normalized amplitudes, are found between 300 m and 1200 m, apart from the terrestrial LAR_MSGL, which is discussed in section 3.4. While the dominant wavelengths span a relatively narrow range, amplitudes reaching several times the background amplitude occur in many settings for λ > 2000 m and for λ < 100 m (Figures 2, 3b, 4a, 4c, and S7). The latter is toward the resolvable limit of short wavelength topography. However, these broad-and fine-scale corrugations have significantly weaker amplitudes than those at wavelengths between approximately 300 and 1200 m. Topographies without clear MSGLs, and even landscapes dominated by streamlined bedrock and crag-and-tail features, do not have dominant wavelengths characterized by normalized amplitudes as high as those from the MSGL landscapes (nowhere do they exceed 8) and have amplitude peaks that often span a wider range of wavelengths than the MSGLs (Figure 4d and Table 1).

Amplitudes of MSGLs and Other Topographies
Normalization of each orientation's spectrum by the background spectra renders the most periodic features, such as MSGLs, immediately apparent. However, normalization presents the amplitudes of periodic wavelengths only in a relative sense, obscuring the true height of MSGLs. Within the detrended data for PI_MSGL_1 (Figure 2b), trough-to-crest heights of some individual MSGLs exceed 10 m. This trough-to-crest measurement is twice the height of our "amplitudes," which are defined in the sense common to physics and signal processing as one half the trough-to-crest range. However, the amplitudes of many MSGLs in PI_MSGL_1 are clearly lower. For example, the average amplitude of the dominant λ = 690 m MSGLs, as calculated by our Fourier analysis, is 1.0 m (Figure 3a). This amplitude is 42 times greater than the median, background amplitude of 0.023 m at the same wavelength (Figures 3a and 3b).
Examination of the unnormalized spectra for each MSGL dataset reveals that, within their characteristic 300-1200 m wavelengths, MSGL amplitudes show a tendency to increase, with wavelength and average about 0.2 m at λ = 300 m and 1.0 m at λ = 1200 m ( Figure 5c). Thus, at the long end of the wavelength spectrum, these MSGL amplitudes are consistent with the 1-3 m average trough-to-crest ranges reported previously [Spagnolo et al., 2014].
The amplitude spectra of MSGLs share similar features to those of other landforms and produce a power law spectral signature, with A = a λ b , where A is amplitude and a and b are parameters fit to the data [Sayles and Thomas, 1978;Huang and Turcotte, 1990]. Our use here of true topographic amplitudes (Fourier coefficients), rather than power, differs slightly from convention, where power law spectra are shown through power spectral densities. Power in power spectral densities is related to the square of amplitude, and thus, the b values in our amplitude spectra are half as large as those in power spectra [cf. Perron et al., 2008]. The non-MSGL terrains analyzed here follow this general pattern of increasing amplitude with increasing wavelength, with b varying slightly among topographies between 1.5 and 2.3 ( Figure 5a). The spectra of MSGL topographies, along their dominant orientations (Figure 5c), are generally bracketed within a similar range of A(λ), from a few decimeters to a few meters, as apparent from their lying between the same black, power law lines drawn in all plots of Figure 5. However, the median spectra of MSGL background topographies (i.e., topography not oriented with the MSGLs themselves) are strikingly different from both the non-MSGL topographies and the dominant MSGL orientations. These background MSGL spectra have amplitudes depressed by a factor of 10 or more (at λ = 500 m) and slopes also depressed over a broad range of λ, from 100 m (below the typical λ of MSGLs) up to about 1300 m (greater than the maximum dominant λ of MSGLs) (Figure 5b). Thus, where MSGLs are present, the vertical relief of landscape features not oriented with MSGLs is considerably subdued relative to landscapes without MSGLs.
These results apply to all analyzed MSGL settings, including the extant Rutford Ice Stream bed. The exceptions to these general patterns are the three LAR datasets, from the terrestrial landscape of northern Alberta, and the iceberg furrows from the Antarctic marine settings (both PI_furrows and MT_furrows). The LAR landscapes have considerably smaller amplitudes than the other topographies (smaller a values), but with similar spectral slopes (b values of about 2). The two furrowed topographies have lower spectral slopes (b values of about 0.5), similar to the background spectra of MSGL datasets. These two exceptions are discussed in sections 3.4 and 4.2.2, respectively.

Downstream Variations in Amplitude and Wavelength Along the Bed of a Former Ice Stream
The extent of the bathymetry from the Getz paleo ice stream bed (Amundsen Sea Embayment, Antarctica, Figure 1a) allows us to identify how topographic spectral properties vary along over 50 km of former ice stream bed. This setting presents a morphological configuration that is typical of many Antarctic paleo ice stream beds, with upstream bedrock outcrops progressively transitioning downstream into crag-and-tail and MSGL features [Wellner et al., 2001]. The results of our approach show that the streamlined bedrock and/or crag-and-tail topography of the upstream portions (ASE_bedrock_1 and 2) consist of above-average median-normalized amplitudes, at many different orientations and wavelengths (Figures 6d and 6e). Specifically, the number of orientations with above-average median-normalized amplitudes spans up to 45°, a much greater value than the typical 20°for MSGL landscapes, while normalized amplitudes (up to 8) are much lower than those of the MSGLs (Figure 4c). The dataset furthest upstream (ASE_bedrock_1), largely composed of streamlined bedrock with a few MSGLs within a region to the NE, has highest normalized amplitudes at wavelengths between 1000 and 2000 m, as well as 200-500 m (Figure 6e). Downstream, the topography is dominated by crag-and-tail features (ASE_bedrock_2) and characterized  (Figure 6c). Very large median-normalized amplitudes are also present in the MSGL datasets further downstream but at clearly different wavelengths. In ASE_MSGL_2, we find peak amplitudes at~350,~600, and~1000 m wavelengths, while in ASE_MSGL_3, the furthest downstream dataset, only wavelengths between 450 and 600 m have the highest amplitudes (Figures 6a and 6b). Overall, a striking difference is found between the Fourier signatures of subglacially modified bedrock-cored landforms and MSGLs, with a clear tendency toward more strongly peaked wavelengths in the MSGLs downstream. Within the MSGL datasets, we find no consistent wavelengths even among datasets spaced just 10 km apart.
To evaluate the extent to which these dominant wavelengths vary along-flow, we resampled the ASE topography (Figure 6f) with 75% overlap between 10 km samples, spaced evenly every 2.5 km. These resampled data are identified as ASE_MSGL_0.25, _0.50, etc., from the transition out of the bedrock-dominated topography, down to ASE_MSGL_3. The orientations with the peak-normalized amplitudes are shown in Figure 7a. We find that, with this greater spatial resolution, individual dominant wavelengths correlate up and downstream and trends among these peak wavelengths are identifiable. We identify four trends in peak wavelength that appear to lengthen over 5, 8, 6, and 4 successive samples (red lines in Figure 7a). We also find weaker evidence of shortening or steady wavelengths (blue lines in Figure 7a). To identify the magnitude of these trends, we measure the wavelength of each local amplitude maximum and fit linear models to the base 10 logarithm of the wavelengths, as a function of position along the paleo ice stream. The slopes of these lines, bounded with 95% confidence intervals and identified from large to small wavelengths, reveal that these MSGL wavelengths increase at the rates of 2:3 þ1:8 À1:7 % km À1 , 0.95 ± 0.52% km À1 , 2.7 ± 1.0% km À1 , and 2:7 þ3:4 À3:3 % km À1 , respectively. Linear regression of wavelength on position yielded more widely varying estimates of the growth rate. To identify the most robust average relationship among the four trend lines, we employ a generalized linear mixed effect model, fit with a random effect of trend line. This modeling framework allows us to estimate a shared, average rate of downstream wavelength change, while allowing for different initial wavelengths. The shared average rate of wavelength increase is 1.7 ± 0.52% km À1 . We also identify two trends, spanning shorter portions of the MSGL field, that reveal steady wavelengths or wavelengths that decrease downstream. The steady wavelengths are found among those DTMs closest to the bedrock at λ = 1200 m, similar to the wavelengths of sculpted bedrock features. Toward the downstream end of the analyzed section of seafloor, we also identify a set of peak wavelengths that appear to decrease from 1500 to 1200 m over 5 km (three DTMs).
Among this set of resampled DTMs, we find that the unnormalized amplitudes ( Figure 7b) differ between those DTMs containing the largest proportion of bedrock and those DTMs dominated by MSGLs. ASE_bedrock_2, ASE_MSGL_0.25, and ASE_MSGL_0.50 reach amplitudes approximately half an order of magnitude greater than the MSGL DTMs. Among the MSGL DTMs, no along-flow variation in amplitude for λ < 600 m is found, while for λ > 600 m, there is a tendency toward reduced amplitudes downstream.

Present-Day Active and Paleo Terrestrial MSGLs
The unnormalized spectra of the Rutford Ice Stream are qualitatively similar to the spectra of the paleo ice streams (Figure 5c). In terms of their amplitudes and their spectral slopes, these extant features are indistinguishable from MSGLs formed by the paleo Pine Island, Marguerite Trough and Getz ice streams.
Most of the analyzed paleo ice stream beds belong to marine settings and are thus well preserved (Table 1). However, the Lower Athabasca dataset is characterized by a terrestrial paleo ice stream bed. The LAR_MSGL dataset consists of an MSGL field modified by terrestrial processes and exhibits a subdued Fourier spectral signature when compared with offshore counterparts. The peak amplitude for LAR_MSGL is only 3.9 times greater than background, lower than the typical normalized amplitudes of 10 or higher. It also occurs at an especially large wavelength (1600 m) (Figure 4). This is possibly due to poor preservation related to fluvial and hillslope erosion and deposition following MSGL formation. Alternatively, the LAR MSGLs may have formed under different boundary conditions than in the marine settings or had insufficient time to develop.

The Distinctive, Fourier Spectral Signature of MSGLs
The 2-D Fourier analysis has quantitatively demonstrated that MSGLs are characterized by (i) several dominant wavelengths whose spectral amplitude rises sharply above the background amplitude spectra, bounded by wavelength gaps with very low amplitudes; (ii) amplitudes of decimeters to a few meters; (iii) median-normalized amplitudes of MSGLs that are usually an order of magnitude larger than the amplitudes in all other orientations and up to 5 times larger than other, non-MSGL topographies; (iv) a narrow range of dominant orientations (<20°complete range); and (v) a limited interval of wavelengths with greatest amplitudes (~300-1200 m).
Some of these results are consistent with estimates from previous qualitative and quantitative studies. For example, MSGL heights of a few meters have been reported before [Spagnolo et al., 2014], although here we have revealed periodic features with average amplitudes of down to the scale of decimeters. The reported spacing between adjacent MSGLs from various settings [e.g., Ottesen et al., 2005] also generally fit within the 300-1200 m interval highlighted here. A previous study of glacial lineations across terrestrial Canadian Arctic sites, where 1-D Fourier analysis was applied to transects drawn across satellite images, also revealed comparable wavelength ranges (150-750 m) [Fowahuh and Clark, 1995]. However, this study is the first time a completely automated 2-D Fourier analysis is applied and to a wide variety of MSGL settings. Numerical ice PERIODIC TOPOGRAPHY OF ICE STREAM BEDS stream models that consider bed topography could adopt (i)-(v) as input data; numerical models attempting to evolve a sedimentary bed into landforms can use these results as tests for their outputs.
Characteristics (i)-(v) are found among all analyzed MSGL datasets and are distinct from those identified for other landforms analyzed in this study. Such a distinct signature offers some promise for the development of tools for the automated identification of MSGLs and the extraction of metrics without the need to map individual features. This is a significant advance, as mapping MSGLs has proved to be difficult, given the spatial continuity of these highly elongated corrugations and the care necessary to judge what constitutes a mappable feature [Spagnolo et al., 2014;Piasecka et al., 2016]. The automated identification of MSGLs based on their spectral signature in the frequency domain could also be extremely useful in the study of buried ice stream beds, such as those imaged from 3-D seismic data across glaciated continental margins [e.g., Dowdeswell et al., 2006;Graham et al., 2007;Andreassen and Winsborrow, 2009;Piasecka et al., 2016].

MSGLs As A Patterned Phenomenon and Implications for Formational Theories
The occurrence of a small number of dominant wavelengths within a narrow range of orientations demonstrates, quantitatively, that MSGLs represent a patterned, periodic topography, with regularly spaced landforms of similar shape and size. These characteristics are typical of spatially self-organized landscapes [Hallet, 1990].
A number of hypotheses have been proposed for MSGL formation. Proponents of subglacial megafloods suggest that MSGLs are formed by erosional vortices within turbulent sheet-floods capable of decoupling the ice sheet from its bed [Shaw et al., 2008]. The rilling instability hypothesis proposes that MSGLs are erosional landforms that emerge at the ice-water-bed interfaces from positive feedbacks in the unstable, coupled flow of ice, water, and sediment [Fowler, 2010]. MSGLs could also be the product of sediment deformation [Boulton, 1987;Clark, 1993]. The pressure-gradient hypothesis suggests that MSGLs are formed through gradients in ice pressure over deforming sediment, with MSGLs emerging from an uneven bed as a result of stoss side erosion and lee side deposition into cavities [Barchyn et al., 2016]. Finally, the groove-plowing hypothesis invokes erosion via basal ice keels plowing through soft sediment [Tulaczyk et al., 2001;Clark et al., 2003]. Although sedimentological evidence from MSGLs exposed on land in Arctic Canada suggests that some lineations have an (at least partially) erosional origin , other studies have shown that MSGLs are constructed by the continuous accretion of sediment . From our analysis, it is not possible to determine whether the MSGLs presented here are erosional, depositional, or a combination of both. However, our data lead to five considerations (detailed below), which are relevant to our understanding of subglacial processes and the formation of MSGLs.

Comparison Between MSGLs and Bedrock-Related Glacial Landforms
The 2-D Fourier spectral analysis has revealed that the erosion of streamlined bedrock by flowing ice does not produce the regular topographic signature that we see where MSGLs are present. Similarly, deposition in the lee of bedrock knobs, producing crag-and-tail features, fails to generate 2-D Fourier spectra similar to that of MSGLs. This is certainly the case for the Getz paleo ice stream bed, where MSGLs and bedrock-related features are only a few kilometers apart and yet have distinct spectra (Figures 6c and 6d). Deposition or erosion along an ice stream bed, when anchored to bedrock outcrops, does not generate the regularly spaced topography of MSGLs. For MSGLs to generate periodic topography, sediment must be able to move freely without fixed anchor points.
The streamlined bedrock and crag-and-tail datasets show some above average median amplitude values (colored lines in Figure 6d). However, these values are found across a wide range of orientations (up to 45°) and wavelengths (spanning an order of magnitude or more). Both these orientation and wavelength ranges far exceed those of MSGLs ( Figure 6). Bedrock-related topographies could reflect wavelengths and orientations not wholly dependent on the erosive action of ice flow. For example, they could represent the geomorphological expression of structures in the sedimentary bedrock or rock fracture patterns [Krabbendam and Bradwell, 2011]. However, the spread of large amplitude orientations for bedrockdominated terrains, roughly aligned with the overall S-N ice-flow direction, is likely to reflect the tortuosity of basal ice flow in the face of large-scale bedrock obstacles. These basal flow directions have the potential to offer a new constraint on the plasticity of basal ice.

10.1002/2016JF004154
Spectral analysis is also capable of resolving the signature of multiple bedforms within a single DTM. Within ASE_bedrock_1 (Figure 6e), we find features with normalized amplitudes near 5 spanning a wide range of wavelengths from 800-3000 m, reflecting the eroded bedrock. We also identify enhanced amplitudes over a 200-500 m wavelength range with a narrower range of orientations (20°, rather than the 50°range of orientations present at 800-3000 m), which are consistent with the spectral signature of MSGLs. These MSGLs are apparent in the NE portion of ASE_bedrock_1 (Figure 6f). This demonstrates that Fourier analysis can identify and distinguish two different landforms within the same setting, as well as their orientation ranges and roughness.

Downstream Variability of MSGL Wavelengths and Amplitudes
The groove-plowing hypothesis suggests that MSGLs may be formed by the plowing action of basal ice keels through soft sediment [Tulaczyk et al., 2001;Clark et al., 2003]. The keels are hypothesized to form either by convergent flow or from differential ice molding when the basal ice flows over outcropping bedrock. In this latter case, one would expect MSGLs found downstream of outcropping bedrock to have wavelengths comparable to those of the bedrock upstream [Clark et al., 2003]. The groove-plowing hypothesis also suggests that keels should progressively melt downstream because of frictional heating, thus broadening the MSGLs and reducing their amplitude, while their wavelength should remain constant [Clark et al., 2003]. These aspects were tested along the Getz paleo ice stream bed. Contrary to groove-plowing prediction, the large number of significant wavelengths, associated with modest (lower than 8 times the median) amplitudes of the streamlined bedrock and crag-and-tail features, does not match the few dominant wavelengths typical of the MSGLs located immediately downstream (Figures 6  and 7).
Rather than preserving the spacing imprint of upstream bedrock, we find that most of the multiple, distinct MSGL wavelength peaks within a DTM evolve toward longer wavelengths downstream. Over the approximately 20 km between ASE_MSGL_1 and ASE_MSGL_3, a given peak wavelength increases by an average of 40%. New peak wavelengths emerge at intervals downstream (at wavelengths between 200 and 300 m, the lower end of the MSGL range) and increase until they reach their maximum cutoff, near 1200 m (Figure 7). At the base of ice streams, both sediment and ice will move downflow over time.
Thus, through the assumption that downstream bedforms reflect temporal evolution, we conclude that individual MSGL spectra evolve continuously, even under conceivably steady conditions, just as a steady, turbulent, jet continuously casts off eddies that grow in scale prior to dissipating in a surrounding viscous fluid.
Two potential exceptions to the trend of downstream-increasing wavelength exist. One is found close to the upstream end of the MSGL field, near the outcropping bedrock, where a similar, peak wavelength near 1200 m is identified across four consecutive datasets. A similar dominant wavelength is found in the bedrock-rich dataset further upstream (ASE_bedrock_1) ( Figure 6). Thus, this exception might represent the only evidence of bedrock-related groove-plowing in our study. The other is toward the downstream end of the MSGL field, where wavelengths decrease over 5 km (across three consecutive datasets). This short-distance trend is found at the upper limit of the wavelength range of MSGLs and might be related to non-MSGL topography.
The amplitudes reveal different patterns for different wavelengths. Decreased downstream amplitudes, observed for wavelengths >600 m, are in agreement with predictions of the groove-plowing theory. However, we find no consistent change in MSGL amplitudes for wavelengths <600 m. Overall, while more MSGL settings should be analyzed to reach broader conclusions, most of the results presented here appear inconsistent with the groove-plowing theory, at least in its present form. Instead, these results provide support for theories that invoke the growth of dominant wavelengths via positive feedbacks occurring at an unstable ice-bed interface.

Instability As A Potential Ingredient in the Formation of MSGLs
Instability theories predict that there is a wavelength that grows most rapidly and tends to naturally dominate over others, with the system evolving toward a periodic landscape [e.g., Baas, 2002]. Such models have been developed for a variety of glacial bedforms, including flutes [Schoof, 2007], drumlins [Hindmarsh, 1999], and MSGLs [Fowler, 2010;Fowler and Chapwanya, 2014]. The distinct MSGL spectral signature, multiple high-amplitude wavelengths within the MSGL fields, and the downstream coarsening of MSGL wavelengths are consistent with hypotheses that see instability driven positive feedbacks playing a Journal of Geophysical Research: Earth Surface 10.1002/2016JF004154 key role in the formation and evolution of these sedimentary subglacial bedforms. However, no current theories predict multiple peak wavelengths within a single patch of ice stream bed. We suggest that the downstream increase in wavelength and the observation of multiple peak wavelengths are consistent with a continuously evolving pattern of bedforms. Thus, the spectral characteristics of an ice stream bed may not be strictly deterministic, even under steady conditions. An observed MSGL pattern could be created by many different ice stream characteristics. This challenges the assertion that MSGL characteristics may, at some point, be inverted for paleo ice-flow characteristics [Barchyn et al., 2016]. MSGL evolution may be affected by pattern coarsening, a phenomenon already suggested for the formation of dunes [e.g., Fourrière et al., 2010]. Instabilities emerging at small wavelengths may progressively grow in scale, increasing their wavelengths until the large wavelength limit (~1200 m), at which point they dissipate under the influence of other basal processes. In this view, the presence of multiple dominant MSGL wavelengths ( Figure 6) could indicate that new MSGLs with short wavelengths are continuously formed while others dissipate.

Relative Amplitudes of MSGLs and Their Background Topography
The results show that MSGLs are characterized by modest amplitudes, from a few decimeters to a few meters ( Figure 5c). Amplitudes of non-MSGL topographies are slightly higher, but mostly within the same order of magnitude ( Figure 5a). However, MSGL spectra are unique within our study due to the large difference between their amplitudes at dominant orientation and those at other orientations. This difference has the potential to arise through either the addition or subtraction of roughness (or a combination of the two) within the range of MSGL wavelengths. In the roughness addition hypothesis, MSGLs are visually apparent (Figures 1 and 2) because they have been constructed by erosion and/or deposition on an initially smooth topography. In the roughness subtraction hypothesis, MSGLs are visually apparent because, on an initially rough topography, roughness at all but a narrow range of orientations is removed.
The different non-normalized spectra can be interpreted through either lens. More analyses of current ice stream beds are required to determine whether MSGLs emerge through the addition or subtraction of roughness or through both processes. Regardless, the difference in ice stream bed amplitudes between ice-flow parallel and background (i.e., other) orientations (Figure 4c), along with the low spectral slopes (Figure 5b), quantitatively demonstrates that ice stream beds rich in MSGLs are characterized by reduced roughness, when compared to other topographies. MSGL terrains are so smooth that, even when the MSGLs are completely obliterated by the destructive action of iceberg keels (as in the two furrowed datasets), the original signature of low background roughness remains largely unaltered (Figure 5a). If MSGLs are an "inevitable" emergent feature of coupled ice-(water?)-till flow, then their relatively low-roughness configuration might represent the balance between MSGL's tendency to grow and the ice stream's tendency to reduce basal drag. This study quantifies amplitudes along and across the ice-flow direction for dominant wavelengths. As such, it provides a better understanding of form drag imposed by subglacial landforms and might be of relevance to inversions of ice-flow speed for basal properties [Joughin et al., 2004;Arthern et al., 2015].

MSGLs Differ From Typical Fractal Topographies
Topography is often described as a fractal phenomenon, in that the morphology of landforms maintains its scale invariant proportions [Mandelbrot, 1982]. This property is present in some of the non-MSGL landscapes that we analyzed, which show amplitudes increasing as a power law function of wavelength (Figure 5a). The exponents that best fit these non-MSGL amplitudes (1.5 < b < 2.3) overlap the range appropriate for fractal surfaces (1 < b < 2). The consistency of the scaling relationship over the nearly 2 orders of magnitude is also in line with the concept of fractals [Perron et al., 2008].
In contrast, deviation from this simple fractal scaling is found in all the analyzed MSGL datasets (Figure 5b), similar to the modification of amplitude spectra by periodic stream valleys [Perron et al., 2008]. In our MSGL datasets, the consistent power law relationship of the background topography is clearly interrupted within the specific wavelength interval of the MSGLs (300-1200 m), where amplitudes grow relatively slowly with wavelength ( Figure 5b). The pace at which MSGL spectral signals grow (i.e., the slope of the lines (b) in Figure 5b) is~0.5, significantly lower than the values of 1.5-2.3 found outside the MSGL wavelength interval and in all other non-MSGL settings we examined (Figures 6a and 6c). MSGL topography spectra also reveal a kink near 300 m, with lower b values for λ > 300 m and greater b values for λ < 300 m. These properties are inconsistent with those of fractals [Huang and Turcotte, 1990;Perron et al., 2008].

Conclusions
We have applied 2-D Fourier spectral analysis to 22 different topographies, mostly associated with ice streams. These datasets include areas rich in MSGLs and, in some cases, their bedrock-dominated areas upstream. Our technique allows for the objective quantification of landscape metrics without the need for landforms to be subjectively or individually mapped. The key conclusions are as follows: 1. MSGL fields consist of long corrugations with decimeter to meter amplitudes periodically spaced every 300-1200 m. The periodic topography of MSGLs consists of the superposition of multiple, similarly oriented corrugations with distinct wavelengths. MSGL patterns are spatially coherent, with regularly spaced landforms of similar shape and size and produce a distinct spectral signature that is not found in the other landforms analyzed in our study. For such a periodic topography to evolve, the anchoring effect of unevenly distributed bedrock knolls must play a minimal role in the formation of MSGLs. 2. A formation hypothesis involving positive feedbacks occurring in an unstable subglacial system (instability theory) is supported by the observation of distinct, periodic MSGL wavelengths. Furthermore, along the Getz paleo ice stream bed in Antarctica, we find (a) a downstream increase in wavelength; (b) the appearance of new, short-wavelength MSGLs downstream; and (c) a lack of relationship between wavelengths of MSGL spectra and the spectra of upstream bedrock. Together, these properties are compatible with an instability-driven formation of MSGLs and challenge theories that require the plowing of ice keels downstream of bedrock outcrops. 3. MSGL spectral amplitudes are comparable to the amplitudes of topographies without MSGLs. However, the presence of MSGLs is coincident with an order-of-magnitude lowering of topographic amplitudes in directions not oriented with MSGLs. This indicates that the roughness of MSGL topographies is lower than that of any other terrains we analyzed.
The distinct signature of MSGL topographies enables their automated identification. Our quantitative results may serve as a test for ice-water-bed interaction models or as input parameters for other numerical models seeking to understand the effect of topography on ice flow or attempting to evolve sedimentary bed landforms.