Self-afﬁne subglacial roughness: consequences for radar scattering and basal water discrimination in northern Greenland

. Subglacial roughness can be determined at a variety of length scales from radio-echo sounding (RES) data either via statistical analysis of topography or inferred from basal radar scattering. Past studies have demonstrated that subglacial terrain exhibits self-afﬁne (power law) roughness scaling behaviour, but existing radar scattering models do not take this into account. Here, using RES data from northern Greenland, we introduce a self-afﬁne statistical framework that enables a consistent integration of topographic-scale roughness with the electromagnetic theory of radar scattering. We demonstrate that the degree of radar scattering, quantiﬁed using the waveform abruptness (pulse peaki-ness), is topographically controlled by the Hurst (roughness power law) exponent. Notably, specular bed reﬂections are associated with a lower Hurst exponent, with diffuse scattering associated with a higher Hurst exponent. Abrupt waveforms (specular reﬂections) have previously been used as a RES diagnostic for basal water, and to test this assumption we compare our radar scattering map with a recent prediction for the basal thermal state. We demonstrate that the majority of thawed regions (above pressure melting point) exhibit a diffuse scattering signature, which is in contradiction to the prior approach. Self-afﬁne statistics provide a generalised model for subglacial terrain and can improve our understanding of the relationship between basal properties and ice-sheet dynamics.


Introduction
With the development of the newest generation of thermomechanical ice-sheet models, there has been a growing awareness that better constraining the physical properties of the glacier bed is essential for improving their predictive capability (e.g.Price et al., 2011;Seroussi et al., 2013;Nowicki et al., 2013;Shannon et al., 2013;Sergienko et al., 2014;Ritz et al., 2015;Cornford et al., 2015).Notably, the basal traction parameterisation -which encapsulates the thermal state, basal roughness, and lithology -is potentially the largest single geophysical uncertainty in projections of the response of ice sheets to climate change (Ritz et al., 2015).Distinction between frozen and thawed regions of the glacier bed is particularly important in constraining ice dynamics, since appreciable basal motion can only occur in regions where the glacier bed is wet (Weertman, 1957;Nye, 1970;Mac-Gregor et al., 2016).Airborne radio-echo sounding (RES) is the only existing remote sensing technique that can acquire bed data with sufficient spatial coverage to enable subglacial information to be obtained across the ice sheets (refer to Pritchard (2014) and Bamber et al. (2013a) for recent Antarctic and Greenland coverage maps).Often, however, there is great ambiguity in RES-derived subglacial information (Matsuoka, 2011), or RES-derived information is suboptimal for direct applicability in ice-sheet models (Wilkens et al., 2015).Subsequently, data analysis methods which seek to improve the clarity and glaciological utility of RES-derived subglacial information are undergoing a period of rapid development (e.g.Oswald and Gogineni, 2008;Li et al., 2010;Fujita et al., Published by Copernicus Publications on behalf of the European Geosciences Union. 2012; Wolovick et al., 2013;Schroeder et al., 2013Schroeder et al., , 2016;;Jordan et al., 2016).
RES data analysis methods for determining subglacial physical properties can be categorised in two ways: those which determine bulk properties (including the discrimination of basal water) and those which determine interfacial properties (subglacial roughness).Bulk material properties of the glacier bed can, in principle, be determined using the basal reflectivity (Bogorodsky et al., 1983;Peters et al., 2005;Jacobel et al., 2009;Schroeder et al., 2016).Performing basal reflectivity analysis on ice-sheet-wide scale is, however, greatly limited by uncertainty and spatial variation in englacial radar attenuation (Matsuoka, 2011;Matsuoka et al., 2012;MacGregor et al., 2012MacGregor et al., , 2015;;Jordan et al., 2016).In contrast to bulk properties, subglacial roughness analysis methods are (nearly) independent of radar attenuation.Subglacial roughness can be determined either via statistical analysis of topography (typically spectral analysis) (Taylor et al., 2004;Siegert et al., 2005;Bingham and Siegert, 2009;Li et al., 2010;Rippin, 2013) or inferred from the electromagnetic scattering properties of the radar pulse (Oswald and Gogineni, 2008;Schroeder et al., 2014;Young et al., 2016).Spectral analysis can provide valuable insight toward aspects past ice dynamics and landscape formation (Siegert et al., 2005;Bingham and Siegert, 2009;Rippin et al., 2014).However, since the technique is limited to investigating length scales greater than the horizontal resolution (typically ∼ 30 m or greater), the relevance of the method of informing contemporary basal sliding physics -which requires metre-scale roughness information (Weertman, 1957;Nye, 1970;Hubbard et al., 2000;Fowler, 2011) -remains unclear.Radar scattering is sensitive to the length scale of the electromagnetic wave (Shepard and Campbell, 1999) (∼ 1-5 m in ice for the majority of airborne sounders) and can potentially reveal finer-scale roughness information, including the geometry of subglacial hydrological systems (Oswald and Gogineni, 2008;Schroeder et al., 2013Schroeder et al., , 2015;;Young et al., 2016).High reflection specularity, such as occurs from deep (> 10 m) subglacial lakes (Oswald and Robin, 1973;Gorman and Siegert, 1999;Palmer et al., 2013), has been proposed as a RES diagnostic for basal water (Oswald andGogineni, 2008, 2012).
Degrees of radar scattering can be mapped either using the waveform properties of the bed echo -e.g. the waveform abruptness (pulse-peakiness) (Oswald andGogineni, 2008, 2012) -or by constraining the angular distribution of scattered energy -e.g. the specularity content (Schroeder et al., 2013;Young et al., 2016).Maps of both scattering parameters indicate defined spatial patterns but, to date, have not been integrated with topographic-scale roughness analysis (horizontal length scales ∼ 10 s of metres and upwards).As such, there is a knowledge gap regarding the topographic control upon radar scattering.Observations indicate that subglacial roughness exhibits self-affine (fractal) scaling behaviour over length scales from ∼ 10 −3 to ∼ 10 2 m (Hub-bard et al., 2000;MacGregor et al., 2013).Self-affine scaling corresponds to when the vertical roughness increases at a fixed slower rate than the horizontal length scale, following a power-law relationship that is parameterised by the Hurst exponent (Malinverno, 1990;Shepard et al., 2001).It is observed for a wide variety of natural terrain (Smith, 2014), including the surface of Mars (Orosei et al., 2003), volcanic lava (Morris et al., 2008), and alluvial channels (Robert, 1988).If widely present, the self-affinity of subglacial roughness poses a challenge for integrating topographic roughness with existing glacial radar scattering models (Berry, 1973;Peters et al., 2005;MacGregor et al., 2013;Schroeder et al., 2015).This is because these are statistically stationary models which assume that roughness is independent of horizontal length scale, and hence an artificial scale separation between high-frequency roughness and low-frequency topography is present (Berry, 1973).Radar scattering models with nonstationary, self-affine statistics naturally incorporate the multiscale dependence of roughness and are in widespread use in other fields of radar geophysics (e.g.Shepard and Campbell, 1999;Franceschetti et al., 1999;Campbell and Shepard, 2003;Oleschko et al., 2003).
In this study, we explore the connection between selfaffine subglacial roughness and radar scattering using recent airborne Operation IceBridge (OIB) RES data from the north-western Greenland Ice Sheet (GrIS).Firstly we review the theory of self-affine roughness statistics, using examples from ice-penetrating radargrams and bed elevation profiles to demonstrate its applicability to subglacial terrain (Sect.2).We then outline analysis methods that enable topographic roughness and radar scattering (quantified using the waveform abruptness) to be extracted from RES flighttrack data (Sect.3).A self-affine radar scattering model, adapted from planetary radar sounding (Shepard and Campbell, 1999;Campbell and Shepard, 2003), is then used to predict the relationship between the Hurst exponent and waveform abruptness (Sect.4).We then present maps of the RESderived roughness and scattering data for the northern GrIS and compare the spatial distribution with bed topography (Bamber et al., 2013a) and a recent prediction for the basal thermal state (MacGregor et al., 2016) (Sect. 5.1).The radar scattering model is then used to quantify self-affine topographic control upon radar scattering, via the Hurst exponent (Sect.5.2).The statistics of the RES-derived data in predicted thawed and frozen regions of the glacier bed are then analysed (Sect.5.3), with the purpose of testing the basal water discrimination method by Oswald andGogineni (2008, 2012), which assumes a specular scattering signature is present.Finally, we discuss the wider consequences of our study, including subglacial landscape classification, the relationship between bed properties and ice-sheet dynamics, basal thaw/water discrimination, and radar scattering theory applied to RES (Sect.6).

