Interactive comment on “ Morphological properties of tunnel valleys of the southern sector of the Laurentide Ice Sheet and implications for their formation

This paper is a contribution to the ongoing research on tunnel valleys, which certainly is a very important topic in paleoglaciology that needs to be deepened if the behavior of past ice sheets (stability, meltwater drainage, landforming processes) is to be better understood. Therefore, this study is a timely contribution that surely would not pass unnoticed and I fully sympathize with the authors in their effort to learn more about these fascinating landforms that have attracted a lot of attention and generated yet more controversies. The aim of this paper is to illuminate the origin of the tunnel valleys in the southern sector of the Laurentide Ice Sheet, in particular to test the catastrophic formation theory against the gradual formation theory. In order to do so around 2000 tunnel valleys have been mapped yielding valuable data on their geometrical characteristics


Introduction
Incised into bedrock or sediment, tunnel valleys are elongate depressions up to several kilometres wide, often with undulating long-profiles, tens of kilometres long and tens to hundreds of metres deep.They are observed in many formerly glaciated landscapes around the world, and tend to be orientated parallel to the direction of former ice flow (e.g.Ussing, 1903;Wright, 1973;Attig et al., 1989;Wingfield, 1990;Piotrowski, 1994;Patterson, 1997;Huuse and Lykke-Anderson, 2000;Jørgensen and Sandersen, 2006).Features with similar dimensions have also been described beneath current ice masses (e.g. Rose et al., 2014).Since first being described (Gottsche, 1897) and then attributed to erosion by subglacial meltwater (Ussing, 1903), such phenomena have attracted a variety of names, which vary according to interpretations about how they form: Tunnel Channels (implies whole depressions occupied and cut by water); Linear Incision (purely descriptive term); Tunnel Valley sensu stricto (taken by some to imply a large linear depression created by activity of a smaller channel which occupied part of it, as in a river channel producing a valley); or the rather vague term Palaeo-valley.In this paper we stick to the most widely used term -Tunnel Valley -but use it in its broadest sense (sensu lato) to include depressions that could actually be tunnel channels.In this manner we initially treat all linear depressions of the appropriate scale and morphological charac- teristics as the same thing, whilst recognising that in detail this might be a grouping of a number of types.
Tunnel valley formation is typically attributed to subglacial meltwater erosion at the base of ice sheets (cf.Ó Cofaigh, 1996;Kehew et al., 2012;van der Vegt et al., 2012), and they are considered an important component of the subglacial hydrological system, providing drainage routeways for large volumes of water and sediment.Understanding their genesis is relevant for reconstructing former ice sheets, elucidating basal processes and exploiting the geomorphological record in a way that is useful for modelling subglacial hydrology.However, despite being debated for over 100 years, there is still uncertainty over how tunnel valleys form.This debate is focused around two genetic models: "outburst" formation and "gradual or steady-state" formation (Fig. 1) (cf.Ó Cofaigh, 1996;Kehew et al., 2012;van der Vegt et al., 2012).
The "outburst" hypothesis (Fig. 1a) ascribes the erosion of tunnel valleys to rapid drainage of sub-or supraglacially stored meltwater.Contemporary observations from the Antarctic and Greenland ice sheets demonstrate the efficacy of meltwater storage and drainage in sub-and supraglacial environments (Zwally et al., 2002;Wingham et al., 2006;Fricker et al., 2007;Das et al., 2008) and it is reasonable to expect that the Laurentide Ice Sheet experienced similar events.In addition, the impoundment of meltwater behind an ice margin frozen to its bed has been linked to tunnel valley formation, for example, along the southern terrestrial margins of the former Laurentide and European ice sheets where permafrost was prevalent (e.g.Piotrowski, 1994;Cutler et al., 2002;Hooke and Jennings, 2006).Genesis is typically thought to occur via repeated low to moderate magnitude floods that may be at or below bankfull flow (e.g.Wright, 1973;Boyd et al., 1988;Wingfield, 1990;Piotrowski, 1994;Cutler et al., 2002;Hooke and Jennings, 2006;Jørgensen and Sandersen, 2006).Catastrophic erosion of entire tunnel val-ley networks by massive sheet floods has also been proposed (e.g.Shaw and Gilbert, 1990;Brennand and Shaw;1994;Shaw, 2002), but has been considered less likely given the very large volumes of stored water required (e.g.Ó Cofaigh et al., 1996;Clarke, 2005).
The "gradual" or "steady-state" hypothesis (Fig. 1b) typically invokes erosion of soft-sediment beds in low-pressure subglacial channels (Boulton and Hindmarsh, 1987;Mooers, 1989;Praeg, 2003;Boulton et al., 2009).In this model, high water pressures transmitted through the substrate to the ice-sheet terminus initiates failure and headward erosion of a conduit (by piping) (Shoemaker, 1986;Boulton and Hindmarsh, 1987;Hooke and Jennings, 2006;Boulton et al., 2009).As the fluid pressure of the conduit is lower than the surrounding substrate, meltwater flows towards the conduit, the walls are enlarged by sapping (i.e.undermining and headward recession of a scarp) and the sediments are mobilised and transported away by the resulting subglacial stream (Boulton and Hindmarsh, 1987).In general, enlargement is suggested to occur via steady-state Darcian flow of water into the conduit (e.g.Boulton and Hindmarsh, 1987;Boulton et al., 2007aBoulton et al., , b, 2009)).Hooke and Jennings (2006) adapted this hypothesis, suggesting that initial headward erosion by piping was followed by more rapid enlargement when the conduit tapped into a subglacial lake, thereby combining both scenarios in Fig. 1.Ravier et al. (2014) emphasised the potential influence of localised high porewater pressures in promoting efficient erosion by hydrofracturing and brecciation of the preglacial bed, while Mooers (1989) considered supraglacial drainage to the bed rather than basal meltwater as the dominant source for gradual tunnel valley erosion.
A range of approaches can be applied to the investigation of tunnel valleys including theoretical, sedimentological and morphological.Thus far, most efforts have used a combination of these approaches, with a lot of data, description and detail, but for a small number of tunnel valleys (see Sect. 2).From these studies it is difficult to extract representative information of the population of tunnel valleys or to gain an understanding of the broader-scale distribution of landforms.Based on previous studies and the availability of digital elevation models (DEMs), we are now able to undertake a systematic and large-scale mapping campaign of the size, shape, spatial arrangement and distribution of tunnel valleys to better understand the spatial properties of this phenomenon, noting that it is useful to know more precisely what it is that requires explanation (e.g.Dunlop and Clark, 2006, for ribbed moraine).In doing so we will address the following questions: (1) how do we define a tunnel valley and how can they be distinguished in a glaciated landscape?(2) What are their morphological characteristics?(3) Is there a characteristic distribution and arrangement?(4) What can associations with other landforms tell us about tunnel valley formation?The southern sector of the Laurentide Ice Sheet was selected because it contains thousands of these landforms and the distinctive geometry of the ice lobes provides information on the water drainage pathways (Fig. 2).Our mapping builds on, and replicates in many places, comprehensive local and regional studies, which include sedimentological details that we draw on.The main purpose of this paper is to present systematic mapping of a large sample of tunnel valleys and to provide basic metrics on their variation in size, morphological characteristics and distribution.We do so to advance knowledge about these landforms and provide representative data that modellers of subglacial hydrology might find useful.We then use our morphological observations and data, along with published field observations to assess existing theories on tunnel valley formation.An obvious limitation to using DEMs of the current land surface is that any post-formational deposition (infilling) of tunnel valleys will mask or modify aspects of the morphology we can see.That it is easy to identify thousands of tunnel valleys in DEMs however, shows that infilling is often only partial and allows us to assess their presence, distribution, width and minimum length, and some properties about undulating long profiles.
Tunnel valleys are commonly observed on the bed of the southern margin of the Laurentide Ice Sheet and typically occur as distinct radiating sets of regularly spaced depressions associated with eskers and terminal or recessional moraines (cf.Kehew et al., 2012).At the bed of the Saginaw Lobe, for instance, valleys are typically spaced at 6-10 km intervals (Fisher et al., 2005;Kehew et al., 2013).Tunnel valleys are incised into glacial sediments up to a depth of 25 m and extend for < 50 km (e.g.Jennings, 2006).However, tunnel valleys up to 150 km long have been documented in the Superior Lobe, Minnesota (Wright, 1973), and valleys are eroded up to 200 m into the bedrock floors of Lake Superior and Lake Michigan (Regis et al., 2003;Jennings, 2006).
Although tunnel valleys are typically sub-parallel, they are also observed to join, split and even cross-cut each other (e.g.Wright, 1973;Mooers, 1989;Kehew et al., 1999Kehew et al., , 2005;;Fisher et al., 2005;Kehew and Kozlowski, 2007).Cross-cutting relationships, both between tunnel valleys and with other glacial landforms (e.g.drumlins, outwash fans, moraines), record a palimpsest signature of tunnel valley erosion.In the Saginaw Lobe, Kehew et al. (1999Kehew et al. ( , 2005) ) and Kehew and Kozlowski (2007) identified a series of palimpsest associations in which partially buried tunnel valleys pass beneath terminal moraines, diamicton and surficial outwash associated with later advances.This palimpsest style is interpreted to result from the collapse of ice and debris into the valley, which becomes (partially) buried by sediment during a re-advance and then re-emerges as the ice melts out (e.g.Kehew and Kozlowski, 2007).
Tunnel valley cross-sectional morphology ranges from sharply defined with constant or downstream increasing dimensions (e.g.Mooers, 1989), to indistinct valleys often associated with hummocky terrain and characterised by beaded or crenulated planforms, or as a series of aligned depressions (e.g.Kehew et al., 1999;Sjogren et al., 2002).Indistinct valleys may be due to partial burial during re-advance events or by melt out of debris rich ice obscuring them (Kehew et al., 1999).Sjogren et al. (2002) also identified indistinct valleys in Michigan that are eroded into the hummocky terrain.
In Wisconsin, Michigan and Minnesota, bands of hills are observed to occur upstream of tunnel valleys (Johnson, 1999).These are interpreted as erosional remnants of an anastomosing subglacial meltwater system that drained through the inter-hill depressions.At their downstream end, tunnel valleys often terminate at outwash fans some of which contain coarse boulder-gravel material (e.g.Attig et al., 1989;Mooers, 1989;Patterson, 1994;Clayton et al., 1999;Johnson, 1999;Kehew et al., 1999;Cutler et al., 2002;Derouin, 2008).Palaeo-discharge estimates from the boulder deposits imply large discharges, and this has been used as evidence for outburst flood events (Cutler et al., 2002).

