Interactive comment on “ Subglacial drainage patterns of Devon Island , Canada : Detailed comparison of river and tunnel valleys ”

I would like to thank Stephen Livingstone for this very detailed and constructive revision, which we believe will greatly enrich the manuscript. My coauthors and I agree on all the general revisions you proposed, and after considering the changes and literature sources you suggested, our plan is to modify the manuscript in the following way: (1) Change the term ‘tunnel valleys’ to ‘subglacial meltwater channels’ throughout the manuscript (2) Addition of literature relevant to the study


Introduction
Tunnel valleys, often referred to as tunnel channels, N-channels, and snake coils in the literature, are the erosional expression of turbulent flows in pressurized subglacial channels.Together with subglacial channels incised in the ice (R-channels), they modulate meltwater accumulation at the base of ice sheets and serve as highly efficient drainage pathways carrying meltwater to the ice terminus (e.g., Röthlisberger, 1972;Weertman, 1972;Nye, 1976;Kehew et al., 2012).Indeed, the transition from distributed to channelized drainage, and the characteristics of the resulting subglacial channel network leads to a reduction in ice flow rates, affecting ice loss and ice surging (e.g., Schoof, 2010).In spite of their importance, some outstanding questions remain: What are the typical lengthscales that characterize subglacial drainage systems (i.e., what is the drainage area, how many individual valleys form, how many tributaries do they have, etc.)?How do tunnel valleys differ geometrically from river valleys?How can we reliably identify tunnel valleys only with remote sensing techniques, including imagery and topographic data?We use a detailed geomorphological study of exposed tunnel valleys on Devon Island (Canadian Arctic Archipelago) to build understanding of these questions.Our study of exposed tunnel valleys represents a complementary approach to in situ and modelling studies of active subglacial drainage, allowing for a network wide characterization of subglacial drainage which adds context to point measures of discharge or effective pressure, and ground truths models of subglacial channel formation.
Exposures of tunnel valleys are rare.During glacial recession, meltwater released from the ice sheet accumulates at the ice marginal area and erodes the channel, with post-glacial sediment accumulation causing burial or partial burial (Le Heron et al., 2009).Vegetation overprint and fluvial incision makes the detailed study of channel geometry and morphology difficult (e.g., Walder and Hallet, 1979).Exceptions include areas with polar desert climate in the Antarctica Dry Valleys (e.g., Sugden et al., 1991) and the Canadian High Arctic (e.g., Dyke, 1999).The reduced rainfall conditions of these sites, recent ice retreat, and null or minimal vegetation cover are key for the preservation of tunnel valleys.The morphology of tunnel valleys at our field study area on Devon Island (Fig. 1), the second-largest of the Queen Elizabeth Islands in the Canadian Arctic Archipelago, is consequently well-preserved, and the retreat of the Devon Island ice cap offers a unique opportunity to compare recently exposed tunnel valleys with tunnel valleys incised during the Younger Drias Innuitian de-glaciation.
The main goal of this study is to characterize tunnel valley morphology using remote sensing data.Although tunnel valleys have been described in northern Europe (e.g., Piotrowski et al., 2006;Jørgensen and Sandersen, 2006), the Dry Valleys in The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-236Manuscript under review for journal The Cryosphere Discussion started: 23 November 2017 c Author(s) 2017.CC BY 4.0 License.self-explanatory.These guidelines are qualitative, and a topographic and imagery based approach to tunnel valley identification is yet to be developed.(e.g., Kehew et al., 2012).
To address this objective we conducted fieldwork on central Devon Island, and complemented it with airborne imagery data of central to eastern Devon Island acquired in a helicopter campaign.In section 2 we provide a detailed description of the field data acquisition and processing for river and tunnel valleys.We acquired GPS borne channel longitudinal profiles, stereo imagery, and derived photogrammetry digital elevation models.We also used, for the first time, a kinematic mobile LiDAR scanning (KLS) system, a novel method of measuring ultra-high resolution topography (<2cm/pixel Digital Elevation Models (DEMs)), which is described in more detail in Section 2.1.2.Section 3 includes the results from field observations regarding channel longitudinal profiles, the Digital Surface Models produced with photogrammetric and KLS datasets, and a comparison of rivers and tunnel valleys within the same geological units.The compiled data allows us to discuss and revisit the guidelines reviewed by Kehew et al. (2012) in the context of remote sensing of previously glaciated areas.
1.1 Field site: Devon Island Devon Island was covered by an extensive Innuitian ice sheet that reached its maximum extent during the last glacial maximum (e.g., England, 1987;Dyke, 1999;England et al., 2006).Shortly after the Younger Dryas, around 10 ka BP, the margin of this ice sheet was under recession to the current coast line, and the final remnants in central Devon Island vanished around 8.8 ka BP (Dyke, 1999), leaving a landscape of plateaus, fiords and deeply incised canyons.We refer to the work by Dyke (1999) and England et al. (2006) for a detailed discussion on the glacial history of the island during and since the Innuitian ice sheet.Since deglaciation, the landscape evolution is mainly the result of periglacial processes and erosion by ephemeral seasonal streams (e.g., McCann et al., 1972;Dyke, 1999).Fluvial incision represents only a small and highly localized contribution to the overall landscape evolution as a consequence of the island's polar desert climate conditions.(e.g., French, 2013).Aerial photography obtained from the National Air Photo Library (National Resources Canada) of Devon Island reveals exposed tunnel valleys (Dyke, 1999) incised into the otherwise flat plateaus that comprise the majority of the topography on the island (typical slopes are 1 to 3 • ).These channels typically drain into deeply incised canyon systems that are believed to pre-date Innuitian glaciation (Dyke, 1999).
The study site comprises an area of approximately 15 km 2 to the E-SE of the Haughton impact structure in central northern Devon Island (Fig. 1).This is a well-preserved 23 km diameter (Osinski and Spray, 2005) 23 Myr.old (Young et al., 2013) meteorite impact structure, which is a well established Mars analogue terrain, and has been the focus of numerous planetary analogue studies including crater morphology and erosion, periglacial landscape evolution on Earth and Mars, and evolution of ancient lake beds, among other activities (Lee and Osinski, 2005).The location provides access to both exposed tunnel valleys and river valleys, allowing for systematic comparisons of their geometry and longitudinal profiles.Geologically, the study area lies entirely within carbonate strata of the Upper Ordovician Allen Bay Formation, specifically the Lower Member, which comprises a uniform succession of medium bedded to massive limestone with dolomitic labyrinthine mottling (Osinski and Spray, 2005;Thorsteinsson and Mayr, 1987), and is overlain by quaternary glacial till (e.g., Dyke, 1999;Osinski and Spray, 2005).Constructional landforms such as eskers, and glacial deposits including glacier moraines and striations are rare on the The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-236Manuscript under review for journal The Cryosphere Discussion started: 23 November 2017 c Author(s) 2017.CC BY 4.0 License.plateau surface of Devon Island (Roots et al., 1963) and only occur sporadically within the Haughton impact structure (Osinski and Spray, 2005).Dyke (1999) found evidence for subglacial meltwater channels.We located individual tunnel valleys with criteria (1) and (2) (i.e., channel direction and undulating longitudinal profiles) using high resolution WorldView imagery of the site with a resolution of 2m/ pixel, a 0.15 arc-sec CDEM which corresponds to 20 by 36 m/ pixel at a latitude of 75 • N, obtained from the Natural Resources Canada website (http://geogratis.gc.ca/site/eng/extraction) and the recently released Arctic DEM (release 5) at a resolution of 2m/pixel available for free at https://www.pgc.umn.edu/guides/arcticdem/distribution/.
We selected and visited 20 drainage systems for detailed, in situ characterizations, of which 14 are tunnel valleys and 6 are rivers for comparison (see Fig. 1).For this study, we selected only first order tributaries and studied them from the origin until the first junction.Downstream of the first junction, meltwater accumulates into streams and fluvial incision is apparent in the profiles and cross sections.In addition to in situ data, we acquired helicopter airborne imagery of sites located in central and eastern Devon Island, and identified tunnel valleys as they are exposed by the current retreat of the ice cap at its terminus.

Longitudinal profile data
To detect topographic undulations in the target drainage systems, we obtained longitudinal profiles (i.e., elevation vs. distance data) of 20 tunnel and river valleys using a GARMIN gpsmap 64s with a horizontal resolution on the range 1-3m depending on polar satellite availability, and a vertical barometric resolution of 3-6 m.We acquired the data by walking or driving along each tunnel or river valley from the head until the first junction, and averaging the multiple profiles acquired in two to three runs.
To minimize the effects of the variable GPS resolution in processing the data, we grouped the channels into 5 groups that correspond to sites visited on the same day during the season.We also recorded the speed of the traverse during the collection of the longitudinal profiles, which varied between walking speed (v walk ∼ 3-5km/h) and All Terrain Vehicle speed (ATV) (v drive ∼ 10-15km/h).We then use this information to filter the GPS data in each group of tunnel and river valleys accordingly (Fig. 1 and Fig. 4).
We processed the GPS raw data by high pass filtering the signal with an upper frequency of twice the Nyquist sampling size, which corresponds to 1 GPS point per 3 seconds in the hiking traverses and 1 point per second in the ATV traverses, and a lower frequency corresponding to the inverse of the time required to drive/ walk along the channel length, that is, 1/v drive in ATV traverses and 1/v walk in hiking traverses.This filtering operation removes the data spikes related to avoiding obstacles on traverses including stream paths, large boulders, and snow patches.

Airborne imagery and photogrammetry
We complemented our in-situ channel characterizations with an extensive collection of high resolution aerial photography of over 50 tunnel valleys throughout the island, from which we derived Digital Surface Models (DSM).DSM generation The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-236Manuscript under review for journal The Cryosphere Discussion started: 23 November 2017 c Author(s) 2017.CC BY 4.0 License.through stereo-photogrammetry processing of image data involves the reconstruction of a three-dimensional body employing measurements in two or more overlapping images, acquired from different positions.Accurate reconstructions require an overlap of more than 50% between each image on a basis of at least 10 images per location, common features identifiable in different images for reference, and detailed spatial coordinates for each site.
For this purpose, we acquired over 1000 helicopter airborne images to capture the topography of multiple tunnel valleys.
To build this image database, we used a GPS-referenced CANON EOS 6D with an image resolution of 72ppp (5472 by 3648 pixels) (e.g., Smith et al., 2009).The built-in GPS has a horizontal spatial resolution of ∼ 10m and a vertical resolution of ∼ 5m.Although the camera GPS resolution also depends on polar satellite availability, resolution variations are minimal given the very small time lapse of image acquisition of all helicopter data.
To construct a Digital Surface Model (DSM), we obtained geo-referenced helicopter borne images for more than 50 channels including the centre-east of the island and the margin of the Devon Island ice cap (black arrows in Fig. 1).Significantly, this survey includes imagery of tunnel valleys currently emerging under the active Devon Island ice cap margin (Fig. 2(d)), which enabled us to ground truth our identification scheme.Figure 2 includes examples of these images acquired at different points in the island.
We process the data in several steps.For each site, we first upload the images into AGISOFT software (e.g., Tonkin et al., 2014), together with the camera-generated EXchangeable Image Format (EXIF) files that include the geo-reference information.The software automatically aligns the imagery using the overlap existent between images.We improve the initial automated alignment with manual alignment of the images by selecting and matching common features (control points).Next, we produce a dense point cloud, a meshed surface model, and a surface model with ground texture.At this stage, we use the recently released Arctic DEM (available for free download at http://www.agic.umn.edu/arcticdem) to manually introduce markers in the model with known coordinates and elevation.This step improves the resolution of the final product by an order of magnitude.With this improved 3D model, we produce an orthophoto and the Digital Surface Model (DSM).In turn, at the final stages of processing we manually crop the DSM to remove noisy areas.
The DSM model reconstructions range in resolution between 0.4 m/pixel to 10 m/pixel depending on helicopter elevation and speed, number of images captured at each site and their overlap, the number of manually introduced control points, and other factors.All products are available upon request from the authors in point cloud (.LAS) and Geotiff (.tif) formats.

Kinematic LiDAR Scan acquisition
We used a novel kinematic backpack LiDAR scanning (KLS) system to capture the the detail of tunnel valley topography from the ground.This study is the first time the KLS system has been deployed in the Arctic, and the first time it has been applied to detailed tunnel valley morphometry measurements.The goal of the survey was to reproduce the surface topography at cm resolution, but also to make a proof of concept for the capabilities of kinematic LiDAR.The system consists of a LiDAR scanner, a GNSS/GPS positioning system, and an inertial measurement unit (IMU) which are mounted within a backpack frame, allowing the user to make ultra-high-resolution LiDAR point clouds (>5,000 points/m 2 ) of features traversed by the user.KLS enables the reproduction of surface topography at cm resolution (2 cm/pixel; Fig. 4    2(e) shows the cross section of a deeply incised canyon emerging from under the ice cap. 2 (f) shows the cross section of a tunnel valley located 89.6 • W and 75.17River and tunnel valley topography data was collected by traversing the centreline of the channel by ATV, with the operator carrying the LiDAR system on his/her back.Continuous scanning was done using 150 Hz profiling, and 500 kHz pulse repetition frequencies.With these settings the maximum range was about 200 m from the scanner (i.e. a 400m-wide channel could be completely scanned), and the along-track line spacing about 1 cm with angular resolution of 1.8 mrad.
For accurate trajectory computations we set up a GNSS base station (Trimble R10) at the Haughton river valley base (75 • 22.42' N, 89 • 31.89'W) constituting of about 5 km base line length to the target channels.The raw GNSS observables were recorded at 5 Hz frequency, as were those at the KLS mapping system.The altitude data for the KLS system were recorded at 200 Hz data rate, and the positions and altitude trajectory were computed in a post-processing step for the point cloud generation.For post-processing we used the base station position using PPP method and the tightly coupled KLS trajectory Waypoint Inertial Explorer 8.60 (NovAtel, Canada).
To produce an elevation model the raw point cloud data was further processed: the points resulting from multiple reflections were removed as well as points with weak return signal (intensity less than 800 in the scanner scale).Some remaining points resulting from the laser beam hitting the rear of the ATV during the capture were manually cleaned out of the data using CloudCompare software.Post-processed data was exported as .LAS files.
The final DEMs (see example in Fig. 4 (c) and (d)) from the georeferenced LiDAR point clouds were created using .LAS Dataset tools in ArcGIS (LAS Dataset to Raster: Bin, Avg, Simple).The effective pixel resolution of 2cm/pixel represents in the DEMs represent the average value of the point cloud within a 2 cm bin (typically 5-10 LiDAR points, depending on proximity to the scanner).Due to the very high spatial coverage of the LiDAR points, minimal interpolation was needed, except in areas of LiDAR shadow.These areas were interpolated using the simple interpolation method outlined in the ArcGIS help section.

Tunnel and river valley longitudinal profiles
Reliably identifying tunnel valleys on the basis of the criteria reviewed in Kehew et al. (2012), and in particular recognizing sections where the subglacial flow eroded up against topographic gradients requires quantitative field longitudinal profile observations.Figure 3  (see Fig. 1 and Table 1 for location reference).From these data we calculate the ratio of the elevation gain in an undulation, defined as the elevation difference between the local minima at the beginning of a section with positive topographic gradient h min , and the local maxima that follows it to the total topographic loss h max : ψ = (h max − h min )/∆h total .We refer to this magnitude as the magnitude of undulation ψ and we use it herein to quantify differences between the fluvial and subglacial longitudinal profiles.Table 1 shows the magnitude of undulation ψ of the different systems in Fig. 3, together with their detailed coordinates.
Channels in group 1 show magnitudes of undulation ψ corresponding to 24%, 3% and 27% of the total topographic loss, which is equivalent to 6.6 m, 1.5 m and 6.5 m respectively for j201, j202 and j204.Channels j201 and j204 display the largest undulations that we analyzed.Group 2 corresponds to longitudinal profiles of currently active river valleys, acquired and analyzed for comparison.In all the rivers in group 2 no undulations are detected above the GPS confidence level.River profiles show a steady decrease of elevation with distance, much more consistent with other examples of fluvial channels in the literature (e.g., Whipple and Tucker, 1999).Group 3 shows two profiles corresponding to active rivers (j231 and j232) and two corresponding to tunnel valleys (j233 and j234).The first two streams show ψ = 4% and ψ = 0%, and the last two show ψ = 3% and ψ = 0% respectively.Tunnel valleys in group 3 and 4 were identified on the basis of their similar morphology and orientation to tunnel valleys in other groups, but do not present significant undulations.Within this group 4, channels j251 and j253 display undulations of ψ = 7%, j252 has ψ = 1%, and j254 displays ψ = 3%.Finally, group 5 displays 5 channels with different levels of undulation, respectively ψ = 0%, ψ = 13%, ψ = 4%, ψ = 4% and ψ = 0%.Based on similar morphology and proximity to other channels with large ψ within the same group, and recurrent N − N W to S − SE orientation, we conclude that all channels in group 5 to be originated in the subglacial regime.
Additionally, Fig. 3 offers a comparison of LiDAR (solid orange line) and GPS-acquired (dashed line) longitudinal profiles.
This was a useful indicator of the reliability of the hand-held GPS profile dataset within the GPS resolution range.We discuss the LiDAR results in more detail below.We also performed an additional comparison of our data (LiDAR and GPS) with corresponding longitudinal profiles extracted from the Arctic DEM at 5m/pixel resolution (Fig. 3 green lines).In most of the profiles the agree is excellent and well within the GPS precision margin.However, profiles j201, j211, j213, j214, j263 and j264 show significant deviations.Profiles j201, j202, j203, and j253 are also significantly noisier than the rest of Arctic DEM derived data.We attribute the discrepancy, in particular the difference in concavity in profiles j263 and j264 with our data, to the presence of snow covering the tunnel and river valleys during the acquisition of the photogrammetry data used to derive the Arctic DEM.We did not observe in the field any of the spikes present in the Arctic DEM profiles j201, j202, j253, and j254, and therefore we argue that they are DEM artefacts.

LiDAR observations
Kinematic LiDAR Scanning (KLS) was acquired in 5 tunnel valleys and one river valley.The LiDAR dataset provides high resolution topography data, which adds robustness to GPS-based undulation observations, targeting criteria (2) in Kehew et al. (2012).Furthermore, KLS highlights a difference in cross sectional shape, scale, and downstream evolution that has not yet been considered as a distinctive characteristic of subglacial erosion, and is not appreciable from GPS profile data.Surface models acquired through photogrammetry show the transition between the tunnel valley erosional regime to the river valley erosional regime (black arrows) in terms of cross sectional shape and deepening of the channels.The orthoimage in panel (b) highlights the similar width of all tributaries, which starts increasing after the junction (see black arrow).Snow and ice packs were a common view in some channels, particularly closer to the Devon ice cap.This was an issue at processing the DSM and textured image, although in some channels the thickness of the snow pack could be estimated.The high resolution product highlights the size of the first order tributaries, that are 10−20m in width from the origin, with no small scale channels or tributaries visible (see DEM and orthoimage in panel (a)).
The photogrammetry dataset and observations also highlight the drainage characteristics of tunnel valley networks.Notice, for example, how a cluster of 12 tunnel valleys merges into a single, meltwater fed river in Fig. 5 (a).Order 1 tunnel valleys throughout the island are generally parallel to each other during the first few hundred m, although we also observed braiding in some tributaries (see Fig. 5 (a) but also 4 (d)).

Morphometric comparison of tunnel valleys and river valleys
Morphological differences between river and tunnel valleys are apparent even before evaluating the criteria exposed before (Kehew et al., 2012), including correlation to former ice margins, direction consistent with estimated ice flow lines, and analysis of the presence of undulations in their longitudinal profiles.On remote sensing data including satellite imagery and topography of the field site and surrounding area, tunnel valleys appear in groups of 5 to 10, parallel to each other and consistently in the N-NW to S-SE direction.Moving to the east of the island, channel directions change on average from W to E, remaining oriented radially towards the current day Devon Island ice cap.
Characteristic tunnel valley lengths are ∼ 1km throughout the distinct targeted channel groups, all incised into the lower member of the Allen Bay formation.The cross section is shallow, flat bottomed, and steep-sided, with widths that are ∼ 5−10m approximately constant from the origin until at least the first junction, and depths of ∼ 3 − 5m, deepening in the downstream direction.In comparison, the geometry of river valleys on the island displays major differences.In particular, river valley cross sectional widths vary continuously downstream by one to two orders of magnitude from the origin (∼ 5 − 10 cm) until the first junction (1 − 5m) (i.e., Fig. 4(c) and (d)), consistent with observations elsewhere (Parker, 1978a, b, e.g.,).River valleys deepen up to ∼ 20 − 50m from the origin until the first junction and form deeply incised canyons with V-shaped cross sections (see Fig. 2(f) and Fig. 4).Another geometrical distinction between both erosional regimes is the width of first order tributaries (Grau Galofre and Jellinek, 2017).Even at the tip of the channel, tunnel valley widths amount up to tens of metres (consistent with arguments in Weertman (1972)), as opposed to widths of first order river channels which are typically sub-metre in scale.
We observed in different sites on the island how the melting of snow packs in tunnel valleys leads to the formation of ephemeral streams, which merge at the channel junctions into larger streams.This transition is often associated with a gradual change in cross section, from the sallow flat-bottomed characteristic of tunnel valleys to a deeply incised V-shape.3(a), (d), and (e).The branching pattern of subglacial channels often differs from the tree-like typical of rivers, also apparent from Fig. 2. Channels also intersect with and diverge from each other over relatively short lengthscales (∼ 100m, see tunnel valleys in Fig. 1 as an example of this anastomosing pattern).Generally, tunnel valleys also show very few or no tributaries (notice this pattern in Fig. 2(a), (b), 5(a)), rarely achieving a stream order higher than 3, much lower than typical values for river valleys.

An additional metric for tunnel valley identification
In section 3.1 we introduce the magnitude of undulation ψ and discuss its role in identifying tunnel valleys.We argue here that differences among channel direction and local topographic gradients are also indicative of subglacial erosion in areas where the ice erosion rate by sliding is lower than the meltwater erosion rate (Weertman, 1972;Paterson, 1994).
We consider the confined flow of pressurized water in a subglacial channel at the base of an ice sheet to follow the x direction, with y perpendicular to ice flow and z perpendicular to the ground surface, so that z b and z i are the bed and ice surface elevations respectively.At steady-state, water flow at the base of the ice is driven by the water pressure potential gradient ∇φ: ∇φ = −ρ i g∇z i − ∆ρg∇z b + ∇N. (1) Here ρ i is the ice density, ∆ρ = ρ w − ρ i is the density difference, g is the gravity, and N = p i − p w is the effective pressure, where p i = ρ i g(z i − z b ) and p w is the water pressure.The topography of Devon Island's plateaus is mostly flat, which implies that the controls in effective pressure gradient arise mostly from ice surface slopes ∇z i and not from surface topographic gradients ∇z b .This picture is true generally if ice surface slope is more important than bed topography, such that ρ i g∇z i /∆ρg∇z b >> 1.Although ice surface slope is correlated with topography at a regional scale it can depart from it the scale of individual channels, driving both channelized and distributed meltwater accordingly.We thus expect tunnel valleys to potentially deviate from local topographic gradients.In this case, neither metric will identify channels as subglacial, and their characterization will depend on other observations such as their individual and drainage geometries and directionality with respect to former ice lines.Criteria (4) in Kehew et al. (2012) can also provide a guide for interpretation in the last case, with an example of that in Fig. 4 group 5.In this group of tunnel valleys, 2 of the channels display no undulations above the GPS detection limit, whereas the other 3 channels do (see Table 1).According with criteria (4) in Kehew et al. (2012), and given the geometrical similarity of all channels in that same cluster and their recurrent direction, it is reasonable to propose that if three channels in the area were incised under the ice sheet, then the other 2 must have a similar subglacial origin.
The advantage of using channel directionality as opposed to undulations lies in its straightforward application, by contrasting the topographic contour lines with channel direction.It is, however, limited by the assumption that channels formed under slow sliding or cold based ice body.

Identification of tunnel valleys from remote sensing data
We distinguish tunnel valleys on the basis of four properties: (1) consistent N-NW to S-SE direction, radial to the paleo-ice margins (Dyke, 1999, Fig. 8) near our field area, changing to W to E near the ice cap margin; (2) Topographic undulations in the longitudinal profiles (Fig. 4 and Table 1), and channel incision with orientations not parallel to the local topographic gradient (Fig. 2(d)); (3) Cross-section size and shape, i.e., shallow wide flat-floored for tunnel valleys and deeply incised V shape for river valleys; and (4) large 1st order channel widths on the order of ∼ 10 m.Not all channels we identify as tunnel valleys from their morphology and direction have undulations in their longitudinal profiles.However, none of the river valley profiles show any detectable undulation.We conclude that the magnitude of the undulation index ψ > 0 unequivocally distinguishes tunnel valleys, but ψ = 0 does not necessarily preclude subglacial incision.
In revisiting the guidelines reviewed in Kehew et al. (2012) and their application to remote sensing data, we find that we can target criteria (2) undulations in the longitudinal profile, with high resolution topographic data (i.e., GPS tracks in Fig. 4, and KLS LiDAR data in figure 5), and criteria (1) consistent directionality relative to active or former ice flow lines, with imagery and topography (i.e., consistent and parallel channel direction, which does not follow topography in Fig. 1).The remote sensing data applications of criteria (3) are limited, as they require additional geological constrains on the location and margins of former ice bodies.We assessed criteria (4) with the morphological similarity among tunnel valleys within the same group.
To these well-established criteria, we introduce three useful indicators of subglacial channel erosion from remote sensing data.The first is the difference between tunnel valley direction and topographic gradients, as exemplified in Fig. 2 (d) and discussed in section 4.1.This criteria is straightforward to evaluate, and requires only a measure of the angle between the topographic gradient direction (i.e., perpendicular to contour lines) and channel direction.The second indicator is the large cross-section widths at the origin of order 1 tunnel valleys (i.e., Fig. 5) which are orders of magnitude larger than in river

Conclusions
In this study we describe a series of tunnel valleys (N-channels) exposed on Devon Island, Canadian Arctic Archipelago, that were emplaced during the retreat of the Devon Island ice cap to its present day location, 8-9.5kyr ago.In particular, we discuss the use of remote sensing techniques to distinguish between systems of river and tunnel valleys, which serves as a complement to existing field based methods.We provide detailed field descriptions of 20 tunnel valleys, including their longitudinal profile characteristics, cross sectional geometry, channel directionality and drainage network morphology, to then revisit and expand the identification methods reviewed by Kehew et al. (2012).Our field observations include GPS mapping of subglacial and fluvial incised channels longitudinal profiles, photogrammetry, kinematic mobile LiDAR data (KLS), and aerial imagery, allowing for both a qualitative and quantitative description.
We find that a quantitative measure of undulation ψ, defined as the topographic loss (local minima to local maxima of the undulation) at an upstream section to the total topographic loss, reliably distinguishes tunnel and river valley profiles (Fig. 4), although the lack of undulations does not rule out subglacial erosion.We also argue that the departure in channel direction from local topographic gradients also reflects subglacial erosion.We then discuss the limitations of both these metrics in identifying tunnel valleys.If both metrics fail, other morphological observations such as channel direction, clustering in groups of parallel, closely spaced tunnel valleys, and spatial correlation with other subglacial features, i.e., other tunnel valleys, are key to identify the channels.
With our data and observations, we revisit the guidelines reviewed by Kehew et al. (2012) to improve identification of tunnel valleys from remote sensing data.We conclude with the following target characteristics of interest: (1) undulations in the longitudinal profile and changes in channel direction with respect to local topographic gradients, (2) consistent channel direction radial from present or former ice sheets and close to the ice margin, (3) order 1 channel widths on the 5-10 m scale that show minimal variation downstream, with wide, flat-bottomed cross sections, and (4) presence of other sublacial features in the area.

Figure 1 .
Figure 1.(a) satellite imagery of Devon Island within the Arctic Archipelago (white box).(b) satellite image of Devon Island, with a white box indicating the selected field site.The map also shows the Innuitian ice sheet terminus lines digitalized from Dyke (1999), with age reference in the legend.(c) Field site (UTM zone 16), with black arrows showing the mapped groups of tunnel and river valleys.

Figure 1 4
Figure 1(a) shows a Digital Globe satellite image of the Canadian Arctic including Devon Island and the rest of the Arctic Archipelago for context.This panel also shows the Innuitian ice sheet termini as digitized from the work by Dyke (1999).5

Figure 1
Figure 1(b) shows a high resolution satellite view of our selected field area.Our target locations consist on 4 distinct groups of tunnel valleys and 1 group of fluvial valleys.To identify the areas incised by subglacial erosion, we used criteria (1) and (4) in Kehew et al. (2012) (i.e., parallel to former ice flow lines and close to other subglacial features) in areas easily accessible from the Haughton structure, and also based on the locations where Dyke (1999) found evidence for subglacial ), which is a higher level of 6 The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-236Manuscript under review for journal The Cryosphere Discussion started: 23 November 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 2 .
Figure 2. Aerial and field imagery of tunnel and river valleys.2(a) corresponds to helicopter imagery of a group of tunnel valleys (89.08 • W,75.28 • N). 2(b) corresponds to a groups of tunnel valleys located at 89.37 • W, 75.18 • N. Notice the recurrent parallel direction and scale.
2(c) corresponds to a group of tunnel valleys emerging underneath the Devon Island ice cap, notice the similar morphology to 2(a) and 2(b).
2(d) shows a tunnel valley located at 85.18 • W, 75.17 • N following a direction perpendicular to the local topographic gradient (black arrow).
• N.accuracy compared to e.g.airborne LiDAR data, which would also be expensive to obtain in remote areas such as the high 7The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-236Manuscript under review for journal The Cryosphere Discussion started: 23 November 2017 c Author(s) 2017.CC BY 4.0 License.Arctic.The KLS dataset and derived surface models improve cross sectional and longitudinal profile analysis and comparison with river valleys, and are a ground-truth for hand-held GPS topographic data acquisition.The LiDAR data was acquired using the AkhkaR3 kinematic backpack LiDAR developed at the Finnish Geospatial Research Institute, which is an updated version from those presented inKukko et al. (2012) andLiang et al. (2015).The system is based on GNSS-IMU (Global Navigation Satellite System-Inertial Measurement Unit) positioning system, consisting on Novatel SPAN: Flexpak6 receiver, UIMU-LCI inertia measurement unit and 702GG antenna, a 360 degrees of field of view crosstrack profiling laser scanner (Riegl VUX-1HA) synchronized to the positioning and operated by a tablet computer (Panasonic Toughpad ZF-G1).
shows the longitudinal profiles of the different channels in each group visited within the field area The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-236Manuscript under review for journal The Cryosphere Discussion started: 23 November 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 3 .
Figure3.Longitudinal profiles of tunnel and river valleys normalized to total topographic loss and length along the channel, with context satellite imagery for each channel group.Group 1 consists of 3 tunnel valleys.Group 2 includes 4 rivers.Group 3 is in the same location than group 1 (see Fig.1) but consists on 2 rivers (j231, j232) and 2 tunnel valleys (j233, j234).Group number 4 includes 4 tunnel valleys.Group 5 consists of 5 tunnel valleys.In the longitudinal profiles, blue crosses represent the raw GPS data for each channel, blue dashed lines are the data after filtering, and orange solid lines represent the LiDAR sections that overlap GPS data for comparison.The profiles obtianed with the Arctic DEM at 5m resolution are shown in green color.

Figure 4
Figure 4 shows the results of using the kinematic backpack LiDAR approach to imaging the topography of a channel.The first panel 4(a) shows the point cloud files produced when investigating the channel cross section.The data is coloured by back scattered intensity at a laser wavelength of 1550 nm used in the KLS LiDAR system, resulting in darker values for wet snow and ice as seen in the image.The ATV (roughly 1.75m long) gives a scale for the illustrations.Panel 4(b) shows the trajectory of the KLS user in blue lines, with the ATV again for scale.Panel 4(c) shows the raster derived from the point cloud files for a river valley (corresponding to j231), and 4(d) shows the raster for the tunnel valleys in group 1.

Figure 4 .
Figure 4. KLS LiDAR observations.Panels (a) and (b) show the point cloud files color coded according to back scattered intensity (dark is low return).Panel (b) shows the trajectory of the KLS user overlapped to the point cloud product for a reference of coverage.Panels (c) and (d) show the raster produced using the point clouds, at a resolution of 9cm/pixel.The black arrows in panels (a) and (b) point at the ATV for scale.Point spacing in 4(a) and 4(b) corresponds to 6 cm, with a total point count of 117,147,558 points for the river valley and for the tunnel valleys in group 1. Raster resolution corresponds to 9 cm in the channels in group 1, and 10cm in the river valley corresponding to j231.

Figure 5 .
Figure 5. Stereo-photogrammetry products derived from helicopter borne photography.The top panels (a) and (b) show the digital elevation model (DEM) at a resolution of 0.48 and 0.56 m/pixel respectively, with the colorbar indicating the elevation of the model surfaces.The images underlying the panels correspond to the textured orthoimages in both locations.
Fig. 2(e)    and 5 exemplify this morphological transition from a group of tunnel valleys to a single meltwater fed river.At a network scale, tunnel valleys occur in clusters of 5-10 channels with finger-like patterns including a large curved channel surrounding the group of tunnel valleys, and several smaller parallel internal systems.Examples of this pattern are in Fig. 2(b), The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-236Manuscript under review for journal The Cryosphere Discussion started: 23 November 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 6 .
Figure 6.Cartoon representing an ice body sliding over its bed, showing the axis notation for reference, and the ice and topographic surfaces zi and z b .

The
Cryosphere Discuss., https://doi.org/10.5194/tc-2017-236Manuscript under review for journal The Cryosphere Discussion started: 23 November 2017 c Author(s) 2017.CC BY 4.0 License.systems.The third and last criteria is the minimal variation in width downstream, from the beginning of the channel until the first junction, as visible from Fig. 5, comparison of panels (c) and (d).

The
Cryosphere Discuss., https://doi.org/10.5194/tc-2017-236Manuscript under review for journal The Cryosphere Discussion started: 23 November 2017 c Author(s) 2017.CC BY 4.0 License.Science Foundation awards 1043681, 1559691, and 1542736.LiDAR studies were made possible by financial funding from the Academy of Finland for "Multi-spectral personal laser scanning for automated environment characterization (300066)" and "Centre of Excellence in Laser Scanning Research (CoE-LaSR) (272195)".We thank the Polar Continental Shelf Project (PCSP) for their logistical support without which this work would not have been possible.17 The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-236Manuscript under review for journal The Cryosphere Discussion started: 23 November 2017 c Author(s) 2017.CC BY 4.0 License.