Overview
Statistical methods to calculate the Hurst exponent, and thus to quantify self-affine scaling behaviour, are well established in the earth and planetary science literature (Malinverno, 1990;Shepard et al., 2001;Kulatilake et al., 1998;Orosei et al., 2003).These space-domain methods extract the Hurst exponent using the variogram (roughness verses profile length) and deviogram (roughness versus horizontal lag).Our motivation for use of these methods, rather than the spectral (frequency-domain) methods previously applied in studies of subglacial roughness (Taylor et al., 2004;Siegert et al., 2005;Bingham and Siegert, 2009;Li et al., 2010;Rippin, 2013), is that they better reveal self-affine scaling behaviour (Turcotte, 1992;Shepard et al., 1995Shepard et al., , 2001)).Since the theory of self-affine roughness and related space-domain methods are not widely discussed in the glaciological literaturethe only example being MacGregor et al. ( 2013) -we now provide a review of the key concepts.

Interfacial roughness parameters
Topographic roughness can be measured by means of statistical parameters that are, in general, a function of horizontal length scale (Shepard et al., 2001;Smith, 2014).Two different interfacial roughness parameters -the root mean square (rms) height and rms deviation -are typically employed in self-affine roughness statistics (Shepard et al., 2001).The rms height is given by where N is the number of sample points within the profile window of length L, z(x i ) is the bed elevation at point x i , and z is the mean bed elevation of the profile.ξ represents the standard deviation in bed elevation about a mean surface and models the topographic roughness as a Gaussian-distributed random variable (Orosei et al., 2003).The rms deviation is given by where x is the horizontal step size (lag).ν has a particular significance in the parameterisation of radar scattering models with self-affine statistics (Shepard and Campbell, 1999;Campbell and Shepard, 2003), and we focus upon this roughness parameter when integrating topographic-scale roughness with radar scattering data.The rms slope, which is proportional to the rms deviation, is also widely used in selfaffine statistics, but we do not do so here.

Self-affine scaling behaviour and the role of the Hurst exponent
Self-affine scaling is a subclass of fractal scaling behaviour and can be parameterised using the Hurst exponent, H (Malinverno, 1990;Shepard et al., 1995Shepard et al., , 2001)).H quantifies the rate at which roughness in the vertical direction increases relative to the horizontal length scale (and is defined for 0 ≤ H ≤ 1).For a self-affine interface the following powerlaw relationships hold: and where L 0 is a reference profile length and x 0 is a reference horizontal lag (Shepard and Campbell, 1999;Shepard et al., 2001).Three limiting cases of self-affine scaling are typically discussed (Shepard and Campbell, 1999).Terrain with H = 1 (where the roughness in the vertical direction increases at the same rate as the horizontal length scale) is referred to as "self-similar".Terrain with H = 0.5 (where the roughness in the vertical direction increases with the square root of horizontal length scale) is referred to as "Brownian".
Terrain with H = 0 (where the roughness in the vertical direction is independent of horizontal length scale) is referred to as "stationary".For a stationary (H = 0) interface it follows from Eqs. (3) and (4) that ξ and ν are independent of L and x respectively (i.e. the roughness parameters are independent of horizontal length scale).
We will later demonstrate that subglacial terrain exhibits near-ubiquitous self-affine scaling behaviour with pronounced spatial structure and variation for H . Examples of OIB ice-penetrating radargrams (Z scopes) (Paden, 2015) and associated bed elevation profiles for terrain with different H are shown in Fig. 1.Clear differences are apparent between the different terrain examples.The black (H ≈ 0.9) terrain (Fig. 1a) and red (H ≈ 0.7) terrain (Fig. 1b) are between Brownian and self-similar scaling behaviour.This terrain exhibits "persistent trends", where neighbouring measurements tend to follow a general trend of increasing or decreasing elevation (refer to Shepard and Campbell (1999) for a full discussion).A feature of terrain with higher H is that it tends to appear relatively rough at larger length scales (low frequency) and smooth at smaller length scales (high frequency).By contrast, the green (H ≈ 0.3) terrain (Fig. 1d) is in the sub-Brownian scaling regime and exhibits "anti-persistent trends", where neighbouring measurements tend to alternate between increasing and decreasing elevation.A feature of lower H terrain such as this example is that it tends to have similar roughness across length scales.The blue (H ≈ 0.5) terrain (Fig. 1c) is close to an ideal Brownian surface and exhibits no overall persistence (with some www.the-cryosphere.net/11/1247/2017/The Cryosphere, 11, 1247-1264, 2017 sections of the profile following an increasing/decreasing elevation trend and other sections alternating).The 10 km profile windows in Fig. 1 represent the length of flight-track data over which H is calculated (see Sect. 3.2).