Data sets and mapping
For this study, we used the National Elevation Data set (NED) (http://nationalmap.gov/elevation.html),utilising DEMs with a resolution of 1/3 arcsec (∼ 10 m) across the entire study area, and 1/9 arcsec (∼ 3 m) in some locations.Surficial and bedrock geology maps (e.g.Fullerton et al., 2003;Soller et al., 2011) were also used to aid identification and interpretation.Glacial landforms were digitised directly into a Geographical Information System (GIS).Polylines were used to map tunnel valleys sides and thalwegs, eskers and moraines.Polygons were used to map hill-hole pairs, outwash fans and dissected hills.

How do we distinguish tunnel valleys in a glaciated landscape?
Apart from tunnel valleys, large elongate depressions with similar dimensions may also form by fluvial erosion (river valleys), proglacial meltwater erosion (spillways), subglacial abrasion and/or plucking (overdeepenings), or arise from geological structures (e.g.fault lines).These phenomena are readily observed today and the formative mechanisms are reasonably well known.In contrast, tunnel valleys have yet to be observed actively forming beneath, or at the margins of, modern day ice sheets, and so their genesis and properties are more difficult to discern.In glacial landscapes they have been distinguished by their large size and characteristics such as their orientation parallel to inferred ice flow, undulating long profiles and associations with subglacial bedforms and eskers; all pointing to a subglacial origin.In particular, undulating thalwegs and their association with eskers and outwash fans, permit them to be distinguished from subaerial rivers.However, negative evidence (e.g.no esker found in a valley) does not necessarily preclude a subglacial origin.Linear incisions similar to tunnel valleys but of much smaller size (tens of metres in width) and generally called subglacial meltwater (or Nye) channels are also common in glaciated landscapes (e.g.Greenwood et al., 2007) but it is generally presumed that these are not part of the same population as tunnel valleys, that they are different landforms distinguished by size but perhaps also by process.
For the purposes of this study we restrict our definition of a tunnel valley to subglacially eroded channel-forms and use the term non-genetically in reference to both tunnel valleys and tunnel channels.Tunnel valleys that could clearly be differentiated as being eroded into bedrock were not mapped as their formation is more difficult to decipher from geological structures or glacial overdeepenings and valleys abraded and plucked by overlying ice.All valley forms that potentially could be interpreted as tunnel valleys or tunnel channels were mapped, and then each was tested to see if it could be shown to have been formed subglacially, and thus, be interpreted to be a tunnel valley or tunnel channel.One way to strengthen a subglacial interpretation would be to demonstrate that the longitudinal profile slopes upward towards an associated ice margin or that the profile undulates.To determine whether the valley bottom is undulating the number of negative and positive slope segments over 100 m length scale were calculated (see later with regard to problem of valley infills contaminating these assessments).Each valley was then assigned a confidence level from one to three, with one being the most certain and three the least (Fig. 3).Channels lacking undulations and that do not contain subglacial bedforms are difficult to differentiate from proglacial or postglacial channel systems and were therefore given a confidence of three.
Valleys with an undulating long-profile, which contain eskers or terminate in outwash fans were classified as "certain" tunnel valleys and given a confidence level of one (Fig. 3a-d).
Only those tunnel valleys with a confidence level of one or two were used in the spatial and morphological analyses.In panels (a), (c) and (d) the tunnel valley clusters terminate at a moraine position.The large valley in panel (e) (Superior) is assigned a confidence level of 3 as it does not contain any subglacial bedforms and exhibits a gradual and consistent change in bed slope consistent with a proglacial spillway.However, the smaller NW-SE valleys that it bisects is given a confidence of 2 as they have undulating thalwegs that cut across moraines.The dendritic valley cluster in panel (f) (North Dakota) is given a confidence of 3 as it is not associated with any subglacial bedforms and has a consistent bed slope indicating water flow towards the south.A braided channel morphology and a widening reach towards the south allows us to interpret this valley system as a proglacial spillway (fed by tunnel valleys emanating from under the ice to the north).

Tunnel valley measurements
Tunnel valley length was computed in the GIS.These measurements are best treated as minimum bounds, because if some valleys have a complete infilling in their upstream reaches we would not be able to recognise this part of the valley and therefore underestimate its length.Where two or more tributaries coalesce, the longest routeway was used to determine length.Tunnel valley width (distance between mapped valley sides) was measured from cross-profile transects positioned at 1 km intervals along the course of each tunnel valley.The relationship between local elevation gra- dient (G loc ) on along valley changes in width (W loc ) was calculated at each 1 km interval (j ) using Eqs.( 1) and (2): where E = elevation, W = width and D = downstream distance from the head of the tunnel valley.To calculate tunnel valley spacing it was necessary to restrict our analysis to clusters comprising distinct populations of similar orientation, which were likely formed during a similar time period (although they may not all have been operating at the same time).We calculated the spacing of 966 tunnel valleys organised in 24 discrete clusters.Spacing (S) was calculated from cross-profile transects orientated perpendicular to the direction of the cluster and positioned at 5 km intervals along each long profile.A median spacing value and the standard deviation (σ ) was calculated for each tunnel valley cluster.To provide an indication of regularity per cluster the coefficient of variation (σ/ mean spacing), expressed as a percentage (σ %), was also calculated (Hovius, 1996;Talling et al., 1997); tunnel valley clusters with a low σ % exhibit low variability in spacing.
To investigate drainage evolution during deglaciation, a subset of meltwater features in Wisconsin were grouped into "drainage-sets", defined as a collection of features interpreted as having formed during the same drainage phase.This was based on cross-cutting relationships (e.g. between channels, outwash fans and moraines) to reconstruct a relative history of drainage activity as the ice sheet retreated across the area.Cross-cutting relationships between tunnel valleys and moraines were classified according to whether the tunnel valley (1) terminates at a moraine at its downstream end and therefore formed contemporaneously with it; (2) is overlain by moraines along its length, thus suggesting that the tunnel valley was no longer active when the moraines were deposited; or (3) breaches moraines along its length, thereby indicating that the tunnel valley continued to drain water, either destroying pre-existing moraines or preventing them from forming during retreat.

Limitations
The partial or complete burial of tunnel valleys by glacial and post-glacial infill (cf.Kehew et al., 2012;van der Vegt et al., 2012) limits the mapping and measurement of tunnel valleys from DEMs.Significantly, we do not extract information on tunnel valley depth for this reason.Width measured at the valley edges is not affected.As noted earlier, our measurements of length are minimum bounds if upstream infilling sufficiently hides the valleys.The long profiles were used to determine whether the valley formed subglacially (i.e.whether there are undulations).As infilling is likely to selectively occur in hollows smoothing the long profile, this is likely to result in some tunnel valleys being missed because we could not ascertain if they had truly undulating long profiles and so our classification of linear incisions as tunnel valleys is therefore a minimum bound.
Completely buried tunnel valleys with no surface expression cannot be identified by the mapping and could therefore present a problem of bias in assessing spatial distributions.However, unlike elsewhere in the world few buried tunnel valleys have been reported along the southern margin of the Laurentide Ice Sheet, and most of the completely buried tunnel valley networks in Europe relate to pre-Weichselian ice advances (e.g.Jørgensen and Sandersen, 2006;Kristensen et al., 2007Kristensen et al., , 2008;;Stewart and Lonergan, 2011), so we suppose this problem to be minimal.Moreover, unless there is a systematic bias predisposing burial in some locations over others then the mapped pattern and distribution of tunnel valleys is likely to be informative.In analysing this large data set of tunnel valleys along the southern margin of the Laurentide Ice Sheet we make the initial presumption that tunnel valleys have a common genesis and then search for circumstances and data that challenge this.This allows us to focus on possible relationships between form, distribution and process.

Results
4.1 Is there a characteristic spatial distribution and arrangement?

Distribution
Figure 4 shows the distribution of all 1931 tunnel valleys (1694 of which have a confidence of 1 or 2) mapped beneath the terrestrial southern sector of the Laurentide Sheet.We estimate that ∼ 80 % of these tunnel valleys have been previously identified and mapped during more localised investigations (e.g.Wright, 1973;Attig et al., 1989;Patterson, 1997;Fisher et al., 2005).The map reveals a tendency for tunnel valleys to cluster together with large intervening areas where no or very few valleys occur.Indeed, while clusters of tunnel valleys are prevalent along much of the Wadena, Itasca, Superior, Chippewa and Saginaw ice lobes, they are rarer and more dispersed or isolated at the southernmost (LGM) margins of the James, Des Moines, Lake Michigan and Huron-Erie ice lobes (Fig. 4).Those that do occur in these ice lobes tend to be positioned up-ice, either at the lateral margins of the LGM lobes (e.g. Green Bay Ice Lobe) or at recessional moraines (e.g.Des Moines Ice Lobe).Tunnel valley clusters often occur down ice-flow of depressions in the landscape (Fig. 4).For example, the Saginaw Lobe tunnel valleys emanate from an arm of the present-day Lake Huron Basin, the Langlade and Chippewa tunnel valleys are downstream of the present-day Lake Superior Basin, and tunnel valleys occur downstream of the low-relief trough of the Des Moines Lobe.Based on modelled hydraulic potential surfaces, Livingstone et al. (2013) predicted that the Lake Superior Basin and NE sector of the Lake Michigan Basin were sites of several subglacial lakes during the last glacial (marked in Fig. 4).There appears to be no clear link between these potential lakes and tunnel valleys.On the other hand, subglacial lakes may also have been present elsewhere in the Great Lake basins and it is noteworthy that tunnel valleys are commonly found downstream of these basins.

Spatial arrangement
The overall shape of tunnel valley clusters varies (Fig. 4), with both broad clusters composed of many short valleys (e.g. Green Bay, James and SE edge of Superior), and narrow clusters composed of long valleys (e.g.Superior, Huron-Erie and Langlade).Cross-cutting of tunnel valleys occurs both between and within clusters.
Tunnel valley spacing (Fig. 5) displays a positively skewed, unimodal distribution with a median spacing of 4.5 km and standard deviation of 4.6 km (σ % = 81).However, the median spacing of individual tunnel valley clusters ranges from 1.9 to 9.1 km.Tunnel valleys in the Green Bay (median: 2.9 km), Superior (median: 3.7 km) and Huron-Erie (median: 1.9 km) lobes are closely spaced.Conversely, tunnel valley clusters in the large Saginaw (median: 5.7 km), Michigan (median: 5.5 km) and Des Moines (median: 5.4 km) lobes and in North Dakota (median: 5.1 km) have a wider than average spacing.In all of the measured clusters the standard deviation of the spacing is less than the mean spacing, and 9 of the 24 clusters are < 60 %.There is no significant correlation between the number of tunnel valleys within a cluster (ranging from 7 to 169) and the standard deviation, but the standard deviation increases as the mean and median cluster spacing increases, hence the use of the coefficient of variation (σ %).

What are the morphological characteristics of tunnel valleys?
The lengths of mapped tunnel valleys display a unimodal, positively skewed distribution, which is approximately lognormal (Fig. 6a).Lengths range from 200 m to 65 km, with a mode of 7-9 km, median of 6.4 km and standard deviation of 8 km.Long and short tunnel valleys occur in most places, although long valleys are less common in the Green Bay and Huron Erie lobes, and dominate in the Superior, Langlade, Wadena, Michigan and Saginaw lobes.The widths of mapped tunnel valleys display a unimodal distribution with a positive skew, which approximates normal when log-transformed (Fig. 6b).Tunnel valley widths vary considerably across the study area, ranging from 15 m to 6.7 km, with a mode of 600-800 m, median of 560 m and standard deviation of 660 m.The Chippewa, Langlade and Michigan valleys are consistently wide (typically > 600 m), while the Huron-Erie, Superior, Green Bay and Des Moines valleys are narrow (< 600 m).Other clusters, in the Saginaw, Superior and Wadena lobes, comprise a mix of wide and narrow valleys.There is a tendency for longer tunnel valleys to be wider (power law function, r 2 = 0.34, p value =< 0.001) (Fig. 7).
Tunnel valley planform shape varies across the study area (Fig. 8).The majority consist of a single valley "thread"; more than two orders of "stream ordering" are rare and tributaries tend to be restricted towards valley heads (Figs. 3,4,8).Valley margins range from sharp to indistinct and from crenulated to straight.Straight margins are more typical of long, thin tunnel valleys (Fig. 8a, d, f).However, many margins are crenulated, with bulbous and abrupt angular morphologies that result in large down-valley changes in width (Fig. 8a-f).Figure 9 demonstrates a weak relationship between tunnel valley width and distance downstream.Valleys both widen and narrow downstream with considerable and abrupt variations in width.The variation in tunnel valley width bears no relation to the local elevation gradient (Fig. 10).Local along-valley elevation gradients are relatively low (typically < ±1.5 • ) and valleys widen and narrow on both reverse and normal slopes.
Tunnel valleys and tunnel valley segments often start and end abruptly and can appear fragmented or contain bulbous depressions (Fig. 8).The gaps between segments of tunnel valleys may show no evidence of modification (Fig. 8e, f); are partially incised by narrower and more discontinuous val-leys or sets of parallel valleys (Fig. 8e); or consist of a series of depressions and hummocks with indistinct valley planform (Fig. 8a, b, d).The up-glacier ends of tunnel valleys range from rounded heads with steep sides (amphitheatre) (Fig. 8a, f) to open or indistinct (Fig. 3c-d).In Fig. 8e-f, tunnel valleys comprise parallel tracks of two or more tightly spaced (< 1 km) valleys.

What can associations with other landforms tell us
about tunnel valley formation?

Moraines
The association between moraines and tunnel valleys varies with some valleys cutting through moraines (Fig. 11a); while in other locations moraines are superimposed on the valley or the valley terminates at a moraine (Fig. 11b).In Fig. 11a, in North Dakota, tunnel valleys cutting through an end moraine are observed to narrow and then trend down-glacier into esker and outwash fan deposits.Up-glacier of the end moraine are low relief (1-2 m) and regularly spaced transverse ridges ("washboard" moraine).They have a cuspate geometry with the horns pointing up-glacier and converging on tunnel valley positions (see also Stewart et al., 1988;Cline et al., 2015).
Figure 11b shows examples of tunnel valleys terminating at, cutting through and overlain by recessional moraines.The tunnel valleys here do not show a consistent pattern, with neighbouring valleys exhibiting different moraine associations.Some valleys are continuous or semi-continuous, with a single outwash fan at, or just down-glacier from the terminus, and a series of on-lapping recessional moraines upglacier.Elsewhere, valleys contain multiple outwash fans deposited at successive moraine positions.

Hill-hole pairs
We mapped 12 hill-hole pairs (Bluemle and Clayton, 1984), 11 of which are found in North Dakota.Typically, hill-hole pairs comprise isolated features, but 4 of them are associated with tunnel valleys (e.g.Figs. 3c,12).These seem to occur at the down-glacier end of the valleys, with smaller channels and eskers emanating from and diverging around the ice-thrust hill (Fig. 12a, b).In Fig. 12a, an esker emanating from one of the hill-hole pairs trends into another tunnel valley segment further down-glacier.

Outwash fans
We mapped 187 outwash fans across the study area, predominantly at the downstream end of, but also within and between segments of tunnel valleys at moraine positions (Fig. 11b).Outwash fans are particularly common in the Chippewa, Wisconsin, Langlade, Green Bay, Superior and Wadena lobes (see also Attig et al., 1989;Clayton et al., 1999;Cutler et al., 2002;Fisher and Taylor, 2002), and at the downstream end of the large bedrock tunnel valleys in Lake Superior (Regis

Giant current ripples
In Minnesota the floor of one tunnel valley is shown to contain regularly spaced sinusoidal bedforms orientated roughly perpendicular to the valley axis (Fig. 13).The bedforms are 0.2-1.9m high (H), 10-60 m long (L) and their crests are straight to slightly sinuous.Our data show that longer bedforms tend to be higher (linear regression, r 2 = 0.5), and that the H : L ratio is ∼ 0.02.The tunnel valley that the bedforms are constrained within is partially incised into underlying drumlins orientated obliquely to the valley long axis.An esker running NW-SE is overprinted on the sinusoidal bedforms.The southern end of the valley is bisected by a large (1 km diameter) circular incision with an intact central island.
The dimensions and shape of the transverse sinusoidal bedforms, the tendency for longer bedforms to be higher and their context in the base of the tunnel valley is consistent with giant current ripples (e.g.Bretz et al., 1956;Carling, 1996;Rudoy, 2005).Given the undulating valley thalweg and superimposition of an esker on top of the ripple-forms, we suggest the simplest explanation is that the valley, circular incision and ripples were formed subglacially by water flowing down this tunnel valley.

Southern sector of the former Laurentide Ice Sheet
There is spatial variation in the large-scale distribution of tunnel valleys along the southern sector of the former Laurentide Ice Sheet (Fig. 4), which likely reflects their sensitivity to regional conditions such as basal thermal regime, ice and bed topography, timing and climate.For instance, the James and Des Moines lobes, which do not contain many tunnel valleys, are younger, underlain by clay-rich tills and rapidly surged to and then retreated from their maximum positions in a warmer climate (Clayton and Moran, 1982;Clayton et al., 1985).Indeed, there is a general trend of fewer tunnel valleys in the more southerly ice lobes (James, Des Moines, Lake Michigan and Huron-Erie), which may be associated with less extensive permafrost (e.g.Johnson, 1990; .In (d), there is an example of two cross cutting tunnel valleys that formed at different times during eastward ice retreat (see also Mooers, 1989).The valley trending E-W terminates at an outwash fan, which must mark the position of the ice margin when it was formed.The valley cross-cutting it can be traced several tens of kilometres further to the west and therefore must have formed during an earlier phase when ice was more extensive.Mickelson and Colgan, 2003).Moreover, the very low icesurface slopes (∼ 0.001 m km −1 ; Clark, 1992) reconstructed for the James and Des Moines ice lobes would have resulted in low subglacial hydraulic gradients not conducive to channel formation (e.g.Hewitt, 2011).Conversely, older or more northerly ice lobes with steeper ice-surface slopes, more extensive permafrost zones and sandier sediments (e.g.Chippewa, Superior and Green Bay - Wright, 1973;Attig et al., 1989;Johnson, 1990;Clark, 1992;Colgan and Mickelson, 1997;Clayton et al., 2001) have a greater occurrence and density of tunnel valleys.The locations of ice lobes along the southern sector of the Laurentide Ice Sheet are topographically controlled, mostly by upstream basins, and ice producing these lobes has been inferred to have surged or have been fast-flowing for at least part of their history (e.g.Mickelson and Colgan, 2003; Mar- showing the expected relationship between width and distance downstream for a fluvial river (Leopold and Maddock, 1953;Leopold et al., 1964) and single flood event (e.g.Lamb and Fonstad, 2010).Note that the measured tunnel valley width variations conform to neither of these expectations, but instead show variations in width greatly exceeding any possible systematic trends.gold et al., 2015).Fast ice-flow is often thought to be promoted by thermomechanical feedbacks which enhance basal meltwater production thereby lubricating the bed (cf.Winsborrow et al., 2010 and references therein).It is therefore interesting that tunnel valleys are preferentially found downglacier of depressions, where the greatest volumes of basal meltwater were likely routed.We explored the hypothesis that some of the tunnel valleys might have been fed from subglacial lakes.No clear links were found to predicted subglacial lake locations or with their obvious drainage corridors except for the Langlade Lobe tunnel valley cluster (Fig. 4) (Livingstone et al., 2013).If the predictions of lake locations of Livingstone et al. (2013) are correct, it suggests that the drainage of subglacially stored water was not the main control on tunnel valley formation.However, it is likely that the modelling underestimates the true extent of subglacial lakes (i.e. the prediction in Fig. 4 is a minimum distribution) as the predictions do not account for the possibility of water ponding behind frozen margins as suggested by Cutler et al. (2002) and Hooke and Jennings (2006).
Measurements of tunnel valley spacing reveal an overall median spacing of 4.5 km with some degree of intra-cluster regularity (Fig. 5).Inter-cluster variation is greater, with median values ranging from 1.9 to 9.1 km for individual clus-Figure 10.Scatter plot compiled to investigate if downstream variation in channel width was controlled by variations in downstream slope gradient (see text for details).That the data are centred on zero and spread fairly evenly around this demonstrates that there is no systematic relationship between elevation gradient (i.e., whether it is a normal or reverse gradient slope) and width (i.e., whether the tunnel valley is narrowing or widening).
ters across the study area.The spacing metrics are within the range of previously reported values for tunnel valleys (Praeg, 2003;Jørgensen and Sandersen, 2006;Stackebrandt, 2009;Moreau et al., 2012;Kehew et al., 2013) but smaller than the average spacing of eskers (Storrar et al., 2014a, and references therein).Theory suggests that the spacing of subglacial conduits is controlled by substrate properties, basal melt rate and the hydraulic potential gradient (e.g.Boulton et al., 2007aBoulton et al., , b, 2009;;Hewitt, 2011).According to such theory the spacing between adjacent tunnel valleys should be wider if (i) bed transmissivity is larger; (ii) melt rate and/or discharge is lower; and/or (iii) the subglacial hydraulic gradient is smaller.Thus the wider than average spacing towards the terminus of major ice lobes where ice surface slopes and thus hydraulic gradients are inferred to be shallower (e.g.Des Moines and Saginaw lobes -Clark, 1992), and a smaller spacing along narrow ice lobes characterised by steeper icesurface and hydraulic gradients (e.g. Green Bay and Superior lobes -Clark, 1992) is consistent with theory.However, cross-cutting relationships indicate that not all tunnel valleys were acting synchronously, even within a drainage cluster (Fig. 11b), which might explain the large variations in spacing.

Wider geographical distribution of tunnel valleys during the last glaciation
Figure 14 displays the geographical distribution of tunnel valleys reported in the Northern Hemisphere and attributed to the last glaciation.It appears that tunnel valleys tend to be associated with the flat southern margins of terrestrial or formerly terrestrial (e.g.North Sea) palaeo-ice sheets.They also tend to occur towards the maximum limit of glaciation and are often found downstream of large basins such as the Witch Ground in the North Sea, Baltic Depression along the south- ern limit of the European Ice Sheet, and Great Lake basins along the southern limit of the Laurentide Ice Sheet.
The tendency for tunnel valleys to form on beds of low relief and gradient near southernmost ice limits implies a genetic association.It might be that melt volumes were sufficiently high at the warm southern extremities of the ice sheets to overcome the ability of the subglacial system to export the water by other means (e.g.groundwater), making tunnel valleys more common.Perhaps it is only in low relief settings where water flow is uninhibited by topography, and can therefore organise itself into a few large catchments, in which tunnel valley forms can arise.The prevalence of tunnel valleys close behind terrestrial margins hints at an important role of permafrost in their formation (e.g.Wright, 1973;Piotrowski, 1994Piotrowski, , 1997;;Cutler et al., 2002;Jørgensen and Sandersen, 2006).It has been proposed that the development of a toe frozen to its bed along the fringe of an ice sheet acted as a barrier to water flow facilitating tunnel valley formation by subglacial ponding and outburst cycles (e.g.Wingfield, 1990;Piotrowski, 1994Piotrowski, , 1997;;Cutler et al., 2002).Moreover, freezing of sediment deposited in channels under the thin fringe of the ice sheet during winter months may have helped to prevent creep-closure of incipient tunnel valleys, thereby stabilising and preserving their forms from year to year.There is abundant evidence for well-developed permafrost conditions in the southern sector of the Laurentide Ice Sheet during and after the LGM (cf.French and Millar, 2014 and references therein), and it has been associated with glacial land systems comprising hummocky moraine, tunnel valleys and hill-hole pairs (e.g.Wright, 1973;Clayton and Moran, 1974;Bluemle and Clayton, 1984;Attig et al., 1989;Ham and Attig, 1996;Clayton et al., 1999Clayton et al., , 2001;;Colgan et al., 2003).The width of the frozen toe is likely to decrease during retreat because adjustment of the thermal structure of the toe will lag considerably behind adjustment of the margin position to an ameliorating climate.Thus, decrease in tunnel valley occurrence away from the maximum ice limit (Fig. 14) may be indicative of a change to temperate glacier conditions.

Morphology of tunnel valleys
The tunnel valleys extend for up to 65 km, although the majority (90 %) are < 17 km long and the median is 6.4 km (Fig. 6a).In comparison, reported tunnel valley lengths from the North Sea range from a few kilometres to around 100 km, with the length of individual segments not normally exceeding 20-30 km (e.g.Huuse and Lykke-Andersen, 2000).Although very wide tunnel valleys were found (maximum width ∼ 6.7 km), the majority (90 %) are 500-3000 m (Fig. 6b).This is similar to tunnel valley widths (500-5000 m) reported in Europe and elsewhere in North America (e.g.Brennand and Shaw, 1994;Huuse and Lykke-Andersen, 2000;Jørgensen and Sandersen, 2006;Kristensen et al., 2007).
Tunnel valley length and width display log-normal distributions (Fig. 6), which is common of other glacial landforms (Fowler et al., 2013;Hillier et al., 2013;Spagnolo et al., 2014;Storrar et al., 2014a).Log-normal distributions are thought to typically emerge from many independent random events in which incremental growth or fragmentation occurs (e.g.Limpert et al., 2001).For drumlins and megascale glacial lineations (MSGLs) a log-normal distribution has been used to suggest a growing phenomenon that occurs randomly, for random durations, or under random conditions (Hillier et al., 2013;Spagnolo et al., 2014), while for eskers it is thought to reflect ridge fragmentation (Storrar et al., 2014a).Examples of aligned tunnel valley segments characterised by abrupt start and end points implies at least some tunnel valley fragmentation, and this may occur due to partial burial during re-advance events or the melt out of debris-rich ice (Kehew et al., 1999).However, this fragmented appearance may also arise by differential erosion along the length of a drainage route (Fig. 8a, e, f; see also Sjogren et al., 2002) or due to water cutting up into the ice as well as down into the sediment (e.g.Fowler, 2011;Livingstone et al., 2016).In other cases, aligned tunnel valley segments could indicate a time-transgressive origin (e.g.Mooers, 1989;Patterson, 1994;Jørgensen and Sandersen, 2006;Janszen et al., 2013).This is particularly apparent where the valley segments terminate in outwash fans, and/or where segments cross-cut each other (Fig. 11 and see also Mooers, 1989).The positive relationship between tunnel valley length and width (Fig. 7) is consistent with a growing phenomenon (e.g. by headward expansion) or continuous flow (e.g. a river).In contrast, the length and width of valleys formed by floods are likely to be independent of each other; length is related to the distance that the stored water body is from the ice margin, while width is a function of the magnitude and/or frequency of drainage.
In fluvial geomorphology, channel width in an equilibrium system increases downstream (Fig. 8f) and has classically been related to discharge, and hence drainage area (Leopold and Maddock Jr., 1953;Leopold et al., 1964).This may be complicated locally by the erodibility of the bed substrate and channel slope (e.g.Finnegan et al., 2005).In contrast, large single source flood events (as may occur during a subglacial or supraglacial lake drainage event), will produce a relatively constant channel width (e.g.Lamb and Fonstad, 2010) (Fig. 8f).The downstream width of tunnel valleys in our data set varies considerably and there is no systematic downstream trend in valley form, although general increases and decreases in width do occur (Figs. 7,.Thus, there is no observable signature of catastrophic (constant width) or stable, bankfull drainage (steady widening).Moreover, the downstream variation in widths is also inconsistent with subglacial drainage channels fed by multiple supraglacial lake inputs (e.g.Palmer et al., 2011), which we would expect to produce a downstream increase in width concomitant with water input.When the rate of retreat is slow there is more time for valleys to grow and widen, whereas more rapid retreat will produce narrower segments.(c) Undulatory conduit erosion.In this theory the width of the channel eroded into sediment depends upon the competition between erosion down into the sediment (canals) vs. melting up into the ice (R-channel) (see Fowler, 2011;Livingstone et al., 2016).Note that each of the conduits (i-iii) have roughly the same area, but that in (ii) no channel forms and in (iii) the channel width is roughly half that of (i).(d) Tunnel valley formation is modulated by the basal thermal regime (modified from Hughes, 1995).Channels are able to develop more easily across warm sediment patches, and the mosaic of cold and warm sediment patches results in variations in width.

Figure 8a-e indicates that local variations in tunnel valley
width are generally more pronounced than any downstream trend.These widenings could arise from basal conditions at the time of formation (e.g.thermal regime), catastrophic drainage (e.g.Sjogren et al., 2002), or a laterally migrating stream at the base of the valley floor.Laterally migrating streams are unlikely as we do not observe terraces, bars or incised braided or meandering channels within the broader tunnel valleys, although this may partially be due to ice and post-glacial modification.The crenulated margins, circular incisions, residual hills, hummocky terrain and valley discontinuities are all analogous to features eroded during large floods by macroturbulent flow (e.g.Sjogren et al., 2002), although these are typically associated with bedrock channels (Baker, 2009 and references therein).However, we see little evidence for other characteristic features of floods, such as irregular anabranching channels (although they are observed elsewhere, e.g.Boyd, 1988;Brennand and Shaw, 1994), inner channels, furrows and large bars (e.g.Channelled Scablands: Bretz, 1923), while residual hills are not typically streamlined.
The alternative to the catastrophic hypothesis is that variations in width are strongly controlled by regional basal and hydrological conditions.There is greater similarity between tunnel valleys from the same cluster (e.g. in form, size and association with other landforms) compared to different clusters, which hints at the importance of regional conditions.Although there is no clear association with bed slope (Fig. 10) or geology (Fig. 4), the strength and therefore stability of tunnel valleys sides could have been strongly modulated by variations in basal thermal regime, substrate properties, water flow and or ice behaviour during glaciation.We therefore propose four ideas that could produce these variations in width, and which we hope will motivate physical modelling studies or field investigations (Fig. 15).Firstly, the variations in tunnel valley width may be a consequence of the very flat beds on which they form (Fig. 14).Water flow in such a landscape will be very sensitive to small changes in bed relief and variations in discharge.Coupled with sluggish water flow due to the low hydraulic gradients, we therefore envisage tunnel valleys, which display large variations in width, as a series of interconnected swampy regions (Fig. 15a).This is analogous to lakes and or swampy ground connected by overspill channels, or wide flood plains comprising dynamic river channels observed in fluvial systems flowing across similarly flat landscapes.However, not all widenings occur in bed lows (e.g.Fig. 8a, f) and so cannot account for all the variation.Secondly, tunnel valley width could relate to the rate of ice retreat, with relatively wide segments developing over longer durations when the ice is either retreating slowly or stable, and narrower segments developing during more rapid retreat (Fig. 15b).This idea is predicated on the assumption that tunnel valleys primarily form and grow close to the margin.If this were the case and width is related to time we might expect the widest segments of tunnel valleys to be associated with still-stands, and as this is not the case (e.g.Fig. 11) we therefore consider this idea unlikely.Thirdly, a basal thermal regime consisting of a mosaic  11b).(b) Reconstructed history of valley formation behind a back-stepping ice margin.Note that some valleys were long-lived during deglaciation and some abandoned shortly after their formation.The relative age relations help explain the variation in lengths between long continuous tunnel valleys and those comprising short fragments.
of cold-and warm-based sediment patches (e.g.Kleman and Glasser, 2007) would locally influence how easily widening could happen (Fig. 15d).Frozen patches would inhibit channel formation and may even result in ponding of meltwater, while warm based patches would be more susceptible to erosion.Finally, the conduit carrying water can cut down into the bed (typically forming tunnel valleys or N-channels), up into the ice (R-channels) or some mixture of the two (R and N channels) (Fig. 15c).Given that the controls on which case occurs are likely to vary over time (e.g.water discharge) and space (e.g.varying basal conditions), we propose that the large variations in tunnel valley width might be the record of how high or low the conduit was positioned in relation to the bed.Consider a conduit with an undulatory long profile cutting deeply and widely into the bed in some places and then rising back into the ice such that the cut channel in the bed narrows and pinches out and then disappears altogether where the conduit becomes entirely englacial.Control of the conduits position in relation to the bed is likely to vary with, for example, the effective pressure and relative strength of ice versus the sediment and has been explored in Fowler (2011) and Livingstone et al. (2016).

Relative timing of tunnel valley formation
Cross-cutting relationships between moraines, outwash fans, and tunnel valleys in Wisconsin have enabled their relative timing of formation to be used to build a history of formation (Figs. 11,16).If a tunnel valley cuts through moraine positions, formation must have occurred during or after the moraine was deposited.These tunnel valleys, and those interrupted by outwash fans mid-way along their length, must therefore have been used as a drainage route either repeatedly or over a long duration during retreat (see Fig. 16b).Conversely, tunnel valleys that are cross-cut by recessional moraines were abandoned as ice retreated.In Fig. 16b these tunnel valleys correspond to the age of a single moraine position, and may have been eroded during a singular "event" (i.e.outburst of a sub-or supra-glacial lake) or been abandoned due to a switch in drainage configuration or supply.

Moraines
The close link between tunnel valleys and moraines throughout the study area (Figs. 4,11; and see also Attig et al., 1989;Mooers, 1989;Patterson, 1997;Smed, 1988;Johnson, 1999;Cutler et al., 2002;Jørgensen and Sandersen, 2006) suggests formation and growth is intimately associated with pauses or slow-downs in ice margin fluctuations and that meltwater drained to the ice margin.The implication is that tunnel valley formation requires a relatively stable ice-sheet configuration to allow headward growth or recharge of source storage areas.It also provides further support for the role of permafrost in tunnel valley formation given that rapid retreat will reduce the width of the frozen toe and consequently reduce the efficacy for water storage.However, whether a reconfiguration of the subglacial hydrological regime via the development of tunnel valleys behind ice margins (moraines) can influence ice retreat (via ice dynamics), for example causing the observed staccato jumps between still-stands (Fig. 16a), remains an open question.
Regularly spaced, low relief transverse ridges in North Dakota and the Des Moines Lobe (e.g.Fig. 11b), termed washboard or corrugation moraine, have been interpreted as (annual) end moraine deposits or as subglacial crevasse fills (Kemmis et al., 1981;Stewart et al., 1988;Patterson, 1997;Jennings, 2006;Cline et al., 2015;Ankerstjerne et al., 2015).The deflection of transverse ridges towards the long axis of tunnel valleys (e.g.Fig. 11b), and buried sand and gravel deposits (see Stewart et al., 1988;Cline et al., 2015), indicates a temporal and possibly genetic relationship.One interpretation is that lower water and pore water pressures in tunnel valley and glaciofluvial deposits respectively, result in slower local ice velocities that cause the pattern of crevasses and thus ridges to be deflected (see Cline et al., 2015).However, high-pressure discharges have also been inferred from coarse-grained outwash fans deposited in front of tunnel valleys (e.g. at the edge of the Green Bay Lobe - Cutler et al., 2002, see Sect. 5.3.4).There may therefore have been multiple modes of meltwater drainage down tunnel valleys (e.g.predominantly low-pressure drainage interrupted by episodic high-pressure outbursts), or regional variation in tunnel valley evolution (e.g.some clusters formed predominantly by high energy drainage events and others in low-pressure channels).

Hill-hole pairs
The formation of tunnel valleys up-glacier from hill-hole pairs in North Dakota (Fig. 12) suggests a causative relationship.Hill-hole-pair formation is believed to require the ice to be strongly coupled to the bed so that it can exert sufficient shear stress to produce failure (Bluemle and Clayton, 1984;Aber et al., 1989).Thus, either the hill-hole pair was produced first, and the tunnel valley grew headward out of the "hole", or once drainage through the tunnel valley had waned, ice re-coupled strongly to the bed and the downstream termination of the valley became the focus of large shear stresses that resulted in failure and formation of the hill-hole pair.We suggest the former is more likely as these tunnel valleys do not terminate at moraine positions as is typical elsewhere, and small channels and eskers emanating from and diverging around the hills appear to record the down-glacier leakage of pressurised water around the obstruction (Fig. 12b).If true, these tunnel valleys appear to be unique in having initiated up-glacier from the margin.The formation of a hill-hole pair may therefore have seeded tunnel valley erosion by providing a low-point to attract water and a pathway for water through a frozen toe.This association highlights the importance of regional variations in controlling tunnel valley formation and morphology; in this case, it is the local geology (Cenozoic and Cretaceous shale and sandstone) and presence of permafrost that likely controlled the initial formation of the hill-hole pair (e.g.Bluemle and Clayton, 1984;Clayton et al., 1985) and which subsequently triggered tunnel valley growth.

Outwash Fans
Outwash fans occur at the down-glacier end of at least 10 % of the tunnel valleys in our study area (e.g.Fig. 11b).The fan sediments at the margin of the Green Bay Lobe include wellrounded pebbles and boulders up to 2 m diameter (Cutler et al., 2002), similar to accumulations documented in front of European tunnel valleys (Piotrowski, 1994;Jørgensen and Sandersen, 2006;Lesemann et al., 2014).The coarse-grained sediments have been interpreted to indicate high-energy discharges and or highly pressured subglacial meltwater flow through the tunnel valleys.Cutler et al. (2002) suggested there was at least one large outburst flood just before the termination of glaciofluvial activity through each tunnel valley in their investigation.These high-energy floods may have been responsible for cutting the valley itself, or the valley could have acted as a preferential drainage route for a flood.The concentration of outwash fans in the Chippewa, Wisconsin, Langlade, Green Bay, Superior and Wadena lobes could indicate a greater or dominant influence of discrete drainage events in these regions, while more gradual processes prevailed in other lobes (e.g.North Dakota, Des Moines, Huron-Erie).

Giant current ripples
The occurrence of giant current ripples stretching across the whole width of a tunnel valley in Minnesota (Fig. 13) implies a large sub-or supra-glacial lake outburst event (e.g.Bretz et al., 1956;Carling, 1996;Rudoy, 2005).This is further supported by the circular incision at the southern end of the tunnel valley, which is similar in form to large potholes generated by macroturbulent eddies.The flood could have cut this particular tunnel valley or the valley pre-existed and became the route of a subglacial flood further modifying and enlarging the valley (Bretz et al., 1956;Carling, 1996;Rudoy, 2005).The unique occurrence of this landform suggests that large floods were rare or such landform signatures are rarely preserved.

Implications for the formation of tunnel valleys
Based on our analysis of the morphological properties of tunnel valleys and associated landforms we provide some insights into their formation.The regular spacing of most of the tunnel valleys (Fig. 5), and that particular clusters have their own characteristic spacing, suggests self-organisation in the basal hydrological system.Individual channels somehow "know of" each other such that spacing can be set and this is most easily interpreted as there being an integrated system of many tunnel valleys operating at roughly the same time.It is difficult, for example, to understand how a collection of separate flood events could produce tunnel valleys that would combine to produce such regularity in spacing, unless the whole cluster of valleys were produced in single large flood events (as suggested for example by Shaw, 2002).
Consistent with the argument of self-organised spacing under steady flow are a series of studies that argue that the spacing may be controlled by a combination of bed transmissivity, meltwater discharge and the hydraulic potential gradient (Piotrowski, 1997;Boulton et al., 2007aBoulton et al., , b, 2009;;Hewitt, 2011).We suggest that reconstructions of drainage history, as demonstrated in Fig. 16, where we show tunnel valleys remaining active and relatively stable over long phases of ice retreat could significantly help advance knowledge of tunnel valley formation, especially when combined with information from field-based investigations and dating.
Recurrent outbursts of stored water to produce incision of whole clusters are appealing where tunnel valleys converge towards up-glacier basins (e.g.Superior and Langlade -Figs.3a, c, 4), where one could infer that subglacial lakes existed.In Evatt et al. (2006) for example they use theory of subglacial drainage to show that lakes should undergo periodic drainage and filling episodes.Perhaps many tunnel valleys are the record of large and repeated drainage events.Against this argument however, many of the clusters are very broad (> 60 km across) and the tunnel valleys relatively parallel (e.g. Green Bay and eastern Superior -Fig.4).To produce these systems would require lakes many tens or even hundreds of kilometres wide.This is difficult to reconcile with mean (< 1 km 2 ) and maximum supraglacial lake areas (up to ∼ 150 km 2 -which equates to a diameter of ∼ 14 km if a circular lake is assumed) on the surface of the presentday Greenland Ice Sheet (e.g.Leeson et al., 2013).Moreover, while very large subglacial lakes do exist beneath the Antarctic Ice Sheet (Wright and Siegert, 2011, e.g.Lake Vostok, > 250 km long by ∼ 80 km wide) and are theorised to have existed in Hudson Bay and the Great Lake Basins (e.g.Shoemaker, 1991Shoemaker, , 1999)), they have neither been predicted by modelling nor identified in the geological record (e.g.Livingstone et al., 2013).
Despite the lack of support for a mega-flood genesis of whole tunnel valley clusters, drainage of stored water down individual valleys almost certainly did happen (after Piotrowski et al., 1994;Cutler et al., 2002;Jørgensen and Sandersen, 2006;Hooke and Jennings, 2006).Not all tunnel valleys formed in clusters or were incised timetransgressively up-glacier (Figs. 4,16), and the simplest explanation for the formation of fans containing boulders (Figs. 3, 11b) (e.g.Piotrowski, 1994;Cutler et al., 2002;Derouin, 2008;Lesemann et al., 2014) and for giant currentripples (Fig. 13) is high discharge (possibly bankfull) events.Indeed, periodic higher energy or pressurised meltwater events (e.g. during penetration of surface meltwater to the bed during summer months) were probably necessary to prevent armouring of the valley sides by coarse sediment, while bedrock tunnel valleys (e.g. in Lake Superior) are difficult to reconcile solely by gradual formation.Evidence for seasonal surface meltwater reaching the bed and draining along tunnel valleys is proffered by Mooers (1989), who identified short esker segments that frequently start at moulin kames and ter-minate at outwash fans at the bed of tunnel valleys in the Superior Lobe.We therefore contend that floods from suband supra-glacial lakes, and by injections of surface meltwater down moulins did occur, contributing to the formation of some tunnel valleys either by eroding new valleys or enlarging existing ones.However, we note that most tunnel valley lengths (Fig. 6a) are an order of magnitude less than the distance up-glacier (tens to hundreds of kilometres) that supraglacial and subglacial lakes are commonly documented in Greenland and Antarctica (e.g.Selmes et al., 2011;Wright and Siegert, 2011).
We suggest that the majority of tunnel valleys along the southern sector of the Laurentide Ice Sheet were initiated at the ice margin and then typically (although not exclusively) eroded gradually up-glacier.Tunnel valley length and width display log-normal distributions and are positively correlated, indicative of a growing phenomenon (cf.Fowler et al., 2013;Hillier et al., 2013).Their strong association with moraine positions (Fig. 3) suggests that formation is time dependent (i.e. they require time to grow), while our drainage history reconstruction (Fig. 16) demonstrates that many of the features remained active for extended periods during ice margin retreat.Growth likely proceeded up-ice from the margin rather than down-ice from a stored water body because tunnel valleys preferentially terminate at ice-margin positions irrespective of their size (e.g.very small tunnel valleys along the southern margin of the Green Bay Lobe, Fig. 4).Further support is provided by some tunnel valleys in Dakota that grew headwards out of the "hole" of a hill-hole pair (Fig. 12).We suggest that when retreat is slow or a stable position is reached (allowing formation of a moraine), tunnel valleys have time to grow up-glacier and to widen and deepen as more water is discharged through them (Fig. 17a).A more unstable and/or rapid ice-retreat will limit the time for growth (headward and lateral) or may even produce a segmented tunnel valley if retreat overtakes headwards incision (Fig. 17b).
Our data indicate that the formation and morphology of tunnel valleys was strongly controlled by regional variations in basal thermal regime, bed and ice topography, timing and climate.At the broad scale, tunnel valleys tend to form on beds of low relief near southern terrestrial ice sheet margins.The paucity of tunnel valleys in the James and Des Moines lobes may be a result of the very low ice-surface slopes inhibiting channel formation and because of their relatively late advance and southerly positions that would have resulted in a less extensive zone of permafrost (Clayton and Moran, 1982;Clayton et al., 1985).Indeed, there is a general trend of fewer tunnel valleys in the more southerly ice lobes (James, Des Moines, Lake Michigan and Huron-Erie), where permafrost was reconstructed as less extensive (e.g.Johnson, 1990;Mickelson and Colgan, 2003).Regionally, we observe large inter-cluster variation in tunnel valley spacing and morphology (form and size), and their association with other glacial landforms (e.g.outwash fans, hill-hole pairs), while down-valley variations in width suggests that incision was sensitive to local conditions (e.g.Fig. 15).
Despite finding evidence for both gradual formation and high discharge events (floods) down tunnel valleys sensu lato, those that could be identified as from floods and defined as tunnel channels are not founds to be morphologically distinct from tunnel valleys sensu stricto and are therefore considered as equifinal landforms.This is unfortunate as it would have been useful to find a clear distinction.For instance, we found that the spacing, width and length of potential tunnel channels, i.e. those that terminated at an outwash fan or contained giant current ripples, or from clusters thought to have experienced large drainage events (e.g. Green Bay Lobe - Cutler et al., 2002) were similar to the overall morphology of tunnel valleys sensu lato (e.g.Fig. 6).They were also similar to tunnel valleys sensu stricto (e.g.North Dakota, where tunnel valleys grew out of hill-hole pairs).Rather, the distinction between outburst flood (tunnel channels) and gradual (tunnel valleys sensu stricto) origins in this and other studies (e.g.Cutler et al., 2002), is based on their association with other glacial landforms such as out-Figure 18.For the southern Laurentide region we consider gradual headward erosion as the usual mechanism, but with some floods down selected valleys -note the potential for stored water to cut their own valleys (e.g.supraglacial lake drainage example) or to drain along pre-existing corridors that may have tapped into a reservoir (e.g.subglacial lake example).wash fans, moraines, hill-hole pairs and giant current ripples.An important next step is to use these landform associations, where they occur, to learn more precisely about the morphological characteristics that define tunnel valleys and tunnel channels and to see if unique forms can be identified.Although we have grouped landforms of supposed different origins, the large-scale distribution and arrangement of tunnel valleys sensu lato suggests some commonality of process (e.g.Hooke and Jennings, 2006), and it may be, for example, that all tunnel valleys grow gradually, but that some experience occasional high-discharge, bankfull events.

Summary and conclusions
To provide new information on the morphological characteristics of tunnel valleys we undertook a large-scale mapping campaign to document the distribution and morphology of about 2000 tunnel valleys and associated bedforms on the bed of the former Laurentide Ice Sheet.Our maps and analyses show that tunnel valleys are semi-regularly spaced (median of 4.5 km) and tend to cluster together.The distribution of tunnel valleys varies across the study area, with clusters of tunnel valleys common across much of the Wadena, Itasca, Superior, Chippewa and Saginaw ice lobes, but much rarer along the more southerly lobes such as James, Des Moines, Lake Michigan and Huron-Erie.The wider geographical distribution suggests that tunnel valleys tend to form on flat, terrestrial beds close to the former southern LGM extent.They are typically < 20 km long and 0.5-3 km wide and longer valleys tend to be wider.The planform edges of tunnel valleys varies considerably across the study region, ranging from straight to crenulated and sharp to indistinct, and while there is no systematic downstream trend in valley form there are pronounced changes in width.There is a close link between tunnel valleys and moraines, while outwash fans occur at the down-glacier end of at least 10 % of valleys in our study area.We also observed one tunnel valley with giant current ripples on its bed, and rare cases where tunnel valleys appear to have grown out of hill-hole pairs.
There have traditionally been two explanations for the formation of tunnel valleys: (1) outburst formation by rapid drainage of sub-and/or supraglacially stored meltwater; and (2) gradual formation by headward sapping in low pressure subglacial channels (Fig. 1) (cf.Ó Cofaigh, 1996;Kehew et al., 2012;van der Vegt et al., 2012).What does our mapping and analyses say about these two hypotheses?
Given our previous work on subglacial lakes beneath the Laurentide Ice Sheet (e.g.Livingstone et al., 2013Livingstone et al., , 2016) ) we specifically explored tunnel valleys with an expectation that they might mostly be the geomorphological record of outburst floods.However, to the contrary, the morphological evidence suggests that most of the tunnel valleys underwent gradual formation, but notably with some contributions from floods from stored water (Fig. 18).In particular, our findings indicate that tunnel valleys comprise organised clusters of regularly spaced (1.9-9.1 km) valleys that formed incrementally during ice retreat.This is a strong argument for selforganising hydrological networks that mostly operated at the same time.We find that tunnel valleys preferentially terminate at moraines (irrespective of their size), which suggests that growth was initiated at and then progressed headwards from stable ice-margin positions.The concept of a growing phenomenon is further supported by log-normally distributed metrics (width, length) of valley size, the positive correlation between length and width and their initiation and growth out of hill-hole pairs.Although we favour gradual headward formation as the primary process, our results also show examples where outburst of supraglacial and or subglacial lakes have incised and/or drained down valleys.Evidence includes giant current ripples and outwash fans with large boulders (Cutler et al., 2002), and that some valleys were only occupied for brief periods during deglaciation suggestive perhaps of a short-lived event.Indeed, our reconstructed drainage history (Fig. 16) demonstrates a time-transgressive origin for many tunnel valleys, with individual clusters forming within the same time frame but individual valleys evolving over different spans involving multiple discrete flow events.Intercluster variation in tunnel valley spacing and morphology (form and size), and their association with other landforms highlights the importance of regional conditions in controlling tunnel valley formation.In particular, the presence of permafrost seems to have played a key role in determining whether tunnel valleys were produced.
Many of our observations are consistent with previous findings (e.g.Kehew et al., 2012 and references therein) and we are not the first to suggest a polygenetic origin (e.g.Hooke and Jennings, 2006).However, whilst geomorphological and sedimentological investigations in certain areas have generally advocated either an outburst or gradual genesis for tunnel valleys (Fig. 1), when their morphology, distribution and association with other glacial landforms are considered at a regional-scale it suggests that both tunnel channels and tunnel valleys sensu stricto can occur and that they appear to be equifinal (Fig. 18).

