Thermokarst lakes formed in buried glacier ice: Observations from Bylot Island, eastern Canadian Arctic

In formerly glaciated permafrost regions, extensive areas are still underlain by a considerable amount of glacier ice 15 buried by glacigenic sediments. Although the extent and volume of undisturbed relict glacier ice are unknown, these ice bodies are predicted to melt with climate warming but their impact on landscape evolution remains poorly studied. The spatial distribution of buried glacier ice can play a significant role in reshaping periglacial landscapes, in particular thermokarst aquatic systems. This study focuses on lake initiation and development in response to the melting of buried glacier ice on Bylot Island, Nunavut. We studied a lake-rich area using lake-sediment cores, detailed bathymetric data, remotely sensed data and 20 observations of buried glacier ice exposures. Our results suggest that initiation of deeper thermokarst lakes was triggered by the melting of buried glacier ice. They have subsequently enlarged through thermal and mechanical shoreline erosion, as well as vertically through thaw consolidation and subsidence, and they later coalesced with neighbouring water bodies to form larger lakes. Thus, these lakes now evolve as “classic” thermokarst lakes that expand in area and volume as a result of the melting of intrasedimental ground ice in the surrounding material and the underlying glaciofluvial and till material. It is 25 expected that the deepening of thaw bulbs (taliks) and the enlargement of Arctic lakes in response to global warming will reach undisturbed buried glacier ice, if any, which in turn will substantially alter lake bathymetry, geochemistry and greenhouse gas emissions from Arctic lowlands.

subsidence and ponding (van Everdingen, 1988;Kokelj and Jorgenson, 2013). In flat-lying terrains, thermokarst processes often result in the formation of numerous wetlands, ponds and lakes, thereby creating or modifying existing 'limnoscapes' (lake-rich landscapes) through thermal and mechanical erosional processes as well as thaw consolidation and subsidence beneath waterbodies (Bouchard et al., 2020;Grosse et al., 2013;Shur et al., 2012). The formation and growth of these lacustrine ecosystems have important effects on the hydrology, ecology, biogeochemistry and geomorphology of affected landscapes. 5 Shoreline erosion may affect key biogeochemical processes within these lakes, such as the burial of organic matter in sediments, and its degradation and release as greenhouse gases (GHG; CO2 and CH4) to the atmosphere (Matveev et al., 2016;Vonk et al., 2015;Wik et al., 2016). Lake basin morphology also influence GHG flux patterns during the open-water season by affecting the mixing regime (Prėskienis et al., 2021;Hughes-Allen et al., 2021). 10 The extent to which permafrost degradation occurs is dependent on the distribution and abundance of ground ice. In formerly glaciated permafrost regions, extensive areas still contain a considerable amount of glacial ice buried beneath glacigenic sediments (Belova, 2015;Coulombe et al., 2019;French and Harry, 1990;Ingólfsson and Lokrantz, 2003;Kanevskiy et al., 2013;Swanger, 2017). Remnants of buried glacier ice remain stable as long as the ground temperature is below freezing, and the active layer thickness (i.e. depth of annual thawing) does not exceed the depth to the massive ice body (Shur, 1988). The 15 transition from partially ice-cored moraines with isolated buried ice blocks to ice-free landscapes may post-date the major glacial retreat by thousands to millions of years (Bibby et al., 2016;Coulombe et al., 2019;Lacelle et al., 2007;Swanger, 2017). The persistence of thick beds of buried Pleistocene glacier ice in contemporary permafrost environments indicates that deglaciation is still incomplete (Astakhov and Isayeva, 1988;Everest and Bradwell, 2003). The broad distribution and the substantial amount of ground ice in permafrost-preserved glacial landscapes make it highly vulnerable to disturbances, such 20 as thermokarst, under the ongoing climate warming. Some of these landscapes are now experiencing climate-driven renewed deglaciation leading to post-glacial landscape evolution that transformed sub-Arctic environments millennia ago. For example, on hillslopes, the thawing of permafrost terrain underlain by remnants of glacial ice triggered mass wasting processes, such as retrogressive thaw slump and active layer detachment slides Rudy et al., 2017). In flat or very gently sloping terrain, formation and evolution of ponds and lakes are typically associated with the melting of intrasedimental ice, 25 such as ice wedges and segregation ice (Bouchard et al., 2017;Grosse et al., 2013). These lakes tend to be shallow, with deeper central pools (~ 2-5 m) and shallow littoral shelves (~ 1 m), or shallow flat-bottomed basins (French, 2010;Grosse et al., 2013;Hinkel et al., 2012). It is generally recognized that numerous Arctic lakes were formed during deglaciation in depressions left by in-situ melting of stagnant blocks of glacier ice (i.e. kettle lakes). However, very few studies have linked lake inception to the thawing of sediments containing glacier ice that had been buried and preserved in permafrost for decades to millennia 30 (Dyke and Savelle, 2000;Henriksen et al., 2003;Jorgenson and Osterkamp, 2005;Worsley, 1999). As a result, there is little information on the spatial distribution and abundance and evolution in modern paraglacial and periglacial environments on these glacial thermokarst lakes (Jorgenson and Osterkamp, 2005).
The Quaternary geology of the eastern Canadian Arctic records several glaciations by ice sheets and local mountain glaciers, which means that the landscape stores vast amount of buried glacial ice and there is potential for significant post-glacial landscape change associated with the ablation of buried ice. This study builds on the findings of Coulombe et al. (2019) conducted on Bylot Island (Nunavut), where blocks of stagnant ice became separated from an ice stream flowing from the Foxe Dome of the Laurentide Ice Sheet and subsequently buried by aggradation of glaciofluvial sands and gravels at the 5 margins of the receding glacier. Subsequent neoglacial cooling resulted in widespread permafrost aggradation and preservation of this glacial ice. Here, we investigate whether the melting of relict glacial ice preserved in the permafrost has triggered the inception of thermokarst lakes on Bylot Island over the Holocene. We hypothesized that remnants of buried glacier ice in lowlands slowly melted during the Holocene, which created depressions that formed glacial thermokarst lakes. The objectives of the present study were therefore (1) to examine the origin of twenty-one lakes in the lower reach of a glacial valley on Bylot 10 Island by investigating the geomorphological properties, stratigraphic profiles, and spatial distribution of this group of lakes and (2) to develop a conceptual model of lake inception and evolution triggered by the delayed melting of buried glacier ice.

