Automated mapping of near bed radio-echo layer disruptions in the Greenland Ice Sheet

One of the key processes for modulating ice ﬂow is the interaction between the ice and the bed, but direct observations of the subglacial environment are sparse and diﬃcult to obtain. In this study we use information from an extensive radio-echo sounding dataset to identify areas of the Greenland Ice Sheet where internal layers have been inﬂuenced by near-bed processes. Based on an automatic algorithm for calculating the slope of the internal radio-echo layers, we identify areas with disrupted layer stratigraphy. We ﬁnd that large parts of the northern portion of the ice sheet are inﬂuenced by locally conﬁned mechanisms that produce up-warping or folds in the layer stratigraphy inconsistent with the surface and bed topography. This is particularly evident at the onset of ice streams, although less dynamic areas close to the ice divide also contain imprints of layer disturbances. Our results show that the disturbances are found in many different ﬂow and thermal regimes, and underscore the need to understand the mechanisms responsible for creating them.


Introduction
During the past several decades, the Greenland Ice Sheet (GrIS) has been extensively surveyed by ice-penetrating radar (also referred to as radio-echo sounding: RES) (e.g., Gudmandsen, 1975;Gogineni et al., 2001).This has enabled improved ice-volume estimates through the years (e.g., Bamber et al., 2001Bamber et al., , 2013a)), and detailed maps of regional and local bed topography (e.g., Bamber et al., 2013b;Gogineni et al., 2014).In addition to measurements of ice thickness, the echograms often also reveal internal reflecting horizons in the ice.These features, here referred to as internal layers, represent a powerful source of information on ice-sheet history and dynamics.The internal layers are the dielectric signature of changes in density, impurity content, acidity and/or crystal orientation (e.g., Millar, 1981;Fujita et al., 1999;Eisen et al., 2007).Several studies have demonstrated that the stratigraphy of internal layers can be used to reconstruct past ice flow, surface mass balance and basal melt rates if the layers are assumed to be isochrones, i.e., of the same age (e.g., Fahnestock et al., 2001;Conway et al., 2002;Ng and Conway, 2004;Siegert et al., 2005;Buchardt and Dahl-Jensen, 2007;Leysinger Vieli et al., 2011).Further, theoretical studies (Parrenin et al., 2006) have demonstrated that the layer slopes also contain information on ice flow, and Ng and King (2011) constructed a mathematical framework to in-vestigate how slope information propagates in near-surface radar data.
Recently, an extensive study by MacGregor et al. (2015) presented a dataset of traced layers covering most of the echograms available from the GrIS, and linked the traced internal layers to icecore sites, where the age of the ice is known.Thus age-horizons have been constructed for most of the ice sheet, providing an important tool for constraining the past dynamics of the GrIS.The combined tracing of internal layers and construction of radiostratigraphy is a laborious and time-consuming task.Therefore, alternative approaches to avoid directly tracing the layers have been developed in several previous studies.This includes the development of semi-automatic methods for picking layers in near-surface RES data (Onana et al., 2015), for automatically identifying subsurface signals in RES data from Mars radar sounders (Ferro andBruzzone, 2012, 2013;Freeman et al., 2010), and for tracing a distinct pattern of internal layers by fitting a shape function (Karlsson et al., 2013).Yet again, other studies have tried to circumvent the layer tracing problem by instead quantifying the continuity of the layers (Karlsson et al., 2012) or calculating the slope of the layers (Sime et al., 2011;Panton, 2014).
Finally, RES data have also been used to directly or indirectly extract information about the basal conditions of ice sheets (Fahnestock et al., 2001;Carter et al., 2007;Buchardt and Dahl-Jensen, 2007;Oswald and Gogineni, 2012).This is aided by the recent improvements to RES equipment and processing (Rodriguez-Morales et al., 2014), making it possible to sound the deepest internal structure of the GrIS with great clarity.Some radio-echograms show substantial disturbances to the internal layers (often in the form of up-welling and/or folding), and in some cases the disturbances take up more than half of the ice sheet column, with a horizontal extent on the order of ten kilometers.Bell et al. (2014) proposed that these features are correlated with the occurrence of refrozen meltwater.In contrast, the stratigraphically warped bottom section of the NEEM ice core was explained by folding due to the differences in rheological properties of the Eemian and glacial ice (Dahl-Jensen et al., 2013).And most recently, Wolovick et al. (2014) suggested that the features could be formed by non-uniform and time-variable basal sliding.
In short, internal layers of ice sheets have been observed to be influenced by processes acting near the bed, resulting in a stratigraphy that cannot be explained by surface or bed topography alone.MacGregor et al. (2015) used the term 'disrupted radiostratigraphy' to describe the RES observations of disturbances in the internal layer stratigraphy apparently induced near the bed.We extend on this term and use 'unit of disrupted radiostratigraphy', or UDR, in order to avoid implying a specific formation process and to indicate that the disruptions have a finite physical extent.
Like MacGregor et al. (2015), we utilize information from the internal layers in the GrIS to map the signature of UDRs.In contrast to methods where the basal conditions are observed by the quantitative properties of the echogram, e.g., where the presence of liquid water is detected using the strength and specularity of the RES signal (Carter et al., 2007;Oswald and Gogineni, 2012;Schroeder et al., 2013Schroeder et al., , 2015)), we instead study the problem based on interpretation of internal layer slopes.
Our analysis is based on an automated method for extracting information on the internal layer slopes (Panton, 2014), and on a simple mathematical framework demonstrating how changes in the basal environment can leave an imprint on the internal layers.