Calculation of the Hurst exponent using the variogram and deviogram
In order to calculate H , and identify the scale regime over which glacial terrain exhibits self-affine behaviour, ξ and ν are plotted as functions of L and x, respectively, on doublelogarithmic-scale plots, referred to as the variogram and deviogram (Kulatilake et al., 1998;Shepard et al., 2001).Variogram and deviogram plots for ξ(L) and ν( x) for the four terrain examples in Fig. 1 are shown in Fig. 2a and b respectively.It follows from Eqs. (3) and (4) that, upon this doublelogarithmic scale, a straight line relationship is predicted for glacial terrain that is self-affine with the gradient equal to H .In practice, a single self-affine relationship only holds over a limited scale regime and a "break-point" transition is often observed (Shepard et al., 2001).We describe how we assess the break points for glacial terrain in Sect.3.2, along with further details regarding the application of the variogram and deviogram to along-track RES data.Figure 2 clearly demonstrates the significance of the Hurst exponent and horizontal length scale when assessing the relative roughness of different terrain.For example, the black (H ≈ 0.9) terrain is rougher than the red (H ≈ 0.7) terrain at larger length scales but is smoother at smaller length scales.The space-domain variogram and deviogram have an approximate correspondence to the frequency-domain power spectrum (Turcotte, 1992;Shepard et al., 1995Shepard et al., , 2001)).In The plots correspond to subglacial terrain profiles in Fig. 1.The Hurst exponent is estimated from the linear gradient of the first five data points (indicated by dashed lines).These space-domain plots are (approximate) equivalents to frequency-domain roughness power spectra, and smaller length scales correspond to higher frequencies.frequency space, self-affine scaling occurs when the power spectrum, S, has a relationship of the form S(k) ∝ k −β , where k is the spatial frequency and −β is the spectral slope.The relationship between β and H is dimensionally dependent and for along-track data is given by H = 1 2 (β − 1) (Turcotte, 1992).Despite this correspondence, the space-domain methods are recommended to calculate H as they are less noisy and less likely to bias slope estimates than the power spectrum method (Shepard et al., 1995).The study by Hubbard et al. (2000) observed self-affine scaling in the roughness power spectrum over length scales from ∼ 10 −3 to ∼ 10 m for different sites across recently deglaciated terrain in the immediate foreground of Tsanfleuron glacier, Switzerland.Their range for measured values of β corresponds to 2.27 < β < 2.48, which implies H ≈ 0.7.Paden, 2015).For the flight lines considered, the alongtrack resolution after synthetic aperture radar (SAR) processing and multi-looking is ∼ 30 m with an along-track-sample spacing of ∼ 15 m (Gogineni et al., 2014).The 2011 and 2014 field seasons were used since they have a higher alongtrack resolution than other recent field seasons and hence enable a clearer connection to be made between radar scattering and topographic-scale roughness.
The study focused on flight-track data from north-western Greenland and encompassed measurements close to three deep ice cores: Camp Century, NEEM, and NorthGRIP (Fig. 3).The first reason for selection of this region is that the data coverage for the 2011 and 2014 field seasons is of high density relative to most other regions of the ice sheet.The second reason is that confidence regarding the basal thermal state is high near to the ice cores and thus enables the validity of the basal water RES analysis by Oswald andGogineni (2008, 2012) to be tested.
Measurements from MCoRDS are supplied as data products with different levels of additional processing (Paden, 2015).Level 2 data correspond to ice thickness, ice surface, and bed elevation data and are used to calculate topographicscale roughness and the Hurst exponent (Sect.3.2).Details regarding the semi-manual picking procedure are described by Paden (2015), and only the highest-quality picks were used.Level 1B data correspond to radar-echo strength profiles and are used to extract the waveform abruptness parameter from the bed echo (Sect.3.3).Basal reflection values can also be extracted from Level 1B data, but we do not do this here because we do not wish to bias our interpretation due to uncertainty in radar attenuation.The preprocessing of the combined channel Level 1B data is also described by Paden (2015).Sequentially this involves channel compensation between each of the antenna phase centres, pulse compression (using a 20 % Tukey window in the time domain), coherent-averaging of the channels, SAR processing with along-track frequency window, channel combination, and waveform combination.

Determination of topographic roughness and Hurst exponent from Level 2 data
The along-track spacing (∼ 15 m) of the Level 2 data is half the horizontal resolution (∼ 30 m), which represents the spacing at which bed elevation measurements are considered as independent.Therefore, to remove local correlation bias, the Level 2 data were down-sampled, considering every second data point (corresponding to a ∼ 30 m along-track spacing).Each flight track was then divided into 10 km alongtrack profile windows, as shown in the examples in Fig. 1a.

T. M. Jordan et al.: Self-affine subglacial roughness
The windows overlap with a sample spacing of 1 km, with the centre of each window defined to be the point to which H and the roughness parameters are geolocated.This "moving window" approach was employed as it enables greater continuity in the estimates for H . Prior to estimating H , ξ(L) and ν( x) were computed following Eqs.( 1) and ( 2) respectively.These calculations used the "interleaving" sampling method described in Shepard et al. (2001), which enables all of the data points to be sampled effectively.The windowing method is similar to that described in Orosei et al. (2003) for the self-affine characterisation of Martian topography, where a non-overlapping 30 km window was assumed.The choice of 10 km for the profile window and 1 km for the effective resolution represents a good tradeoff between resolution and the smoothness of the derived data fields.
In this study we are interested in calculating H at the length scale of the Fresnel zone (∼ 100 m), since this enables the most accurate parameterisation of the radar scattering model described in Sect. 4. Additionally, due to the break point transitions that occur at larger length scales, the focus on smaller length scales is a robust approach to calculate H (Shepard et al., 2001).For the data we consider, the lower bounds of the horizontal length scales are ∼ 90 m for ξ(L) (since three elevation measurements are the minimum required to calculate ξ(L) using Eq. 1) and ∼ 30 m for ν( x).ν( x) therefore better enables the estimation of H at smaller length scales and we primarily focused on the deviogram method (Fig. 2b).Additionally, as suggested in Fig. 2, the relationships for ν( x) are, in general, significantly smoother than ξ(L).The upper length scales in the deviogram and variogram were set to be x = 1 km and L = 1 km respectively, which follows from the recommendation by Shepard et al. (2001) that at least 10 independent sections of track are used in the calculations.As shown in Fig. 2, the gradients (H ) were calculated using the first five data points (which, for the deviogram, are over the range x ∼ 30-150 m).Self-affine scaling behaviour often extends beyond these smaller length scales and we estimated the break points for ξ(L) and ν( x) using a segmented linear regression procedure.Briefly, this involved firstly calculating the gradient (H ) for the first five data points.Additional data points at increasing length scales were then added into each linear regression model, and the gradient was recalculated.Finally, break points in the linear relationship were identified by testing if the new gradient exceeded a specified tolerance from the original estimate.

Determination of waveform abruptness from Level 1B data
The post-processing of the Level 1B data (analysis of the basal waveform) uses the procedure described in Jordan et al. (2016), which, in turn, is largely based upon Oswald and Gogineni (2008).Firstly, this involved performing an alongtrack average of the basal waveform, where adjacent basal waveforms are stacked about their peak power values and arithmetically averaged.This averaging approach is phaseincoherent and acts to smooth power fluctuations due to electromagnetic interference (Oswald and Gogineni, 2008).The size of the averaging window varies as a function of Fresnel zone radius, and subsequently each along-track averaged waveform corresponds to approximately a separately illuminated region of the glacier bed (see Jordan et al., 2016, for details).The degree of radar scattering is quantified using the waveform abruptness where P peak is the peak power of the bed echo and P agg is the aggregated power, which is calculated by a discrete summation of the bed-echo power measurements in each depth range bin.P agg was introduced by Oswald and Gogineni (2008) since, based upon energy conservation arguments, it is argued to be more directly related to the predicted (specular) reflection coefficients than equivalent peak power values.In radar altimetry, the waveform abruptness is commonly called "pulse peakiness" (e.g.Peacock and Laxon, 2004;Zygmuntowska et al., 2013).Observed values of A range from ∼ 0.03 to 0.60, and in Sect.4.3 we theoretically constrain the maximum value to be ∼ 0.65.Three examples of basal waveforms, along with their corresponding A values, are shown in Fig. 4. Higher A values are associated with specular reflections from smoother regions of the glacier bed (e.g. the blue waveform), whilst lower A values are associated with diffuse reflections from rougher regions (e.g. the green waveform) (Oswald and Gogineni, 2008).The positions of the peak power were established by firstly using Level 2 data picks, then applying a local re-tracker to centre over the peak power.When calculating the summation for P agg (both fore and aft of the peak power so as to best capture the energy contained in the echo envelope), a signal-noise-ratio threshold was implemented by testing for decay of the peak power to specified percentage above the noise floor.Thresholds of 1, 2, and 5 % were considered and 2 % was found to give the best coverage, whilst excluding obvious anomalies.Due to this qualityfiltering step there are therefore sometimes small gaps in the along-track A data.
As RES over ice employs a nadir-facing sounder, the scattering contribution toward the waveform abruptness is mainly from coherent reflection (as opposed to side-looking SAR instruments which would be mainly diffuse scattering).Whether coherent pre-processing (either coherent presumming of Doppler focusing) of the raw data acts to increase or decrease the value of A depends upon the exact character and roughness of the surface.As a first example, if the specular/nadir component of the echo is assumed to be coherent, whilst the diffuse/off-nadir component is assumed to be incoherent (e.g.Grima et al., 2014), then coherent processing would cause the specular component of the signal to increase with coherent gain but not the diffuse (incoherent)