Study area
The study area is in the Qarlikturvik Valley (73˚09' N, 79˚57' W) on the southwest plain of Bylot Island, in the Canadian Arctic Archipelago (Fig. 1). The landscape was glaciated several times in the late Quaternary by both local mountain glaciers 15 and the Laurentide Ice Sheet (LIS; Klassen, 1993). The study area was most likely a confluence zone between LIS ice and local alpine glaciers with glacier ice flowing out of major valleys but Laurentide ice flowing into the southern plain and up the valleys (Dyke and Hooper, 2001;Lacelle et al., 2018). The maximum extent of the LIS is clearly outlined by the Eclipse moraine, a major moraine system across the outer coastal mountains of Bylot Island and parts of adjacent Baffin Island (Klassen and Fisher, 1988). Today, Bylot Island remains 40% glacierized as numerous valleys and piedmont glaciers still flow 20 from the local ice cap and terminate in lowlands underlain by sedimentary rock of Cretaceous-Tertiary ages (Dowdeswell et al., 2007). The Qarlikturvik Valley is one of the many U-shaped glacial valleys with ice-rich sediments dating back to the Late Pleistocene and Holocene, which are highly susceptible to thermokarst (Fortier and Allard, 2004;Bouchard et al., 2020). The valley contains abundant and diverse water bodies, including a proglacial river, lakes, through and polygon ponds, small streams, and gullies (Godin et al., 2014;Muster et al., 2017;Prėskienis et al., 2021). This lake-rich valley is an ideal location 25 to study buried glacier ice dynamics in thermokarst-affected and glaciated permafrost landscapes. With glaciers ending within the continuous permafrost zone, this typical glacial valley represents a small geosystem shaped by proglacial, paraglacial and periglacial processes (Church and Ryder, 1972).
In the Qarlikturvik Valley, mounds of reworked till and ice-contact stratified sediments mark former positions of the glacier 30 front (Fig. 1c). The earliest postglacial radiocarbon date from marine shells retrieved from marine clays is 11,331 cal yr BP (IntCal20), suggesting that the valley was partially ice-free by this time (Allard, 1996). About 2-3 meters of ice-rich Quaternary https://doi.org /10.5194/tc-2021-302 Preprint.  silt and sand derived from aeolian deposition, interstratified with peat, overlies ice-poor glaciofluvial outwash deposits (Fig. 1c;Fortier and Allard, 2004). Syngenetic ice wedge growth has created an extensive polygonal patterned ground. Thermokarst is an active landscape change mechanism operating in the valley, as demonstrated by the abundance of lakes, thermo-erosional gullies and thaw slumps within the study area (Bouchard et al., 2020;Fortier et al., 2007;Godin et al., 2014). Previous work in the area have examined various aspects of thermokarst lake dynamics such as GHG exchanges (Bouchard et al., 2015a;5 Prėskienis et al., 2021), photochemical and microbial decomposition of organic matter , microbial diversity (Negandhi et al., 2014), nutrient inputs from the goose colony (Côté et al., 2010) and methylmercury , as well as lake development in syngenetic ice-wedge polygon terrain (Bouchard et al., 2020).
The mean annual air temperature at Pond Inlet for the 1981-2010 normal is −14.6°C, which is 0.5°C higher than the previous 10 1971-2000record (Environment Canada, 2021. The mean annual precipitation for the 1981-2010 period was 189 mm per year, with rainfall representing 91 mm. Bylot Island is located within the continuous permafrost zone. The regional thickness of permafrost was estimated to be at least 100-500 m extrapolated from thermal measurements on nearby Somerset and Devon Island (Smith and Burgess, 2002). On average, the active-layer thickness varies between 0.4-0.7 m in peaty and silty soils, to more than 1 m in drained unvegetated sands and gravels (Allard et al., 2016). Thawing and freezing indices averaged  2010 period) 473 degree-days above 0 °C and 5736 degree-days below 0 °C, respectively (Environment Canada, 2021).

Materials and methods
Two spatial scales were used to investigate the role of buried glacier ice in the formation and evolution of thermokarst lakes.
First, we focused on the Qarlikturvik Valley (~ 75 km 2 ), where buried Pleistocene glacier ice has been found in permafrost (Fig. 1;Coulombe et al., 2019). We examined the morphology and conducted bathymetric surveys of 21 lakes and analysed 20 lake sediment cores from three of these lakes to infer probable lake origin. Then we examined the spatial distribution of lakes on the broader coastal plain of Bylot Island to link the extension of former local and regional glaciations to lake distribution.

Landforms, surficial deposits and lake mapping
We used contemporary high-resolution GeoEye satellite imagery (2010, pixel = 0.5 m), historical aerial photos (1961,1982) and ArcticDEM data (pixel = 2 m) to map glacier frontal positions, Quaternary surficial deposits, and landforms in the 25 Qarlikturvik Valley. A Sentinel-2 image mosaic (2016, pixel = 10 m) of the southern plain of Bylot Island served as the basis for mapping the water bodies outside the valley (Copernicus, 2016). We also used the Google Earth Engine Timelapse dataset (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019) to visually assess terrain change and sediment accumulation at the glacier terminus based on Tasseled cap trend analysis calculated from a stack of Landsat satellite images (Fraser et al., 2012;Gorelick et al., 2017). Data processing and analyses were performed using QGIS (v.3.16;QGIS Development Team, 2021). To extract all water bodies, we used the 30 reflectance properties of water in the Green and NIR bands (McFeeters, 1996). Because water bodies have high 'Normalized Difference Water Index' (NDWI) values, a simple thresholding technique was used to isolate most water features. Lake shorelines were extracted as vector data and converted to polygon topology. Lakes smaller than 1000 m 2 were automatically removed from the analysis to exclude ice-wedge polygon ponds and collapsed ice-wedge troughs filled with water.