Radio-echo sounding data
We rely on a publicly available dataset acquired over the years 1999-2014 by the Center for Remote Sensing of Ice Sheets (CReSIS; University of Kansas).The data acquisition campaigns have covered almost all of the GrIS allowing for the construction of, for example, detailed maps of bed topography (Bamber et al., 2013a) and internal layer depths (MacGregor et al., 2015;Karlsson et al., 2013).In the dataset we are using, the bedrock and surface have already been picked by CReSIS.
The data for this study were acquired using two different and continuously developed radar systems.The radar system operating from 1999 to 2009 had a centre frequency of 150 MHz, and a bandwidth of 17 MHz with a peak transmit power of 200 W. From 2010 onwards, the centre frequency was changed to 195 MHz with a bandwidth of 30 MHz and with increasing peak transmit power up to 1200 W.More information on the radar systems is provided by Gogineni et al. (2001) and Rodriguez-Morales et al. (2014).

Methods
Our automatic detection of UDRs is based on the observation that the disturbances exhibit increasing layer slopes with depth.To further develop this hypothesis, we discuss the evolution of sloping internal layers using a simple theoretical framework in the section below.Then we briefly present the automatic layer slope algorithm first published in Panton (2014) and finally we show how the layer slope field can be used to detect the signature of UDRs.Our methods described below are comparable to what was developed by MacGregor et al. (2015) and the differences between these will be discussed later.

