Increased channelization of subglacial drainage during deglaciation of the Laurentide Ice Sheet

Recent estimates have shown that mass loss from the Greenland Ice Sheet and West Antarctic Ice Sheet is accelerating (Rignot et al., 2011). This has been attributed to surface melt and dynamic losses caused by the acceleration, thinning, and retreat of outlet glaciers (van den Broeke et al., 2009). These two processes are not mutually exclusive in that some workers have suggested that surface meltwater may drain to the bed and cause dynamic flow acceleration (Zwally et al., 2002; Parizek and Alley, 2004). However, others have argued that an increased melt flux may be counteracted by the evolution of the drainage system, which responds by generating large channels that efficiently evacuate excess surface meltwater (Bartholomew et al., 2010; Schoof, 2010; Sundal et al., 2011; Shannon et al., 2013), similar to the seasonal evolution inferred beneath alpine glaciers (Hubbard and Nienow, 1997; Iken and Truffer, 1997). These potential switches in subglacial drainage introduce uncertainty when attempting to predict the future behavior of ice sheets: Could increased surface melt lead to more widespread dynamic instabilities, or will subglacial drainage systems become more efficient? Moreover, observations of subglacial drainage systems are difficult, often indirect, and generally limited to just a few years (Zwally et al., 2002; Bartholomew et al., 2010; Sundal et al., 2011; Cowton et al., 2012). Thus, our understanding of how drainage systems might evolve over time scales of centuries to millennia, the impact this may have on ice dynamics, and our ability to model future ice sheet behavior remains limited.

One approach to understanding the evolution of subglacial drainage systems at centennial to millennial time scales is to examine the geological record of former ice sheets that have subsequently deglaciated. In particular, evidence of channelized drainage beneath ice sheets is typically preserved in the form of eskers: elongate ridges of glaciofluvial sand and gravel deposited in subglacial, englacial, or supraglacial drainage channels (Brennand, 2000). Here, we use eskers to investigate the evolution of channelized drainage beneath the Laurentide Ice Sheet (LIS) during 6000 yr of deglaciation from 13 to 7 kyr B.P. We quantify the spatial pattern of over 20,000 eskers and show that the ice sheet was progressively drained by an increasing number of subglacial channels as it deglaciated. Our hypothesis is that the subglacial drainage system evolved to cope with an increasing amount of meltwater during deglaciation (Carlson et al., 2008, 2009), in a manner that is reminiscent of the seasonal evolution of drainage under alpine glaciers, but which occurred over much longer time scales. An important implication is that the area of the ice sheet bed subjected to a less efficient “distributed” system would have decreased, perhaps limiting the potential for flow instabilities (surges or ice streaming) (Kamb, 1987). This suggests that the final demise of the ice sheet was intimately and predictably linked to its surface mass balance.

Eskers were mapped from Landsat Enhanced Thematic Mapper Plus (ETM+) imagery of Canada (Figs. 1A and 1B), and the data are published as a map in Storrar et al. (2013), the only quantitative map of Canadian eskers from a consistent data source (cf. Aylsworth et al., 2012). They were identified using the criteria of Margold and Jansson (2012) and are conspicuous as narrow, sinuous ridges with distinct (highly reflective) spectral signatures related to the presence of glaciofluvial sediments. Their ridge crest lines were mapped as line features directly into a GIS. Given the ice-sheet scale of the mapping, smaller eskers (i.e., below the image resolution) and associated landforms were not included in the analysis but, given that the focus is on the evolution of large-scale patterns, this is unlikely to impact on the results. Esker ridge lengths, spacing, and number of tributaries were extracted from the GIS, and we analyzed their location in relation to geomorphologically and chronometrically constrained ice margin positions during deglaciation. Specifically, we calculated the density of eskers (number of eskers per 100 km of ice margin) that intersected ice margin positions from the best available ice margin chronology, published by Dyke et al. (2003). The chronology is based on ∼4000 dates, most of which are from radiocarbon dating, which has an associated uncertainty not shown here (see Dyke et al., 2003). We also added intermediate margin positions by interpolating between dated margins, to enable additional measurements. Where interpolations were made, the age uncertainties are shown by error bars (e.g., Fig. 2). It is important to note, however, that our focus is on the relative changes in eskers through time as the ice margin retreated, such that absolute ages (and their uncertainties) are less important. Tributaries were identified where two eskers meet and confluences were then counted between ice margins. We assume that the mapped eskers reflect the location of the majority of large meltwater channels beneath the ice sheet. The eskers in our database exhibit lengths of up to several hundred kilometers, and, because there is little evidence that eskers of such great length form synchronously, we assume that they formed time-transgressively at a retreating ice margin (e.g., Aylsworth and Shilts, 1989; Clark and Walder, 1994), although some degree of synchronous deposition is possible over short (1–10 km) length scales (cf. Brennand, 2000). Thus, long eskers may intersect more than one margin.