Distribution of lakes in the valley and the southern plain of Bylot Island
We analyzed the spatial distribution of lakes, and their position in relation to past positions of glaciers on the island. To map 5 the density of lakes in the Qarlikturvik Valley and the broader southern plain of Bylot Island, a kernel density estimation was performed using the 'spatstat' package in R (v. 3.5.3;Baddeley et al., 2019;R Core Team, 2021). Input for kernel density came from lake centroids obtained from the vector polygon, which were calculated automatically in R as the geometric center of the lakes. We defined the extent as all areas of the Qarlikturvik Valley and the broader southern coastal plain of Bylot Island, excluding the bedrock outcrops, steep slopes (> 5°), glaciers and outwash plains. To analyse lake distribution, we also 10 performed a clustering analysis using the inhomogeneous pairwise correlation function (ginhom), which accounts for spatial inhomogeneity in lake locations (quadra test; p = 0.001). This technique allows distinguishing between dispersed (R > 1) and clustered (R < 1) spatial patterns by comparing the observed point patterns against the expectation for a randomly distributed sample population (CSR model), which assumes that the objects can be distributed anywhere in the region of interest.

Lake morphology in the Qarlikturvik Valley 15
Detailed bathymetric data were collected for 21 lakes across the valley using a Humminbird 859XD Sonar with a built-in global positioning system. Lake bathymetric surveys were conducted with an inflatable boat when the lake was free of ice (June to August). Geographic location and water depth were recorded each second along transect lines that were spaced at approximately 5 to 25 m (depending on lake size) to entirely cover the lake. Some uncontrolled conditions have degraded the accuracy of the survey, such as the presence of littoral vegetation and waves, especially in shallow lakes. Depth and location 20 data were imported into QGIS for visualization and additional processing. Initial processing included the removal of spurious data points (outliers) such as single-point depths located substantially above or below the general depth of lake-bottom. We used a spline algorithm to generate an interpolated surface from the individual depths. Ground penetrating radar (GPR) surveys with 50 Hz antenna were conducted across frozen lakes to investigate the lake-bottom morphology (see Supplementary material S2 for further details). For each lake, we also calculated the area, perimeter, elongation ratio (long axis/short axis), 25 and shoreline development or DL from the digitized shoreline polygons in order to compare lake metrics and determine if they can be used to discriminate between thermokarst lakes formed in ice-wedge polygon terrain and deeper thermokarst lakes formed by the melting of buried glacier ice. The shoreline development ratio (DL) is a standard measure of the complexity of the shoreline, which is the ratio of the length of the shoreline of a lake (i.e. its perimeter) to the circumference of a circle of area equal to that of the lake (Equation 1; Hutchinson, 1957).  DL for a perfect circle is 1.0, and its value increases as the shape of the lake surface deviates from that of a circle, indicating the shoreline is more dendritic or irregular. With their circular shape, glacial thermokarst lakes should have low complexity values (~1) whereas thermokarst lakes formed through thermal and mechanical shoreline erosion should be more irregular and have values >1. A highly indented shoreline may also indicate coalescent lakes formed by shoreline expansion. The elongation ratio can be used to separate round features (perfect circle = 1) from elongated ones (>1). Correlation between shoreline 5 morphology variables and basin morphometry (maximum depth) were tested using the non-parametric Kendall tau rank correlation for non-normally distributed data. All statistical tests were run in the open-source software R (R Core Team, 2021).

Lake sediments and vertical structure
We selected three nearby lakes (G, K and L) exhibiting different morphometry to compare their stratigraphic profiles and vertical profiles of temperature and dissolved oxygen. According to the bathymetric surveys, lakes K and L are the deepest 10 lake in the valley (max. depth = 12.2 m and 11.7 respectively), but lake K is smaller than lake L. We also sampled lake G (max. depth = 4.1 m) as lake bottom imagery revealed submerged ice-wedge polygons (~1 m depth) and degraded ice-wedge troughs in shallower lake G, which confirmed that this lake (G) is evolving through the melting of permafrost intrasedimental ice and ice wedges as revealed by subsiding and eroding shores (see video supplement in Bouchard et al., 2020). Two sediment cores of 109 cm and 114 cm were collected in spring 2015 from lakes G and K, respectively. The cores were collected through 15 a 2-m thick ice cover using a 7-cm diameter handheld percussion corer (Aquatic Research Instruments), sealed, and returned to the laboratory where they were stored in the dark at 4˚C. Coring occasionally caused minor deformations to the sediments owing to friction and pressure along coring tubes. Both cores were observed under X-ray computed tomography (CT), allowing to visualize and reconstruct the internal structure (2D and 3D) of the cores. Details on the CT scanning procedure are provided in supplementary material S1. Facies were identified based on visual inspection and physical properties, including sedimentary 20 structures, grain size, colour, and density. Percentage dry weight was determined for all samples (drying overnight at 105°C).
Organic matter content was determined by weight loss (loss-on-ignition, LOI), following a combustion of dried samples at 550°C for 4 h (Heiri et al., 2001). Sediment grain size was measured in triplicates using a Malvern Mastersizer 2000 and Hydro2000G liquid handling unit. Bulk sediment and fossil plant fragments were radiocarbon-dated by accelerator mass spectrometry (AMS) at Keck Carbon Cycle AMS Facility (University of California, Irvine, CA, USA). Calibrated ages (cal yr 25 BP) were calculated using "CALIB 8.2" (Stuiver et al., 2021; IntCal20 dataset, Reimer et al., 2020). In the case of lake G, facies are described in more details based on other proxies, such as organic content and fossil diatoms (Bouchard et al., 2020;Figs. 4 and 5).
We also profiled the same two lakes (G, K) as well as lake L in late winter under the ice cover (early June 2015) and during 30 the ice-free period (July and August) to examine differences in water temperature and dissolved oxygen (DO) between lake types, and discuss the effects of water depth and vertical structure on GHG emissions and fish habitats. Discrete profiles were measured manually from lakes G and K with a ProODO profiler (YSI Inc.), while submersible temperature loggers (Vemco Minilog-II-T installed at 2, 4, 6, 8 and 10 m depth) and DO loggers (PME MiniDOT; 2 and 9 m depth) were installed in lake L to record annual cycles of stratification (1h frequency), from which profiles were selected to match the discrete profiles obtained from the other two lakes. Sensor specifications can be found in Prėskienis et al. (2021).