Theoretical framework
The depth of an internal layer is determined by the ice thickness, surface accumulation, basal melting, spatial gradients in ice flow velocity, and ice sheet evolution.The conditions at the bed may affect the internal layers through variations in basal melt rate and changes in ice flow mode.In this context the basal melt rate influences the vertical velocity by dragging the layers downwards in the ice column for increasing melt rates (e.g., Fahnestock et al., 2001;Buchardt and Dahl-Jensen, 2007).Changes in ice flow mode are exemplified by the "Weertman effect" (Weertman, 1976) where the internal layers move downwards when there is a change from no sliding to sliding, or upwards when the change is from sliding to no sliding.The latter case has also been referred to as "sticky spots" (Alley, 1993), and has been in observed numerous places particularly in West Antarctica (e.g., Anandakrishnan and Alley, 1994;Joughin et al., 2006).
We calculate the depth of an internal layer by considering the motion of particles in a two-dimensional velocity field {u(x, z, t), w(x, z, t)}, where x is the horizontal coordinate, z is the vertical coordinate and t is time.We define the direction of horizontal ice flow to be along the x-axis.A line connecting two particles P 1 and P 2 will form an angle α with the horizontal, defined by tan α = dz dx , where dx = x P 2 − x P 1 is the horizontal distance between the particles and dz = z P 2 − z P 1 is the vertical distance (cf.Fig. 1A, inset).If the particles are of the same age, the line connecting them will by definition be an isochrone.Generally, two particles will remain at the same vertical distance dz relative to each other if they experience the same vertical velocity w.Thus, if the vertical velocity field is constant in time and does not change horizontally, dz will always be zero, if the particles start at the same z-position, and the angle α will remain zero regardless of spatial changes in the horizontal velocity.In other words, there must be a horizontal gradient in the vertical velocity in order to induce a slope in the isochrones.
A horizontal gradient in the vertical velocity may occur due to different processes.The full expression for the vertical velocity includes surface accumulation rate, basal melt rate, horizontal velocity, horizontal gradients in surface elevation, horizontal gradients in bed elevation, and so-called velocity shape functions that describe how the velocities vary with depth (cf.Waddington et al., 2007).All of these processes may induce a horizontal gradient in the vertical velocity, however, over the spatial scale of a typical radio-echogram (10 1 -10 2 km) some of the processes can be neglected.For example, in the interior of GrIS the surface accumulation rate varies smoothly over distances typically ≥ 10 2 km (e.g., Burgess et al., 2010), and thus have a small horizontal gradient over the extent of a typical radio-echogram.The same is true for the surface slope and the horizontal velocity.Our method of detecting increasing layer slopes with depth will therefore primarily map areas with changes in basal melt rates, velocity shape function, and bed topography.
We will begin by discussing the impact of changing basal melt rates.Assume a constant and uniform accumulation rate and ice thickness, and small horizontal changes in u (i.e.∂u/∂x ∂ w/∂ z).
At the surface the vertical velocities will be the same everywhere because w at the surface is controlled by the accumulation rate.At the bed the vertical velocity is controlled by the basal melt rate, thus if the basal melt rate varies horizontally along the direction of flow, the change in w with depth will differ along flow.This introduces a horizontal gradient in w and this gradient will increase with depth as the vertical velocity approaches the local basal melt rate.Thus, the difference in vertical velocity between two points will increase with depth, leading to an increase in dz and a corresponding increase in layer slope α with depth.This is illustrated in Fig. 1.A change in velocity shape function can be caused by different processes: 1) The difference in ice flow at an ice divide (longitudinal extension and vertical compression) compared to the flow at the flanks of the divide (dominated by shearing).This process may influence the internal layers through so-called "Raymond bumps" (Raymond, 1983) which have been observed at multiple sites in Antarctica (e.g., Conway et al., 1999;Drews et al., 2015) but never in Greenland.2) A transition from deformational ice flow to ice flowing by sliding over the bed (also termed plug flow or streaming flow).This process is typically linked to changes at the ice-bed interface when the temperature rises above the pressure melting point (e.g., Pattyn, 2010).Even if the process is not caused by changing basal melt rates, we would still expect the change in velocity shape profile to be accompanied by a change in basal conditions (e.g., change in till strength at the bed).This is the kind of near-bed process we wish to map. 3) Ice flowing over obstacles at the bed, a process that is associated with and probably dwarfed by changes in bed topography, which is discussed below.Finally, we note that the velocity shape function may also be altered by changes in the ice deformational properties such as such refrozen meltwater or folded Eemian ice.However, since both refreezing and folding indicate processes acting near the bed, we still wish to map their occurrence and therefore do not attempt to distinguish between deformation and a simple change in basal melt rate.
Distinguishing the signal caused by the bed topography from the signal of changing velocity shape function or melt rate is the biggest challenge.Depending on the spatial scale of the variations in bed topography, the isochrones will either override the topography or drape it (Gudmundsson, 2003;Hindmarsh et al., 2006).As discussed later, for this reason we only consider cases where the layers cannot be explained by the bed topography.
It is important to note that once the layer slope has been induced, it will be maintained if there are no further changes in w and the horizontal velocity stays constant.However, in practice this never occurs.In almost all cases the horizontal velocity increases in the flow direction, and there will be an increase in dx and thereby a decrease in layer slope.This is also shown in Fig. 1.This implies that the detection of UDRs becomes problematic in fast flow areas.We will return to this in the discussion.
In the paragraphs above we have for the sake of simplicity restricted the discussion to two dimensions.However, the theoretical framework is also applicable in three dimensions as long as the vertical ice flow velocity changes along the direction of the radioechogram.

