Origin, burial and preservation of late Pleistocene-age glacier ice in Arctic permafrost (Bylot Island, NU, Canada)

Over the past decades, observations of buried glacier ice exposed in coastal bluffs and headwalls of retrogressive thaw slumps of the Arctic have indicated that considerable amounts of late Pleistocene glacier ice survived the deglaciation and are still preserved in permafrost. In exposures, relict glacier ice and intrasedimental ice often coexist and look alike but their genesis is strikingly different. This paper aims to present a detailed description and infer the origin of a massive ice body preserved in the permafrost of Bylot Island (Nunavut). The massive ice exposure and core samples were described according to the cryostratigraphic approach, combining the analysis of permafrost cryofacies and cryostructures, ice crystallography, stable O-H isotopes and cation contents. The ice was clear to whitish in appearance with large crystals (cm) and small gas inclusions (mm) at crystal intersections, similar to observations of englacial ice facies commonly found on contemporary glaciers and ice sheets. However, the δ¹⁸O composition ($-\mathrm{34.0}\pm \mathrm{0.4}$ ‰) of the massive ice was markedly lower than contemporary glacier ice and was consistent with the late Pleistocene age ice in the Barnes Ice Cap. This ice predates the aggradation of the surrounding permafrost and can be used as an archive to infer palaeo-environmental conditions at the study site. As most of the glaciated Arctic landscapes are still strongly determined by their glacial legacy, the melting of these large ice bodies could lead to extensive slope failures and settlement of the ground surface, with significant impact on permafrost geosystem landscape dynamics, terrestrial and aquatic ecosystems and infrastructure.

1 Introduction

In the Arctic, extensive areas of ridged and hummocky moraines are underlain by buried glacier ice (Dyke and Savelle, 2000; Kokelj et al., 2017; Solomatin, 1986). Massive-ice bodies within these landscapes are important indicators of past glacial, hydrologic and hydrogeologic conditions and are used to reconstruct regional palaeo-environments and palaeo-climates (Jorgenson and Shur, 2008; Murton et al., 2005; Lacelle et al., 2018). Areas with buried glacier ice are also becoming increasingly vulnerable to climate warming (Kokelj et al., 2017). Glacier ice is the most common type of buried ice in permafrost and its occurrence was reported in Russia (Astakhov, 1986; Belova et al., 2008; Ingólfsson and Lokrantz, 2003; Kaplanskaya and Tarnogradskiy, 1986; Solomatin, 1986), the Canadian Arctic (Dallimore and Wolfe, 1988; French and Harry, 1990), Alaska (Jorgenson and Shur, 2008; Kanevskiy et al., 2013) and Antarctica (Sugden et al., 1995; Swanger, 2017). Buried glacier ice has been commonly observed in the proglacial zone of contemporary glaciers and can be preserved in formerly glaciated areas (the paraglacial to periglacial zones) within large moraine belts, hummocky till and glaciofluvial deposits (Everest and Bradwell, 2003; Tonkin et al., 2016). In the Canadian Arctic, buried glacier ice has been mainly reported in the western and central regions, especially in the Mackenzie Delta region (French and Harry, 1990), Tuktoyaktuk Coastlands (Murton et al., 2005), Herschel Island (Fritz et al., 2011), Yukon (Lacelle et al., 2007), Banks Island (Lakeman and England, 2012) and Victoria Island (Dyke and Savelle, 2000; Lorrain and Demeur, 1985). While the permafrost of the eastern Canadian Arctic is expected to contain remnants of Pleistocene ice sheets or glaciers, very few have been reported so far.

Distinguishing between buried glacier ice and other types of massive ground ice in the permafrost is usually based on cryostratigraphy combined with detailed studies of physical, geochemical and isotopic properties of the ice that may also include analyses of occluded gases (Cardyn et al., 2007; Fritz et al., 2011; Ingólfsson and Lokrantz, 2003). When classifying glacier ice facies, a distinction is made between englacial and basal ice facies (Fortier et al., 2012; Lawson, 1979). The englacial (firn-derived) ice facies is formed by the gradual snow compaction and recrystallization, a process called firnification, and has a low debris content (Benn and Evans, 2010). Basal ice has distinctive physical and chemical characteristics and has a much higher debris content than the overlying englacial ice as a result of subglacial processes operating at the glacier bed: regelation, glaciohydraulic supercooling accretion, net basal adfreezing, incorporation of ice and sediments by overriding ice during glacier advance, glacio-tectonics, ice lens aggradation and downward propagation of cold temperature in the sediment at the glacier bed (congelation) (Alley et al., 1998; Fortier et al., 2012; Hubbard and Sharp, 1995; Knight, 1997; Lawson, 1979; Sharp et al., 1994). Both types of glacier ice may experience burial but basal ice is probably the most common form of buried glacier ice according to reports from various circumpolar permafrost regions (Belova et al., 2008; Fritz et al., 2011; Murton et al., 2005; St-Onge and McMartin, 1999). The process of burial of glacier ice has been described by Shur (1988), Solomatin (1986), Harris and Murton (2005) and citations therein. Burial of glacier ice may occur as a result of (1) accumulation of fluvial, lacustrine, aeolian or slope sediments on top of the ice; (2) glaciotectonic processes or (3) formation of insulating blanket of supraglacial melt-out till. With the latter scenario, sediment-rich basal ice has a greater potential to persist in a buried state than englacial ice with little debris. Buried glacier ice remains stable for a long period of time only if the soil temperature is below freezing, and the active-layer thickness does not exceed a depth to the massive ice body (Shur, 1988).