Spatial distribution of lakes on the southern plain of Bylot Island 5
Using remote sensing classification, we detected 845 lakes larger than 1000 m 2 within the study area (total lake area reaching 14 km 2 over ~ 1700 km 2 , or 0.8% of the area), of which 189 lakes (totalling 1.6 km 2 ) are in the Qarlikturvik Valley (~ 122 km 2 ; 1.3% of the area). The spatial distribution of the lakes showed a significant aggregation pattern in both Qarlikturvik Valley and the southern coastal plain of Bylot Island ( Fig. 6c and d). Patterns of distribution emerge on this level with higher densities, with 10 to 85 lakes per km 2 , detected nearby mounds of ice-contact deposits or in areas of unvegetated moraine in 10 front of glaciers C-79 and C-93 (Fig. 6c). A third group of lakes is also observed on the plateau bordering glacier C-93.
According to the point pattern analyses, the lakes in the valley show significant clustering in short distance (far above the 95% confidence envelope; r < 0.85 km) and a regular distribution further away (r > 0.85 km; Fig. A2). On the southern plain of Bylot Island, the highest densities occur directly in front of contemporary glaciers and within the extent of local mountain glaciations, with up to 40 lakes per km 2 (Fig. 6d). The observed points (lake centroids) show considerable clustering at small 15 distances less than about 3.3 km, but show regularity beyond about 4 km (Fig. A2). In addition, we analyzed the formation of lakes in front of glaciers C-93 and C-79 combining a detailed lake inventory with records of glacier retreat since the Little Ice Age (LIA; 140 cal yr BP; 120 14 C BP; Klassen, 1993). By 2020 the glacier front had retreated ca. 2 km since the end of the LIA, and 383 new lakes have developed in the ice-free zone over a period of about two centuries (Fig. 6a). The TC trend analysis method reveals sediment accumulation at the front of the receding glaciers between 2000 and 2019 as represented by 20 red colours (drier and unvegetated areas) on TC images. This shows the active burial of glacier ice at the front of glaciers C93 and C79, thereby providing a modern analogue for the past burial of ice when the glacier was several km further down-valley.

Morphology of lakes in the Qarlikturvik Valley
In the valley, we identified two groups of lakes according to their depth range and lake-floor morphometry (Fig. 2). The first group of lakes (n=8) stands out by their greater depths and sizes, and in some cases, the presence of multiple sub-basins. The 25 maximum measured depths recorded in these lakes range from 8.4 to 12.2 m. Most of these lakes are characterized by a relatively deep central lake basin surrounded by shallower areas, ranging between 0.5 and 1.0 m (mean depth= 0.6 ± 0.4 m; supporting information Table S3). Three lakes (L, N, U) have two or three steep-sided and confined sub-basins that are surrounded by a relatively shallow marginal platform. The GPR profiles indicated that these deeper lakes usually have smoother lake bottom, whereas lake L also exhibits an irregular lake floor morphology in the shallowest areas (Fig. S2). The 30 bathymetric map of lake L also revealed a deeper depression in the lake bottom that is aligned with the small lakeside thaw slump exposing buried glacier ice (Figs. 2 and 3). The bottoms of lakes K and L are respectively 5.5 m and 5 m below the current sea level. In the valley (zones 1 and 2), the deeper lakes are located near mounds of ice-contact glaciofluvial deposits.
By contrast, the second group of lakes (n=13) showed markedly different characteristics (Fig. 2). At the lake scale, these shallow lakes have relatively flat and homogeneous lake beds with a deeper central basin surrounded by shallower nearshore zones (< 2 m deep), where submerged polygons are clearly visible (see the video supplement in Bouchard et al., 2020). These 5 lakes have maximum depths ranging between 1 and 4 m, with mean depth reaching 1.4 ± 0.7 m across their platform (Fig. 2).
They are characterized by irregular shorelines, which generally follow the deep troughs caused by the melting of ice wedges from their tops. Despite the above-stated differences, lakes from both subgroups present similar shapes and shoreline characteristics ( Fig. A1) as their morphometric properties (area, perimeter, elongation ratio, complexity) were not significantly different (Mann-Whitney-Wilcoxon test, p > 0.05). In addition, mean or maximum depth did not show any significant 10 correlation with the other morphometric variables (all p > 0.1).

Vertical structure of water column of lakes
Profiles done in late winter under the ice cover of lakes K and L (deep; group 1) and lake G (shallow; group 2) showed an inverse thermal stratification with bottom water temperature reaching 0.6°C (lake G), 1.4°C (lake K) and 2.0 °C (lake L) (Fig.   5). DO was much lower in Lake G (13% saturation below the ice cover, decreasing to ~1% near sediment at 4 m depth in 15 2015) than in lake K (43% below the ice, <1% below 6 m depth in 2015) and lake L (70% below the ice, 26% at 9 m depth in 2016, but down to < 2% by mid-May in 2019; no data in spring 2015). Quickly after the ice cover melted at the beginning of July, the water column became weakly stratified in all three lakes, and warmer in lake G (above 8°C at the surface, as compared to ~5°C at the surface of the deeper lakes). By then, the water column already showed signs of DO depletion in the deeper lakes (K, L) but the summer stratification period was short, lasting for about a month in the larger lake L. At the bottom of 20 lake L, oxygen depletion occurred as soon as stratification established from the beginning of July and decreased down to 84% of saturation by mid-August 2015 (65% in 2018), until the autumnal turnover kicked saturation level up again (unpubl. data).
While only weak hypoxia was encountered at the bottom of lake L in late summer, anoxia was reached in the hypolimnion of lake K (Fig. 5). On its side, lake G was generally well-mixed during the open-water period, but presented weakly stratified periods during warm and calm days (e.g. on 3 August 2016). Early August profiles indicate that the entire water column was 25 above 13°C in the shallow thermokarst lake G, while the surface of lake K and L was slightly colder (respectively ~10°C and 11°C in 2016).