Data availability
Underlying data are available by request to Stephen J. Livingstone.

Figure 1 .
Figure 1.Cartoons depicting the two main models for tunnel valley formation: (a) outburst floods from supraglacial and/or subglacially stored water; and (b) gradual headward growth by sapping.

Figure 3 .
Figure 3. Examples of mapped valleys and the assignment of confidence levels (1 = high confidence to 3 = low confidence) along the southern sector of the former Laurentide Ice Sheet.Valleys in panels (a) Superior, (b) Saginaw, (c) North Dakota and (d) Green Bay are assigned a confidence of 1.The relict valleys contain eskers, are parallel and relatively straight, and do not trend along the regional slope.In panels (a), (c) and (d) the tunnel valley clusters terminate at a moraine position.The large valley in panel (e) (Superior) is assigned a confidence level of 3 as it does not contain any subglacial bedforms and exhibits a gradual and consistent change in bed slope consistent with a proglacial spillway.However, the smaller NW-SE valleys that it bisects is given a confidence of 2 as they have undulating thalwegs that cut across moraines.The dendritic valley cluster in panel (f) (North Dakota) is given a confidence of 3 as it is not associated with any subglacial bedforms and has a consistent bed slope indicating water flow towards the south.A braided channel morphology and a widening reach towards the south allows us to interpret this valley system as a proglacial spillway (fed by tunnel valleys emanating from under the ice to the north).