Bed-echo waveforms
Figure 4. Examples of bed-echo waveforms and their abruptness (pulse peakiness).Observed values for A range from ∼ 0.03 (associated with diffuse scattering) to ∼ 0.60 (associated with specular reflection).For the purpose of comparative plotting, the waveforms are normalised about their peak power values with the sample bin of the peak power set to zero.The sample bin spacing corresponds to a depth-range spacing of ∼ 2.81 m in ice.
signal.Therefore the measured A value would decrease with gain.As a second example, if both the specular/nadir and diffuse/off-nadir components of the echo are assumed to be coherent (e.g.Schroeder et al., 2013Schroeder et al., , 2015)), then for small SAR processing angles (coherent pre-summing) the waveform abruptness should be largely unaffected.However, for larger angles (exceeding the angle spanned by the specular component of the echo in the scattering function) the A value will decrease with coherent pre-processing.
The basal waveform (and hence the calculated values of A) results from a superposition of along-track and cross-track energy (Young et al., 2016).Subsequently, the anisotropy of radar scattering (and inferences regarding the anisotropy of subglacial roughness) is not explicitly revealed by A. Hence, the studies of Oswald andGogineni (2008, 2012) treat A as an isotropic parameter, and we follow this approach here.
4 Radar scattering model for self-affine roughness

Overview
The waveform abruptness has previously been discussed without reference to roughness statistics, and here we do this using a self-affine radar scattering model.Radar scattering models from natural terrain fall into two different categories: "coherent", which incorporates deterministic phase interference, and "incoherent", which incorporates random phase interference (Ulaby et al., 1982;Campbell and Shepard, 2003;Grima et al., 2014).Coherent scattering models are applicable where the reflecting region is orientated nearly perpendicular to the incident pulse (the nadir regime) and the reflecting region is fairly smooth at the scale of the illumi-nating wavelength (Campbell and Shepard, 2003), which is normally assumed to be a good approximation for the RES of glacier beds (e.g.Peters et al., 2005;MacGregor et al., 2013;Schroeder et al., 2015).Volume (Mie) scattering is typically neglected from basal RES scattering analysis and would hypothetically require scatterer dimensions of the order of the radar wavelength (∼ 0.5 to 5 m dependent on the bed dielectric and radar system).This neglection of volume scattering is justified given the ∼ 10 −6 to 10 −3 m scale of water pore radii in typical bed materials (Nimmo, 2004).Moreover, even in the extreme case of planetary ice regoliths (which are colder than terrestrial ice and will therefore sustain larger heterogeneities), scatterer dimensions are ∼ 10 −3 to 10 −2 m and volume scattering losses are small (Aglyamov et al., 2017).
Below we describe and adapt a coherent scattering model, first developed for the nadir regime of planetary radar sounding measurements, which incorporates self-affine roughness statistics (Shepard and Campbell, 1999;Campbell and Shepard, 2003).The model is parameterised using the Hurst exponent values derived from the subglacial topography (Sect.3.2) and thus enables a connection to be made between the topographic roughness and radar scattering.Coherent scattering models can be used to model a decrease in specularly reflected power as a function of rms roughness (Berry, 1973;Peters et al., 2005), and this is the central aspect of the model which we focus upon here.Specifically, we show that, under assumptions of energy conservation, this power decrease can be used to theoretically predict the relationship between the Hurst exponent and waveform abruptness.

Modelling the coherent power
The physical assumptions behind the self-affine scattering model are summarised in Shepard and Campbell (1999).The central assumption that differentiates the model from coherent stationary (H = 0) models (Berry, 1973;Peters et al., 2005;MacGregor et al., 2013;Grima et al., 2014;Schroeder et al., 2015) is that the rms height increases as a function of radius, r, about any given point, following the self-affine relationship where ν λ = ν( x = λ) is the wavelength-scale rms deviation.Equation ( 6) assumes radial isotropy for H and ξ .Since we are focusing upon constraining the (near-) isotropic abruptness parameter, this is a justifiable approximation.The statistical distribution for ξ(r) is assumed to be Gaussian, which is similar to most H = 0 models (but with an additional radial dependence).Via ν λ , the self-affine model is explicitly formulated with respect to the horizontal scale of rms roughness.The radio wavelength of MCoRDS in ice is ∼ 0.87 m, and hence wavelength-scale rms deviation is approximately equivalent to metre-scale rms deviation.An un-avoidable caveat to the parameterisation of the radar scattering model using Eq. ( 6) is that the H values derived from the topography (length scale ∼ 30-150 m) are extrapolated downwards to the wavelength scale.
An expression for the radar backscatter coefficient (radar cross section per unit area) is then derived by considering a phase variation, 4πξ(r)  λ , integrated across the Fresnel zone (Shepard and Campbell, 1999;Campbell and Shepard, 2003).For nadir reflection the radar backscatter coefficient is given by where r = r λ is the wavelength-scaled radius, rmax is the wavelength-scaled radius of the illuminated area (the Fresnel zone), and R e is the reflection coefficient for the electric field (Campbell and Shepard, 2003).The coherent power, P , can then be obtained by dividing Eq. ( 7) by 4π 2 r2 max (a geometric factor which follows from the backscatter coefficient of a flat conducting plate; Ulaby et al., 1982) to obtain For the case where H = 0, ξ(r) in Eq. ( 6) is independent of radius.It follows that ξ 2 = 1 2 ν 2 λ and the exponent in Eq. ( 8) is also independent of radius, which gives Equation ( 9) is the same power decay formula as coherent H = 0 models (Peters et al., 2005;MacGregor et al., 2013;Grima et al., 2014;Schroeder et al., 2015), where it is sometimes multiplied by a first-order Bessel function (which enables some of the incoherent energy contribution to be captured; MacGregor et al., 2013).Thus the stationary limit of the self-affine model is consistent with previous glacial basal scattering models.It is clear that the coherent power for the self-affine model, Eq. ( 8), has two roughness degrees of freedom: H and ν λ , which can be conceptually related to the gradient and the intercept of the deviogram (Fig. 2).This contrasts with the stationary model, Eq. ( 9), which has one degree of freedom: ξ .