Sediment profiles of lakes
Four distinct lithofacies or units, labelled from core bottom to top, were identified in lake K (max. depth: 12.2 m) based on visual analysis of CT-scan images and field description: A) sandy silt and gravel with interspersed peat/organic debris (114-30 silt and fine sand) and gravel with scattered organic material, which was dated near its top (95 cm) at 3531 cal yr BP (3330 14 C BP; Fig. 4). The coring operation did not reach the bottom of this unit, so its total thickness is unknown. Compared to other units, it has a higher mean density (2.0 to 2.5 g cm -3 ), typical of dominantly mineral material. Unit B consists of dry and fibrous organic-rich material (peat) with fine sand and silt. Unit C consists of laminated to massive sandy silt containing very sparse and fine gravels, which is massive in its uppermost 5 cm. This unit also displays sharp lower and upper contacts, and includes 5 some deformation structures, caused by the coring operation (layers bended downward near the coring tube walls). Bulk sediment near the top of this unit (73 cm) has been dated at 4036 cal yr BP (3700 14 C BP ; Fig. 4). The uppermost unit (D) is composed of laminated gyttja that grades upwards into soft and loose gyttja. The CT-scan image along with the LOI profile show that the mineral input steadily decreases towards the middle of unit B. Light-coloured thin laminae of silt (0.3 and 0.9 cm) are common in the upper part of the sequence. The upper section of this unit (50-0 cm depth) is less compact compared 10 to deeper sediments (> 50 cm depth) and it has a high water content that becomes dryer towards the bottom of the unit (moisture content decreases from about 80% at the surface to 15% at the bottom; Fig. 4). A similar trend is also observed for the organic content as it decreases from ~15% at the surface to near 1% in the lower portions of the core (Fig. 4).
Three lithofacies were identified in lake G (max. depth: 4.1 m): A) organic-poor sandy silt (109-80 cm); B) organic-rich sandy 15 silt interstratified with peat (80-10 cm); C) sandy silt gyttja (organic lacustrine mud; 10-0 cm; Fig. 4). Unit A is composed of sand and gravel with scattered, cm-scale peat and organic debris. The lower section of this unit has a higher density (~ 2.0 g cm -3 ) and a low organic content (mean: 7.7 % ± 3.7) compared to the upper sections of the core. The bottom deposit of unit A contains organic matter older than 5507 cal yr BP (4805 14 C BP), based on dating of a wooden fragment at a depth of 108 cm (Bouchard et al., 2020). Unit B consists of medium to dark brown peat, dated to 4567 cal yr BP (35 cm;4070 14 C 20 BP), and interbedded with mm-to cm-thick silt and sand laminations, and with gradational upper and lower boundaries.
Throughout the unit, sand-silt layers (aeolian) are roughly interbedded with layers of organic detritus as recorded by sifts in organic matter contents. Towards the upper boundary, unit B progressively grades into dark brown gyttja (unit C). The upper part of unit C is faintly stratified, as the organic-rich material becomes regularly interspersed with silty material. These silty laminae (0.3 to 0.9 cm) are visually distinguishable by their light-grey colour on the CT-scan image and their higher density. 25 The top sediments have the highest water (70.6 % ± 6.1) and organic matter contents (20.3 % ± 5.2) compared to the bottomed units. Contrary to lake K the boundaries between the units are diffuse. The bottom of this unit (bulk sample collected at 10 cm) yielded an age of around 2061 cal yr BP (2100 14 C BP).