Internal layer slope
In order to find the slope of the internal layers, we used the method described in Panton ( 2014), which we will briefly outline here: First the echograms were aligned to the picked ice surface and de-trended with a high-pass filter.Then the echograms were convolved with slanted two-dimensional Gaussian filters, for n different slant angles.Finally, the angle of the filter that yields the strongest response to the convolution was used as an estimate of the layer slope at each point in the echogram.To account for changes in acquisition techniques and radar setup, we normalized the layer slope by the physical distance between the vertical traces and the radar depth resolution.
We use different parameters and a different interpolation scheme compared to what was originally published in Panton (2014).We varied the filter angle θ in n = 50 steps, with a maximum filter slant θ max = 55 • and the spread of the Gaussian filter σ x , σ z = (16, 0.125).Whereas the original method used weighted splines to interpolate the slope field in areas with poor response (e.g.due to missing layers), we use a simple linear regression, weighed by the response to the convolution.As the bedrock has a particularly strong reflection, the bottom 20% of the convolution response is nulled, in order to avoid over-fitting to the bedrock slope.See the supplementary material for an illustration of interpolating the slope field derived from RES echograms.
The choice of a linear fit was chosen over the second-order polynomial used by MacGregor et al. (2015) for three reasons: 1) For the case shown in Fig. 1, localized melt introduce a gradient in vertical velocity.Theoretical vertical velocities from ice flow models such as the Nye and Dansgaard-Johnsen solutions is almost linear with depth (Dansgaard and Johnsen, 1969), until fairly close to the bed (where we do not interpolate the slope in both methods).Therefore it is reasonable that a small change in the vertical velocity, as captured by the internal layers, can be fitted by linear regression.2) In-situ vertical velocity profiles captured by Kingslake et al. (2014), using phase sensitive radar, show that the vertical velocity in most cases can be fitted by linear regression.And 3) for the areas where the vertical velocity function can not be estimated, we find that a linear fit is a good conservative estimate for extrapolating 20% or more beyond the data coverage.

Detection and segmentation
As previously discussed, we intend to find the locations in the ice sheet where the internal layer stratigraphy is disturbed, and where these disturbances are not caused by the surface or bed topography.
To detect this, we introduce two variables: The slope of the observed bed, S o , which can be found from the pre-picked bed in the echograms, and the reconstructed bed slope, S r , which is the value of the interpolated internal layer slope at the depth of the bedrock.The method and variables are illustrated in Fig. 2. Note that the sign of S o and S r depends on both the slope direction (i.e. up or down) and on flight direction, however, this is corrected for in the subsequent mapping.
Since ice flow usually smooths the bedrock signature in the internal layers, we would expect the reconstructed bed slope S r to be smaller than the observed bed slope: S r ≤ S o (for upwards sloping layers, i.e. for S r > 0), if there is no disturbance in the layers.
On the other hand, if the reconstructed bed slope is larger than the observed bed slope, S r > S o (for S r > 0), we see this as a sign of inconsistency between the slope of the internal layers and the bedrock topography.In other words, the slope of the internal layers near the bed must be greater than the slope of the bed to detect the UDR.From this, we can set up the following boolean criteria to classify whether an area can be considered a UDR, We then apply the classification to all the flightlines (Figs. 4  and 6 are examples of radio-echograms, that contains locations where the Boolean criteria holds true).These locations have a spatial extent in the along track horizontal direction.In Figs. 4 and 6 they are marked with light red boxes.In the following, we will refer to each box as a segment, i.e., the term segment refers to a subset of a radio-echogram where a UDR has been detected.
We also consider along-track length L of the segment and mean slope difference, M, between S o and S r ,