Predicted relationship between the Hurst exponent and waveform abruptness
The utility of the waveform abruptness in quantifying different degrees of scattering rests upon the assumption that the majority of the overall energy is contained within the echo envelope (Oswald and Gogineni, 2008).In other words, it is assumed that, for reflection from the same bulk material, the aggregated/integrated power from a rough interface (ν λ > 0) is equivalent to the peak power from a given smooth interface; i.e.P agg ≈ P (ν λ = 0).This energy equivalence was demonstrated to hold well for the waveform processing procedure and Greenland RES systems by Oswald and Gogineni (2008).It follows from this energy equivalence that the abruptness, A, can be expressed in terms of the coherent power, Eq. ( 8), as where C is a proportionality constant that corresponds to the theoretical maximum abruptness value, which occurs when the radar pulse is specularly reflected and P agg = P peak .For a perfectly specular reflection the pulse is the shape of compressed chirp (absolute value of a sinc function with the width determined by the signal bandwidth).If the depthrange sample spacing of the waveform (Fig. 4) were the same as the depth-range resolution then C would be near unity.However, C can be estimated from the ratio of the sample spacing (∼ 2.8 m) to the range resolution (∼ 4.3 m) to give C ∼ 0.65.Finally, substituting Eq. ( 8) into Eq.( 10) gives As is the case for P in Eqs. ( 8) and (11) A has two roughness degrees of freedom: H and ν λ .Shepard and Campbell (1999) note that the primary dependence for P (and hence A) is upon H , with a weaker secondary dependence upon ν λ .In order to illustrate this dependency, we consider first the relationship between A and H for fixed ν λ (Fig. 5a) and secondly the relationship between A and ν λ for fixed H (Fig. 5b). Figure 5a demonstrates that higher values of ν λ (the black curve) result in negligible A for all but the lowest values of H . Intermediate values of ν λ (the red and blue curves) exhibit a sharp transition from higher to lower values of A as H increases.Low ν λ (the green curve) has high A for all H . Figure 5b demonstrates a monotonic decrease in A with ν λ for each value of H , with the decay length decreasing rapidly with increasing H .It is important to note that the predictions of the self-affine radar scattering model are consistent with the specular RES scattering signature that we would expect from electrically deep subglacial lakes.Under the self-affine roughness framework, a large geometrically flat feature such as a lake would have a negligible value of H and ν λ .This scenario occurs for the low H limit of the green curve in Fig. 5a, where predicted values for A are ∼ 0.65 (corresponding to a perfectly specular reflection).
The physical explanation for the strong dependence of the coherent power upon H , and the relationships which we observe in Fig. 5, is discussed by Shepard and Campbell (1999)  and Campbell and Shepard (2003).It relates to the fact that significant coherent returns can only occur from annular regions where ξ(r) < λ 8 (the Rayleigh criterion).It follows from Eq. ( 6) that high values of H lead to a rapid increase in roughness with radius that rapidly exceeds this threshold.Subsequently, for high H interfaces, the roughness at the wavelength scale, ν λ , must be a couple orders of magnitude smaller than the Rayleigh criterion to enable significant coherent returns (i.e.non-negligible A).The curves in Fig. 5 assume rmax = 100 (corresponding to a Fresnel zone radius ∼ 115 m for the ice wavelength ∼ 0.87 m).In general, the relationships in Fig. 5 are insensitive to this choice of radius.This is because the radii of the coherent annular regions are typically significantly less than the Fresnel zone and thus act as the dominant length scale for the integration limit in Eq. ( 11).

Results
Firstly, we describe maps for the rms deviation and Hurst exponent (topographic-scale roughness) and the waveform abruptness (radar scattering) in the northern Greenland (Sect.5.1).In this analysis we compare the RES-derived data with the Greenland bed digital elevation model (DEM) (Bamber et al., 2013a) and the predicted basal thermal state (MacGregor et al., 2016).Secondly, by comparing the theoretical predictions of the self-affine radar scattering model with the observed relationship between the Hurst exponent and waveform abruptness, we quantitatively assess topographic control upon radar scattering (Sect.5.2).Thirdly, we perform a statistical analysis of the RES-derived data in predicted thawed and frozen regions of the glacier bed (Sect.5.3), which enables us to assess the validity of the basal water discrimination algorithm in Oswald andGogineni (2008, 2012).Finally, we present uncertainty estimates for the RES-derived data (Sect.5.4).

Maps for self-affine roughness and radar scattering in northern Greenland
In Fig. 6 flight-track maps for the RES-derived roughness and scattering data are compared with the Greenland bed DEM (Bamber et al., 2013a) and the predicted basal thermal state (Fig. 11 in MacGregor et al., 2016).The flighttrack maps all demonstrate a high degree of spatial structure, with some notable correlations present (both between each other and the DEM).There is a clear inverse relationship between the rms deviation, ν (shown at two different length scales in Fig. 6a and b), and the waveform abruptness, A (Fig. 6c), with higher abruptness (specular reflections) present in smoother regions of the ice-sheet bed and lower abruptness (diffuse scattering) present in rougher regions.For example, smoother regions (lower ν, higher A) occur for flight tracks in the region inland from the settlement of Qaanaaq and around Camp Century (including the green profile, Fig. 1d), around the NorthGRIP ice core, and a region ∼ 150 km ENE of the NEEM ice core.Whilst these smoother regions are at a range of bed elevations (ranging from ∼ 800 m NE of Qaanaaq to around sea level in the interior), they are all spatially correlated with flatter bed topography (Fig. 6e).Correspondingly, many rougher regions (higher ν, lower A) are spatially correlated with more complex topography -e.g. the region of the ice sheet inland from the Melville Bugt coast (including the red profile, Fig. 1b).However, some rougher regions of the bed have a less obvious correlation with higher contour gradients -e.g. the flatter regions inland from the Humboldt glacier.Pronounced spatial variation in the Hurst exponent, H , is evident in Fig. 6d.H also has a inverse relationship with A and spatially correlates with the bed topography in a similar manner to ν.In other words, lower H is associated with higher A and flatter regions of the bed -e.g.near Camp www.the-cryosphere.net/11/1247/2017/The Cryosphere, 11, 1247-1264, 2017 Century -whilst higher H is associated with lower A and generally more complex bed topography -e.g.inland from the Melville Bugt coast and inland from Ryder glacier (including the black profile, Fig. 1a).In Sect.5.2 a quantitative assessment of this relationship is made using the radar scattering model.The simple notion that, at the topographic scale, rougher regions of the bed correspond to higher H can be related back to the power-law scaling relationship in the deviogram (Fig. 2b).The length scales for the rms deviation maps ν( x = 30 m) in Fig. 6a and ν( x = 150 m) in Fig. 6b are chosen as they are the lower and upper bounds in the deviogram calculation for H .It is notable that, despite the clear spatial variation in H in Fig. 6d, the overall spatial distributions for ν( x = 30 m) and ν( x = 150 m) are remarkably similar.Thus, from a purely visual inspection of ν at differ-ent length scales, the pronounced spatial variation in H is not immediately apparent.
The basal thermal state prediction by MacGregor et al. ( 2016) (Fig. 6f) represents an up-to-date best estimate for the GrIS at a 5 km resolution.It is based upon a trinary classification: likely thawed/above pressure melting point (red), likely frozen/below pressure melting point (blue), and uncertain (grey).The mask was determined using four independent methods: thermomechanical modelling of basal temperature, basal melting inferred from radiostratigraphy, surface velocity, and surface texture.The mask is therefore independent of our RES-derived data fields.There are some obvious correlations between the basal thermal state prediction and the RES-derived roughness and scattering data.For example, many predicted thawed regions toward the margins -e.g. the region of the ice sheet inland from the Melville Bugt coastcorrespond to rougher terrain (higher H and ν) and diffuse scattering (lower A).However, there are regions of predicted thaw that demonstrate the opposite behaviour (lower H and ν and higher A) -for example, the two interior regions previously identified as smooth around the NorthGRIP ice core and the region ENE of NEEM.The scattering signature of predicted thawed regions is therefore non-distinct and can be either specular or diffuse.Predicted frozen regions tend to be smoother with specular reflections (higher A), although it is clear that spatial variation is present with some regions exhibiting more diffuse scattering (lower A).Section 5.3 provides a more detailed statistical analysis.
There are some clear discontinuities in the flight-track maps for ν and H in Fig. 6.These can be explained by either roughness anisotropy or the self-affine terrain model breaking down in certain regions (e.g. a sharp terrain discontinuity such as a subglacial cliff).By contrast the map for A is smoother, which is consistent with its interpretation as an isotropic scattering parameter.

