Buried remnants of the Laurentide Ice Sheet and connections to its surface elevation

The Laurentide Ice Sheet (LIS) occupied a large part of North-America during the late Pleistocene. Determining the proper surface geometry and elevation of the LIS is of central importance to estimate global changes in sea-level and atmospheric circulation patterns during the late Pleistocene and Holocene. Despite largely disappearing from the landscape during the late Holocene, LIS remnants are found in the Penny and Barnes ice caps on Baffin Island (Canada) and ongoing permafrost degradation has been exposing relics of the LIS buried along its northern margin since the late Pleistocene. Here, we use the δ18O records of six LIS remnants and the late Pleistocene δ18O-elevation relation to establish ice elevation in their source area during the last glacial maximum (LGM). Contrary to some modeled reconstructions, our findings indicate an asymmetric LIS topography with higher ice on Keewatin Dome (~3200 m) and thinner ice in the prairies along the Plains divide (1700–2100 m) during LGM. The resiliency of icy permafrost to past warm intervals preserved relics of the LIS; these ice-marginal landscapes, now poised for thaw, should uncover more valuable clues about the conditions of the last major ice sheet on Earth.

spatial extent of grounded LIS ice, yield shear stresses and the ice frozen to the bed of Hudson Bay at LGM 2 ; no independent assessment of the reconstructed elevations using the δ 18 O-elevation relation was made.
In addition to the LIS remnants in Penny and Barnes ice caps, climate warming in the Canadian Arctic has been increasing thermokarst activity and buried relics of the LIS are being exposed along headwall of thaw slumps and other natural exposures 24 . Such sites have been discovered near the northern maximum limit of the LIS, including: (1) Peel Plateau, NWT, Canada; (2) Richards Island, NWT; (3) north-west Victoria Island, NU, Canada; (4) south-west Bylot Island, NU, Canada (Fig. 1). By analogy to the ice found along the margin of Greenland Ice Sheet 25, 26 , the buried LIS ice near its maximum northern extent must have originated from their local accumulation centers and their δ 18 O records likely preserve details about paleo-ice elevation that can be extracted using the δ 18   Sites of remnants of LIS ice and its reconstructed elevation during the last glacial maximum. The extent of the LIS (including ice-shelves) at last glacial maximum is derived from ref. 1 . The surface elevation of the LIS is derived from the steady-state model of ref. 2 which is based on the empirical margins of the ice sheet (minimum-concept of ice margins and excluding ice-shelves), a simple plastic ice rheology and assumes hardbed conditions in the Hudson Bay sector. Surface elevations are in 100′s of meters above present-day sea-level (errors are 5-7%). The thick dashed black line is the boundary between deformable beds in the Prairies and Great Lakes regions and hard beds for interior and eastern regions. The dashed red lines are inferred ice flow and source area for the four buried LIS sites. The underlying topography is from GTOPO30 digital elevation data (https://lta.cr.usgs.gov/). LGM] = 12-15‰) greatly exceeds that expected due to changing air temperatures only (a shift of 3-10‰ is expected from a zonally-averaged model over latitudes of 60-80°N that includes the effect of changes in sea-ice front 21 ; Fig. 2).
Buried LIS ice in ice-cored hummocky permafrost terrain. The four buried LIS sites consist of ice that was stranded or buried shortly after the LIS reached its maximum extent at ~21-17 ka BP 1,8 , with the LIS receding from these regions by 10-14 ka BP 1 . This timing of advance and retreat of the LIS at the sites offers a relative chronology for its burial and preservation in permafrost with the receding ages providing a minimum age for the ice because of the time needed to transport ice from the accumulation center to the margin (the ice is likely LGM-age, 29-19 ka BP) 27 .
Peel Plateau site, NWT, Canada (67°37.1′N; 135°33.4′W; 350 m a.s.l.). On the Peel Plateau ( Fig. 1), buried LIS ice was identified in ice-cored hummocky terrain near the maximum westward extent of the LIS (reached near 18 k cal yr BP) 28 and within the Tutsieta Lake Phase recessional limit dated to 15 k cal yr BP 29 . The ice sheet had receded from the region by 13.5 k cal yr BP based on 14 C ages obtained from remains of steppe bison near Tsigehtchic 30 . In this region, a major ice lobe flowed north-northwest along the Mackenzie Valley with the Mackenzie Through and Bear Lake ice streams likely bringing LIS ice from the more southerly Plains ice divide 2,14 .
The buried ice was exposed over a 70-m wide and ca. 5 to 6 m high headwall. Leaf fragments in the soils overlying the buried ice were dated to 9534 cal yr BP. Along the headwall, four types of ice were identified: (1) clear ice with mm-size spherical gas inclusions; (2) sub-vertically banded clear ice and fine sediments; (3) bubble-poor blue ice; (4) white ice rich in spherical gas inclusions, similar to the bubble-rich white ice band on Barnes Ice Cap. The four ice types were characterized by a Ca-Na-Cl geochemical facies (Fig. 3), with concentrations being lowest in the white ice (less than 5 mg L −1 ). The white ice also has the lowest δ 18 O values (−31.2 ± 0.6‰; ranged from -32.0 to −29.5‰), whereas the other three ice types have slightly higher δ 18 O values, ranging from −31.2 to −27.3‰ (Fig. 3). The D-excess values are in the 5.0 to 9.8‰ range for the white ice (within the range of that on Penny and Barnes ice caps; Fig. S1), and slightly lower in the other ice types (in the −0.5 to 7.2‰ range). The physical, geochemical and δ 18 O and D-excess characteristics of the ice and its stratigraphic location suggest that it corresponds to late Wisconsinan age LIS ice.  Along the shore of YaYa Lake, massive ice was identified by radar and borehole drilling beneath 7-m of glaciofluvial sand 33  Near Loch Point on Prince Albert Peninsula ( Fig. 1), ref. 34 identified buried LIS ice beneath 1-m thick till extending 60-450 m laterally. The englacial ice, which contained a small amount of dispersed sediments, had a sub-vertical banding with alternating bubble-rich and bubble-free ice and its δ 18 O averaged −30.1 ± 1.4‰ (ranged from −32 to −27.5‰; Fig. 3).

Bylot Island site. (73°09′N; 79°57′W; 25 m a.s.l.). Buried LIS ice was discovered in the Qarlikturvik
Valley on southwest Bylot Island 37 . The late Wisconsinan LIS limit was first positioned south of Bylot Island and it was suggested that alpine glaciers, similar in size to those present today, occupied Bylot Island during that time 38 . However, recent mapping by ref. 39 provided a re-interpretation of LIS extent and it is now suggested that late Wisconsinan LIS ice filled Navy Board Inlet (as evidenced from mega-scale glacial lineations on the sea-bed) with ice flowing northward from the Lancaster Sound ice stream 40 . Radiocarbon age of shells in nearby marine sediments indicate the region was deglaciated near 9860 ± 140 14 C BP 41 .
The exposed buried ice in the valley was discovered beneath an ice-contact stratified drift 37 . The ice had a clear to whitish appearance, similar to the late Pleistocene-age white ice band on Barnes Ice Cap. The buried ice contained many small spherical and coalescing gas inclusions. The ice was dominated by Ca 2+ , Na + , Mg 2+ cations and their concentration was low and within the range of ice sampled from the nearby C-93 glacier (Fig. 3). The δ 18 O of the buried ice ranged from −34.4 to −33.4‰ (average of −34.0 ± 0.4‰; Fig. 2) and D-excess averaged 6.6 ± 2.5‰, both similar to values obtained from Barnes Ice Cap (Fig. S1).  42 ). Further, the ~4‰ difference between the δ 18 O record of the LIS remnants cannot solely be explained by regional temperatures or major moisture source differences because the δ 18 O composition of modern precipitation is strongly zonal and varies by less than 1‰ along latitudes of 67-72°; late Pleistocene δ 18 O precipitation records from 154 empirical proxy data showed a similar zonal distribution 43 . However, like Penny and Barnes ice caps 20 , the difference in Δ 18 O [modern-late Wisconsinan] for the buried LIS sites is higher than expected for their latitudes (Fig. 2). Therefore, a more likely explanation of the differing δ 18 O records between the LIS remnants, including their higher than expected δ 18 O shift between late Wisconsinan and modern conditions, is differences in ice elevation in their source areas. On ice sheets and ice caps, precipitation accumulates in the catchment area above the equilibrium line and is transferred to lower elevation along flow lines; this occurs while maintaining the δ 18 O stratigraphy. This was demonstrated on the Greenland Ice Sheet where a horizontal δ 18 O transect at the edge of the ice sheet was nearly identical (with a ~1-2‰ variation) to a deep ice core at the head of its flow line near the center of the ice sheet (Pakitsoq site and GISP2 ice core) 25,26 ; such ice transport and preservation of δ 18 O record was also incorporated in 3D isotope stratigraphy models 44 (Fig. S1). Therefore, the δ 18 O records of the buried LIS ice near its maximum northern extent likely preserve details about ice surface elevation in their source areas and this information can be extracted from the δ 18 O-elevation relation.
The δ 18 O records from Penny and Barnes ice caps have a sufficient chronological control, such that we can infer ice surface elevations at these sites for the late Wisconsinan and LGM with good confidence. The buried LIS ice in the ice-cored permafrost terrains near the maximum extent of the LIS lack a robust chronological control due to challenges of obtaining reliable ages from ice itself; however, the timing of advance and retreat of the LIS at these localities indicates that the ice is of late Wisconsinan-age and most likely LGM-age. To test the effect of random sampling of the buried LIS ice and if it may represent LGM conditions (29-19 ka BP) 27 , we performed a frequency distribution analysis of the late Pleistocene age δ 18 O records of GISP2, Penny and Barnes ice caps (Fig. 5); the likelihood of randomly sampling ice with δ 18 O values in the range of LGM-age ice is high (28-36%). Additionally, the variance in δ 18 O for the late Wisconsinan ice was assessed for the GISP2 and Penny δ 18 O records and standard deviations of 1.1 to 1.4‰ were calculated. As such, random sampling of buried late Wisconsinan ice does not have δ 18 O values that would deviate greatly from LGM-age ice (these would translate into errors of <200 m in our paleo-elevations reconstruction). Increasing the uncertainties in δ 18 O to 3‰ would translate into an elevation error of 500 m.
Together, the paleo-elevation reconstruction derived from the δ 18 O records of six LIS remnants sites indicate an asymmetric LIS during the LGM, with ice surface elevations near 2200-2400 m on Foxe Dome, 3200 ± 100 m on the Keewatin Dome and 1700-2100 m along the Plains ice divide (Fig. 1). Overall, our reconstruction of ice surface elevation from northern sectors of the LIS near LGM is consistent with the Fisher-1985 model, using both minimum and maximum grounded ice margin extents (Figs 1 and S2), and the Tarasov-2012 model 5 54 . The ICE5G, ICE6G and NAICE models generally predicted too much ice in the Keewatin Dome and/or Plains ice divide. Given some uncertainties in our approach and in geophysical models, we highlight differences in the order of 1000 m between these various paleo-elevation reconstructions at LGM. For example, the ICE5G model 55 has similar ice elevations for the Foxe Dome and along the Mc'Clintock and Plains ice divides, but it has a much thicker Keewatin Dome (~4000 m). Conversely, the ICE6G model 6 has similar ice elevations for the Keewatin and Foxe domes, but increased ice elevation along the Plains ice divide (2800-3200 m). The NAICE model 4 also has much higher ice elevation along the Plains ice divide (~2500-3000 m).
The concept of using δ 18    volume at LGM of ~21.1 × 10 6 km 3 (ref. 2 ; increases to ~25.9 × 10 6 km 3 if the maximum-concept of ice margins at LGM is used). Some of the latest modeled glacio-geophysical LIS elevation reconstruction are generally found to be over-estimating its surface elevation near the LGM on the Keewatin Dome and along the Plains ice divide, although these may also relate to uncertainties in our ice-age estimates and the dynamic nature of the LIS. There are probably other remnants near former ice-marginal extents of past ice sheets in North-America and Siberia that contain information that would help reconstruct their interior geometry and global-scale impacts. Given that the ongoing climate warming and associated thermokarst is exposing new sites that potentially contain LIS relics in permafrost 24 , these should be found and sampled before they melt away as they contain key information about the conditions prevailing for the last major ice sheet on Earth (i.e., trapped gases, impurities).