Sr
if So < 0 and Sr > So Sr if So > 0 and Sr < So Sr − So elsewhere (2) where Sr and So is the mean value of S r and S o along the seg- ment.The two first cases handle the few instances where the bed might slope downwards, but the internal layers would slope upwards.In those cases we calculate M as if the observed bed was flat, S o = 0, in order not to get a misleadingly large M due to the bed topography.Finally, for each detection, we calculate the rise R = L • M, which is mean "area" of the layer inconsistency along the flight path.It is the value of the rise, R, that we use to quantify the size of the UDR.The reason for using R rather than the length of the segment, L, or the mean slope difference, M, is that we wish to avoid emphasizing horizontally elongated discrepancies or very narrow, vertically extensive features.A high value of R indicates both a large horizontal and vertical extent of a UDR.
While the classification was computed for all flightlines, some detections were discarded due to lack of layers in the segment, or if the segment is small L < 25 m, which corresponds to less than three traces in the echogram.Finally, we set a threshold value for R such that small and large values are discarded: R / ∈ [30 m, 2 km].A rise above the latter limit (2 km) is considered to be an outlier since that is approaching or exceeding the typical thickness of the ice column.

Results
All the processed flightlines can be seen in Fig. 3, along with the locations of the UDRs.Northern Greenland has by far the highest density of detected UDRs.The largest features can be seen at the onset of the Petermann glacier and west of the NEEM drill site on the ice divide, parallel to the Northeast Greenland Ice Stream (NEGIS).Like MacGregor et al. (2015), we consistently find the UDRs in repeat, cross-over, parallel and neighboring flightlines.In addition, we find a similar repeatability in the directionality of the slopes of the UDRs.The UDRs span a significant spatial area, typically on the range of 10 2 km 2 to 10 4 km 2 .Furthermore, the signal also returns identical magnitude and orientation showing that the method is robust across different generations of radar systems.In contrast, very few of such features were identified in central Greenland, where the bulk of the signal is found close to the summit of the GrIS, and along the ice divide north of the summit.Compared to northern Greenland, very few detections are made in this part of the ice sheet.In southern Greenland the method indicates that features are present in the Jakobshavn basin and along the east coast.The apparent high density of UDRs near the Jacobshavn Isbrae in western Greenland is partly caused by the high density of flightlines in the region, and the some of the detections in south-east Greenland are due to decreased layer quality.Our results are similar to, but have a higher along-track resolution (individual traces) than MacGregor et al. (2015) (fixed 5 km intervals).However, we do not recover as high a density in near Jacobshavn Isbrae and south-east Greenland as they do, which we discuss later.
We found the imprint of UDRs in 2.3% of the distance covered by echograms in Greenland.This metric increased to 4.3% if we only considered northern Greenland above 75 • north, and 6.4% of the flightline distance, if we confine the survey to echograms acquired over at least 1000 m of ice thickness over northern Greenland.
An example of a mapped feature can be seen in Fig. 4. It is located on the ice divide, 51 km southeast of the NEEM ice core drill site.Comparing the echogram (panel A) with the identified segments of UDRs (panel B) reveals an excellent agreement between the two.Importantly, the spatial extent of the feature is clearly seen across multiple flightlines (panel C), and there is very good agreement between the location of the identified slopes and their orientation and signal magnitude.
Another example is shown in Fig. 6.This shows the onset of Petermann Glacier where numerous UDRs have been identified.Again, visual inspection of the echogram (panel A) shows excellent agreement with the automated detection method (panel B).The densely sampled RES grid reveals that the up-warping UDRs at Petermann form elongated features oriented along the direction of ice flow.The direction of the features are indicated with pink lines.Note that the UDRs show up in the main flow towards Petermann, but they disappear once the surface velocity reaches 100 m/a.
We identify UDRs in both slow flowing and fast flowing areas (Fig. 5).In Northeast Greenland the features appear outside of the main NEGIS flow, and almost vanish once the surface velocity exceeds 100 m/a.Overall, however, we find no statistical correlation between surface velocities and R-values (ρ V surf ,R = −0.03).6.9% of all detections are found in areas where the velocity is less than 10 m/a, and 30% of detections are in areas where the velocity is less than 40 m/a.There is an indication that fast flow obscures or obliterates some of the features, since less than 25% of detections are in areas of fast flow, in this context defined as above 100 m/a.However, this could be due to the fact that the quality of the echograms is also decreasing due to the fast flow.Note that it is necessary to normalize the percentages with respect to the On the ice divide and in the area between the ice divide and NEGIS, the features are primarily elongated perpendicularly to the direction of ice flow.As the surface velocity increases, the orientation of the features change and it becomes aligned with the flow direction.This alignment is clearly evident east of the ice divide (marked with a black arrows in Fig. 5), where the orientation of the features (marked with pink lines in Fig. 5) is changing down flow.Finally, we note that both at Petermann and NEGIS our method clearly recovers the structures reported by Bell et al. (2014), which was based on manual picks, but we also report numerous detections not mentioned in that study.