Statistics for topographic control upon radar scattering and comparison with radar scattering model
Before we consider a quantitative comparison between the predictions of the radar scattering model and the RESderived data, we first summarise the statistics for the Hurst exponent, H .The total frequency distribution for H , corresponding to the flight-track data in Fig. 6d, is shown in Fig. 7a.The distribution is divided into three categories: (i) H > 0.75 ("high" H ), (ii) 0.5 < H ≤ 0.75 ("medium H "), and (iii) H ≤ 0.5 ("low" H ), which we later use to compare with the radar scattering model predictions.These categories correspond to approximately 30, 50, and 20 % of the total data respectively.Approximately 0.1 % of the H estimates are > 1 and none of the H estimates are < 0, representing near-ubiquitous self-affine scaling behaviour (0 < H < 1).An overall negative skew for the distribution of H is observed with a mean value of 0.65, indicating that the majority of the subglacial terrain along the flight tracks lies between Brownian (H = 0.5) and self-similar (H = 1) scaling regimes.The spatial coverage of the radar flight tracks in Fig. 6d is, however, more comprehensive in regions of higher H . Thus the mean value and skew of H in Fig. 7a are likely overestimates and underestimates of true (equal area) averaged values for the region of the northern GrIS in Fig. 3.The self-affine coherent scattering model (Sect.4) predicts that there are two roughness degrees of freedom that control A: H (the primary control) and ν λ (the secondary control).At metre scale, ν λ is significantly smaller than the alongtrack resolution (∼ 30 m) and therefore cannot be observed directly.Additionally, given the theoretically predicted primary dependence of A upon H , a natural starting point is to compare with the observed relationship between A and H (Fig. 6).Based upon the assumption that ν λ varies spatially, a statistically distributed inverse relationship between H and A is predicted which corresponds to the family of predicted curves in H -A space in Fig. 5.This approach assumes a downward extrapolation of H from the topographic scale to the wavelength scale in the radar scattering model.
In order to test this prediction, we considered the statistics of three separate A distributions for each H category, which are shown for high H in Fig. 7b, medium H in Fig. 7c, and low H in 7d.A nearest-neighbour interpolation was used to pair each A value (∼ 100-150 m along-track spacing) with each H value (1 km along-track spacing).The lowest mean value, smallest variance, and strongest positive skew are observed for the high-H category.This supports the general prediction in Fig. 5 that higher A values (specular reflections) are suppressed in regions of higher H , with lower A values (diffuse scattering) being more probable.The highest mean value, greatest variance, and weakest positive skew are observed for the low-H category.Again, this supports the prediction in Fig. 5 that A is less constrained in regions of lower H , with a tendency toward higher values (specular reflections).As would be expected, the A-distribution statistics for the medium H category lie between the high-H and low-H categories with intermediate mean values, variance, and skewness.Finally, the observed values of A in Fig. 7 range from ∼ 0.03 to 0.60, which is in agreement with the theoretically constrained maximum value of 0.65.

Statistics in thawed and frozen regions
Here we summarise the statistics of the RES-derived roughness and scattering data in predicted thawed and frozen regions of the glacier bed, with an overall purpose of testing the basal water discrimination algorithm by Oswald andGogineni (2008, 2012).Conceptually, their approach assumes that water in thawed regions has a similar RES signature to deep subglacial lakes which exhibit brighter and more specular reflections than surrounding regions (e.g.Oswald and Robin, 1973;Gorman and Siegert, 1999;Palmer et al., 2013).In their algorithm wet regions are discriminated if (i) the relative bed reflectivity is above a threshold (using an attenuation model where the attenuation rate has an inverse relationship with surface elevation) and (ii) the abruptness is also above a threshold (around 0.3).Thus, in their approach, high abruptness (specular reflections) is a necessary, but not sufficient, criterion for identifying basal water.A further feature of their approach is that spatial continuity for water is imposed, i.e.only larger-scale regions (∼ 100s of km 2 and upwards) are considered.
The distributions for all RES-derived data exhibit pronounced statistical differences between thawed and frozen regions (Fig. 8).The mean value for H in thawed regions is 0.74 with a strong negative skew (Fig. 8a), whereas the mean value for H in frozen regions is 0.54 with a weak negative skew (Fig. 8b).The mean value for ν( x = 30 m) in www.the-cryosphere.net/11/1247/2017/The Cryosphere, 11, 1247-1264, 2017 Distribution for Hurst exponent thawed regions is 6.36 m, which is over double the mean value of 2.80 m in frozen regions.A qualitatively similar distinction between thawed and frozen regions is also present for ν( x = 150 m), with a mean value of 21.7 m in thawed regions and 7.2 m in frozen regions (not shown).The thawed distribution for A is similar to the high-H category in Fig. 7b, with a mean A value of 0.165 and strong positive skew.
The frozen distribution is similar to the low-H category in Fig. 7d with a mean A value of 0.264 and a weak positive skew.These statistics demonstrate a contradiction with the basal water discrimination algorithm of Oswald andGogineni (2008, 2012).Lower abruptness (diffuse scattering) is more common in thawed regions where basal water is likely to be present.Moreover, the necessary high-abruptness (specular reflections) condition for water is generally not satisfied (particularly at the larger spatial scales that were considered by Oswald andGogineni (2008, 2012) when mapping basal water).

Uncertainty and consistency of RES-derived data
In RES data analysis, cross-over distributions at flight-track intersections can give an indication of uncertainty based upon internal consistency (e.g.MacGregor et al., 2015;Jordan et al., 2016).However, due to the anisotropy in Fig. 6d, cross-over analysis for H (ν( x)) cannot be applied directly.Hence repeat estimates were made using the variogram to calculate H (ξ(L)) (i.e.calculating H using rms height).The map for H (ξ(L)) (not shown) has a similar spatial distribution as Fig. 6d but with greater high-frequency noise ap-parent.Differencing the estimates as H (ν( x))-H (ξ(L)) and performing cross-over analysis gives a mean bias of −0.026 and a standard deviation of 0.10 (10 % of the parameter range).The small mean bias is potentially explained by the variogram estimates being at a slightly larger length scale (L ∼ 90-210 m).Additional cross-over analysis using different profile window sizes (e.g. 15 km) confirms that 0.10 serves a reasonable estimate for the uncertainty of H . Since A is assumed to be isotropic, the uncertainty can be estimated via cross-over analysis of flight-track intersections.This gives a cross-over standard of ∼ 0.05 (again ∼ 10 % the parameter range).
As part of the analysis we also considered estimation of the breakpoint transitions for H (ξ(L)) and H (ν( x)) using the segmented linear regression procedure described Sect.3.2.The exact values of the breakpoints depend upon how strict the stopping criterion is, so here we just discuss some general trends.Firstly, the self-affine scaling relationships often extend over a much greater length scale than the upper length scale used in the calculation of H (often over 500 m as occurs in Fig. 2).Secondly, the breakpoints for H (ν( x)) generally occur at greater length scales than for H (ξ(L)). Thirdly, the break points for both H (ν( x)) and H (ξ(L)) tend to be greater toward the ice-sheet margins where H is higher.