Glacial origin of lakes 30
The presence of deep lakes and numerous thaw slumps in Qarlikturvik Valley indicates the delayed melting of several bodies of buried glacier ice as compared to the Holocene glacier retreat. The ice-free zones of Bylot Island are therefore still strongly influenced by its glacial legacy given the presence of late Pleistocene-age glacier ice buried in the permafrost in Qarlikturvik Valley (Coulombe et al., 2019) and in other valleys and coastal plains of the island (Klassen, 1993;Moorman and Michel, 2000).These ice-cored landforms have been adjusting to non-glacial conditions and their evolution is strongly linked with climate, geomorphological processes and local topography. Contemporary thawing of ice-rich permafrost is now triggering rapid and strong cascading effects on Arctic landscapes, including reshaping local topography and mobilizing previously 5 frozen materials to downstream aquatic ecosystems such as rivers and lakes (Kokelj et al., 2020;Vonk et al., 2015).
The bathymetric data revealed the coexistence of two groups of lakes with different morphological characteristics. We also found that different sedimentary facies were present in the cores collected from each group, suggesting different origins and evolutionary conditions. Sixty-two percent of the lakes in Qarlikturvik Valley are shallow (~2-4 m) and relatively flat at their 10 bottom, with a central deeper pool. This group of shallow lakes displays maximum depths very similar to those in 'classic' thermokarst lakes (~1-4 m deep), excluding lakes formed in Yedoma-type permafrost (Hinkel et al., 2012;Kanevskiy et al., 2014). Their depth is controlled by the depth of the syngenetic ice wedges, and by the amount and distribution of ground ice in the substrate (Grosse et al., 2013). Their maximum depths are also in accordance with the thickness (~2-3 m) of the peaty silt sequence forming the surrounding material, which developed during the Late Holocene (Fortier et al., 2006). Subsequent 15 thermokarst evolution in those basins is not likely to result in substantial subsidence of the lake or basin floor, which can be explained by the stratigraphy of the study area. In lake G, the lowermost unit (A) represents the top of a glaciofluvial unit. The lake has been slowly expanding in the frozen silt-peat terrace, and thawing has reached the underlying glaciofluvial sand (Bouchard et al., 2020). The intermediate unit (B) includes a layer (~35-55 cm) of convoluted horizons, which is absent in lake K, and likely originates from collapsed bank material in response to thermo-mechanical erosion processes or broken horizon 20 due to lake bottom subsidence after ground ice melting. The sediment profile from lake G is very similar to those found in lakes initiated by the degradation of ice-wedge and intrasedimental ice. These lakes typically have a transitional organic-rich layer containing peat derived from permafrost thawing and subsidence, underlying a layer of laminated organic-rich lacustrine mud (Biskaborn et al., 2013;Bouchard et al., 2017;Farquharson et al., 2016;Murton, 1996). Such an interpretation is further supported by the fossil diatom record investigated in lake G ('Gull Lake'; Bouchard et al.,  The other lakes (38%) stand out by their notably deeper basins (~ 5-12 m), and in some cases the presence of multiple subbasins (e.g., lake L). Owing to the size of the lake depressions and their location just next to mounds of ice-contact deposits, 30 these deeper lakes were primarily formed by the melting of buried glacier ice. This interpretation is supported by the presence of two exposures of glacier ice revealed by small lakeside slumps, which also indicates that the shoreline of lake L is still icecored by glacier ice (Fig. 3). Côté et al. (2010)  greater than 5 m (Fig. 6d; depth range: 5-21 m). The stratigraphic sequence in Qarlikturvik valley suggests that the amount of intrasedimental ice in the sediment layers is not sufficient to create lakes with depths reaching up to 12 m. The sediments underlying the silt and peat sequence are glaciofluvial sands with a gravimetric content of 20 to 25%, which is too low to allow for significant subsidence. Given the gravimetric water content measured in this unit (71.5 ± 12.7%; Fortier et al., 2004) and typical porosity (0.3) and density (2.0 g/cm 3 ), the sand and gravel unit is likely undersaturated in ice. Therefore, there should 5 be little or no further subsidence. The highest estimates for potential thaw settlement varied from 1.2 m in gravel and sandy soils to 2.4 m in very ice-rich silty and peaty soils. In addition, the deepest sections of some of these lakes, for example lakes K and L, are ~ 5 m below current sea level, indicating burial in a glaciomarine/glaciofluvial environment followed by isostatic uplift.

10
The glacial origin of several lakes in the valley is further corroborated by the analysis of the sediment core collected in lake K (12.2 m), which differs significantly from the one obtained in lake G (Fig. 4). In lake K, four depositional stages were inferred from the sediment profiles. The inception of lake K began with the collapse of supraglacial material during melt-out of stagnant glacier ice, which resulted in re-sedimentation of sand and gravel from glaciofluvial and mass movement deposits into a forming basin. Inclusions of fibrous peat in unit A and the ~10 cm darker layer (Unit B) of organic debris and inorganic 15 material (mostly silt and sand) were probably derived from surficial vegetation on upland surfaces washed in as the lake basin developed. This basin was then filled by sandy mud (Unit C) deposited by the combined action of meltwater streams and aeolian activity. The most prominent features of this core from lake K are the sharp boundaries (Units B and C), indicating a shift in the depositional conditions, which are not reflecting gradual deposition within a stable lake floor (Henriksen et al. 2003). Similarly to the upper core section in lake G, the upper unit D includes recent lacustrine sediments composed of organic-20 rich mud (gyttja). This material was deposited within deeper, calmer waters where fine material can settle down. The laminations reflect variability in minerogenic inputs and are likely related to 1) terrestrial runoff; 2) aeolian activity, or 3) discrete episodes of bank collapse or slumping where amounts of silt were released to the lake. Similar stages of sedimentation were identified from glacial thermokarst lake basins, including one lake with multiple sub-basins (up to > 14 m deep) in the continuous zone of permafrost in northern European Russia, where buried glacier ice has survived for ca. 80 000 years from 25 the last glaciation (Henriksen et al., 2003). Reversal of ages in the core suggests that older organic matter was also washed into the lake during the mid-Holocene, causing abnormally old dates in basal core sediments, a common dating problem in high-latitude lakes (Bouchard et al., 2017;Wolfe et al., 2004).
Moreover, our spatial analysis demonstrated that lake distribution is strongly linked to the maximum and recessional positions 30 of local mountain glaciers and the LIS in both Qarlikturvik Valley and the southern plain of Bylot Island. This is particularly evident in the valley, as shown by the three well-defined lake clusters (Fig. 6c). For most of the southern lowland coastal plains, moderate to high point densities were also encountered within the extent of the LIS. There is also a notable increase in lake density close to the contemporary ice margin, which we interpret to be the result of a relatively recent and continuous deglaciation process. We found that, even after accounting for landscape heterogeneity (i.e. high slope gradients, bedrock exposures), the lakes are still far more clustered when compared to a random spatial distribution. As a result, additional clustering mechanisms control the location and development of lakes, such as the presence of patches of buried glacier ice.

Conceptual model of thermokarst lake development in buried glacier ice
Several studies have described the stages of thermokarst lake development in different permafrost environments, such as in 5 Yedoma environments (Morgenstern et al., 2011;Shur et al., 2012), ice-rich cryogenic mounds (Calmels et al., 2008), or icewedge and intrasedimental ice (Bouchard et al., 2020). Based on the geomorpholology of the deeper lakes and lake sediment profiles of lake K, we developed a four-stage model of glacial thermokarst lake formation and evolution in the specific stratigraphic context of buried glacier ice within the study area (Fig. 7).