In this study, we describe the occurrence of massive ice preserved in the permafrost of the Qarlikturvik Valley, south-western Bylot Island (NU, Arctic Canada). We investigated the physical and geochemical properties of a recently exposed body of massive ice and compared them with those of other ice types in the region (snow, glacier ice, wedge ice, segregated ice) to infer its origin. A cryostratigraphic approach was used to delineate cryostratigraphic units on the basis of their cryostructures, physical properties and thaw unconformities (French and Shur, 2010; Gilbert et al., 2016), combined with crystallography and geochemical analyses of the different ground-ice types. The origin of the massive ice is discussed.

Figure 1 (a) Location of Bylot Island, (b) location of the Qarlikturvik valley, (c) location of the study area within the Qarlikturvik Valley (valley of glaciers C-79 and C-93), Bylot Island, Nunavut, Canada. Wedge ice, segregated ice and snow were also sampled near (<1 km) the massive ice exposure. The thaw depths were measured with a steel probe at every 10 m along a 150 m transect.

2 Regional setting

With ice-capped summits dominating the central highlands and glaciated valleys that extend near the coast, the mountainous central section forms a striking contrast with the relatively flat coastal lowlands. The Byam Martin Mountains (∼1400 m a.s.l.) consist primarily of Archean–Aphebian crystalline igneous rocks and Proterozoic metasedimentary and metamorphic rocks (Jackson and Davidson, 1975). Klassen (1993) suggested that alpine glaciers, larger but similar in size to those present today, occupied Bylot Island during the late Wisconsinan. At the Last Glacial Maximum (LGM), ice streams of the Laurentide Ice Sheet (LIS) flowed in adjacent marine channels (Lancaster Sound, Navy Board Inlet and Eclipse Sound) and reached Bylot Island (De Angelis and Kleman, 2007). The study area is situated in the Qarlikturvik Valley (73^∘09^′ N, 79^∘57^′ W, 25 m a.s.l.) on south-west Bylot Island at about 80 km north-west of Mittimatalik (Pond Inlet; Fig. 1b). This valley was eroded through a Cretaceous–Tertiary sequence of poorly consolidated sandstone and shale (Jackson and Davidson, 1975).

The Qarlikturvik Valley is a typical U-shaped glacial valley with surface sediments reflecting the complex history of the valley: presence of unconsolidated glacial, colluvial, alluvial, marine, aeolian and organic deposits dating back to the Late Pleistocene and Holocene (Allard, 1996; Fortier and Allard, 2004). The valley comprises low-lying ice-wedge polygon terraces bordering a proglacial braided river running in a glaciofluvial outwash plain and forming a delta in Navy Board Inlet. During the Holocene glaciers C-93 and C-79 retreated up-valley and today they are located about 14 km from the coast. Following glacial retreat, the valley became partially submerged between 11 335 to 6020 cal yr BP as a result of a marine transgression (Allard, 1996). Alternating layers of peat and aeolian sands and silts (∼2–3 m) cover the glaciofluvial terraces of the valley, where an extensive network of syngenetic ice-wedge polygons have developed after 6000 cal yr BP (Fortier and Allard, 2004). Mounds of reworked till and ice-contact-stratified sediments mark a former position of the glacier in the valley.

Bylot Island belongs to the Arctic Cordillera and the northern Arctic terrestrial ecozones. The mean annual air temperature between 1971 and 2000 in Pond Inlet (NU, Canada) was $-\mathrm{15.1}\pm \mathrm{5.1}$ ^∘C, increasing slightly to $-\mathrm{14.6}\pm \mathrm{4.9}$ ^∘C between the 1981 and 2010 (Environment Canada, 2015). No significant trends in precipitation have been observed over the last decades, with a mean annual precipitation of 189 mm, much of it falling as snow. Thawing and freezing indices averaged (1981–2010) 473 degree days above 0 ^∘C and 5736.1 degree days below 0 ^∘C respectively (Environment Canada, 2015). Vegetation in the valley is typical of Arctic tundra environments and is largely determined by soil moisture and slope (Duclos et al., 2006). Wetlands occur in low-lying areas, commonly with low-centred polygons, and are typically dominated by grasses and sedges. Mesic tundra is found in areas characterized by better-drained soils (i.e. uplands and hillslopes). The climatic and vegetation conditions determine the presence of continuous permafrost in the ice-free areas of Bylot Island. Locally, permafrost thickness >400 m has been detected from thermal measurements on the nearby Somerset and Devon islands (Smith and Burgess, 2002). Active-layer thickness varies from ∼1 m in drained unvegetated sands and gravels to ∼0.3–0.7 m in peaty and silty soils (Allard et al., 2016).

Figure 2 A photograph and a schematic cross section showing generalized stratigraphy of the massive ice exposure and the overlying sediments. The lower and lateral contacts of the massive ice have not been reached. The thaw depths measured in late July (2013) at the headwall reaches 55 cm. The red star indicates the sampling location of the organic material and the red dots show the sampling points for stable O-H isotope and hydrochemistry.

3 Material and methods