Discussion
Our results demonstrate that self-affine scaling behaviour is a near-ubiquitous property of the subglacial topography of northern Greenland.Moreover, there is both spatial structure and variability in the Hurst exponent, which can range from being near-self similar (H ≈ 1) to sub-Brownian (H < 0.5).The Hurst exponent is valuable as it provides a way to integrate maps of topographic-scale roughness metrics (e.g.rms height and rms deviation) and maps of radar scattering parameters (e.g. the waveform abruptness), which provide finer-scale roughness information.Notably, theoretical predictions and observations both demonstrate that higher values of the abruptness (specular reflections) are suppressed in rougher regions of the bed with a higher Hurst exponent.Additionally, extended continuous regions of higher abruptness are generally limited to occur in smoother regions with a lower Hurst exponent.This finding implies that maps of radar scattering information -including the waveform abruptness parameter in this study and in Oswald andGogineni (2008, 2012) and the specularity content in Schroeder et al. (2013) and Young et al. (2016) -will benefit from analysis that incorporates self-affine topographic control.
The Hurst exponent provides information about the relationship that exists between vertical roughness and the horizontal length scale.Whilst it is related to the slope of the roughness power spectrum, past spectral analysis of glaciological terrain tends to obscure this information (since an integrated "total roughness" metric is typically used) (Taylor et al., 2004;Siegert et al., 2005;Bingham and Siegert, 2009;Li et al., 2010;Rippin, 2013).Subsequently, the Hurst exponent represents new subglacial roughness information that could potentially be utilised much more widely than our current application in constraining radar scattering.For example, planetary scientists have previous employed the Hurst exponent in a geostatistical classification of Martian terrain (Orosei et al., 2003).Interestingly, the spatial distribution of the Hurst exponent for the Martian surface has a similar level of spatial variation and coherence to what we observe for glacial terrain.Additionally, the distribution of H for Martian terrain is skewed toward higher, self-similar values with nearcontinuous regions of lower H limited to mid-latitude plains.For Greenland, this self-affine statistical landscape classification could be integrated with existing knowledge of geology (e.g.Henriksen, 2008) and larger-scale landscape features including subglacial drainage networks (Cooper et al., 2016;Chu et al., 2016;Livingstone et al., 2017) and palaeofluvial canyons (such the "mega canyon" feature observed www.the-cryosphere.net/11/1247/2017/The Cryosphere, 11, 1247-1264, 2017 in Fig. 6e, which has Petermann glacier as its modern-day terminus) (Bamber et al., 2013b).
The Hurst exponent has previously been shown to play a dynamical role in the flow resistance of alluvial channels (Robert, 1988).Whilst basal sliding is clearly a different physical phenomena -modulated by enhanced plastic flow and regelation (Weertman, 1957;Nye, 1970;Hubbard et al., 2000;Fowler, 2011) -it is possible that the Hurst exponent may provide a useful radar-derived parameter for our understanding of geometric control upon this process.In Sect.5.1 and 5.3 we observed that toward the ice-sheet margins, such as inland from the Melville Bugt coast, predicted thawed regions are characterised by higher (often near self-similar) values of H .One could therefore speculate that the persistent behaviour associated with high-H interfaces (neighbouring points follow a similar elevation trend; Sect.2.3) could act to promote basal sliding.However, as is widely acknowledged, attributing a direct link between subglacial roughness and contemporary ice dynamics is a complex topic (Siegert et al., 2005;Bingham and Siegert, 2009;Rippin et al., 2014).Therefore, as with other measures or basal roughness, the spatial variation in the Hurst exponent is likely to also originate from different glaciological processes at a variety of spatial scales, including erosion and deposition.Additionally, we recommend that future works which investigate the connection between the Hurst exponent and glaciological processes should be discussed with reference to anisotropy and flow direction.
The statistical analysis of the waveform abruptness in predicted frozen and thawed regions (Sect.5.3) demonstrates that, overall, very different RES scattering signatures are present than assumed by Oswald andGogineni (2008, 2012).Firstly, the majority of the predicted thawed regions have lower abruptness (diffuse scattering).In their algorithm, this would correspond to false-negative detection of basal water (since the necessary high abruptness condition is not satisfied).Secondly, high abruptness is often present in predicted frozen regions, many of which are interpreted as wet by Oswald and Gogineni (2008Gogineni ( , 2012) (e.g.some of the region of higher abruptness near to the Camp Century ice core, which at high bed elevation is likely to correspond to harder bedrock).It is, however, important to note that some of the smoother regions discriminated as wet by Oswald andGogineni (2008, 2012) are consistent with basal thermal state prediction by MacGregor et al. ( 2016) (e.g.near NorthGRIP).Radar bed reflectivity was also used by Oswald andGogineni (2008, 2012) in their discrimination of thawed beds.However, since these original studies, the role that uncertainty in radar attenuation plays in biasing the spatial distribution of radar bed reflectivity has become much better understood (Matsuoka, 2011;MacGregor et al., 2012;Jordan et al., 2016).For example, if an attenuation model has a constant systematic bias in attenuation rate, then there will be an ice-thickness-correlated bias in estimated bed reflectivity (Jordan et al., 2016).Thus, spatially correlated bias in the attenuation model is one explanation for why elevated reflectivity was observed in some predicted frozen regions.Additionally, geological transitions, between less-reflective sediment and more-reflective bedrock (see Bogorodsky et al. (1983) and Peters et al. (2005) for reflectivity values) could also play a role in complicating the analysis.
Subglacial hydrological systems are understood to produce more complex and variable scattering signatures than the specular lake-like reflection assumed by Oswald andGogineni (2008, 2012).For example, concentrated hydrological channels act as an anisotropic rough surface capable of orientation-dependent scattering (Schroeder et al., 2013;Young et al., 2016).Additionally, due to scattering from the lake bottom and related interference effects, shallower (depth < 10 m) subglacial lakes can produce diffuse scattering (Gorman and Siegert, 1999).Whilst the majority of the thawed regions have lower abruptness, there are some smaller, localised patches of higher abruptness present in Fig. 6c.These regions are consistent with the presence of deep lake-like water (in the sense that specular reflections are observed in a region predicted to be above pressure melting point).However, because the frozen abruptness distribution in Fig. 8f indicates that basal water is not required to produce highly specular reflections, it is not possible to confirm this without additional analysis.This is because the frozen abruptness distribution in Fig. 8f indicates that basal water is not required to produce highly specular reflections, and thus smooth regions of bedrock may be responsible for the high abruptness.The presence of at least some localised patches of high abruptness in thawed regions is consistent with the recent discovery of two small subglacial lakes in north-western Greenland of ∼ 8 and ∼ 10 km 2 in extent (Palmer et al., 2013).More generally, however, the relative rarity of high abruptness in thawed regions is in agreement with hydrological potential analysis (Livingstone et al., 2013), which predicts that deep subglacial lakes are both rare and small in the north-west of the GrIS.Instead, channelised drainage networks -such as the system recently identified beneath Humboldt glacier (Livingstone et al., 2017) -are likely to be common in thawed regions (and are consistent with the generally diffuse scattering signature that we observe).
The anisotropy of the Hurst exponent was not considered in the radar scattering model, which was justifiable because we were interested in understanding how the Hurst exponent relates to the (near-) isotropic waveform abruptness.However, in certain regions of the ice sheets, basal radar scattering is known to be highly anisotropic, as revealed by maps of the specularity content for Thwaites glacier (Schroeder et al., 2013) and Byrd glacier (Young et al., 2016).Thus a clear direction of future research would be to modify the self-affine radar scattering model (Sect.4) to take into account anisotropy in H and then to compare this model with maps for the specularity content.The pronounced spatial heterogeneity for H implies that estimation of roughness statistics from H = 0 radar scattering models (Eq.9; Berry, 1973;Ulaby et al., 1982;Peters et al., 2005;Grima et al., 2014) may give erroneous results, particularly when comparing the overall spatial distribution between regions with different H values.Additionally, the radar scattering model is formulated with respect to wavelength-scale (approximately metre-scale) roughness and thus provides a way to estimate metre-scale roughness (i.e.given A and H obtain an estimate for ν λ in accordance with the curves in Fig. 5).This could have important glaciological consequences, since the physical processes which influence basal sliding operate at the metre scale (Weertman, 1957;Nye, 1970;Hubbard et al., 2000;Fowler, 2011).
Finally, geostatistically based interpolation methods which employ aspects of self-affine statistics (Goff and Jordan, 1988) have found recent application in generating synthetic subglacial topography (Goff et al., 2014).The selfaffine characterisation of subglacial topography described here informs such techniques and, in turn, could be used to inform the ice-sheet-wide interpolation of future Greenland (Bamber et al., 2013a;Morlighem et al., 2014) and Antarctic (Fretwell et al., 2013) subglacial digital elevation models.

Summary and conclusions
In this study we used recent OIB RES data to demonstrate that subglacial roughness in northern Greenland exhibits self-affine scaling behaviour, with pronounced spatial variation in the Hurst (roughness power law) exponent.We modified a planetary radar scattering model to predict how the Hurst exponent exerts control upon the degree of scattering, which we parameterised using the waveform abruptness.We then demonstrated an agreement between the predictions of the radar scattering model and the statistically distributed inverse relationship that is observed between the Hurst exponent and waveform abruptness.This enables us to conclude that self-affine statistics provide a valuable framework in understanding the topographic control which influences ice-penetrating radar scattering from glacier beds.Self-affine statistics also provide a generalised model for subglacial terrain and in the future could be used to further explore the relationship between bed properties, ice-sheet dynamics, and landscape formation.
An additional glaciological motivation behind our study was to establish whether the waveform abruptness could be used to aid in the discrimination of basal water (and to test the prior assumption that subglacial hydrological systems in Greenland produce abrupt bed echoes; Oswald andGogineni, 2008, 2012).To do this we compared our RES-derived data fields with a recent basal thermal state prediction for northern Greenland (MacGregor et al., 2016).The analysis demonstrated that thawed regions of the glacier bed have statistically lower values of the waveform abruptness than frozen regions (more diffuse scattering).The simple explanation is that many thawed regions are relatively rough with a higher Hurst exponent, whilst many frozen regions are relatively smooth with a lower Hurst exponent.This finding should not be viewed as a new RES diagnostic for basal water (since deep subglacial lakes do have the specular signature proposed by Oswald andGogineni, 2008, 2012).However, it indicates that the diagnostic in Oswald andGogineni (2008, 2012) is likely to yield both false negatives (failing to identify water in rougher regions and where hydrological systems have more complex scattering signatures) and false positives (identifying some smoother frozen regions as wet).

TFigure 1 .
Figure 1.Example radargrams (top panel) and 10 km bed elevation profiles (bottom panel) for subglacial terrain with different Hurst exponent, H : (a) H ≈ 0.9 (near self-similar), (b) H ≈ 0.7 (between Brownian and self-similar), (c) H ≈ 0.5 (Brownian), and (d) H ≈ 0.3 (sub-Brownian).The location of the profiles are shown in Fig. 3. Evident in the radargrams are the surface reflection (pink line), the bed reflection (red line), and reflections from internal layers in ice.The bed elevation profiles are linearly detrended about zero with horizontal resolution ∼ 30 m.The horizontal-vertical aspect ratio of the bottom panels differs between (a), (b) and (c), (d) by a factor of ∼ 10.

Figure 2 .
Figure 2. (a) Variogram for rms height, ξ , versus profile length L (log-log scale).(b) Deviogram for rms deviation, ν, versus horizontal lag, x (log-log scale).The plots correspond to subglacial terrain profiles in Fig.1.The Hurst exponent is estimated from the linear gradient of the first five data points (indicated by dashed lines).These space-domain plots are (approximate) equivalents to frequency-domain roughness power spectra, and smaller length scales correspond to higher frequencies.

3
Analysis of RES data 3.1 Ice-penetrating radar system and coverage region The airborne RES data used in this study were collected by the Center for Remote Sensing of Ice Sheets (CReSIS) within the OIB project, over the months March-May in years 2011 and 2014.For all measurements the radar instrument, the Multichannel Coherent Radar Depth Sounder (MCoRDS), was installed upon a NASA P-3B Orion aircraft.The sounder has a frequency range from 180 to 210 MHz, corresponding to a centre wavelength ∼ 0.87 m in ice.After accounting for pulse shaping and windowing, this results in a depth-range resolution in ice of ∼ 4.3 m (Rodriguez-Morales et al., 2014;

Figure 3 .
Figure 3. Data coverage map for OIB flight tracks and region of interest.The locations of the Camp Century, NEEM, and North-GRIP ice cores are indicated, along with the terrain profile sections in Fig. 1.

Figure 5 .
Figure 5. Parametric dependence of the self-affine radar scattering model.(a) Abruptness, A, as a function of the Hurst exponent, H , for sections of constant wavelength-scale rms deviation, ν λ .(b) A as a function of ν λ for sections of constant of H .The plots illustrate primary dependence for A upon H and secondary dependence for A upon ν λ .High A is suppressed for high H except in the case of exceptionally small ν λ .

Figure 7 .
Figure 7. Relationship between Hurst exponent, H , and waveform abruptness, A (corresponding to flight-track data in Fig. 6).(a) Total distribution for Hurst exponent.(b) Abruptness distribution for high H , (H > 0.75).(c) Abruptness distribution for medium H , (0.5 < H 0.75).(d) Abruptness distribution for low H , (H ≤ 0.5).The observed distributions in (b), (c), and (d) confirm the theoretical prediction of the self-affine radar scattering model that a statistically distributed inverse relationship exists between H and A.

Figure 8 .
Figure 8. Distributions from basal RES analysis in thawed and frozen regions of the northern GrIS (corresponding to flight-track data in Fig. 6): (a) Hurst exponent, H , in thawed regions; (b) H in frozen regions; (c) rms deviation, ν( x = 30 m), in thawed regions; (d) ν( x = 30 m) in frozen regions; (e) abruptness, A, in thawed regions; (f) A in frozen regions.The data subsets correspond to the red (thawed) and blue (frozen) regions of the map in Fig. 6f.