Discussion
As shown in Fig. 3, we find occurrences of UDRs across the entire northern GrIS and along the ice divide towards the summit.The UDRs that dominate our results are large folds or plumes that push the stratigraphy upwards.We find that the largest continuous areas of UDRs are located near the fast flowing ice streams, Petermann and NEGIS, however, when we correct for sampling density, we find that the presence of UDRs is not correlated to the surface velocities less than 100 m/a.
Based on our results we assert that nearly all of the UDRs seen in the northern and central part of GrIS are locally confined in the form of up-warping or folds.We also observe that the UDRs fall into two groups; In areas with surface velocities less than 10 m/a, the UDRs take the form of bumps, orthogonal to the ice flow.In areas with higher surface velocities, the orientation of the UDRs become aligned with the flow.A clear example of the latter can be seen in the Petermann ice stream (Fig. 6).
Infrequently, we find areas with down-warping stratigraphy, indicative a locally increased basal melt rate.In Fig. 4 we see such an example, where the deep stratigraphy unexpectedly warps downwards, then after a few kilometers downstream a bump is formed.While we do not recover the large scale melt signatures on the scale of 10 2 km, such as those found by Fahnestock et al. (2001), we do capture the rapid changes within the melt area on the scale of a few kilometers.This is because the slope change for these signatures are insignificant compared to the roughness of the bed and such analysis becomes impossible over long distances.Traced internal layers are a better tool for such an analysis.
While our method does not directly reveal the origin or mechanism behind the UDRs, their location and orientation might pro-vide some insights into their development.We hypothesize that there is a transport mechanism in place, that moves and enhances the folded layers downstream through horizontal advection, and that the areas of slow flow could be inception points.An observation that supports this hypothesis is marked with a green arrow in Fig. 5.Here the radar coverage completely maps the extent of a feature, and the layer slopes can be seen to increase in the echograms when going from flightline to flightline, and then suddenly disappear.Thus, the layer slopes do not steadily decrease again, but instead stops abruptly.This indicates that the feature could be growing from an inception point roughly in the direction of ice flow.
A recent modeling study by Wolovick et al. (2014) demonstrated that large overturned structures can be formed in response to changing basal friction.However, in the north-central part of the GrIS (covered by Fig. 5), where we find the UDRs transverse to the ice flow, there is little indication of substantial local variation in melt rates, and we do not find a consistent signature of small hot spots upstream from the UDRs.Based on visual inspection of the echograms, we observe that ice tends to fold at the rheological contrast known to exist between glacial and Eemian ice (e.g., Dahl-Jensen et al., 2013).Based on this, we do not expect refrozen meltwater to be the primary source of their formation, but rather caused by reduced bed friction, further enhanced by overturning at the rheological boundaries.
In contrast, the areas where the along-flow UDRs occur could be at the pressure melting point, but there is a significant diversity in such prediction (e.g., Greve, 2005;Aschwanden et al., 2012;Seroussi et al., 2013), We found that the surface velocity over all the along-flow UDRs exceeds 10 m/a.At Petermann these features take up half the ice thickness in some locations, forming UDRs along the flow.As seen in Fig. 6, a trough (marked with a blue line) is formed between the UDRs, where the internal layers are located significantly deeper in the ice than on the other side of the UDRs, which is indicative of basal melt in the trough upstream from the echogram.The direction of the UDRs and their proximity to a probable meltwater source supports the hypothesis of Bell et al. (2014), where the formation of basal units is explained by the refreezing of meltwater.If this is the case, then the latent heat is absorbed by the surrounding glacial ice, which must lead to an overall warming of the base of the ice sheet.This could be the reason why field studies have found evidence of temperate basal layers in ice streams (Lüthi et al., 2002).
We find most of the UDRs in northern Greenland, east of the ice divide, which is the same overall pattern that was found by a reflectivity-based basal melt assessment by Oswald and Gogineni (2012), where they place most of the basal melt in the same area.However, the local agreement between detection of basal melt and our detection of UDRs is poor as we only detect UDRs based on rapid slope changes, which might not be present in areas with uniform melt rates.
It is important to note that the determination of the extent and origin of UDRs requires a high density of flightlines, which has only recently been made available in Greenland during the NASA Operation IceBridge campaigns.We observe UDRs disappear between two flightlines in the flow direction, where the line spacing is on the order of 10 km.This emphasizes the need for high density sampling, in order to understand the origin, extent and evolution of UDRs -and by extension the basal conditions of the GrIS.
The results presented in this study have been generated using a completely automated method, with several differences and improvements to what has been presented by MacGregor et al. (2015): Our method does not rely on the internal layers to be traced, but rather uses the layer slope, and our analysis is made using a higher resolution data product that respects the slope orientation of the UDRs.The method is transferable across RES datasets and the time needed to obtain information on UDRs is substantially reduced, compared to manual inspection of RES data.Finally, we would like to highlight that the mapping from our analysis generally matches the map produced by MacGregor et al. (2015), manual inspection by Bell et al. (2014) and our own manual assessment (e.g.Fig. 4 and Fig. 6).
A limiting factor for our method is the quality of the imaged layers.The quality of internal layers in some areas can be very poor.For the GrIS, the south-western part has few layers (MacGregor et al., 2015), probably due to the large amount of surface melt each summer (Abdalati and Steffen, 2001).Also in the glaciers and ice streams, the layer quality decreases due to a combination of the fast flow, increasing ice temperature and increasing surface clutter.This may also explain drop in detections of UDRs when the surface velocity exceeds 100 m/a.
The scattered detections of UDRs near Jacobshavn Isbrae in western Greenland, does not have a good enough correlation to nearby flightlines to show the physical extent of the features despite high density of flightlines in the area.Manual inspection of echograms from this region revealed that in areas where layer quality was good, UDRs were correctly identified, but they were difficult to constrain in the echograms, due to the reduced layer quality by steeply sloping layers, rapid shifts in ice thickness and artifacts in the echograms (see supplementary material for an example).As some artifacts from the acquisition process resemble layers, this have led to some cases of misclassification by the automated algorithm.
MacGregor et al. ( 2015) detects more UDRs near the Jacobshavn Isbrae and to the southeast, than we do.This is most likely due to differences in detection methodology.We produced similar results prior to extending the mean slope difference, M, to include the special cases where S r and S o had opposite signs and diverged.The bedrock under the southeastern ice sheet is particularly rough (Rippin, 2013), which can explain the increased number and size of UDR detections in that area, which we consider to be overemphasized.