To avoid the complexity associated with retreat from a predominantly “soft” sedimentary bed (where eskers are rare) to the relatively “hard” crystalline bedrock of the Canadian Shield (where eskers are far more common) (Clark and Walder, 1994; Storrar et al., 2013), we focus on deglaciation across the Canadian Shield, from ca. 12.7 kyr B.P. to ca. 7.45 kyr B.P. To account for the diminishing size of the ice sheet (which would theoretically produce fewer eskers than a larger ice sheet), the number of eskers at the ice sheet margin was normalized per 100 km of margin length. Mean ice sheet retreat rate was calculated from 20 transects from the youngest to oldest dated margins in each area (Fig. 1C). Standard deviation values are given as error bars in Figure 2.

Our database includes a total of 20,186 eskers, and individual ridges are up to 97.5 km in length (mean length of 3.5 km). Systems of eskers can be traced for up to 760 km, if small gaps (less than a few kilometers) created by postglacial modification (or breaks in deposition) are taken into account. Some appear as single ridges, whereas others possess up to fourth-order streams (using the Strahler method). The mean number of eskers per 100 km of ice margin is 1.24, but this ranges from 0.05 to 3.31. The most obvious pattern is their abundance in branching dendritic systems on the Canadian Shield, which emanate away from the positions of the major ice divides in Keewatin and Labrador (Fig. 1). We note, however, that eskers are conspicuously absent beneath the final location of these ice divides.

Retreat rates in Keewatin were generally between 100 and 200 m yr^-1 from 13 to 9.5 kyr B.P., but increased rapidly between 9.5 and 9 kyr B.P., followed by a sharp decrease to between 8.5 and 8 kyr B.P. (see Fig. 2). The density of eskers at the ice sheet margin broadly matches the retreat rate in both sectors, increasing from 0.7 to 2.9 eskers per 100 km of ice margin between 12 and 9 kyr B.P. in Keewatin, then abruptly decreasing from 2.9 to 0.8 between 9 and 8.5 kyr B.P., after which it decreased more gradually. In contrast, the number of tributaries per 100 km of eskers decreased steadily from 25.8 to 0.1 from 13 to 7 kyr B.P. To the east, the Labrador sector of the ice sheet shows a similar pattern, esker density increasing from 0.1 to 1.8 eskers per 100 km of ice margin between 13 and 7.6 kyr B.P. and then decreasing to 0.7 by 7 kyr B.P. (Fig. 2B). The number of tributaries per 100 km of eskers decreased from 1.9 to 0.3 between 11.4 and 7 kyr B.P.

The most obvious pattern from our data is a clear increase in esker density during deglaciation, which we summarize in Figure 3. Esker formation has been linked to variations in sediment supply (Aylsworth and Shilts, 1989) and substrate lithological controls (Clark and Walder, 1994), but we rule these out due to the relatively uniform lithology of the study area (Wheeler et al., 1996). It has also been shown that esker sediments are derived primarily from the local substrate (Bolduc, 1992), with no obvious variation in lithology in this region that might explain the observed increase in eskers during deglaciation and their absence beneath ice divides (Figs. 1 and 3). The decreasing number of tributaries during deglaciation indicates that the increase in esker density is not simply due to increased branching upstream.