Figure 4 .
Figure 4. Distribution of mapped tunnel valleys and moraines along the southern sector of the Laurentide Ice Sheet.Likely subglacial lake locations are predictions from Livingstone et al. (2013).The Last Glacial Maximum extent is fromDyke et al. (2004) and surficial deposits are fromFullerton et al. (2003).

Figure 5 .
Figure 5. Frequency histogram of the spacing of 966 tunnel valleys from 24 discrete clusters across the southern sector of the former Laurentide Ice Sheet.

Figure 6 .
Figure 6.Frequency histogram of tunnel valley length and width (for confidence levels 1 and 2).Line is the log-normal distribution for comparison.Width values were extracted at 1 km intervals along the centre-line of each tunnel valley.

Figure 7 .
Figure7.Relationship between tunnel valley length and average width (for single thread valleys with a confidence level of 1 and 2, N = 1135).Note, there is a tendency for longer tunnel valleys to be wider.

Figure 8 .
Figure 8. Examples of tunnel valley morphology of tunnel valleys with a confidence of 1 or 2. (a) Superior Lobe (note the amphitheatre heads of some valleys); (b) Wadena Lobe (note the large downstream changes in tunnel valley width); (c) Langlade Lobe; (d) SaginawLobe; and (e) Wadena (note the parallel valleys) and (f) Huron-Erie (note the abrupt start and end points of the tunnel valleys and parallel organisation).In (d), there is an example of two cross cutting tunnel valleys that formed at different times during eastward ice retreat (see alsoMooers, 1989).The valley trending E-W terminates at an outwash fan, which must mark the position of the ice margin when it was formed.The valley cross-cutting it can be traced several tens of kilometres further to the west and therefore must have formed during an earlier phase when ice was more extensive.

