Glacial geomorphology of the northern Kivalliq region, Nunavut, Canada, with an emphasis on meltwater drainage systems

ABSTRACT This paper presents a glacial geomorphological map of glacial lineations, ribbed terrain, moraines, meltwater channels (subglacial and ice-marginal/proglacial), eskers, glaciofluvial deposits, ice-contact outwash fans and deltas and abandoned shorelines on the bed of the former Laurentide Ice Sheet in northern Canada. Mapping was compiled from satellite imagery and digital elevation data and landforms were digitised directly into a Geographical Information System. The map reveals a complex glacial history characterised by multiple ice-flow events, including fast-flowing ice streams. Moraines record a series of pauses or re-advances during overall SE retreat towards the Keewatin Ice Divide. The distribution of subglacial meltwater landforms indicates that several distinctive scales and modes of drainage system operated beneath the retreating ice sheet. This includes a large (>100 km) integrated network of meltwater channels, eskers, ice-contact outwash fans and deltas and glaciofluvial deposits that originates at the northern edge of Aberdeen Lake. The map comprises zone 66 of the Canadian National Topographic System, which encompasses an area of 160,000 km2. It is presented at a scale of 1:500,000 and is designed to be printed at A0 size.


Introduction
The Laurentide Ice Sheet covered a substantial proportion of North America during the Late Wisconsinan glaciation (Dyke et al., 2002) and its subsequent demise was closely coupled to dramatic climatic changes (e.g. Barber et al., 1999;Clark, 1994;Clark, Alley, & Pollard, 1999). Understanding the fluctuations of ice sheets and the processes which take place beneath them, and how they interact with climate, is important for predicting how ice sheets will respond to contemporary climate change. High resolution climatic records extending over the last glacial cycle (and beyond) are widely available (e.g. Shakun et al., 2012), providing a wealth of data which can be compared with the palaeo-ice sheet record. Previous reconstructions have focused principally on ice-margin fluctuations (e.g. Dyke & Prest, 1987;Dyke, Moore, & Robertson, 2003) and the activity of ice streams (e.g. De Angelis & Kleman, 2005;Margold, Stokes, & Clark, 2015;Ó Cofaigh, Evans, & Smith, 2010;Ross, Campbell, Parent, & Adams, 2009;Stokes & Clark, 2001;Stokes, Clark, & Storrar, 2009;Stokes, Margold, Clark, & Tarasov, 2016), which rapidly drain large portions of ice sheets.
Recent work on contemporary ice sheets has drawn attention to the importance of meltwater in controlling ice dynamics (e.g. Bartholomew et al., 2010;Schoof, 2010;Zwally et al., 2002). However, whilst meltwater was abundant during the deglaciation of the Laurentide Ice Sheet (e.g. Carlson et al., 2009;Storrar, Stokes, & Evans, 2014a) and produced a large geomorphological imprint (e.g. Brennand, 2000;Brennand & Shaw, 1994;Mullins & Hinchey, 1989;Prest, Grant, & Rampton, 1968;Storrar, Stokes, & Evans, 2013), the geomorphological record of meltwater has to date played only a relatively minor role in ice sheet reconstructions. This is typically limited to providing supplementary information about ice geometries and flow direction (Kleman & Borgström, 1996), rather than to reconstruct basal and hydrological processes in space and time (see review by Greenwood, Clason, Helanow, & Margold, 2016). The long-term effect is that the role of meltwater in governing or modulating ice sheet dynamics is relatively underexplored (although see Stokes, Tarasov, & Dyke, 2012;Tarasov & Peltier, 2004. This paper aims to provide the geomorphological mapping basis for a first attempt to incorporate meltwater landforms into a larger ice dynamic reconstruction, using a study area in northern Canada. meridians, encompassing an area of approximately 160,000 km 2 (Figure 1). Despite its large size, the study area is predominantly within the relatively flat Back Lowland and Thelon Plain (Dyke & Dredge, 1989), with a maximum elevation of 361 m a.s.l. Queen Maud Gulf lies to the north and the area contains myriad small lakes and several larger lakes, including Aberdeen Lake (1106 km 2 ), Schultz Lake (410 km 2 ) and the Garry Lakes (776 km 2 ). The area lies on the Precambrian Shield and is composed mainly of granitioid rocks, with belts of greenstone and metamorphosed sediments in the north and sedimentary basins in the south (Wheeler et al., 1996).
Early mapping of areas like, and including, zone 66 (e.g. Craig, 1961Craig, , 1964 by the Geological Survey of Canada (GSC) resulted in the publication of the Glacial Map of Canada (Prest et al., 1968), which shows the generalised glacial geomorphology. Since the 1980s, the GSC have produced several more detailed surficial geology map sheets covering subsets of the area at scales ranging from 1:100,000 to 1:1,000,000, the spatial distribution of which is shown in Figure 2 (Aylsworth, 1990;Aylsworth & Clarke, 1989;Aylsworth & Shilts, 1989a;Aylsworth, Cunningham, & Shilts, 1990;Helie, 1984;McMartin, Dredge, & Aylsworth, 2008;McMartin, Dredge, & Robertson, 2005;St-Onge & Kerr, 2013Thomas, 1981aThomas, , 1981b. This mapping is based on the interpretation of aerial photographs and field observations and contains different collections of landforms, depending on the mapper. A summary of the landforms recorded on each map is provided in Table 1. Whilst some landforms, such as eskers, are presented in a high level of detail, others, notably glacial lineations, are necessarily generalised and individual bedforms are often not mapped. Moreover, Figure 2 shows that several portions of the study   area are restricted to mapping at small scales of 1:500,000 or 1:1,000,000. De Angelis (2007) produced a map of the glacial geomorphology of the east-central Canadian Arctic, which is based predominantly on the interpretation of Landsat Enhanced Thematic Mapper (ETM)+ imagery. This map portrays glacial landforms required for a glacial inversion model (Kleman et al., 2007). Like the Glacial Map of Canada (Prest et al., 1968), some features (e.g. glacial lineations) are displayed in a generalised pattern, rather than each individual landform being mapped. The focus of previous mapping has been predominantly on ice sheet dynamics, with the result that important meltwater features, including meltwater channels, have not previously been mapped in this region. Further mapping is therefore required to provide a higher level of detail for reconstructions, and to interpret the meltwater history of this part of the Laurentide Ice Sheet as it deglaciated.

A brief glacial history
The study area is likely to have been ice covered for much of the last glacial cycle (Stokes et al., 2012) and the landform record attests to three distinctive signatures of ice sheet history. First, prior to deglaciation the study area was located in a core region of the Keewatin sector of the Laurentide Ice Sheet, and acted as a centre of ice dispersal (Dyke et al., 2002;McMartin & Henderson, 2004). A series of ice divides developed over the area, which reconfigured over time, having a dramatic influence on ice-flow direction as recorded by cross-cutting glacial lineations and striae (Boulton & Clark, 1990;Hodder, Ross, & Menzies, 2016;McMartin & Henderson, 2004).
Second, the Dubawnt Lake Ice Stream flowed across the southern part of the area towards the NW/WNW, during deglaciation from approximately 10.2 to 9 ka, producing numerous mega-scale glacial lineations (MSGLs) (Stokes & Clark, 2003a, 2003b. It has also been suggested that, prior to deglaciation, another ice stream operated in the northern and central parts of the area, which flowed to the north-east and fed the Gulf of Boothia Ice Stream (Hodder et al., 2016).
Third, the final deglaciation of the area occurred between ca. 10 and 8 ka (Dyke et al., 2003), as ice retreated to the east-south-east towards the final location of the Keewatin Ice Divide (McMartin & Henderson, 2004). The deglacial signature is characterised by prominent moraines and eskers (Aylsworth & Shilts, 1989a, 1989bPrest et al., 1968).

Methods and data sources
Glacial landforms were identified from satellite imagery and digital elevation data and were digitised directly into a Geographical Information System (GIS). Landsat 7 ETM+ and 8 OLI imagery has a spatial resolution of 30 m (and 15 m in the panchromatic band) and was downloaded from the USGS (available at: earthexplorer.usgs.gov). Individual images were viewed at a variety of scales, typically between 1:50,000 and 1:500,000, as false colour composites in different band combinations: principally 4,3,2 and 7,5,2 (R,G,B). Landforms were digitised at approximately 1:50,000. The 1:50,000 scale Canadian Digital Elevation Dataset (CDED: available from geogratis.gc.ca) was used to produce a Digital Elevation Model (DEM) with a horizontal resolution of 15.5 m. In order to avoid azimuth biasing (Smith & Clark, 2005;Smith & Wise, 2007), the DEM was visualised in three different ways: (i) as two relief-shaded models with different azimuths (315°and 45°); (ii) as a slope gradient-shaded model and (iii) as a simple elevation-shaded model.
Although both sets of imagery were used to identify landforms, some features were more readily identifiable in one particular set. For example, low relief corridors of glaciofluvial deposits were identifiable only by their distinctive spectral signature in the Landsat ETM+ imagery. In contrast, meltwater channels are characterised by a topographic signature and were typically easier to identify in the CDED data. All mapping was cross-checked with previous investigations (see Table 1) to ensure consistency.
Glacial lineations, ribbed terrain, esker crestlines, meltwater channel thalwegs, abandoned shorelines and moraine crestlines were digitised as polylines. Glaciofluvial deposits (deposits of sand and/or gravel characterised by high reflectance in Landsat imagery) and ice-contact outwash fans and deltas were digitised as polygons. Lake outlines used on the map were downloaded from the North American Atlas (available from: nationalatlas.gov).
Most of the above landforms are clearly distinguishable in different sets of remotely sensed imagery. However, meltwater channels (erosional products of meltwater flow that may range from metre to kilometre scale) are often ambiguous. They may form in subglacial, proglacial or ice-marginal positions and often occupy topographic lows, leading to potential confusion with modern or postglacial drainage. Indeed, postglacial drainage patterns on the Canadian Shield are characteristically deranged as a result of the effect of repeated glaciations on reshaping the surface, with insufficient time for graded fluvial systems to develop (Dyke & Dredge, 1989). Many modern channels therefore occupy former meltwater channels, making them difficult to objectively separate. All potential meltwater channels were therefore mapped and individually assessed to determine whether they formed subglacially, ice-marginally or proglacially (after Greenwood, Clark, & Hughes, 2007). The criteria for distinguishing subglacial meltwater channels include: (i) an association with other glacial landforms such as moraines and eskers; (ii) an undulating channel thalweg and (iii) an oblique channel orientation relative to the regional drainage slope. Proglacial meltwater channels flow down slope, often start abruptly at former icemargin positions and meander downstream. Ice-marginal meltwater channels document the ice margin at the time of their formation and are classically related to parallel flights of lateral channels perched on valley sides (e.g. Dyke, 1993;Greenwood et al., 2007;Maag, 1969;Margold, Jansson, Kleman, & Stroeven, 2011).

Mapped landforms
A new map of the glacial geomorphology of NTS zone 66 is presented in the supplementary data (Main Map). The mapping reveals a rich and complex suite of glacial landforms and provides a consistent level of detail across the study area, in places exceeding that of previous work (Figure 2). Eight different types of landform were mapped and these are briefly described below.
The map contains 28,061 glacial lineations (Figure 3), which range in size from relatively small drumlins (Figure 3(A,B)) a few 10s to 100s of m long, to spectacular MSGLs (Clark, 1993) of the Dubawnt Lake Ice Stream bed (Stokes & Clark, 2003a) that reach up to 23 km long (Figure 3(C,D)). Glacial lineations reflect different ice-flow phases, with principal flows aligned E-W (in the southern part) and S-N (in the east). Smaller lineation clusters depict a complex arrangement of orientations, including cross-cutting, suggesting that ice flow changed direction during the last glacial. Some of the large (10s km long) S-N trending lineations in the south of the study area have previously been identified as 'mega-scale transverse bedforms' akin to giant ribbed terrain (Greenwood & Kleman, 2010). These landforms have more recently been interpreted by Hodder et al. (2016) as relict MSGL. This interpretation is supported by parallel striations, which would be difficult to explain if the features formed transverse to ice flow. For this reason we follow Hodder et al. (2016) and classify these features as glacial lineations, though we note some uncertainty in this interpretation.
Moraine ridges ( Figure 5) demarcate former icemargin positions. Most of the moraine ridges that were mapped are consistent with deglacial ice margins reflecting the gradual retreat of the ice sheet from the NW to the SE (Dyke et al., 2003;Dyke & Prest, 1987). In some locations, several moraine ridges are aligned in close proximity, whereas elsewhere they occur as single isolated ridges.
Eskers ( Figure 6) are straight-to-sinuous ridges composed of glaciofluvial sand and gravel, which form by deposition of sediment predominantly in R-channels, which are aligned roughly normal to the ice margin (Brennand, 2000). Numerous eskers were mapped in the area and they record the meltwater drainage pattern and ice-margin positions during final deglaciation over the area, approximately from the NNW to SSE towards the final location of the Keewatin Ice Divide   (Dyke & Prest, 1987) shown in (A) CDED and (B) Landsat ETM+ (R,G,B 4,3,2) imagery. Location is given in Figure 1. (Aylsworth & Shilts, 1989b). Most eskers are up to 10 km long, although some reach lengths up to 46 km. Eskers are quasi-regularly spaced approximately 12 km apart across the study area (Storrar, Stokes, & Evans, 2014b).
Meltwater channels (Figure 7) are up to 102 km long and may be incised into either bedrock or sediment. Of the 404 'definite' and 'probable' channels mapped in the study area, 306 were interpreted as being subglacial in origin (76% of the total population) and 98 either proglacial or ice-marginal channels (24% of the total population). Ice-marginal meltwater channels are up to 30 km long and tend to be orientated transverse to glacial lineations and subglacial meltwater products in association with moraine still-stand positions. The majority of subglacial meltwater channels occur between Aberdeen and Garry lakes along a roughly SE-NW to S-N axis. The channels are up to 1.5 km wide and extend for long distances (up to 100 km), consistent with the dimensions of tunnel valleys (e.g. Kehew, Piotrowski, & Jørgensen, 2012;Livingstone & Clark, 2016;Ó Cofaigh, 1996;Wright, 1973). Many appear to initiate on the northern edge of Aberdeen Lake, whilst a smaller number originate further south and trend through it. Subglacial meltwater channels are often associated with eskers, which either occur within the channels or form a continuation with them, producing integrated drainage networks. In contrast to other locations (e.g. Livingstone & Clark, 2016), about half of the subglacial meltwater channels throughout the study area do not terminate at moraine positions, despite the prevalence of moraines in the area. Some meltwater channels are aligned conformably with lineation trends, whereas in some places meltwater channels cut obliquely to the surrounding lineations. North of the Garry lakes, there is a series of long (up to 100 km) relatively straight channels that may have formed subglacially. However, they typically trend down normal slopes and are not associated with other glacial landforms, like eskers and moraines, and therefore have been classified as 'possible' subglacial channels and omitted from the map.
Ice-contact outwash fans and deltas (Figure 8) comprise approximately triangular mounds and flat-topped hills up to a few kilometres in diameter, composed of glaciofluvial sand and gravel. They are typically found at the terminus of esker segments and/or at moraine positions (see also Helie, 1984;St-Onge & Kerr, 2013. In contrast, the meltwater channels are rarely associated with fans and deltas.  Glaciofluvial deposits (Figure 9) form wide (up to 6 km) corridor-like tracts of glaciofluvial sand and gravel, often located within erosional corridors (cf. Rampton, 2000;Utting, Ward, & Little, 2009). These deposits are commonly associated with aligned eskers, which may suggest that they form in subglacial positions (St-Onge, 1984). They occur exclusively north of Aberdeen Lake and are particularly extensive in the area between Aberdeen and Garry Lakes, where they are aligned with a relatively strong N-S orientation.
Abandoned shorelines (Figure 10) are present in a restricted band between approximately 80 and 280 m a.s.l., surrounding Aberdeen Lake and Shultz Lake (Aylsworth, 1990). Raised shorelines are also present around the Garry lakes, a few metres above present lake levels. They reflect a series of higher lake levels and subsequent marine incursion (Dyke et al., 2003;McMartin et al., 2008;Prest et al., 1968). Lakes likely formed in proglacial positions, dammed by ice to the south and east.

Discussion and implications
This paper presents a map of the glacial geomorphology of zone 66 of the NTS system of Canada. Mapping from Landsat imagery and CDED elevation data reveal a large number of glacial landforms, representing a complex glacial history. Ubiquitous glacial lineations include the well-documented Dubawnt Lake ice stream bed in the southern part of the area, which records rapid ice flow roughly from east to west (Stokes & Clark, 2003a, 2003b. A further major flow is recorded along a south-north axis across much of the eastern part of the area, which may have fed the Boothia Ice Stream to the north-east (see also Hodder et al., 2016;Kleman et al., 2007;Margold et al., 2015). In addition to these large flow events, the landform record also reveals several smaller events which, in some locations, are recorded by complex clusters of lineations cross-cutting in up to three different directions.
Meltwater landforms, including eskers, meltwater channels, ice-contact outwash fans and deltas and glaciofluvial deposits, record the hydrological component of the ice sheet during deglaciation as it retreated from the north-west towards the final location of the Keewatin Ice Divide to the east-south-east of the area (McMartin & Henderson, 2004). We identified several different arrangements, scales and types of meltwater landform, which suggests that the mode of subglacial drainage varied in space and time. In particular, tracts of glaciofluvial deposits and larger meltwater channels (tunnel valleys) form a distinctive spatial pattern between Aberdeen Lake and the Garry Lakes, and this indicates the intriguing possibility that they may reflect the transportation of meltwater between the two lakes. It is possible that the depressions in which these lakes now sit were occupied by subglacial lakes (Livingstone, Clark, & Tarasov, 2013) prior to the development of proglacial lakes and then marine inundation during deglaciation (cf. Dyke, 2004;Dyke et al., 2003;Prest et al., 1968;Stokes & Clark, 2004). Further analysis is required to determine the precise nature of meltwater drainage in this area but this highlights the new insights that can be gleaned by incorporating meltwater dynamics into palaeo-ice sheet reconstructions.
As well as providing a detailed database that may be used for ice-margin and ice dynamic reconstructions, extensive glaciofluvial landforms, including landforms not previously mapped detail such as corridors of glaciofluvial deposits and meltwater channels, permit the evolving hydrology of the Laurentide Ice Sheet to be reconstructed. It is therefore anticipated that the data contained within this map will be used to produce a detailed reconstruction of (i) ice-margin retreat; (ii) evolving ice dynamics and (iii) evolving meltwater systems of this sector of the Laurentide Ice Sheet during its deglaciation.
Software Esri ArcGIS 10.3 was used to produce the map and figures. manuscript and map. We thank John Abraham for his cartographic review which improved the presentation of the map.

Disclosure statement
No potential conflict of interest was reported by the authors.