A large body of massive ice was found exposed beneath ∼1.7 m of sediments along the headwall of a thaw slump in the Qarlikturvik Valley, ∼10 km from the terminus of C-93 and C-79 glaciers (Fig. 1c). Since a great amount of slump material covered the ice, the exposure was cleaned and enlarged to allow a better description and sampling. The excavated section along the headwall of the slump attained a height exceeding 7 m and extended laterally over 10 m. The lower and lateral contacts between the ice body and the surrounding sediments could not be delineated. The exposure was subdivided into three units from the bottom: (a) massive ice, (b) sand and gravel and (c) muddy sand diamicton overlaid by peat (Fig. 2). The thaw depths were measured with a steel probe at every 10 m along a 150 m transect that started at the upper headwall of the thaw slumps. Samples from the massive ice (unit A) were collected using an axe and a portable core drill equipped with an 8 cm diamond carbide core barrel. Between 2011 and 2013, ice cores were extracted every 10 cm or less from depths ranging from 1 to 3 m below the surface (Fig. 2) without reaching the bottom of the ice body. For comparison, modern glacier ice (englacial ice) was sampled from freshly collapsed ice blocks at the margin of glacier C-93 located a few kilometres up the valley (Fig. 1c). Wedge ice, segregated ice and snow were also sampled for geochemical analysis at the site located near (<1 km) the massive ice exposure. Samples (n=5) from the sediments (units B and C) overlying the massive ice were collected and characterized. All samples were kept frozen, shipped and stored in our laboratory for further analysis.

A cryostratigraphic approach was used to describe the massive-ice body and the overlying sediments (Fortier et al., 2012; French and Shur, 2010; Gilbert et al., 2016; Murton and French, 1994). Cryostratigraphic units were delineated based on cryostructures and cryofacies. Cryofacies are bodies of frozen sediment that are visually distinct from adjacent frozen sediments based on their cryostructures, volumetric ice content, ice-crystal size and sediment texture. Cryostructures refer to the shape, amount and distribution of ice within the frozen sediment. Gas inclusions visible within the ice and the deformation structures in the ice and sediments were also described (Hambrey and Lawson, 2000). To further investigate the cryostratigraphic characteristics of the massive ice body, all samples were observed under X-ray-computed tomography (CT) scanning (Siemens SOMATOM Sensation 64), a technique that allows visualizing and reconstructing the internal structure (2-D and 3-D) of permafrost samples at <1 mm resolution (Calmels et al., 2010; Dillon et al., 2008; Fortier et al., 2012). The complete CT-scanning procedure used in this study is presented in Appendix A. Crystallographic analysis of the massive ice and modern glacier ice was conducted to measure their crystal size and shape (surface area, long axis and circularity ratio) as these parameters contain information about the conditions under which the ice was formed (French and Shur, 2010). The crystalline structure was investigated through thin sections of ice mounted on a glass slide under cross-polarized light. Thin sections were made by cutting the ice sample vertically and/or horizontally into thin slices (thickness: 0.2 to 0.4 mm). C-axis orientations of the crystals have not been measured since the horizontal orientation of the ice samples was not preserved following the sampling. Fiji image analysis software was used to measure the area, long axis and circularity ratio of each crystal (Schindelin et al., 2012). Differences in the crystal shape (area, long axis, circularity ratio) of horizontal and vertical thin sections were tested with the Mann–Whitney U test using R, which is a programming language software for data analysis and graphics (R Development Core Team, 2016).