Figure 9 .
Figure 9. Along-valley plots highlighting normalised (width of tunnel valley at a point/average width of the whole tunnel valley) tunnel valley width variations.(a) Saginaw Lobe; (b) North Dakota; (c) Green Bay Lobe; (d) Superior Lobe; and (e) Wadena Lobe.(f) Cartoonshowing the expected relationship between width and distance downstream for a fluvial river(Leopold and Maddock, 1953;Leopold et al., 1964) and single flood event (e.g.Lamb and Fonstad, 2010).Note that the measured tunnel valley width variations conform to neither of these expectations, but instead show variations in width greatly exceeding any possible systematic trends.

Figure 11 .
Figure 11.The varied cross-cutting associations between moraines, outwash fans and tunnel valleys in: (a) North Dakota -note how the washboard moraines curve up-glacier towards the tunnel valleys; and (b) Wisconsin (Chippewa Lobe).

Figure 12 .
Figure 12.Hill-hole pairs in North Dakota and their association with tunnel valleys.(a) Note the esker downstream of the hill, which trends into an aligned tunnel valley segment.(b) Note the secondary meltwater channels and eskers that diverge around the hill.

Figure 13 .
Figure13.Giant Current Ripples spanning the width of a shallow tunnel valley that is cut into an obliquely-oriented drumlin field (water flow to the south).These sinusoidal bedforms are interpreted as giant current ripples, which formed during a large subglacial flood.Note the undulating thalwegs and esker in the valley that indicates subglacial deposition, and the circular incision (with a remnant island in its centre) in the south of the valley that may have been formed by a large eddy during high-energy turbulent flow.