10
Stage 0: Initial conditions. During the last glacial maximum, the LIS and local ice caps covered much of the Qarlikturvik Valley, and many outlet glaciers were channelled in major valleys and terminate in the lowlands.
Stage 1: Burial of glacier ice. Beyond the active margins of local glaciers or the LIS, wide areas of glacier ice were likely buried in situ by glacigenic sediments transported and reworked on top of an active or stagnant glacier margin by mass 15 movement and meltwater (Fig.7a). The burial of glacial ice can still be seen today at the margins of many glaciers on the island (Appendix A3). Progressively, stagnant ice blocks became isolated from the upper active flowing ice. On Bylot Island, bodies of glacier ice were preserved at various places in the outwash plain, in mounds of ice-contact stratified drift and moraines (Coulombe et al., 2019;Klassen, 1993). Interpretations of the sedimentary sequence overlying the buried ice studied in the valley indicated that the burial of the ice involved glaciofluvial deposition directly on the ice, which was followed by plant 20 colonization. In some cases, glaciofluvial sands and gravels were also covered by colluvial sediments as debris were transferred away from topographic highs by mass movements and meltwater (Coulombe et al., 2019). Preservation of the ice for several millennia was possible because the sediment cover became sufficiently stable and reached or exceeded the active layer thickness, but also because the neoglacial climatic conditions during the second half of the Holocene were conducive to syngenetic permafrost aggradation following glacier retreat (Fortier et al., 2006;Fortier and Allard, 2004;Bouchard et al., 25 2020).
Stage 2: Initiation of buried ice melting. Melting of the upper glacier ice and formation of a depression begins when the deepening active layer reaches the glacier ice (Fig. 7b). Local factors such as topography, thickness of the insulating layer, snow accumulation and water pooling in pre-existing depressions, as well as thermal properties of soil all play a role in 30 initiating ice melting. We suggest three scenarios of initial pond inception following the burial of glacier ice: (1) ponds may have formed quickly in the proglacial environment or later during deglaciation as water pooled in these pre-existing topographic depressions caused by the uneven ablation of the glacier surface (i.e. differential melting); (2) the ice was first https://doi.org /10.5194/tc-2021-302 Preprint.  buried under a thin cover of glacigenic sediment, which was close to the active layer thickness. The sediment cover was thick enough to slow down the melting rate of underlying ice without completely preserving it, which allowed a gradual melting of the ice creating small depressions; (3) the ice was buried under thick sediment cover, acting as a barrier to heat transfer, and preserving the ice in the long term. However, unusual warm and wet conditions have periodically caused the active layer to deepen considerably, initiating the thawing of the underlying ice and creating new depressions. Following the burial of the ice, 5 the uneven ablation of the glacier surface produced an irregular topography of ridges and mounds. Since buried glacier ice is still present in the study area, thaw slump activity is thought to have been a fundamental driver of its degradation by exposing the ice and accelerating its melting (Coulombe et al., 2019).
Stage 3: Lake inception and syngenetic ice growth. Episodes of warming sufficient to cause degradation of the existing 10 permafrost and the buried glacier ice in the valley triggered thermokarst lake initiation (Fig. 7c). Summer-melt layers from the Agassiz Ice Cap and Greenland ice sheet provide a robust record of warmer decades since the mid-Holocene (Fisher et al., 1995;Westhoff et al., 2021). This situation applied to thin debris cover or topographic depressions where snow and water accumulated. This resulted in subsidence of the terrain surface, deeper snow accumulation in winter and ponding of surface water during the warm season, which began to thaw the underlying permafrost. During the first phase of lake development, 15 relict glacier ice can serve as a focal point for the onset of accelerated thermokarst degradation. If exposed, the ice core then undergoes accelerated wastage through the effects of solar radiation or becomes buried again under the slumping material until a new thermal balance can be reached. We cannot deduce absolute timing for the inception of lake K since no reliable basal dates are available. However, we suggest that lake inception of these deeper glacial lakes occurred sometime during the mid-Holocene and precedes that of the shallowest lakes formed by the thawing of ice wedges, which were 14 C dated at around 2100 20 yr BP (Bouchard et al., 2020).
Stage 4: Thermokarst lake inception and lake expansion within permafrost. Once a lake gets deeper than the maximum thickness of the winter ice cover (~2 m ± 20 cm in the valley as measured by Prėskienis et al., 2021), it will continue to grow laterally (thermo-mechanical erosion) and vertically (subsidence) by thermokarst processes each year; the rate of expansion 25 being dependant of the local climatic conditions, ground-ice content and lake bed temperature (Fig.7d). In cases where buried glacier ice remains present beneath the lakebed, the ice will slowly continue to melt, causing lake bottom subsidence. Further ground ice melting and the resulting thaw slumps contribute to lake expansion, as shown by the headscarps located close to the shoreline of lake L (Fig. 3). Other studies have shown that thaw slumping is an important mechanism of lake expansion. (Hinkel et al., 2012;Plug and West, 2009;Kokelj et al., 2009) Our results indicate that these thermokarst glacial 30 lakes also evolved at a later stage as 'classic' thermokarst lakes that are now slowly expanding in area and volume, because of the melting of intrasedimental ground ice and ice wedges in the frozen silt-peat terrace and in the underlying glaciofluvial and till material. The shorelines of glacial lakes, such as kettle lakes, are expected to be very smooth and close to the shape of a perfect circle, but most lakes studied here display slightly irregular shorelines (Fig. S6). The shoreline morphology of the https://doi.org/10.5194/tc-2021-302 Preprint. Discussion started: 20 November 2021 c Author(s) 2021. CC BY 4.0 License. deeper glacial lakes is very similar to the other thermokarst lakes, indicating that all lakes are now laterally expanding in the polygon terrace by thermal and mechanical erosion. Thermokarst is an active landscape change mechanism currently operating in the valley and on the island in general (Bouchard et al., 2020;Fortier et al., 2007;Godin et al., 2014). Today, the lakes can expand by thermal subsidence and different shoreline erosional processes including: (1) the development of thermomechanical erosional niches; 2) the mechanical erosion caused by lake ice pushing against the shore, and 3) the incorporation of adjacent 5 polygonal ponds into the lake (Jones et al., 2011). Eventually, the lakes may cease expanding in the event of evaporation, partial or complete surface/subsurface drainage, which can sometimes be catastrophic (Bouchard et al., 2017 and references therein). This shows the interplay of climatic (external) and local landscape (internal) processes in the inception and evolution of thaw lakes in general, including the ones developed through melting of buried glacier ice.