Rather, we favor a simpler explanation, which is that the patterns are connected to changes in the volume of meltwater during deglaciation. Specifically, the Northern Hemisphere temperature record shows pronounced warming from 12.5 kyr B.P., following the Younger Dryas, to ca. 9 kyr B.P. (Fig. 2). This warming would have led to a negative mass balance, with an associated increase in the production of surface meltwater (Carlson et al., 2009). It may also be manifest as an increase in the rate of ice margin retreat, which does indeed occur, e.g., from 100–200 m yr^-1 to almost 400 m yr^-1 from 11 to 9 kyr B.P. in Keewatin, and from 0–100 to over 400 m yr^-1 from 8.5 to 7.5 kyr B.P. in Labrador (though we note a lag between the peak esker density and ice margin retreat rates in the Keewatin sector). Thus, we suggest that the increase in esker density during deglaciation is a reflection of an increased volume of surface meltwater entering the subglacial system via moulins and englacial channels, similar to that observed in the ablation zone of the Greenland Ice Sheet (Catania and Neumann, 2010) and earlier hypothesized to explain some esker patterns by Brennand (2000). This explanation is consistent with numerical modeling experiments, which predict a decrease in esker spacing (i.e., an increase in esker density) with increasing meltwater supply (Boulton et al., 2009; Hewitt, 2011).

We note that the increase in esker density in Labrador lags the Keewatin sector by ∼1000 yr (Fig. 2), which is consistent with the stabilization of this sector of the ice sheet near the Gulf of St. Lawrence in the early Holocene (Hillaire-Marcel and Occhietti, 1980). Re-equilibration moraines in the St. Lawrence lowlands indicate that margin retreat paused as it switched from being marine- to land-terminating and its mass was redistributed (Hillaire-Marcel and Occhietti, 1980). Following this switch, and the redistribution of mass, the ice margin continued its retreat. This is further supported by geomorphological evidence from Clark et al. (2000), who suggested that the ice sheet in this sector was strongly asymmetric, with the southern sector retreating rapidly compared with the northern sector.

In view of the above, the abrupt decrease in esker density after 8.5–9 kyr B.P. (7.5 kyr B.P. in Labrador) and their complete absence from the location of former ice divides requires explanation because the ice mass was still >200 km from ice divide to margin. We speculate that the LIS was sufficiently thin by this point that much of the remaining ice would have been cold and mostly cold based (Kleman and Hättestrand, 1999; Kleman and Glasser, 2007), potentially precluding the development of subglacial meltwater channels and forcing meltwater to drain from the surface without depositing eskers. In addition, while sediment supply is unlikely to explain the large-scale pattern of eskers, their absence beneath the smaller ice caps closer to final deglaciation may be related to limited sediment availability, in that low subglacial erosion rates are likely to have persisted under the ice divides. Furthermore, if a minimum length of conduit is required to erode the bed (for example, Cowton et al., [2012] documented erosion from a conduit >50 km in length beneath the Leverett glacier, Greenland) and provide sufficient sediment to build eskers, this may also explain the absence of eskers toward the center of ice dispersal.

A key implication of our work is that the final demise of the LIS may have been driven largely by surface melt. Independent support for this is found in numerical modeling of LIS mass balance, which is thought to have been strongly negative during the early Holocene, attributed to enhanced boreal summer insolation (Carlson et al., 2008, 2009). Dynamic mass loss (e.g., through ice streaming or surging) is typically dependent on distributed drainage of meltwater at the bed (Clarke, 1987) and so an important corollary of our findings is that dynamic mass loss may have been less significant during final deglaciation. This is supported by the relative paucity of ice streams hypothesized to have been operational at this time (Kleman and Glasser, 2007; Stokes and Tarasov, 2010), although we note that few paleo–ice streams have been dated with any confidence.

Using the pattern and spacing of >20,000 eskers formed under the LIS, alongside an established ice margin chronology (with intervening interpolated ice margins), we demonstrate that the meltwater drainage systems of ice sheets appear to evolve in response to changes in meltwater input over thousands of years, reminiscent of seasonal changes observed in valley glaciers (Hubbard and Nienow, 1997; Iken and Truffer, 1997). We suggest that increased surface melting in a warming climate following the Younger Dryas led to the formation of additional subglacial drainage channels beneath the LIS, which increased the proportion of the bed drained by efficient drainage networks. As a result, dynamic instabilities, such as ice streaming and surging, may have been less important in contributing to its eventual demise (Carlson et al., 2008, 2009). These findings have implications for numerical models of ice sheets that incorporate basal hydrology. Ideally, the models should account for changes in the structure of meltwater drainage systems in response to increased meltwater supply, and the associated impact on ice sheet dynamics.

We are grateful to Anders Carlson and two anonymous reviewers, whose comments improved the clarity of the manuscript. This work was funded by a Natural Environment Research Council (NERC) Ph.D. studentship awarded to Storrar.