Glacial geomorphology of the central sector of the Cordilleran Ice Sheet, Northern British Columbia, Canada

ABSTRACT Northern British Columbia was repeatedly covered by the Cordilleran Ice Sheet (CIS) during the glacial periods. However, its mountainous terrain and remote location have thus far impeded our understanding of the central sector of the ice sheet. The improved resolution of remotely sensed data provides new opportunities to unravel the glacial history of this inaccessible location. Here, we present a comprehensive map of glacial landforms for the central sector of the CIS (55° to 60° N). Seven landform categories were mapped: ice flow parallel lineations, moraines (CIS outlet glacier moraines, Late Glacial moraines and moraines of unknown origin), meltwater channels (lateral and submarginal, subglacial, proglacial, and meltwater channels of unknown origin), kame terraces, eskers (single ridges and esker complexes), perched deltas and subglacial ribs. Collectively, these landforms provide a record of the extent, thickness and behaviour of the CIS, the direction of its movement and pattern of ice retreat.


Introduction
Mountainous British Columbia in western Canada has been repeatedly covered by the Cordilleran Ice Sheet (CIS) during the Quaternary. At the Last Glacial Maximum (LGM), it attained a volume and area similar to that of the present-day Greenland Ice Sheet. While it has long been recognized that the high mountains in northern British Columbia played an important role in establishing the ice sheet (Clague & Ward, 2011;Davis & Mathews, 1944;Fulton, 1991), the mountainous topography and remote location have thus far impeded our understanding of the central sector of the CIS and we still do not fully comprehend its LGM configuration or pattern of ice retreat. This makes it one of the least understood areas of any ephemeral Pleistocene ice sheet.
Ice flow is well documented in British Columbia. Both Kleman et al. (2010) and Shaw et al. (2010) give generalized ice flow maps of the entire North American Ice Sheet Complex (NAISC), which was made up of the Cordilleran, Laurentide and Innuitian ice sheets at the LGM, and Arnold et al. (2016) provide a compilation of all published streamlined landforms within British Columbia. However, the Glacial Map of Canada produced by Prest et al. (1968) at a 1:5,000,000 scale is the only ice-sheet-wide map of glacial landforms for the NAISC to date. This map is sufficient for deciphering the overall ice retreat pattern of the Laurentide Ice Sheet using the occurrence of eskers, which are thought to be deposited perpendicular to the retreating ice margin, drumlin orientation and ice marginal positions delineated by moraine systems (Dyke, 2004;Dyke et al., 2003;Dyke & Prest, 1987;Margold et al., 2018). However, in northern British Columbia the majority of ice marginal positions are recorded by meltwater landforms, such as lateral meltwater channels, kame terraces, eskers and perched deltas (Lakeman et al., 2008;Margold et al., 2013aMargold et al., , 2013b, which are not adequately recorded on the Glacial Map of Canada. Therefore, more detailed landform mapping is required to understand retreat of the central sector of the CIS. The increased resolution and coverage of digital elevation models (DEMs) and satellite imagery presents new opportunities to study the subglacial bed of paleo-ice sheets in detail (Chandler et al., 2018;Stokes et al., 2015). In order to gain a better understanding of the CIS, we use high resolution remotely sensed data to create a detailed glacial landform map centered on the LGM ice divide region in northern British Columbia. This map covers an area of ∼200,000 km 2 located between 55-60°N and 122.6-135.5°W (Figure 1).

Data
Glacial geomorphological mapping was primarily undertaken in ArcMap 10.6.1 using hillshade imagery derived from the Canadian Digital Elevation Model (DEM) of 0.75 arc second spatial resolution (12-20 m; complete coverage) (Government of Canada, 2019), and high resolution satellite images from Planet Lab (3-5 m resolution) (Planet Team, 2017) and Google Earth. Additionally, the German Aerospace Centre allowed us to download TanDEM-X DEM tiles of 0.4 arc second resolution (12 m) data covering 100,000 km 2 through their science project (see Figure 1 for coverage) (German Aerospace Center, 2018).

Landform mapping
The map is bound by the Yukon Territory-British Columbia border at 60°N, the Rocky Mountain Trench to the east and the Coast Mountains to the west, and the 55°N parallel to the south (Figure 1). The physiographic regions labelled on the Main Map are from Mathews (1986). A repeat-pass method was used to identify each landform through the entire Topographic System tiles that intersect with the mapped area are labelled. The green colour shows the coverage of TanDEM X DEM tiles used in this study. The orange boxes show the coverage of surficial geology maps from the Geological Survey of Canada that contain glacial landforms. The yellow box corresponds to the broad-scale glacial geomorphology map of Margold et al. (2011) and the blue box corresponds to the map area of Mathews et al. (1975). The location of the remaining figures is shown by the black stars. (b) Inset map showing the extent of the western margin of the North American Ice Sheet Complex at 22.1 cal ka BP drawn from Dalton et al. (2020) with the approximate position of the ice divides and the Keewatin ice dome (K) drawn in dark blue and the suture zone between the CIS and Laurentide ice sheet shown by the blue dashed line (Menounos et al., 2017;. The location of the study area is shown by the black box within the Cordilleran Ice Sheet. study area. To increase consistency, mapping was conducted at a variety of scales between 1:30,000 and 1:70,000 depending on landform type. Each landform has been identified from the imagery based on its morphology, spatial arrangement and association with other landforms (Table 1).

Ice flow parallel lineations
Ice flow parallel lineations include drumlins, flutes, grooved bedrock, and crag-and-tails. These landforms represent a variety of depositional and erosional ridges that are elongated in the direction of ice flow (e.g. Boulton & Clark, 1990a, 1990bClark, 1999;King et al., 2009). Ice flow parallel lineations often occur in fields or swams that are made up of hundreds of individual lineations that have regular morphology, spacing and orientation ( Figure 2). Where the stoss and lee side of a lineation was identified, for example, by the tapering sedimentary tail on the lee side of the crag-and-tails, the ice flow direction has been drawn on the Main Map. Complex arrangements of lineations with different orientations reflect different ice flow phases, which can be used to determine how flow changed over time (e.g. Greenwood & Clark, 2009;Jansson et al., 2002;Kleman et al., 1997).

Meltwater channels
Meltwater channels are formed in three main positions within an ice sheet: they are formed by water flowing along the ice margins (lateral and submarginal; Figure 3), by channelized flow at the bed of the ice sheet (subglacial; Figure 4), and by meltwater draining away from the ice sheet terminus (proglacial) (Greenwood et al., 2007(Greenwood et al., , 2016Mannerfelt, 1949;Margold et al., 2011). Each of these types of meltwater channels provide important paleo-glaciological information about the former ice sheet and therefore, where possible, they were distinguished based on their morphology and relationship to the surrounding topography (Table 1). However, sometimes a channel may transport different sources of meltwater at different stages of the ice sheet's evolution and a channel may have been draining water from both a glacially dominated source and from a fluvial or glaciofluvial source at a later stage. Therefore, where the meltwater channels has a glacial source but the type of meltwater channel could not be identified, they have been mapped as meltwater channels of unknown origin.

Kame terraces
Kame terraces are narrow, gently sloping accumulations of glaciofluvial sediment perched on valley sides that are formed by deposition of sediment along the lateral margins of a valley glacier or an ice sheet lobe (Borsellino et al., 2017;Evans, 2005;Turner et al., 2014) (Figure 3). They often occur in a series and are distinguished from river terraces by their high elevation valley margin location and close association with other glaciofluvial landforms, such as lateral and submarginal meltwater channels.

Eskers
Eskers are linear depositional ridges composed of glaciofluvial sand and gravel that are deposited by meltwater flowing through conduits beneath ice masses (Hebrand & Åmark, 1989;Shreve, 1985;Storrar et al., 2014). They exist as individual segments that often align to form networks up to 200 km in length and their morphology can vary from continuous single ridges to multiple ridged eskers to large esker complexes or deltas (Margold et al., 2011;Storrar et al., 2020). Within esker complexes it can be difficult to distinguish the individual esker ridges and where this occurs they are mapped as a polygon ( Figure 5). It is hypothesized that eskers form within Röthlisberger channels (R-channels) that are aligned roughly normal to the ice margin, and, therefore, they record the meltwater drainage pattern and ice marginal positions during ice retreat (Brennand, 2000;Hewitt & Creyts, 2019;Livingstone et al., 2015;Stroeven et al., 2016).

Moraines
Moraines occur as sharp-crested, straight or arcuateshaped ridges that are formed by the deposition or deformation of glacigenic sediment at the margins of active glaciers (Benn & Evans, 2010). A number of previous studies have recognized that moraines of different origins are located within northern British Columbia (Lakeman et al., 2008;Menounos et al., 2017;Ryder & Maynard, 1991). In particular, Lakeman et al. (2008) describe large, regionally extensive, sharp-crested moraines that extend up to 9 km beyond Little Ice Age moraines, which are interpreted to be deposited by a late Pleistocene advance of independent alpine glaciers (Lakeman et al., 2008;Menounos et al., 2017). While morphologically similar moraines have been noted across the Cassiar and Omineca mountains, their distribution has never been mapped.
We have mapped three different types of moraines based on their morphology, orientation and location: (1) CIS outlet glacier moraines; (2) Late Glacial moraines; and (3) moraines of unknown origin. The first type of moraine was deposited by outlet glaciers of the CIS that were flowing along glacial troughs during deglaciation ( Figure 6(c, d)) and they are identified by their close association with other ice marginal landforms, such as kame terraces and lateral and submarginal meltwater channels. The second type of moraine was deposited by glaciers emanating from mountain peaks, and they are identified as sharp-crested ridges that form on cirque and valley floors (Figure 6(a, b)). The final category is used when the origin of the Table 1. Diagnostic criteria for the identification of glacial geomorphological landforms from remotely sensed data (after Greenwood et al., 2007;Margold & Jansson, 2012;Norris et al., 2017). Broadly linear ridges that exhibit variation in tone across their crestline and shadowing due to changes in the relative relief. Exhibit both sharp and broad ridge crests. Associated with other ice contact landforms, such as deltas, kame terraces and meltwater channels. Lateral moraines perched on valley sides occur roughly parallel with the contours.
Landform crest drawn as a single line Late Glacial moraines Broadly linear ridges exhibit variation in tone across their crestline and shadowing due to changes in the relative relief. Terminal moraines are located at the end of a valley or cirque and have an arcuate shape. Often exhibit sharp ridge crests and may occur as multiple ridges. Unknown origin Broadly linear ridges which exhibit variation in tone across their crestline and shadowing due to changes in the relative relief. The origin of the moraine could not be identified. Distinguished from eskers based on their lack of sinuosity and location within a valley.

Perched deltas
Flat-topped accumulations of glaciofluvial sediment with steeply dipping frontal beds. Landform outline drawn as a polygon Subglacial ribs Ridges of sediment lying transverse to the flow direction. Typically occur in groups or swarms. May be drumlanised. Ridge sizes are highly variable. Landform outline mapped as a polygon around the break in slope.
moraine is uncertain and also includes moraines that could have been formed during the ice advance phase.

Perched deltas
Deltas are deposited when sediment that is transported by a river or stream is discharged into a standing body of water (e.g. a lake) and they typically have a flat top and steeply dipping front. As an ice margin retreats it can block the natural water drainage pathways resulting in transient ice-dammed lakes. Sediments deposited into these glacial lakes form deltas that remain perched on valley slopes once the glacial lake drains (Mannerfelt, 1949;Stroeven et al., 2016) ( Figure 3). It is possible that river terraces could be mistaken for perched deltas because they have a similar morphology. In order to reduce misidentification, we have only mapped deltas that are perched above the present-day rivers.

Subglacial ribs
Subglacial ribs are large, regularly spaced transverse ridges that are formed subglacially (Aylsworth & Shilts, 1989;Dunlop & Clark, 2006;Hättestrand & Kleman, 1999;Lundqvist, 1989) (Figure 7). While these landforms have a variety of names within the literature, including Rogen or ribbed moraines, here we use the non-genetic term subglacial ribs following the BRITICE Glacial Map . Although the morphology and size of subglacial ribs is highly variable (Dunlop & Clark, 2006;Stokes et al., 2016), they often have a curved or anastomosing pattern. Individual ridges usually have an asymmetric profile and multiple ribs often occur together forming fields.

Accuracy and comprehensiveness
We believe the distribution in landforms presented on this map reflects the true distribution pattern rather than being related to variable data sources. Although TanDEM-X tiles only cover half of our mapped area, we found that this higher resolution data did not allow us to identify and map new landforms, it rather allowed us to map certain landforms, such as moraines, with more confidence. However, some moraines, eskers and lateral and submarginal meltwater channels are near the resolution limit of the remotely sensed data, and therefore, it is possible that smaller features may have been missed. Landforms that may be misinterpreted, such as eskers and moraines, are listed in Table 1.

Ice flow parallel lineations
In total 52,300 ice flow parallel lineations were mapped during this study. Overall, our lineations are consistent with the generalized ice flow maps of Kleman et al. (2010) and Shaw et al. (2010) and the detailed ice flow indicator compilation by Arnold et al. (2016). The ice flow parallel lineations range in size from tens of meters to 2 km in length but megascale glacial lineations were not identified within the mapped area. Swarms of ice flow parallel lineations occur throughout the mapped area. Their location within the mountainous terrain, either at the bottom of large valleys (e.g. Teslin Trench), on high elevation peaks and plateaus, and on the plains surrounding the mountainous region, e.g. the Liard Lowland (Figure 2(c, d)), likely corresponds to different phases of ice flow. Furthermore, it is common for the orientation of the ice flow parallel lineations at high elevations (e.g. on mountain summits) to differ from the orientation of those located in the adjacent valleys (e.g. Atlin; see Main Map).
On the floor of some valleys ice flow parallel lineations overprint subglacial ribs (e.g. Teslin Trench). At other locations, they form a divergent, fan-shaped pattern (e.g. Liard Lowland; Figure 2(c, d)). Occasionally, the lineations have complex orientations with cross-cutting relationships, implying there were significant differences in the ice flow direction at different stages of ice sheet evolution (Figure 2(a, b)).

Lateral and submarginal meltwater channels
In total 18,420 lateral and submarginal meltwater channels were mapped during this study, making them the most common type of meltwater channel in northern British Columbia. These meltwater channels vary in length from a few hundred meters to 15 km and they generally have a low sinuosity, with sinuosity increasing when the channel becomes submarginal. Lateral and submarginal meltwater channels can occur at any elevation: the highest meltwater channels often cut small notches in ridges or spurs emanating from mountain peaks (up to 1750m asl); and the lowest channels form close to the valley floor (down to 700 m asl).
These channels are the most distinctive glacial landform identified within the Cassiar and Omineca Mountains, where thousands of channels have been mapped, while they are noticeably absent from the Coast and Skeena mountains. Additionally, while the studies of Margold et al. (2011Margold et al. ( , 2013aMargold et al. ( , 2013b determine the overall direction of glacial meltwater flow through generalized mapping of lateral meltwater channels, this study presents the first map of individual meltwater channels for the central sector of the CIS.

Subglacial meltwater channels
In total 245 subglacial meltwater channels were mapped during this study. While they can occur at any topographic position, and sometimes dissect mountain ridges, they are most common on the valley floor. The mapped subglacial channels often contain plunge pools that form present-day ponds within the channel (Figure 4). The majority of large subglacial meltwater channels (widths between 200 and 800 m) are located in the Liard Lowlands (see Main Map). Here some of the subglacial meltwater channels are aligned with eskers and esker complexes and have roughly consistent orientations over tens of km.

Proglacial meltwater channels
In total 9 wide, meandering proglacial meltwater channels (up to 1 km wide) were identified within the map area. These meltwater channels are primarily recognized based on their relationship to other glacial landforms. For example, the large proglacial meltwater channels in the Liard Lowland crosscut esker complexes and subglacial meltwater channels (see Main Map), indicating that large volumes of water flowed through these channels after the ice margin retreated.

Meltwater channels of unknown origin
In total 304 glacial meltwater channels were identified within the map area where the source of glacial meltwater could not be identified. They vary in length from a few hundred meters to 30 km.

Kame terraces
In total 983 kame terraces were mapped in the study area and they are located on the valley walls throughout the Omineca and Cassiar mountains. They vary in Figure 4. Example of subglacial meltwater channels (a) high-resolution Planet Lab satellite imagery, (b) geomorphological mapping and (c) elevation profile along the channel A-A ′ produced using the Canadian DEM. The elevation profile shows the channel has an undulating long profile, which is a unique feature of subglacial meltwater channels (Table 1), and plunge pools can be seen in the satellite imagery. The location of these figures is shown in Figure 1. size from a few hundred meters to 1 km in length and they are up to a few hundred meters wide. They are often associated with other ice marginal landforms, such as lateral and submarginal meltwater channels.

Esker ridges
Esker ridges have been mapped throughout the study area with the exception of the Coast Mountains, Nass Depression and western Skeena Mountains. In total 3574 esker ridges were mapped and they are most commonly located on valley floors, where semi-continuous ridges can be traced for tens of km and anastomosing ridges are common. At some locations esker ridges have a preferred orientation. For example, the esker ridges located in the Cassiar Mountains and the Liard Lowland in the northeast corner of the map are generally oriented in a northeast-southwest direction.

Esker complexes
In total 80 esker complexes were mapped in this study, which range from 1 to 15 km in length. Esker complexes are primarily located on the eastern side of the mapped area, with the largest complexes occurring within the Liard Lowland and the Rocky Mountain Trench (see Main Map). However, they also occur on the floor of some valleys within the Cassiar and Omineca mountains.

CIS outlet glacier moraines
In total 97 ridges have been mapped as CIS outlet glacier moraines during this study. These moraines are made up of both single and multiple ridges and they vary in length between 2 and 4 km, with widths of ∼250 m. All of CIS outlet glacier moraines are located within the Cassiar Mountains, with the exception of 5  moraine ridges that were mapped in the Omineca Mountains.

Late Glacial moraines
In total 1459 ridges have been mapped as Late Glacial moraines, making them the most common type of moraine mapped within the study area. These moraines are made up of single and multiple ridges and they vary in length from a few hundred meters up to a few kilometers. The largest Late Glacial moraine, which is located southeast of Tuya Lake (see Main Map), consists of multiple ridges with a total length of 8 km. The minimum elevation of the Late Glacial moraines is 1100 m asl. They are primarily located in the Cassiar and Omineca mountains, but, they also occur on the Tuku Plateau and Meszah Peak, while they are absent from the Skeena and Coast mountains.

Moraines of unknown origin
In total 378 ridges have been mapped as moraines of unknown origin (20%). Some of these ridges have a subdued topography and may have been formed during the ice sheet advance stage, however, this cannot be determined by remote sensing methods alone.

Perched deltas
In total 223 perched deltas were mapped during this study, which range in size from a few hundred meters to 16 km wide and 20 km long. The largest perched deltas are located in the Liard Lowland and along the Rocky Mountain Trench, while smaller perched deltas are concentrated in the Cassiar and Omineca mountains.

Subglacial ribs
Here we map subglacial ribs in northern British Columbia for the first time with a total of 683 ribs mapped. These subglacial ribs have a variety of sizes and shapes, with the largest ribs being 5 km long and 1 km wide. The ribs are spaced between 500 m and 1 km apart. While subglacial ribs occur on valley floors in a few locations across the mapped area, the majority are located within the Teslin Trench. Ice flow parallel lineations are sometimes superimposed on the subglacial ribs.

Regional implications and conclusions
The Main Map provides a detailed broad-scale glacial landform record of the central sector of the CIS. Overall, the mapped landforms are consistent with the landform record described in the existing literature (e.g. Kleman et al., 2010;Lakeman et al., 2008;Margold et al., 2013aMargold et al., , 2013bMargold et al., , 2014Shaw et al., 2010). For example, our mapped ice flow parallel lineations generally agree with the generalized flow maps of Kleman et al. (2010) and Shaw et al. (2010) and the ice flow indicator map of Arnold et al. (2016). However, our glacial landform map presents the first large-scale, detailed map for many of the landforms (e.g. lateral and submarginal meltwater channels and Late Glacial moraines), and therefore, thousands of previously unrecognized glacial landforms, including eskers, moraines, meltwater channels, kame terraces and subglacial ribs, have been identified and mapped using the high resolution remotely sensed data. The majority of newly mapped landforms are located in the Cassiar and Omineca mountains.
It is clear from the map that there is a marked difference between the glacial landform assemblage within the Cassiar and Omineca mountains to the east, where thousands of lateral meltwater channels, eskers and Late Glacial moraines have been mapped, and the Skeena and Coast mountains to the west, where these glacial landforms are noticeably absent. This glacial geomorphological map can now be used to unravel the complex glacial history beneath the LGM ice divide region of the CIS and help us to better understand the advance and retreat dynamics of this ephemeral Pleistocene ice sheet.

Software
The hillshade surfaces were produced from the DEM data within ESRI ArcMap 10.6.1. On-screen digitizing of landform crestlines and break-in-slopes was conducted in ArcMap 10.6.1 in the ESRI shapefile format and in Google Earth as KML files that were then converted to ESRI shapefiles in ArcMap 10.6.1. The pdf map was exported and the final map was created in Adobe Illustrator 2020.

Availability of data
The ESRI shapefiles produced for each glacial landform type are supplied with this paper.