Implications on Arctic lakes ecosystem dynamics 10
Lake morphometry, specifically depth, plays an important role in regulating lake water temperature and associated biogeochemistry. It will influence the mixing regime and the number of thermal overturn events per year during the openwater period (i.e., if the lake is monomictic, dimictic or polymictic, the latter being more common for Arctic lakes; Rautio et al., 2011). This is intrinsically linked to water column aeration and light regime, and thus exerting a strong control on respiration and primary production (Vincent, 2010). Results indicate that the three studied lakes can be considered as cold 15 polymictic, or dimictic depending on the year for lake K (although a mooring would be needed to validate this). Among the three studied lakes, bottom water of shallow lake G was the coldest by late winter but the warmest by late summer, a pattern directly linked to its mixing regime where meteorological conditions are more likely to influence bottom water temperature and talik formation. We also found that lake morphology influences dissolved oxygen. Lake G showed the lowest oxygen concentration at the end of the winter (< 2 mg L -1 or < 13% saturation in 2015), likely linked to its large sediment area to water 20 volume ratio and its higher organic content at the lake bottom (submerged peat polygons, as opposed to less organic sediments in the deeper lakes K and L) leading to a faster depletion of oxygen (Vincent, 2010;Ward et al., 2017). This can be a significant limiting factor for overwintering fish populations (Leppi et al., 2016). For the deeper lakes (K and L), the difference in water column stability controls bottom oxygen saturation levels during the open-water period, which decreased well below 60% in lake K (reaching anoxia just above sediment), while it always remained above this level in lake L . The stronger gradient in lake 25 K is likely related to its smaller size (smaller fetch) and greater depth. Climate change may therefore not only affect water temperature, mixing regime and oxygen availability through warming and summer lengthening, but also through effects on the evolution of lake morphology from the melting of buried glacier ice.
These differences in the mixing regime and oxygen availability, controlled by lake morphology (size and especially depth), 30 exert a strong control on the timing (seasonal differences) and magnitude of GHG emissions (CH4, CO2, and their relative proportion) from the water column to the atmosphere (Prėskienis et al., 2021;Bouchard et al., 2015a;Hughes-Allen et al., 2021;Matveev et al., 2016). Previous studies showed that lake G generally maintained high GHG fluxes during the open-water https://doi.org /10.5194/tc-2021-302 Preprint.  as compared to the deeper lakes K and L (named kettle lakes in Prėskienis et al., 2021). Once more, the combination of warmer temperatures and higher organic content of lake G likely explains its higher GHG emissions (Prėskienis et al. 2021, see Fig. 3 and 4, and Table 4), averaging 27.1 mmol CO2 equivalent m -2 d -1 , as compared to 10.8 mmol CO2 equivalent m -2 d -1 from lake K (mainly caused by differences in CH4 ebullition; not assessed in lake L). Considering that most GHG are emitted from lake sediment (Bastviken et al., 2004), it is important to underline that the largest sediment area of a lake is in contact with 5 epilimnetic (shallow) waters, and therefore not only bottom water temperature in the deepest pelagic section of a lake needs to be assessed. Due to the importance and diversity of lakes across the circumpolar Arctic, a better knowledge of their bathymetry and landscape variability is necessary to upscale local biogeochemical assessments to regional or continental scales.

Conclusion 10
In proglacial environments, bodies of buried glacier ice are slowly melting simultaneously with the deglaciation of the area but can remain in sedimentary deposits for much longer. As a result, buried glacier ice can play a role in the inception of thaw lakes and determine their long-term evolution in formerly glaciated permafrost landscapes. The origin and growth of numerous thermokarst lakes in the Qarlikturvik Valley, Bylot Island, has been examined using bathymetric and field surveys, highresolution remotely sensed imagery, and lake sediment analysis. Our results suggest that the melting of buried glacier ice has 15 triggered the inception of the deeper lakes (> 5 m, up to 12 m) in the study area, whereas the formation and evolution of the other shallower lakes (< 5 m) are related to the melting of ice wedges and intrasedimental ice. There are four key arguments in support of this: 1) Past and current presence of buried glacier ice in the study area was shown by the presence of late Pleistocene glacier ice preserved within permafrost, as well as numerous stable and active thaw slumps. 20 2) Co-occurrence of shallow and deeper lakes within the same depositional environment, which means that these lakes have been subjected to the same environmental and climatic conditions, and therefore notable depth difference must be related to different ground ice volume or time spanned since inception.
3) The morpho-stratigraphic analysis of two lake sediment cores revealed two distinct basin types in terms of sediment accumulation. This dissimilarity indicates that these basins have a different origin and evolution 25 throughout the Holocene. This suggests that glacial thermokarst lakes initiated by the melting of buried glacier ice have a distinct depositional history and sedimentological signature.
Moreover, analysis of lake morphometry and distribution revealed that lakes are more densely distributed near the most recent ice positions. It suggests a relationship between the formation of lakes and the deglaciation patterns in both Qarlikturvik Valley and the broader southern plain of Bylot Island. Given future climate projections, it is likely that Arctic lowlands with glacier 30 ice buried in permafrost will change dynamically because of surface permafrost degradation and melting of relict glacial ice.
It is expected that the deepening of the active layer and talik development, as well as the enlargement of Arctic lakes in response https://doi.org /10.5194/tc-2021-302 Preprint. Discussion started: 20 November 2021c Author(s) 2021 to global warming, will reach undisturbed buried glacier ice. This could create new aquatic ecosystems and strongly modify existing ones through the lateral expansion of lakes caused by wind-and circulation-driven erosion, thaw slumping and thaw subsidence along lake margins. In turn, this will likely have pervasive effects on geomorphological, hydrological, ecological processes of affected landscapes, including the high-latitude and global carbon budgets, oxythermal quality of fish habitats and human infrastructure. Transition" (ADAPT). The fieldwork benefited from the logistical support provided by Gilles Gauthier and his team (U. Laval).