Surface Expression of Basal and Englacial Features, Properties, and Processes of the Greenland Ice Sheet

Radar‐sounding surveys measuring ice thickness in Greenland have enabled an increasingly “complete” knowledge of basal topography and glaciological processes. Where such observations are spatially limited, bed elevation has been interpolated through mass conservation or kriging. Ordinary kriging fails to resolve anisotropy in bed geometry, however, leaving complex topography misrepresented in elevation models of the ice sheet bed. Here, we demonstrate the potential of new high‐resolution (≤5 m) surface topography data (ArcticDEM) to provide enhanced insight into basal and englacial geometry and processes. Notable surface features, quantified via residual surface elevation, are observed coincident with documented subglacial channels, and new, smaller‐scale tributaries (<2,000 m in width) and valley‐like structures are clearly identified. Residual surface elevation also allows the extent of basal ice units to be mapped, which in conjunction with radar data indicate that they act as “false bottoms,” likely due to a rheological contrast in the ice column.


Introduction
Advances in understanding basal processes of the Greenland Ice Sheet (GrIS) have primarily been driven by the need to refine projections of the ice sheet's response to climate change using numerical models. An increase in observations of ice thickness, obtained using airborne ice-penetrating radar (IPR), has enabled improved representation of Greenland's bed topography and, consequently, the production of progressively accurate and spatially "complete" (without holes) bed digital elevation models (DEMs; e.g., Bamber et al., 2001;Morlighem et al., 2014Morlighem et al., , 2017. IPR data, therefore, provide a necessary input to ice sheet modeling and also allow for the assessment of more specific basal and englacial properties, including subglacial roughness and its anisotropy (e.g., Jordan et al., 2017;Lindbäck & Pettersson, 2015;Rippin, 2013); ice rheology, through the measurement of englacial stratigraphy (e.g., Karlsson et al., 2013;; englacial temperature from radar attenuation (e.g., MacGregor, Li, et al., 2015); and determining the distribution of basal water from bed-echo reflectivity (e.g., Jordan et al., 2018;Oswald et al., 2018). Despite the influence of basal processes upon ice dynamics being theoretically well constrained (Cuffey & Paterson, 2010;van der Veen, 2013), our understanding of spatial variation in subsurface conditions and processes remains restricted by the paucity of IPR data; consequently, observation-led interrogation is required to understand the magnitude of their effects in situ.
While the number of IPR ice thickness observations has roughly tripled in the past decade, there remain large swathes of the GrIS left unobserved (see Figure 1a in Morlighem et al., 2017). Limitations in IPR coverage, and spatial bias due to variation in observation density, often result in the mapping and assessment of basal geometry and processes being reliant on spatial interpolation between flight tracks. Interpolation has primarily been applied to the construction of bed topographies Morlighem et al., 2017), but also the tracing of englacial stratigraphy (Bons et al., 2016;, and the mapping of subglacial roughness (Rippin, 2013). Interpolated bed DEMs have enabled the "discovery" of large-scale landscape systems beneath the GrIS (e.g., Bamber, Siegert, et al., 2013;Cooper et al., 2016;Livingstone et al., 2017) confirming that the identification of more nuanced interrelationships between basal topography and the overlying ice is now possible. In this vein, process-based assessments have attempted to understand the relationship between ice flow and bed topography, with focus on aspects of landscape genesis, evolution, preservation, and inheritance in Antarctica (e.g., Bingham et al., 2017;Rose et al., 2013Rose et al., , 2015Ross et al., 2014) and Greenland (Cooper et al., 2016;Livingstone et al., 2017).
The most recent data sets for Greenland bed topography have employed a "mass conservation" interpolation between flight lines in fast-flowing regions (Morlighem et al., 2014(Morlighem et al., , 2017. However, the majority of the bed in the latest product (BedMachine v3, referred to here as BM3), primarily the slow-flowing ice sheet interior, was interpolated using ordinary (isotropic) kriging (Morlighem et al., 2017). In observation-sparse areas, this approach generates surfaces that fail to represent the "true" geometry of anisotropic basal features (i.e., quasi-linear geometries): Rendered surfaces can appear artificially smooth  or result in "bulls-eye" anomalies manifesting as a series of isolated peaks and/or hollows (Dentith & Mudge, 2014). Although BM3 is provided at a posting of 150 m, the true resolution of the estimated surface is dependent upon the local observation density (i.e., the level of clustering) and sampling regime (i.e., the flight line spacing; Morlighem et al., 2017). Since parallel track spacing over the ice sheet interior is up to 75 km, it is important, therefore, to consider the need for a spatially consistent method through which to assess basal morphology and associated basal and englacial processes without the artifacts associated with interpolation.
The transfer of basal information from "bed-to-surface" has long been appreciated, with knowledge of local ice velocity and flow parameters being used to calculate bed topography along a profile of the surface (Budd, 1970). More recently, surface and subsurface data have been used to better constrain this transfer of basal information, with particular reference to spatial variability in basal traction (De Rydt et al., 2013;Gudmundsson, 2003;Gudmundsson & Raymond, 2008). However, it is widely understood that local variations in flow regime and surface topography are more suitably explained by undulations in subglacial topography (De Rydt et al., 2013). It is yet to be determined whether this modeled basal transfer (as in Gudmundsson, 2003) can be applied generally, especially where ice thickness and speed vary considerably. This is perhaps exemplified through the application of mass conservation interpolation being limited to regions of fast flow (Morlighem et al., 2014(Morlighem et al., , 2017. While recent efforts have been made to characterize basal information transfer across the interior of the GrIS, where velocity is low and basal slip is spatially heterogenous, the success of this is limited to more large-scale patterns, arguably by computational means Ng et al., 2018). Outside theoretical frameworks, and particularly in the interior of Antarctica, numerous observational studies have demonstrated the relationship between surface and basal geometry (Le Brocq et al., 2008;Rémy & Minster, 1997;Rose et al., 2014;Ross et al., 2014). For example, subglacial lakes manifest as flat surfaces in the overlying ice sheet, due to the negligible basal friction above them (Ridley et al., 1993;Siegert & Ridley, 1998). Such research has benefited from ice sheet surface brightness or topography, including analyses of surface curvature, to determine bed features and subglacial landscape (Le Brocq et al., 2008;Rose et al., 2014;Ross et al., 2014).
While both ice surface observation-and model-led approaches reveal important information about subsurface processes in areas of sparse observational data, much of the small-scale topography of the Greenland interior remains uncharted; this is primarily due to the relative coarseness (spatial resolution) of previously available surface imagery (e.g., Moderate Resolution Imaging Spectroradiometer surface brightness, as used in Ross et al., 2014, and others), and the complexities in modeling transfer of information from bed to surface (De Rydt et al., 2013). However, with the recent release of a new high-resolution (≤5 m) pan-Arctic DEM (ArcticDEM; Morin et al., 2016; https://www.pgc.umn.edu/data/arcticdem/), there is an opportunity to improve our knowledge of spatially heterogenous basal conditions using the ice surface topography rather than, or in conjunction with, IPR data.
Here, we explore the potential for ArcticDEM ice surface information-specifically the residual surface elevation (RSE)-to provide a "lens" for subsurface processes in Greenland. From this, we present evidence for previously undocumented complex basal topographies and characterize the influence of ice dynamics on the ice surface. We also provide evidence that some surface expressions are likely a result of englacially discrete, rheological boundaries and improve the mapping of "units of disrupted radiostratigraphy" (UDRs; Bell et al., 2014;Dow et al., 2018;Leysinger Vieli et al., 2018;Wolovick et al., 2014). We compare results inferred from surface data with BM3 and highlight the importance of topographic anisotropy and scale both in the planning of geophysical surveys and spatial interpolation. Finally, to understand the relationship between the surface expression and basal processes we compare with complementary data sets: basal water , geothermal heat flux , basal slip (MacGregor et al., 2016), and ice velocity (Joughin et al., 2016(Joughin et al., , 2017.

Calculating RSE
The ArcticDEM is a newly released collection of Earth surface DEMs produced using high-resolution (~0.5 m) stereopair imagery, collected by DigitalGlobe Inc. optical satellites. The DEM is supplied by the Polar Geospatial Center in two formats: the "ArcticDEM Strip" format (2-m resolution) and tiled "Mosaic" format (5-m resolution). The Mosaic format was used here owing to both its spatial regularity (50 km × 50 km; 2,500-km 2 tiles), increased data coverage, improved absolute accuracy, and fewer artifacts. A total of 36 tiles, covering an area of~90,000 km 2 over northwest Greenland, south of Petermann Glacier, defines the study region (Figures 1 and S1 in the supporting information). This region was chosen due to its relatively high, and spatially dense, coverage of IPR data ( Figure 1a); additionally, studies mapping and understanding of both englacial (e.g., Bell et al., 2014;Chu et al., 2018) and basal features/properties (e.g., Livingstone et al., 2017) have also focused on this region. This facilitates both the comparison and validation of this paper's conclusions, in order to better assess the capability of high-resolution surface information in extracting subsurface information. Ice sheet surface topography is influenced by underlying basal conditions (i.e., geometry) and processes (e.g., Gudmundsson, 2003;Rémy & Minster, 1997). Compared to the dominant signal of an ice sheet's surface slope (which is approximately parabolic, or dome like), however, this influence yields a much smaller change at the ice surface, both in terms of amplitude and wavelength. To help resolve the influence of basal geometry, large-scale surface slopes must be removed. We calculated a third-order polynomial trend surface for each ArcticDEM tile, using least squares regression, to detrend surface topography, leaving the RSE (Figure 1b). The RSE preserves local surface features down to a minimum wavelength of about 1 km, which are then interpreted as manifestations of basal topography or englacial surfaces (section 3). While the ArcticDEM boasts a much finer scale than this, it is probable that smaller-wavelength features are not resolvable due to the dampened transfer of basal information as a result of ice thickness (De Rydt et al., 2013); for such features "direct" observation using IPR data are likely necessary. Additional second-and fourth-order polynomial surfaces were used; however, no marked improvement, in terms of interpretable "features," was observed. By contrast, the third-order surface performed better over ice divides and at ice margins when compared to fitting a second-order polynomial.
Various feature extraction methods have been employed to obtain bed geometry from ice surface topographies, and surface brightness. This includes edge detectors (e.g., Ekholm et al., 1998) and isolating local topographic change through the calculation of profile curvature (of Moderate Resolution Imaging Spectroradiometer imagery, e.g., Ross et al., 2014). While these methods are usually dimensionless (without vertical information) and only reveal relative changes at the bed, spatial patterns in basal geometry can still be interpreted. We repeated these approaches for the new ArcticDEM ( Figures S2a-S2c), alongside a calculation of RSE for a more coarse surface topography (90-m resolution) DEM ( Figure S2d; Greenland Ice Mapping Project, GIMP; Howat et al., 2014). Although a 30-m GIMP DEM is available, the latter comparison was undertaken to explicitly assess whether topography at a greater spatial resolution (i.e., ArcticDEM) can provide enhanced subsurface information. RSE calculated for the 30-m GIMP DEM provides a near-similar level of information provided by the ArcticDEM; however, for the purpose of this paper we rely on the new, and higher-resolution ArcticDEM for further analysis.

Interpreting Surface Expressions
RSE ( Figure 1b) can be interpreted similar to a relative DEM, where surface undulations are attributed to bed undulations or, as we later demonstrate, basal-englacial UDRs. Various geometric features, including lineations, troughs, and peaks, are visible in the RSE (Movie S1). Interpretation of these features is facilitated by the comparison with BM3, (Morlighem et al., 2017) and along-track radargrams from National Aeronautics and Space Administration's Operation IceBridge (Rodriguez-Morales et al., 2014; Figure 1a), as well as with complementary data sets (Figure 2), including interferometric synthetic aperture radar-derived surface ice velocity (Joughin et al., 2016), radar-derived basal water , and magnetic-derived geothermal heat flux . On initial comparison to BM3, the RSE (Figure 1b) replicates fundamental landscape features. Most notably these are large, and previously documented, dendritic channel networks, such as the "mega-canyon" flowing toward Petermann Glacier (Bamber, Siegert, et al., 2013), and the two channels linked with Humbolt Glacier (Livingstone et al., 2017; MC and LC, respectively; Figure 1a). Further, RSE broadly simulates the overall trends in basal topography (acting similarly to a topographic "hillshade," Figure 2a), and in observationsparse regions it delineates previously undocumented morphological features (see sections 3.1 and 3.2 and Figure 1b). This confirms that the ice surface can be used to explore basal geometries, which are poorly constrained by the shortcomings of isotropic interpolation between sparse/transectbased observations.
There are, however, notable signals and features present in the RSE that do not appear to reflect likely basal topography; instead, these features result from ice dynamics and/or UDRs. In the first instance, a marked "rippling" is observed in the RSE coincident with Petermann Glacier, extending southeast ( Figure 1b). While this signal appears linked to a spatial transition in the ice speed, there is no direct relationship to magnitude (Figures 2a and S3). Evidence of streaming ice is identified by the tight spacing of velocity contours, delineating pronounced shear margins ( Figure 2a) and an elevated basal slip ratio (Figure 2c; MacGregor et al., 2016). It is at shear margins, where strain rates and shear stress are high, and within fast-flowing ice, that ice is likely to become decoupled from the bed (Cuffey & Paterson, 2010). Spatial similarity between the observed RSE rippling features, shear margins as evidenced by velocity contours, and modeled driving and basal shear stress is observed (see Figure 3 in Sergienko et al., 2014). It is, therefore, plausible that the coincident response in ice surface topography may result from changes in ice flow parameters (i.e., decoupling) and a more efficient transfer of large-scale basal undulations, resulting from increased slip (De Rydt et al., 2013;Gudmundsson, 2003).
On the other hand, south of Petermann Glacier, in a relatively flat and low-lying area of the bed (Figure 1a), multiple ovular features are observed in the RSE (Figure 1b). UDRs have been documented across this sector of the GrIS in radiostratigraphy Dow et al., 2018;Leysinger Vieli et al., 2018;Wolovick et al., 2014). The established locations of most of these units correspond to the ovular/lobed features observed in the RSE (Figure 2b). The genesis and interrelationship with local ice flow parameters of these UDRs are discussed below.

"Ground Truthing" RSE
To disentangle the primary causes of the surface features (Figure 1b), comparison can be made against radar-derived bed elevation and stratigraphy. Three illustrative profiles are shown in Figure 3. Profile A provides a near-perpendicular cross section of the MC (Bamber, Siegert, et al., 2013), where profiles B and C present near-orthogonal bisectors of a UDR. Although a marked response in RSE is apparent along all three transects, the nature of these differ (Figures 3a-3c). RSE along profile A follows basal topography with a corresponding negative, smoothed, facsimile of bed elevation; however, this response is notably out of phase. This is likely a function of ice flow direction (~65°from profile; Figure S3), as well as ice rheology, thickness, and thermal regime, all of which affect the efficiency of bed-to-surface transmission (De Rydt et al., 2013;Gudmundsson, 2003). This informs us that although additional influences impinge the bed-to-surface transfer, we can still extract information regarding the scale and geometry of the basal feature using RSE. While a large positive response in RSE is observed along profiles B and C, little variation is observed in sampled bed elevation (standard deviation ≤19 m); it is unlikely, therefore, that basal topography is the primary cause.  Figure S4). Englacial layers that correlate with "surface ice" are shown (white lines), where the thickest line represents the first complete layer above the unit of disrupted radiostratigraphy; dashed white line denotes the separation between Holocene ice (above) and Last Glacial Period ice (below). Radiostratigraphy within reflective basal ice are shown in purple, and near-bed reflectors beneath nonreflective ice are blue. The asterisk in (f) is interpreted as an inner ellipse within a sheath fold. Bed echo (red) and ice surface returns are labeled. Red arrows x axes (e and f) mark the location at which (f) orthogonally (~90.1°) bisects the plane of (e). There is a vertical exaggeration by a factor of 19.35 and 11.45, for (a)-(c) and (d)-(f), respectively. RSE = residual surface elevation.
When considering returned power across the MC (profile A; Figure 3d), englacial stratigraphy generally conforms with basal topography. Specifically, some layer disturbance is partially visible toward the base of the ice column, but at this depth the ice has few radio-reflecting horizons. The dashed white line on all profiles indicates the transition in englacial stratigraphy between the Holocene ice above and the "softer" ice of the Last Glacial Period, below (Chu et al., 2018;Dahl-Jensen & Gundestrup, 1987;Huybrechts, 1994). Above the canyon (Figure 3d) this better documents the "dampening" effect in the transfer of basal information in englacial stratigraphy toward the surface (De Rydt et al., 2013), presenting a smoother transfer of the bed geometry in the RSE (Figure 3a).
Radargrams from profiles B and C (Figures 3e, 3f, and S4) depict a significant UDR (Figure 3a, Bell et al., 2014 and Figure 1, Wolovick et al., 2014), with a thickness up to one half of the ice column (~1,000 m). When comparing the observed structure and reflectance characteristics of the UDR (Figures 3e, 3f, and S5) against similar features described in Greenland  and Antarctica (Wrona et al., 2017), we see several similarities. Large-scale cylindrical folding of this unit is also present, with the axis aligned perpendicular to flow (Figure 3e); more notably, a lobular extension from the top of the basal ice, elongated in the direction of surface ice flow, is similar to the "finger structures" observed in Antarctica formed through the entrainment of the softer basal unit into the overlying ice (Wrona et al., 2017). Across flow a central "eyelet" is observed within the UDR (asterisk, Figure 3f), indicative of extension/deformation of the unit parallel to the flow direction, characteristic of sheath folding Wrona et al., 2017). The radiostratigraphy broadly shows the overlying ice flow "up-and-over" the UDR (Figures 3e and 3f) mostly undisturbed; this reflects the likely cause of the positive rise in RSE (Figure 3b), whereby the basal unit is hypothesized to act as a "false topographic bed," over which the upper ice layers conform. It is, therefore, possible that by using RSE we are able to better map or extrapolate the 2-D areal/horizontal extent of UDRs where directly observed in radargrams; however, in observation-sparse regions (away from IPR flight lines), the cause of such ovular surface expressions is likely to be more ambiguous, as internal surfaces and certain bed geometry (e.g., basins and topographic rises) could present similarly, as such "additional" UDRs (those not previously documented) cannot be confidently mapped using RSE. To better answer this would require the implementation of a bed-to-surface transfer framework (as in Gudmundsson, 2003;De Rydt et al., 2013), which is beyond the scope of this paper.
Specifically, at the point of bisection between profiles B and C, we can use radar power anisotropy to make a preliminary assessment of vertical changes in the horizontal dielectric anisotropy of the ice column and related asymmetry of the ice crystal orientation fabric (COF; Fujita et al., 2006;Hargreaves, 1978). Radar power anisotropy is influenced by both birefringent propagation (associated with horizontal COF anisotropy that varies smoothly with ice depth) and anisotropic scattering (associated with sharp vertical transitions in the COF; e.g., Fujita et al., 2006). Specifically, at middepths (~900-1,700 m), we observe oscillations in the across/along-flow power difference as a function of depth ( Figure S5b). The~10-dB amplitude and smoothly varying depth periodicity are seen to be consistent with power modulation dominated by birefringent radio wave propagation through the ice column (Fujita et al., 2006;Matsuoka et al., 2012). This implies smoothly varying horizontal asymmetry of the COF (i.e., the presence of either a vertical girdle fabric distribution or elongated single maximum). In both more shallow and deeper ice, the power difference shows greater isotropy, implying a weaker horizontal asymmetry of the COF. This is likely representative of a near-random fabric in shallow ice and a strengthened single maximum in deeper ice. As previously noted, the Holocene-Last Glacial Period transition is marked by a change in returned power (Chu et al., 2018) and is observed at~900 m (Figures 3d-3f and S5a), approximately coincident with the point at which we infer the development of horizontal asymmetry in the COF.
Multiple hypotheses surrounding the inception of UDRs exist, including bottom-up accretion through supercooling and basal freeze-on Leysinger Vieli et al., 2018), thermomechanical changes resulting from anisotropic ice rheology, or changes in rheology affected by "stick-slip" mechanisms (Bons et al., 2016;Weertman, 1976 ;Wolovick et al., 2014). The latter of these mechanisms, initiated either by velocity changes within the ice column or in basal thermal state, would give rise to large-scale folding and deformation like that seen here (Weertman, 1976;Wolovick et al., 2014). While it has been noted that some UDRs in Greenland do not have a clear relationship to known basal water networks , the distribution of observed features is consistent with radar basal water predictions (Chu et al., 2018;Jordan et al., 2018). Additionally, these are within a region of elevated geothermal heat flux (>58 mW/m 2 ; Martos et al., 2018;Figure 2b), and transitional basal thermal state (Chu et al., 2018;MacGregor et al., 2016). While hydrological modeling suggests minimal basal freeze-on rates (Dow et al., 2018;MacGregor et al., 2016), the correlation with these complementary geophysical properties, as well as evidence of multilobed sheath folding, may be symptomatic of heterogenous basal sliding or altered ice rheology, ultimately leading to the formation of UDRs. This is compounded not only by an enhanced predicted basal motion over a particularly large UDR-feature observed in the RSE (Figure 2c; asterisk, 4; MacGregor et al., 2016) but also through the inferred depth transitions from a near horizontally isotropic to horizontally anisotropic COF ( Figure  S5). Altogether, this suggests that these UDR features may act, dynamically, as "false beds" with localized englacial "decoupling" across the upper boundary.

Mapping Basal Characteristics
Although the presence, and causes, of surface expressions observed in the RSE can be directly attributed to basal topographies and associated subglacial processes using IPR data, much of the true extent of these features cannot be realized due to sparse/absent observations. Figure 4 presents an annotated RSE, including our interpretation for the primary origin of surface expressions. We can map with greatest confidence quasi-linear topographic features, such as channel systems and dendritic fluvial landscape morphologies, though it is apparent that sites of significant basal UDRs, and velocity induced changes, plausibly indicative of glacier-bed decoupling, are also interpretable. We are able to discern prominent bed geometries in the RSE, notably the additional detail in the landscape surrounding the MC and the northeast of the study region. Multiple, sinuous tributaries exist, as well as complex channel network and dendritic valley-like structures, suggestive of headwater catchments with several "minor" channels visible (<2,000 m in width; e.g., HS, Figure 4). Further, we see in tile E4 (Figure 4) that the MC appears to split into two upstream tributaries. It should be noted that where the RSE calculated for the more coarse resolution (90 m) GIMP DEM ( Figure S2d) provides prominent basal channels (i.e., the MC), "fine" scale detail (i.e., HS, Figure 4) is not interpretable from the more "smoothed" information. While many channels are aligned quasi-parallel to contemporary ice flow direction, there are some channels that are aligned perpendicular to flow and, in one example, span the present-day ice divide (DC, Figure 4). Complex landscapes such as these are indicative of past fluvial activity and are often interpreted as inherited or partially preserved ancient fluvial systems, predating the inception of contemporary ice sheet. More "major" channels (in terms of great width and relative depth), and associated channel networks, are seen to influence contemporary ice flow and the routing of contemporary subglacial water (Bamber, Siegert, et al., 2013;Bingham et al., 2017;Cooper et al., 2016;Jordan et al., 2018;Livingstone et al., 2017).

Conclusions
Here, we established the potential of the high-resolution (5 m) ArcticDEM to elucidate basal geometry with previously unobtainable detail. Where current geostatistically derived bed topographies rely on isotropic methods of interpolation (i.e., ordinary kriging), RSE provides a spatially complete proxy to constrain local (horizontal length scales~1-10 km) morphology. This approach is particularly advantageous for anisotropic features within the ice sheet interior. Additionally, via combination with complementary IPR data, we can constrain dynamical basal features (e.g., subglacial drainage pathways) and englacial processes (e.g., the 2-D areal/horizontal extent of UDRs). Additionally, via the analysis of two bisecting radargrams, we have discussed how one such UDR has formed, as well as its dynamic relationship with the ice column. This permits correlation with, and influences of, and upon, surface flow characteristics (e.g., sheath folding, fold elongation, and localized englacial decoupling).
In summary, we illustrate a complex subglacial and englacial system in northwest Greenland, notably improving the delineation of anisotropic geometric features that had previously been possible, or poorly captured, using isotropic interpolation. While qualitative, knowledge regarding the presence of complex basal topography, or indeed basal and englacial processes, may prove valuable not only in facilitating the further understanding of bed-to-surface transfer, but also in working to constrain future interpolation methods, and aiding the design of more directed geophysical campaigns. In the first instance, RSE is able to provide vertical information (in meter units), absent in previous edge detection and curvature methods . Where additionally, the positive identification of quasi-linear channels and dendritic valley networks can provide channel centerlines for physically based, anisotropic interpolation of bed topographies (exemplified in Goff et al., 2014;Herzfeld et al., 2011;Williams et al., 2017).