The massive ice body, along with glacier C-93 ice, ice wedges, intrasedimental ice and snow, were analyzed for their geochemical and isotopic (δ¹⁸O, δD) composition, an approach that can shed light into the origin of ground ice (Lacelle and Vasil'chuk, 2013). Sampling of the massive ice body was done at 10 cm vertical intervals or less, depending on the unit change; prior to sampling, at least 10 cm of the ice surface was removed with an ice axe. In 2013, all samples (n=80) were melted, filtered (0.45 µm diameter filter) and acidified in the laboratory in sealed polyethylene bottles prior to analyses. The concentration of major cations in the ice and snow (Al_tot, Ca²⁺, Fe_tot, K⁺, Mg²⁺, Mn_tot, Na⁺) was determined by inductively coupled plasma optical emission spectrometry (Vista Pro ICP-OES) at the University of Ottawa. Solutes are expressed in milligrams per litre and analytical reproducibility was ±1 %. The stable isotope ratios of oxygen (¹⁸O∕¹⁶O) and hydrogen (D∕H) were determined using a Los Gatos Research high-precision liquid water analyzer coupled to a CTC LC-PAL autosampler. The results are presented using the δ-notation (δ¹⁸O and δD), where δ represents the parts per thousand differences for ¹⁸O∕¹⁶O or D∕H in a sample with respect to Vienna Standard Mean Ocean Water (VSMOW). Analytical reproducibility for δ¹⁸O and δD was ±0.3 ‰ and ±1 ‰, respectively.

The origin of the sediment overlying the massive ice was inferred from particle-size distribution, the clasts and sand-size quartz grain morphoscopy and the geochemical and isotopic (δD-δ¹⁸O) composition of pore water. Particle-size distributions were determined by dry sieving at 1/2φ intervals (size ranges −12 to 4φ). The hydrometer method was used to determine the distribution of the finer particles smaller than 4φ (ASTM Standard D422, 2007). Descriptive statistics (mean grain size, sorting, skewness) and Folk and Ward sediment classes were determined using the RYSGRAN package for R (Gilbert et al., 2014; R Development Core Team, 2016). Fifty in situ pebble- to cobble-sized clasts were randomly collected from the uppermost unit (Unit C). These clasts were analyzed for shape, roundness and lithology using techniques described by Benn (2004). Morphoscopic analyses of small quartz grains (0.5–1.0 mm) were conducted using a binocular microscope (Cailleux and Tricart, 1963). Additionally, the ground ice in the sediments was analyzed for soluble ions (major cations) and δD-δ¹⁸O composition following the method described above. A fragment of poorly decomposed peat sampled in unit B was radiocarbon dated by accelerator mass spectroscopy (AMS) (ULA-6505, Centre for Northern Studies, Université Laval). Radiocarbon age was calculated as −8033ln (F¹⁴C) and reported in ¹⁴C yr BP (BP = AD 1950) (Stuiver and Polach, 1977) and then corrected to calendar years (cal yr BP) using Calib 7.10 and the IntCal13 calibration curve (Reimer et al., 2013; Stuiver et al., 2017).

Figure 3 Photographs (left) and CT scans (right) of the massive ice body. (a) Pure-ice facies. The dotted red line highlights one single crystal; (b) ice-poor sediment (sands and gravels) with suspended and crustal cryostructures forming a band in the massive ice.

Figure 4 Comparison between the massive ice body and modern glacier ice (glacier C-93, Bylot Island). The first column shows unprocessed photographs of the ice taken under a plain light with surficial sediment inclusions highlighting the crystal boundaries. The dotted red line highlights one crystal. The second column shows thin sections of ice sample viewed under direct cross-polarized light. The third column shows the microstructure (crystal boundaries) extracted from the thin-section photograph.

4 Results

4.1 Cryostratigraphy and properties of the massive ice body

The exposed massive ice body (unit A) was >10 m thick and had a clear to milky white appearance due to its high bubble content (Fig. 3a). With a volumetric ice content near 100 %, it refers to the “pure ice” cryofacies described by Murton and French (1994). Occasional thin bands of sediments (sands and gravels) with suspended and crustal cryostructures were cutting across the ice (Fig. 3b). These discrete and planar bands were <2 cm thick, subparallel to one another and showed a dip direction (21 to 31^∘) down the ice body in the south-eastern direction. Crystallographic analysis of thin sections of the massive ice under cross-polarized light showed that the crystals had mostly bluish colours, suggesting that the ice crystals had similar orientations (Figs. 4 and 5a). Coarse-grained ice crystals characterized the massive ice body: long-axis average of 7.97±2.97 mm (3.13–16.58 mm); average surface area of 34.9±25 mm² (5.8 and 153.5 mm²); and average circularity ratio of 0.65±0.09, indicating the crystals were relatively rounded (Fig. 6). No significant differences in the shape properties (surface area, long axis, circularity ratio) were observed between the horizontal and vertical thin sections (Mann–Whitney–Wilcoxon test, p>0.05), indicating that the ice crystals were nearly equiaxial. The glacier C-93 ice displayed ice crystals with varying colours (Fig. 4). The ice crystals were larger (average surface area of 125.29±148 mm²) than those of the massive ice body (Figs. 4 and 6). Glacier C-93 ice consists of relatively rounded crystals as their mean circularity ratio averaged 0.65±0.01. The total volumetric content of gas inclusions varied from 2 % to 10 % for both the massive ice and glacier C-93 ice, with bubbles being mostly located along grain boundaries (Fig. 5a). The massive ice body contained three types of gas inclusions: (A) spherical bubbles, (B) flattened disks and (C) clusters of deformed and coalescent bubbles (Fig. 5b, c, d). The tiny gas bubbles had a mean surface area of 0.85±1.04 mm² with circularity ratio averaging 0.89±0.17. Gas bubbles observed in glacier C-93 were mostly spherical and small, with a few clusters of deformed coalescent bubbles. The gas bubbles had a mean surface area of 0.13±0.41 mm² and an average circularity ratio of 0.89±0.18.

Figure 5 (a) A thin section of the massive ice viewed under cross-polarized light and a transverse cross section from a scan showing the gas inclusions within the ice. (air is black; ice is dark grey). Photos to the right show patterns of gas inclusions; (b) small (submillimetre to millimetre) spherical bubbles (vertical bars are from measuring ruler of microscope stage); (c) small disks up to 6 mm in diameter; (d) coalescent bubbles and small disks all flattened in the same direction.

Major cations in the massive ice body (i.e. Ca²⁺, Na⁺ Mg²⁺ and K⁺) all occurred in low concentrations (mean <1.76 mg L⁻¹; Fig. 7). The cation concentrations of glacier C-93 ice were very similar to those of the massive ice, whereas the ice wedge had slightly higher cation concentrations, with average concentrations ranging from 1.11±0.18 mg L⁻¹ for Mg²⁺, to 3.32±0.44 mg L⁻¹ for Ca²⁺. The δD-δ¹⁸O composition of the massive ice, along with those from glacier C-93 and other types of ground ice present in the valley are shown in Fig. 8b. The δ¹⁸O composition of the massive ice had a narrow range (average δ¹⁸O: $-\mathrm{34.0}\pm \mathrm{0.4}$ ‰) with D-excess (d=δD – 8δ¹⁸O) averaging 6.6±2.5 ‰. The δ¹⁸O composition of the massive ice was much lower than that of the snow (average δ¹⁸O: $-\mathrm{30.4}\pm \mathrm{1.8}$ ‰), ice-wedge ice (average δ¹⁸O: $-\mathrm{25.6}\pm \mathrm{0.95}$ ‰) and glacier C-93 ice (average δ¹⁸O: $-\mathrm{25.0}\pm \mathrm{3.1}$ ‰). The D-excess of snow, ice wedge, glacier C-93 and ground ice (unit B) samples averages 8.9±3.4 ‰; 9.3±7.1 ‰, 5.2±5.8 ‰ and $-\mathrm{52.4}\pm \mathrm{31.4}$ ‰ respectively. In a δD-δ¹⁸O diagram, the samples from the massive ice, glacier C-93, snow and ice wedges were distributed along linear regression slope values of 6.0 (R²=0.44), 7.1 (R²=0.95), 6.6 (R²=0.96) and 5.2 (R²=0.36) respectively.

Figure 6 Box plots comparing the distribution of ice crystal characteristics (area, long axis, circularity ratio) of horizontal thin sections (H) and vertical thin sections (V) obtained from massive ice samples. C-93 represents data obtained from a sample of modern glacier ice sampled from glacier C-93.

Figure 7 Box plots showing soluble cation concentration of the massive ice, ice wedge, glacier C-93 and intrasedimental ice sampled within the sediments layers covering the massive ice unit.

Figure 8 (a) δ¹⁸O and D-excess depth profiles including both the massive ice unit and the intrasedimental ice from the sediment cover; (b) δ¹⁸O-δD diagram of the massive ice and other types of ground ice (ice wedge, intrasedimental ice), snow and modern glacier ice (C-93) sampled on Bylot Island. The solid red line represents the Global Meteoric Water Line (GMWL): d=δD–8δ¹⁸O+10 (Dansgaard, 1964). Also shown are the δ¹⁸O and δD values of precipitation recorded at Pond Inlet (1990–1992) and Resolute Bay (1989–1993; IAEA/WMO, 2016). LMWL is the local meteoric water line for Resolute Bay, Nunavut (δD=7.8δ¹⁸O +5.33; IAEA/WMO, 2016) and Pond Inlet, Nunavut (δD=8.0δ¹⁸O +10.55; IAEA/WMO, 2016).

Figure 9 Sedimentological data from the stratigraphic section. From left to right: a stratigraphic log showing the mean grain size of each unit; gravel, sand and mud percentages; and skewness and sorting.

4.2 Cryostratigraphy and properties of the overlying sediments

The massive ice body was covered by ∼1.7 m of sediments. The thaw depths measured in late July (2013) ranges between 19 and 55 cm (mean: 30±9 cm). A sharp, subplanar and unconformable contact separated the ice from the overlying sediments (units B and C) along the exposed section (Fig. 2). Unit B (∼115 cm thick) directly overlaid the massive ice body and has ice-poor sediment cryofacies with a structureless cryostructure, essentially made of pore ice. Unit B is subdivided into three subunits: B₁ is coarse sandy gravel, B₂ is sandy gravel, and B₃ is stratified gravelly sand. All subunits were texturally similar, consisting mostly of sands and gravels, with silt- and clay-sized particles constituting <3 % of the sediment (Fig. 9). The grain-size fraction >32 mm was not included in profile distribution, but abundant cobbles were observed within the bottom subunits. Sediments were poorly sorted to very poorly sorted (sorting values ranging between 3 and 5) and were coarser at the base and finer near the top (Fig. 9). Subunits B₁ and B₂ were both characterized by unstratified sandy gravel that became finer towards the uppermost subunit B₃, which consisted of stratified gravelly sand with thin subhorizontal laminae (<1 cm). Morphoscopic analysis of subunit B₃ showed that quartz grains were mostly angular (70 %) to subangular (smooth and polished, glazing grains: 25 %). Rounded and frosted grains represented <5 % of the total sand fraction. Plants fragments located at the base of subunit B₃ were dated to 885±15¹⁴C yr BP (1164 cal yr BP; 1σ range: 1058–1204). The sand and gravel sequence of unit B was sharply overlain by ∼55 cm of an unstratified (massive) thawed diamicton, with abundant pebble- to cobble-sized clasts (unit C; Figs. 2 and 9). The matrix (sandy silt) was very poorly sorted and the grain-size distribution tended to be fine skewed. Clasts had a wide range of shapes, with predominance in the subangular and subrounded classes (76 %). Only 24 % of all clasts were rounded or angular. Gneiss was the dominant cobble type with an average of 60 % in the matrix, whereas sedimentary and igneous rocks averaged 36 % and 4 %. A thin, continuous and irregular layer of dark fibrous peat with roots (subunit C2) overlaid the diamicton.

Supernatant water samples (n=5) from the ice-poor sediments in unit B showed markedly higher δ¹⁸O and cation values compared to the underlying massive-ice body (Figs. 7 and 8a). In unit B, the δ¹⁸O values became progressively higher, from −25.7 ‰ to −15.4 ‰, as we moved upward from the ice–sediment contact towards the surface (Fig. 8a). The cations content in the ice-poor sediment was 1 to 2 orders of magnitude higher than in the massive ice and also shifted to higher concentrations towards the surface (Fig. 7).

Table 1 Comparisons of the physical properties of different massive tabular ground ice found in the permafrost.

¹ Parallel to the heat flow direction, indicating that freezing is downward. ² Occurrences of debris bands, boudinage and pinch-and-swell structures, folding and thrust faulting. ³ Depending on the position of the ice within the glacier or ice sheet at the time it becomes buried.

5 Discussion

5.1 Origin of the massive ice

The cryostratigraphic and crystallographic properties of the massive ice along with its isotopic and geochemical composition indicate that the exposed ice consists of relict Pleistocene englacial ice buried and preserved in the permafrost of Bylot Island. Table 1 summarizes the cryostratigraphic and crystallographic properties of different types of tabular massive ice described in the literature: buried glacier ice (basal and englacial), massive segregated-intrusive ice and buried snowbanks. The appearance and structure of the buried massive ice body are similar to those of englacial ice typically observed at the margin of glaciers, ice caps or ice sheets. The buried massive ice body has a whitish appearance owing to its high concentration of air bubbles. Coarse-bubbly ice is the most abundant type (90 %–95 %) of englacial ice found in glaciers (Allen et al., 1960). Our results also show that the cross-sectional area of the crystals of the buried massive ice is smaller than that of neighbouring C-93 glacier ice, but there is no significant difference in their circularity ratio (Mann–Whitney–Wilcoxon test, p=0.89). However, the difference in ice crystal size is not unforeseen since the ones from glacier ice can show variations of the order of a few millimetres to several centimetres in diameter (Gow, 1963; Thorsteinsson et al., 1997). Patterns of preferred crystal orientation combined with the occurrence of deformation features in the form of debris bands suggest that the ice has been subjected to long-continued shear stress caused by the motion of a glacier (Benn and Evans, 2010; Lawson, 1979). The debris bands cross-cutting the buried glacier ice are comparable to those observed in the terminus zone of Stagnation Glacier on Bylot Island, where basal sediments were transported to the glacier surface through shear planes (Moorman and Michel, 2000).

Cations and stable water isotopes measured in the buried massive ice also support its glacial origin. The low cation content in the buried massive ice is statistically similar to that of the ice of glacier C-93. Although a slope of 6.0 was calculated between δD-δ¹⁸O of the buried ice, lower than the Global Meteoric Water Line (slope = 8; Craig, 1961), this was due to the small range of the data (δ¹⁸O: −34.4 ‰ to −33.4 ‰), which prevented the calculation of a reliable regression slope. The δD-δ¹⁸O values of the buried massive ice are similar to those of Pleistocene-age ice from the Barnes Ice Cap on the adjacent Baffin Island (median δ¹⁸O: −33.1 ‰; Zdanowicz et al., 2002). This is supported by the average D-excess of 6.6±2.5 ‰ for the buried massive ice, which is within the range of the values of Barnes Ice Cap (Lacelle et al., 2018) and glacier C-93 (5.2±5.8 ‰) on Bylot Island.

The average δ¹⁸O value of the buried massive ice ($-\mathrm{34.0}\pm \mathrm{0.4}$ ‰) is much lower than that of the other ice types sampled in the study area (Fig. 8b). The δ¹⁸O values of the buried massive ice are also lower than in the ice of Penny Ice Cap from the last glacial period (δ¹⁸O $\sim -\mathrm{31.3}\pm \mathrm{1.1}$ ‰; Fisher et al., 1998), but in the range of the late Pleistocene ice of Barnes Ice Cap (Zdanowicz et al., 2002). The late Pleistocene δ¹⁸O values on Barnes Ice Cap were 6 ‰–10 ‰ lower than the expected δ¹⁸O [modern – last glacial maximum] and interpreted to be the product of two factors: (1) temperature effect – cold climate during the last glacial period that caused a δ¹⁸O shift of 6 ‰–7 ‰ at that latitude (Fisher and Alt, 1985) – and (2) elevation effect – ice sourced from higher elevation on Foxe Dome (∼2200–2400 m), which caused an additional depletion in δ¹⁸O of local precipitation (Hooke and Clausen, 1982; Zdanowicz et al., 2002). Klassen (1993) suggested that the alpine glaciers on Bylot Island were larger but did not change much in elevation during the late Pleistocene; however, the δ¹⁸O values of our buried ice require ice sourced from elevations in the 2150–2350 m range and as such was likely sourced from the Foxe Dome (Lacelle et al., 2018). Mega-scale glacial lineations and streamlined landforms were mapped on the floor of Navy Board Inlet, Eclipse Sound and Lancaster Sound by De Angelis and Kleman (2007) and interpreted as a product of LIS ice streams. We proposed that the study area was most likely an area of confluence of LIS ice and local alpine glaciers during the LGM as an ice stream moved through Navy Board Inlet and onto south-western Bylot as proposed by Dyke and Hooper (2001).

5.2 Burial and preservation of late Pleistocene glacier ice in permafrost

The texture (sand and gravel), stratification and poorly sorted nature of the sediments (unit B) directly overlying the buried glacier ice suggest a glaciofluvial ice-contact deposit. Similar sediment characteristics were obtained by Fortier and Allard (2004) for glaciofluvial sands (angular grains: 75 %) located a few kilometres away from the study site. The abundance of angular grains (70 %) in unit B indicates little abrasion, which is consistent with transport by glaciofluvial water whereby grains were carried over short distances allowing little wear. Furthermore, the sharp and unconformable contact between the buried glacier ice and the overlying sediment suggests that thermal erosion caused by sediment-laden water flow affected the top of the ice of unit A. Our interpretation is similar to other studies that reported the burial of glacier ice by glaciofluvial sedimentation (Dallimore and Wolfe, 1988; Kaplanskaya and Tarnogradskiy, 1986; Ingólfsson and Lokrantz, 2003; Belova et al., 2008). The uppermost unit (C) is a diamicton, which has undergone reworking by non-glacial processes such as gravitational mass wasting (e.g. slumping and solifluction). In ice-contact environments, the sediment cover is subject to several cycles of subsidence and redeposition as the ice undergoes progressive and partial melting (Schomacker, 2008). The shape of the clasts found within this unit provides evidence that it has experienced active glacial transport, which is also supported by the occurrence of erratic clasts (gneiss) derived from distal bedrock located several kilometres from the study site. To summarize, large stagnant ice blocks could have been covered by glaciogenic sediments at or near the ice margins during the advance and stagnation of a glacier within the Qarlikturvik Valley. The burial of the ice occurred as meltwater streams deposited sediments in direct contact with glacier ice followed by the reworking and redeposition of supraglacial sediments, which formed a surficial cover of mud and sand, later affected by cryoturbations as is indicated by the incorporation of organic material dated at 885±15¹⁴C yr BP (1164 cal yr BP; 1σ range: 1058–1204).

A clear discontinuity in the δ¹⁸O and cations profiles is observed at the ice–sediment contact and in the overlying units (Figs. 7 and 8a). Strong contrasting profiles between buried ice and the overlying sediment are related to different sources of water and formation history. The supernatant water from the ground ice in the overlying sediments (units B and C) has much higher average δ¹⁸O values of $-\mathrm{17.6}\pm \mathrm{5.4}$ ‰. These values are comparable to δ¹⁸O values obtained from precipitation at Pond Inlet between 1990 and 1992 (−23.8 ‰; IAEA/WMO, 2016) and modern segregation ice which generally has δ¹⁸O values of −18 ‰ to −22 ‰ (Michel, 2011). The average concentrations of Ca²⁺, Na⁺ Mg²⁺ and K⁺ in units B and C are 10–86 times greater than those in the buried glacier ice. The substantial increase observed in the cation content from the sediment layers can be attributed to the great amount of mineral dissolution before the water freezes (Lacelle and Vasil'chuk, 2013).

The formation of these ice-cored landforms follows a typical sequence of events where the ice is first buried under a sufficient amount of sediments exceeding the active-layer thickness (Lukas et al., 2005; Benn and Evans, 2010). Beyond the active margins of many glaciers and ice sheets, there are broad areas of glacial deposits, till and outwash underlain by glacier ice. Gravitational and glaciofluvial processes are often identified as the most important processes of sediment deposition and reworking that lead to the burial of englacial ice (Benn and Evans, 2010). In basal zones of glaciers and ice sheets, in situ melting of stagnant debris-rich basal ice produces a supraglacial sediment cover (i.e. melt-out till) that accumulates on the ice, which inhibits its ablation. Both stagnant englacial and basal ice becomes progressively isolated from the upper active-flowing ice. These ice-cored landforms adjust to non-glacial conditions and their evolution is strongly linked with climate-driven processes (e.g. active-layer deepening due to warmer atmospheric temperature or active-layer detachment slides following heavy rains). In the Arctic, large blocks of buried late Pleistocene glacier ice were left undisturbed for several millennia owing to cold and dry climatic conditions that favoured permafrost aggradation following deglaciation. Such ice-cored landforms inherited from the partial melting of buried glacier ice are common in formerly glaciated permafrost regions across northern Canada and Siberia but their spatial distribution remains poorly known (Belova et al., 2008; Ingólfsson and Lokrantz, 2003; Kokelj et al., 2017; Dyke and Savelle, 2000; Lakeman and England, 2012). In this study, the deposition of a cover of coarse and well-drained sediment exceeding the average active-layer thickness of the area has probably been the most important factor limiting the melting of the ice. Furthermore, plant colonization and the development of a continuous vegetation cover with organic accumulation have changed the thermal properties of the ground, which reduced heat flow from the atmosphere to the permafrost during the summer and favoured heat loss during winter. This resulted in a thinner active layer (the permafrost table moved upward) and ground ice aggradation. The preservation potential of buried glacier ice on a millennial timescale following the glacial retreat depends on the complex interactions between climate, geomorphology and the physical properties of the sediment cover. Between 1999 and 2016, ground thermal regime monitoring in an intact low-centred polygon near the study site showed maximum active-layer depths varying between 0.3 and 0.7 m, while the buried ice is located >1 m below the ground surface. While the properties of the sediment cover had positive feedback on the long-term preservation of the buried ice, thaw slumping activity, and not atmospheric temperature per say, was a fundamental driver of its degradation by exposing the ice and accelerating its melting.

6 Conclusions

This study confirms that the permafrost of Bylot Island contains remnants of Pleistocene glacier ice that survived the last deglaciation. Geomorphic and cryostratigraphic observations along with the crystallographic properties of the massive ice suggest its englacial origin. Evidence that support the englacial origin of the massive ice are (1) sharp and unconformable upper contact between the ice and the overlying glaciofluvial sediments; (2) clear to whitish ice, with large crystals; (3) bubble-rich ice, with small gas inclusions (∼ mm) mainly located at crystal junctions; (4) occasional debris bands of sand and fine gravel cross-cutting older debris-free ice; (5) geochemical similarities with contemporary glacier ice. The geochemical data show a clear discontinuity at the buried ice–sediment contact as well as the very low cation content similar to that of modern glacier ice. An origin of massive ice from a Pleistocene glaciation is suggested based on the isotope data: the buried englacial ice is strongly depleted in heavy isotopes, similarly to regional Pleistocene glacier ice. Glacier-derived permafrost contains ice that predates the aggradation of the permafrost and represents unique environmental archives to reconstruct palaeo-environmental conditions at the study site. Although stable isotopes cannot yield information on absolute ages, these data show that the glacier ice body originated from an ice stream flowing from the Foxe Dome of the Laurentide Ice Sheet, which was subsequently buried and preserved in the permafrost of Bylot Island for thousands of years. The sedimentological data and interpretations presented in this study demonstrate that the first phase in the burial of the ice involved glaciofluvial deposition directly on the ice, which was followed by mass wasting and plant colonization. Knowledge regarding the occurrence, origin and preservation of buried glacier ice is of a great interest due to its potential impacts on the landscape stability upon melting. Since Bylot Island has experienced several periods of native and foreign Pleistocene glaciations and based on our findings, we propose that buried glacier ice should occur elsewhere on the island and probably on adjacent Baffin Island. In a context of climate change, active-layer deepening and increased activity of slope processes, such as active-layer detachment slides, thaw slumping and thermo-erosional gullying, will very likely expose buried glacier ice and initiate major landscape changes, with cascade effects on the ecosystems.

Appendix A

Figure A1 Grain-size distribution curves of the sedimentary units (two samples per subunit). The gravel fraction was removed prior to analysis.

Appendix B

To further investigate the cryostratigraphic characteristics of the ice, all samples were observed under X-ray-computed tomography (CT) scanning (Siemens SOMATOM Sensation 64). This technique relies on the calculation of the linear attenuation coefficient that measured the density of an object passed through an X-ray beam at different angles. A CT scan produces cross-sectional images (usually 512 by 512 pixels matrix) of an object where each pixel of the image is assigned an X-ray attenuation value (μ), also called a CT number. CT numbers are then standardized using the Hounsfield scale, where the radiodensity of water (μ_water) is arbitrarily defined as 0 HU (Hounsfield units) according to Eq. (1). Different shades of grey are assigned specific CT numbers to create the displayed image using a specific image processing software (Fiji) dedicated to DICOM (Digital Imaging and Communications in Medicine) images. In permafrost samples, unconsolidated sediments and rock (high-density minerals) appear white, as the attenuation of these materials is very high. Gases inclusions and water appear black, and other materials, such as ice, can have various shades of grey depending on their density.

$$\begin{array}{}\text{(B1)}& {\displaystyle }\text{HU}\text{value}={\displaystyle \frac{\mathit{\mu }-{\mathit{\mu }}_{\text{water}}}{{\mathit{\mu }}_{\text{water}}}}\times \mathrm{1000}\end{array}$$

This tool helps to refine cryostratigraphic characterization of permafrost cores as it can reveal characteristics otherwise difficult or even impossible to observe with the naked eye. It allows visualization and characterization of the internal components of the frozen sample, such as the ice, grain-size variations and gas inclusions. From a quantitative perspective, it has been used to segment images into regions of ice, gas and sediment in order to quantify the volumetric content of the sample scanned (Calmels et al., 2010; Dillon et al., 2008).

Appendix C

The typical method of measuring ice texture and fabric requires that an ice core or block is sliced into 1 cm thick vertical and horizontal sections with a band saw and then thinned to the desired thickness (<0.5 mm) following the standard microtoming procedure. Three photographs were taken of each sample at different angles (0^∘, 45^∘, 90^∘) under cross-polarized light in order to show the crystal boundaries as clearly as possible. The analysis is based mostly on qualitative observations, although statistics on crystal dimensions and geometry were obtained from a digital analysis of thin-section photographs. It should be noted that the values obtained from these parameters do not represent exact values and give limited information since a thin section is only a two-dimensional representation of the crystals (sectioning effect). A digital image analysis using ImageJ software allowed us to delineate and measure three geometrical parameters of each crystal. The area (mm²), the long axis (mm) and the circularity ratio determine the crystal texture (size and shape). The long axis is a measurement of the primary axis of the best-fit ellipse to an ice crystal. The circularity ratio (C) is a dimensionless parameter used to quantitatively describe the shape or outline form of a crystal. It is a function of the perimeter P and the area A:

$$\begin{array}{}\text{(C1)}& {\displaystyle }C={\displaystyle \frac{\mathrm{4}\mathit{\pi }A}{P^{\mathrm{2}}}}.\end{array}$$

As C approaches 1, crystals become more circular, and as it approaches 0, crystals become very narrow relative to their length. Therefore, the circularity ratio is a measure of the deviation of a crystal cross-section shape from a perfect circle.