Methods
Identifying buried LIS ice. Distinguishing between buried glacial ice and other massive ice types in permafrost is usually based on cryostratigraphy and combines detailed physical, geochemical and isotopic measurements of the ice and may also include analyses of occluded gases, while its chronology usually relies on the age of the surrounding sediments 28 .
Sampling. The buried ice on the Peel Plateau and Bylot Island were sampled in summer 2016 and 2012, respectively. At both sites, a cryostratigraphic approach was used to describe the exposed ice and the overlying sediments. Samples of ice were collected using either an ice pick or a portable core-drill. All samples were shipped and stored in our laboratories until analyses.
Water isotope analysis. The stable isotope ratios of oxygen ( 18 O/ 16 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 (δ 18  Ice flow in ice sheets. Ice sheet flow models that predicts that the margins of glaciers contain ice that originate from their local accumulation area at much higher elevation 25, 44,56 . The latter was observed on Greenland Ice Sheet where the δ 18 O record of ice sampled along a horizontal transect from the western edge (Pakitsoq site) was remarkably similar to that of the GISP2 δ 18 O ice core record situated near the center of its accumulation area due to ice resurfacing in the ablation zone 25, 26 . The ice-flow on Greenland Ice Sheet was successfully modeled using a 3D transport model fitted with δ 18 O data from ice cores 44 . In fact, both Pakitsoq and Barnes Ice Cap exposed late Pleistocene-age ice near their margins over a horizontal distances of ~200-500 m; with Holocene-age ice above it (Fig. S1). . The steeper slope in Arctic region is attributed to colder-drier conditions and near adiabatic lapse rate with little mixing of air masses as they gain elevation; this is particularly the case on large ice sheets. In the lower latitudes, vertical advective mixing is much stronger and the mixing of air masses reduces the δ 18 O-elevation slope.

Holocene to late
The δ 18 O of precipitation during the late Wisconsinan is largely affected by changes in air temperature, site elevation, seasonal distribution in amount of precipitation and atmospheric circulation, the latter two relate to changes in moisture source. The effect of a change in zonally-averaged air temperature during the late Wisconsinan, including a shift in the position of sea-ice front, was associated with a latitudinal-dependent shift of 3-10‰ in δ 18 O of precipitation, with the higher latitudes experiencing a greater shift 21 . This trend was also observed in 154 late Pleistocene δ 18 O proxy records 43 . The modern δ 18 O-elevation slope cannot be directly transferred to other and different climate periods, such as the cooler late Pleistocene. However, based on empirical data, ref. 43 suggested that the global δ 18 O-elevation relation during the late Pleistocene changed by −0.05‰ per 100 m. We therefore use a late Pleistocene δ 18 O-elevation relation of −0.67 ± 0.03‰ per 100 m for Arctic regions. This late Pleistocene δ 18 O-elevation slope is nearly identical to that observed in the colder Antarctica (−0.68‰ per 100 m; ref. 50 ). Changes in seasonal distribution in amount precipitation and atmospheric circulation can be assessed from the D-excess parameter, which can provide clues into changes in distance to moisture source (increasing values with increasing distance from moisture) 57 . Examining the D-excess record of Penny Ice Cap (Fig. S1)  (and the east coast of the Queen Elizabeth Islands) largely originates from nearby Baffin Bay with the air masses experiencing adiabatic cooling during uplift over the mountain ranges followed by isobaric cooling in the interior (i.e., rainout at constant condensation altitude with Rayleigh-type distillation) 17,49 . However, potential changes in moisture source during the late Wisconsinan in the north-central (i.e., Victoria Island site) and north-western (i.e., Peel Plateau and Richards Island sites) sectors of the LIS have more uncertainty. The D-excess values of the buried white ice on the Peel Plateau are within the range of that on Penny and Barnes ice caps and also of that calculated in modern-day precipitation at nearby sites (Inuvik: 14.9 ± 7.6‰; Cambridge Bay: 9.6 ± 4.8‰) 58 . Therefore, if there was a change in distance to moisture source, it is not reflected in the D-excess values. However, to truly assess source distribution changes, one would have to do a trajectory tracing of ancient climates over the proposed ice sheet geometries. The biggest unknown relates to the age of the sampled buried ices. We ascribed to the four buried LIS ice a near LGM-age (29-19 ka BP) 27 based on the existing chronologies of the timing of advance and retreat of the LIS in these areas. The LIS reached its maximum extent at ~21k-17 ka BP 1,8 and had receded from these regions by 10-14 ka BP 1 . This timing of retreat of the LIS at the sites provides a minimum age constraints for the sampled ice with the actual age at each site being dependent on the ice flow rate in the ice sheet. The uncertainty in the age of the sampled buried LIS sites is the biggest uncertainty when comparing reconstructed paleo-elevations derived from the δ 18 O-elevation relation to those from various glaciological and geophysical models that have reconstruction for specific times because the LIS topography is evolving over time.
Modeling LIS surface elevation with minimum-and maximum-concept ice margins. The reconstruction of late Wisconsinan LIS surface elevation at 18 ka is from ref. 2 . It uses a simple plastic ice model that is insensitive to unknown parameters and uses as inputs the margins of the ice sheet, present-day topography and an assumed yield shear stress. References 59,60 showed that elevation estimates using this model are good to 6% on present-day Greenland, with ridges and domes having smaller errors. No assumptions are made in advance about ice divides and accumulation centers. At the time of its publication in 1985, some uncertainty still existed about northern and northeastern LIS margins, thus ref. 2 presented the LIS surface elevation for the minimum-concept of grounded ice margins. The maximum-concept of ice margins was produced at that time also and is shown in Fig. S2. It uses all the same inputs as the minimum-concept model. Increasing the ice margins allow for slightly higher LIS elevations and for a higher ice volume (25.9 × 10 6 km 3 ). Note that the Cordilleran and Innuitian ice sheets have not been included in these reconstructions so their volumes are not included.