Conclusion
We present an automated method for identifying features of disrupted radiostratigraphy of the Greenland Ice Sheet.The method relies on detecting the signature that the features leave in the internal layer slope, and shows good agreement with results produced by a similar method based on traced internal layers, and another based on the manual inspection of radio-echograms.As the layer slopes are automatically computed without tracing individual reflectors, the method can be readily applied to extensive RES datasets.
Using this method, we find that large parts of the RES data from the northern GrIS exhibit locally confined up-warping of the internal layer stratigraphy, with the largest imprint near ice streams.We find no correlation between the size of the features, the frequency of occurrence and surface velocity.However, we do observe that some of the features appear to have a specific orientation, which appears to be correlated to ice flow.
Based on the spatial distribution of the features, we propose that folds and up-warping could occur regardless of the presence of basal melt, but that the presence of meltwater might change or enhance the large scale structure of the near bed features.While we cannot decisively determine the formation processes behind the present observations, it is likely that more than one formation mechanism must be in place.
Our results highlight the need to understand the mechanisms behind the formation of the disrupted radiostratigraphy, the basal environment of the GrIS and its influence on ice flow.

Fig. 1 .
Fig. 1.Modeled flow line from the ice divide (x = 0 km) with a "sticky spot" (area of low basal melt) located between x > 100 km and x < 175 km (in blue).A) Vertical velocity field (grayscale background) with different particle paths (arrows) and isochrones (black solid lines).The inset shows α the slope of the layers.B) The slope of the isochrones with depth at x = 120 km, red line indicates the linear interpolation of isochrone slope with depth.C) The mean slope change with depth (normalized), along the flow line.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Comparison of the reconstructed bedrock slope (S r ) and the slope of the observed picked bedrock (S o ).To the left, a schematic echogram is shown with surface, internal layers and bedrock.To the right, the corresponding layer slope with depth is shown, at the location of the red line in the echogram.For case A, the shape of the bedrock dominates the internal layers, whereas for case B, the change in internal layer slope cannot be explained by the shape of the bedrock.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Map of Greenland, with the flightlines used for this study in gray.Blue dots indicate the location and size (R) of the UDRs.Yellow and orange colors indicate surface velocities from Joughin et al. (2010), at the 40 m/a and 100 m/a contours respectively.Red boxes correspond to areas presented in detail in Figs. 4, 5 and 6. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. UDRs seen on ice divide, 51 km southeast (upstream) of the NEEM ice core drill site.A: Radio-echogram of the UDRs (frames 20110506_01_011 and 012), and B: the corresponding detections generated by our algorithm in red boxes.The red boxes are the detections and their area is the rise R. The black line represents the value of S o , the red the value of S r .C: Map of the local area, with flightlines in gray and colored dots for the individual detections.The size of the dots represents the size of the UDR (the R-value) and the color represents the directivity of the UDR slope as detected in the echogram, mimicking a shaded relief (corresponding to the color rose in the lower left, where the center of the rose represent a hill).The multiple southeast-northwest flightlines were flown close to and parallel with the ice divide.In the right panel, the black arrow is the location of the echogram and the orange arrows are the direction of ice flow based on surface velocity.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Map of northern Greenland showing the location of UDRs, with the highest density at the onset of the Petermann glacier to the northwest, and parallel to Northeast Greenland Ice Stream to the east.Smaller areas of UDRs can be seen over most of the interior, across multiple flightlines.The pink lines are our manual interpretation of where the peak of up-warping UDRs are located.Several of these are marked with arrows, where the green arrows are showing large areas of correlated UDRs west of the ice divide and black to the east.The size of the dots represents the size of the UDR (the R-value) and the color represents the directivity of the UDR slope as detected in the echogram, mimicking a shaded relief (corresponding to the color rose in the lower right, where the center of the rose represent a hill).Surface velocities shown as from yellow to orange at 10, 20, 40 and 100 m/a intervals.Orange arrows show the ice flow direction, based on surface velocities.Flightlines in gray.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.The onset of the Petermann Glacier, showing the manual interpretation of where the ridge-like UDRs are located (pink lines on the map, white arrows in the echogram).The figure layout is similar to Fig. 4. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)