Figure 14 .
Figure 14.Currently known Northern Hemisphere distribution of tunnel valleys that have been attributed to the last glaciation.The opaque blue shading is the Last Glacial Maximum ice sheet distribution.Black lines are the mapped tunnel valleys from Fig. 3 and black boxes are where tunnel valleys have been identified.The preference for tunnel valleys to occur in flat areas close to southern limits of the ice sheets is striking, but to what extent does this map capture the true extent of tunnel valleys, perhaps they are selectively unreported from other regions?

Figure 15 .
Figure 15.Cartoons showing four theories to explain the downstream variation in tunnel valley width.(a) Swampy ground (blue stipples)and overspill channels (blue lines) associated with water flow across very flat ground.In such a flat landscape tunnel valleys are able to easily expand laterally, in response to small changes in water flux, and there is little impetus for rapid vertical erosion due to shallow hydraulic gradients.(b) Variable retreat rate.When the rate of retreat is slow there is more time for valleys to grow and widen, whereas more rapid retreat will produce narrower segments.(c) Undulatory conduit erosion.In this theory the width of the channel eroded into sediment depends upon the competition between erosion down into the sediment (canals) vs. melting up into the ice (R-channel) (seeFowler, 2011;Livingstone et al., 2016).Note that each of the conduits (i-iii) have roughly the same area, but that in (ii) no channel forms and in (iii) the channel width is roughly half that of (i).(d) Tunnel valley formation is modulated by the basal thermal regime (modified fromHughes, 1995).Channels are able to develop more easily across warm sediment patches, and the mosaic of cold and warm sediment patches results in variations in width.

Figure 16 .
Figure 16.Using cross-cutting relationships to reconstruct tunnel valley evolution during ice margin retreat.(a) Mapping of tunnel valleys and associated glacial bedforms in Wisconsin (Chippewa Lobe) (from Fig.11b).(b) Reconstructed history of valley formation behind a back-stepping ice margin.Note that some valleys were long-lived during deglaciation and some abandoned shortly after their formation.The relative age relations help explain the variation in lengths between long continuous tunnel valleys and those comprising short fragments.

Figure 17 .
Figure 17.Cartoon demonstrating the dependence of tunnel valley evolution (by headward growth) on ice margin retreat rate.(a) If headward growth of a tunnel valley is faster than the rate of ice retreat the valley will be able to extend continuously up-glacier and its length will only be limited by water supply and hydraulic properties of the bed.(b) If, however, headward growth of a tunnel valley is slower than the rate of ice retreat the valley is likely to be discontinuous, only being able to form and extend up-ice during slow-downs or pauses in retreat.