INTRODUCTION
Large tabular bodies of massive ground ice are often located in formerly glaciated terrain near the maximum extent of ice sheets or within their recessional limits (French and Harry, Reference French and Harry1990; Dyke and Savelle, Reference Dyke and Savelle2000; Ingólfsson and Lokrantz, Reference Ingólfsson and Lokrantz2003; Lacelle et al., Reference Lacelle, Fisher, Coulombe, Fortier and Frappier2018; Coulombe et al., Reference Coulombe, Fortier, Lacelle, Kanevskiy and Shur2019). These permafrost terrains are, as a result, vulnerable to thaw due to their high ground ice content (Rudy et al., Reference Rudy, Lamoureux, Kokelj, Smith and England2017). For example, regional mapping of thermokarst terrain in western Arctic Canada revealed that slump-affected locations delineate the maximum and recessional positions of the Laurentide Ice Sheet (Kokelj et al., Reference Kokelj, Lantz, Tunnicliffe, Segal and Lacelle2017). As such, the development and preservation of tabular massive ice bodies reflect a heritage of Late Quaternary glacial history and ground thermal conditions. In western Arctic Canada, the origin of tabular massive ice bodies preserved in permafrost has been associated with either (1) the burial of glacial ice by an insulating cover of ablation till with thickness greater than the active layer (e.g., Dyke and Savelle, Reference Dyke and Savelle2000; Lakeman and England, Reference Lakeman and England2012) or (2) the in situ freezing of glacial meltwater migrating in a talik toward a freezing front during permafrost aggradation (i.e., segregated-intrusive origin) (Dallimore and Wolfe, Reference Dallimore and Wolfe1988; Harry et al., Reference Harry, French and Pollard1988; Mackay and Dallimore, Reference Mackay and Dallimore1992; Lacelle et al., Reference Lacelle, Bjornson, Lauriol, Clark and Troutet2004). Based on evidence from sites in the Mackenzie delta and the Yukon Coastal Plain, Rampton (Reference Rampton1988) hypothesized that the formation of massive ice of segregated-intrusive origin occurred following flow of glacial meltwater under hydrostatic pressure in a supra-permafrost talik that would then freeze during downward permafrost aggradation following thinning and retreat of the ice sheet. Accordingly, this mechanism that requires a warm-based ice sheet would not apply to regions with cold-based ice sheets, like some places in the high Arctic Canada.
Similar to the Mackenzie delta and Yukon Coastal Plain, the Fosheim Peninsula in the Queen Elizabeth Islands (high Arctic Canada) is a region with abundant exposures of thick tabular massive ground ice bodies (Hodgson and Nixon, Reference Hodgson and Nixon1998). However, regional mapping of thermokarst terrain showed that slump-affected terrain does not occur in glaciogenic sediments along the maximum or recessional positions of the Innuitian Ice Sheet (IIS), but instead, is found uniquely in raised Holocene marine-deltaic sediments (Pollard, Reference Pollard, Garneau and Alt2000a; Ward Jones et al., Reference Ward Jones, Pollard and Jones2019). Based on the spatial distribution and elevation of the tabular massive ice exposures and cryostratigraphic investigations, it was suggested that the ice has a segregated origin and formed synchronously with the Holocene marine regression and subsequent permafrost aggradation in the marine sediments (Pollard, Reference Pollard1991, Reference Pollard, Garneau and Alt2000b; Pollard et al., Reference Pollard, Ward and Becker2015). Pollard (Reference Pollard, Garneau and Alt2000b) later hypothesized that groundwater movement in an unconsolidated and weathered sandstone bedrock aquifer would be driven by a hydraulic and thermal gradient to a freezing front during permafrost aggradation following marine regression and land emergence. Overall, this conceptual model shares some similarities with the one proposed by Rampton (Reference Rampton1988) for the western Canadian Arctic (i.e., permafrost aggrading following a receding ice sheet); however, it has major differences. On the Fosheim Peninsula, permafrost aggraded in marine sediments, and time-transgressive permafrost aggradation is constrained by land emergence during isostatic rebound and not the receding ice sheet. Additionally, the source water recharging the aquifer, although mostly of glacial meltwater origin, may receive contribution of brackish seawater. However, the conceptual model remains a source of controversy, as it does not explain the δD-δ18O regression slope values that were near the meteoric water line, thus leaving a buried glacial ice origin a possibility (Pollard, Reference Pollard1991).
This study tests the conceptual model for the development of tabular massive ice in the raised marine-deltaic sediments of the Fosheim Peninsula based on new regional cryostratigraphic data from three sites found between the coastline and the Holocene marine limit (Fig. 1). This objective is reached by describing the cryostructures and measuring the major ions, δD-δ18O, and 14CDOC of the ice. The results are compared with the output of a numerical model (FREEZCH9) of the progressive freezing of an aquifer that is being recharged with two sources of water having different isotopic composition. Modeling allows for an estimate of the source(s) of water recharging the aquifer and an assessment of the value of the δD-δ18O regression slope during open-system freezing. Overall, the study provides insights into the process of massive ice formation in raised marine-deltaic sediments and is relevant to other localities where ice is forming from a reservoir receiving inputs of water of different δD-δ18O composition.
Location of sites sampled in the Eureka Sound Lowlands, high Arctic Canada. (A–C) Maps showing sampling sites around Eureka (Gemini, Station Creek, and Dump sites). (D) Comparison of elevation of field sites with the regional isostatic curves. Site elevation range from near the marine limit (125 m asl) to near the coast (40 m asl), which reflects a subaerial exposure of ca. 8000 BP to ca. 5000 14C yr BP, respectively. Marine regression curves are from Hodgson and Nixon (Reference Hodgson and Nixon1998) and Bell and Hodgson (Reference Bell, Hodgson, Garneau and Alt2000). Recessional glacial limits are from England et al. (Reference England, Atkinson, Bednarski, Dyke, Hodgson and Ó Cofaigh2006). AIC, Agassiz Ice Cap; GE, Gemini site; IIS, Innuitian Ice Sheet; Nu, Nunavut; SC, Station Creek site; DS, Dump site.
STUDY AREA
The Fosheim Peninsula is an intermontane region located in the west-central section of Ellesmere Island and the southeastern section of Axel Heiberg Island (Fig. 1). The region is surrounded by massifs exceeding 1000 m above sea level (m asl) and hosts two main physiographic regions: the uplands that exceed 500 m asl and the Eureka Sound Lowland (ESL) that consists of flat to gently rolling terrain situated at elevations <200 m asl. The region is underlain by the Sverdrup Basin, a folded and faulted Carboniferous to Paleogene sedimentary bedrock composed of sandstone, siltstone, and shale strata (Bell, Reference Bell1996; Hodgson and Nixon, Reference Hodgson and Nixon1998).
The Late Quaternary glacial history of the northernmost cluster of islands in Canada's Arctic archipelago has long been debated with contrasting views (Hodgson, Reference Hodgson1985; Bell, Reference Bell1996; England et al., Reference England, Atkinson, Bednarski, Dyke, Hodgson and Ó Cofaigh2006). The most recent glacial history reconstruction suggests that the IIS covered most of this region and reached its maximum extent at about 15,800 14C yr BP (~19,100 cal yr BP), with the ice being ca. 1.6 km thick in the ESL (Dyke et al., Reference Dyke, Andrews, Clark, England, Miller, Shaw and Veillette2002; England et al., Reference England, Atkinson, Bednarski, Dyke, Hodgson and Ó Cofaigh2006). During the last glacial maximum, the Fosheim Peninsula was an ice saddle, where ice flowing west of the Ellesmere Ice Divide coalesced with ice flowing east from the Axel Heiberg alpine divide. The coalescing ice flow in fjords and interisland channels likely resulted in local warm-based conditions (England et al., Reference England, Atkinson, Bednarski, Dyke, Hodgson and Ó Cofaigh2006). Deglaciation started between 10,300–8800 14C yr BP (~12,100–9800 cal yr BP) (England et al., Reference England, Atkinson, Bednarski, Dyke, Hodgson and Ó Cofaigh2006), and remnants of the ice sheet are still found at the nearby Agassiz Ice Cap. Marine transgression followed until 8800–7800 14C yr BP (~9800–8600 cal yr BP) with the Holocene marine limit found between 139 and 150 m asl (Bell, Reference Bell1996). The surficial sediments on the Fosheim Peninsula reflect the regional glacial history, with weathered bedrock, till, and glaciofluvial sediments found above the marine limit and Early Holocene–age marine and deltaic deposits found below the marine limit interspersed with patches of fluvial and alluvial deposits and small areas of thin peat (Hodgson and Nixon, Reference Hodgson and Nixon1998; Bell and Hodgson, Reference Bell, Hodgson, Garneau and Alt2000. The thickness of marine sediments varies from a thin veneer to thick horizontally bedded sandy clay (rhythmites). Solifluction, rill washing, and slope failures reworked these deposits into colluvial deposits along hillslopes (Bell and Hodgson, Reference Bell, Hodgson, Garneau and Alt2000).
The Fosheim Peninsula is characterized by a polar desert climate, with cold winter and cool summer air temperatures. The 1980–2016 mean annual air temperature recorded at Eureka weather station is −18.5 ± 1.4°C, and the total annual precipitation reaches 77.6 ± 25.8 mm yr−1, making it one of the coldest-driest places in Canada (Environment Canada, 2010). The typical vegetation consists of a Salix–Dryas community occupying up to 50% of hummocky tundra and less than 20% on drier and more exposed slopes (Becker et al., Reference Becker, Davies and Pollard2016). Recently disturbed areas are largely vegetation free and covered by efflorescence salts (Kokelj and Lewkowicz, Reference Kokelj and Lewkowicz1999). Holocene air temperatures were reconstructed from δ18O measurements from the nearby Agassiz Ice Cap. The 25-yr mean annual air temperature record shows a rapid Early Holocene warming, with temperatures being 6–8°C warmer than today at 10,000 yr b2k followed by a gradual cooling to AD 1700 (Lecavalier et al., Reference Lecavalier, Fisher, Milne, Vinther, Tarasov, Huybrechts and Lacelle2017). Presently, air temperatures are rapidly warming and are now at their warmest in the past 6800–7800 yr (Lecavalier et al., Reference Lecavalier, Fisher, Milne, Vinther, Tarasov, Huybrechts and Lacelle2017; Copland et al., Reference Copland, Lacelle, Fisher, Delaney, Thomson, Main and Burgess2020).
The climate and vegetation conditions ensure that permafrost is cold, continuous, and ~500 m thick (Pollard, Reference Pollard, Garneau and Alt2000a). The thickness of the active layer is typically 40–80 cm but varies due to microclimatic conditions. The mean annual ground temperature measured at 15.4 m, the depth of zero annual amplitude, is −16.5°C (Pollard et al., Reference Pollard, Ward and Becker2015). Permafrost is generally ice-rich, as revealed by 152 boreholes drilled to depths of 8 m with average ice content of 53%vol; the marine-deltaic deposits contained the highest amount of ice (average of 64%vol) (Hodgson and Nixon, Reference Hodgson and Nixon1998). The landscape is dominated by ice wedge polygons, and they occupy ~50% of the ESL (Couture and Pollard, Reference Couture and Pollard1998; Bernard-Grand'Maison and Pollard, Reference Bernard-Grand'Maison and Pollard2018). Based on 14CDOC dating, the ice wedges started to develop immediately after marine regression (Campbell-Heaton et al., Reference Campbell-Heaton, Lacelle, Fisher and Pollard2021). The ice wedges, although mostly found within the Holocene marine limit, are also found in glacial sediments; however, the large tabular bodies of ice are found exclusively within the fine-grained Holocene marine-deltaic sediments (Hodgson and Nixon, Reference Hodgson and Nixon1998; Pollard, Reference Pollard, Garneau and Alt2000a, Reference Pollard, Garneau and Alt2000b). Thaw slumps often develop along the ice wedge network, but they also develop in areas underlain by tabular massive ice bodies.
METHODS
Field sampling
On the Fosheim Peninsula, the tabular massive ice exposures are found in two settings: (1) massive ice beneath a thin marine-deltaic sediment cover and (2) massive ice beneath >10-m-thick marine-deltaic deposits (Pollard, Reference Pollard, Garneau and Alt2000a, Reference Pollard, Garneau and Alt2000b). The tabular massive ice bodies were sampled along the headwall of three thaw slumps located between the Holocene marine limit and the coastline: Gemini slump (120–125 m asl), Station Creek slump (43–45 m asl), and Dump slump (40–43 m asl) (Fig. 2). Based on the marine emergence curve, the Gemini site was subaerially exposed at ca. 8850 14C yr BP (~10,000 cal yr BP), while the two lower-elevation sites became exposed at ca. 5000 14C yr BP (~5700 cal yr BP) (Fig. 1).
Field photographs of tabular massive ice exposed in headwalls of thaw slumps in the Eureka Sound Lowland, NU. (A) Dump site (DS). The tabular massive ice is exposed beneath 0.5 to 1 m of laminated fine-grained marine sediments. (B and C) Station Creek (SC) site. The tabular massive ice is found beneath ~10 m of marine sediments containing reticulate ice. (D and E) Gemini site, where the tabular massive ice is found beneath ~1 to 1.5 m of laminated fine-grained marine sediments.
At each site, the stratigraphy of the material exposed in the headwall was described. At Dump and Gemini sites, the massive ice was <2 m below the ground surface, and a CRREL permafrost corer was used to collect 8.9-cm-diameter (3.5-inch-diameter) vertical cores from the permafrost table to a maximum depth of 1.5 m. At Station Creek, the ice was found beneath >2 m of sediments, and vertical cores up to 1.5 m depth were collected directly into the massive ice by constructing a drilling platform on the slump floor. Individual grab samples from the cleaned ablating ice face and overlying icy permafrost were also collected where possible. All frozen samples were wrapped in aluminum foil and stored in a freezer at −15°C in Eureka. In addition, two late-lying snowbanks and water from Blacktop Creek and Slidre Fjord were sampled and stored in 20 ml HDPE vials. All samples were shipped to the CryoLab at the University of Ottawa, where the active layer soil and the water samples were stored at +4°C and the permafrost samples were stored at −20°C.
Laboratory analyses
Analyses of the massive ice samples included a description of cryostructures, volumetric ice contents, δD-δ18O, major ions, dissolved organic carbon (DOC), and 14CDOC on the melted and filtered samples (prerinsed 0.45 μm cellulose nitrate filters). The cores of massive ice were first sliced along their length using a precleaned circular saw with a 0.8-mm-thick cutting blade and then subsampled at 2 cm intervals; samples were allowed to thaw in sealed and sterile 50 ml falcon tubes.
The half-cores were photographed using a digital camera under a light table to investigate cryostructures. The gravimetric water content was determined from the mass of water and the dry mass of the sample (dried at 105°C for 24 hours) (e.g., Van Everdingen, Reference Van Everdingen2005). The volumetric ice content and excess ice contents were calculated following Pollard and French (Reference Pollard and French1980). The concentrations of major cations (Ca2+, K+, Mg2+, Na+) and anions (SO42−, Cl−, NO3−) of the filtered meltwater were measured using an Agilent 4200 inductively coupled plasma atomic emission spectrometer and ion chromatography (Dionex ICS-2100), respectively. Before the analysis of cations, water samples were acidified to pH of 2 with ultrapure HNO3 acid. Analytical reproducibility for solute analysis is ±5%.
The 18O/16O and D/H ratios of the meltwater of the massive ice were determined using a Los Gatos Research liquid water analyzer coupled to a CTC LC-PAL autosampler for simultaneous 18O/16O and D/H ratios measurements of H2O and verified for spectral interference contamination. The results are presented using the δ-notation (δ18O and δD), where δ represents the parts per thousand differences for 18O/16O or D/H in a sample with respect to Vienna Standard Mean Ocean Water (VSMOW). Analytical reproducibility for δ18O and δD is ±0.3‰ and ±1‰, respectively. Based on the δD and δ18O values, the D-excess was calculated according to the equation provided in Dansgaard (Reference Dansgaard1964).
For DOC and 14CDOC analysis, the massive ice samples were first melted in prebaked 500 mL glass beakers covered with plastic wrap and sealed. Samples were then filtered in prerinsed 0.45 μm cellulose nitrate filters, with subsamples transferred in prebaked 1 L glass amber bottles. The DOC concentration of the meltwater was first measured using a wet TOC analyzer interfaced with a Thermo DeltaPlus XP isotope-ratio mass spectrometer using methods described by St-Jean (Reference St-Jean2003) at the Ján Veizer Stable Isotope Laboratory, University Ottawa. The 2σ analytical precision is ±0.5 ppm. Radiocarbon analysis of the DOC was subsequently performed at the A.E. Lalonde Accelerator Mass Spectrometry Laboratory, University of Ottawa. Sample preparation, including the extraction of DOC from waters and graphitization for 14C analysis, is described in Grinter et al. (Reference Grinter, Lacelle, Baranova, Murseli and Clark2019) and Murseli et al. (Reference Murseli, Middlestead, St-Jean, Zhao, Jean, Crann, Kieser and Clark2019). Graphitized samples were analyzed on a 3 MV tandem mass spectrometer (Kieser et al., Reference Kieser, Zhao, Clark, Cornett, Litherland, Klein, Mous and Alary2015). The 14C/12C ratios are expressed as fraction of modern carbon (F14C) and corrected for spectrometer and preparation fractionation using the AMS measured 13C/12C ratio (Crann et al., Reference Crann, Murseli, St-Jean, Zhao, Clark and Kieser2017). Radiocarbon ages are calculated as −8033ln(F14C) and reported in 14C yr BP (BP = AD 1950) (Stuiver and Polach, Reference Stuiver and Polach1977). In the text, radiocarbon ages with a standard error between 50 and 1000 yr are rounded to the nearest 10 yr. Our 14CDOC ages were calibrated using OxCal and the IntCal20 calibration curve (Bronk Ramsey, Reference Bronk Ramsey2009; Reimer et al., Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey and Butzin2020). All other calibrated ages in parentheses in the text are approximate (median age rounded to the nearest century) and were done by the editors of the journal.
RESULTS
Stratigraphy and chemistry of the tabular massive ice exposures
At the three sites, the headwall of slumps exposes 2 to 20 m of tabular massive ice. At the Dump and Gemini sites, the massive ice is found beneath ~1 to 1.5 m of laminated fine-grained marine sediments (silty clay), whereas at the Station Creek site, the ice is found beneath ~10 m of marine sediments containing reticulate ice. The upper contact between the ice and the overlying marine sediments is conformable and gradational (Fig. 2). Marine sediments also occur as thin layers within the upper section of the massive ice. The massive ice shows a horizontal to subhorizontal bedding, concordant with the structure of the overlying sediments (Fig. 2). In all cores, the volumetric ice content reached nearly 95%. When examined under the light table, the massive ice contained small amounts of suspended sediments, giving it a translucid color, with several small pockets of sediments and some ~0.5- to 2-cm-thick horizontal fine-grained sediment bands (Fig. 3).
Photographs of the ice cores under a light table: (A) Dump site (DS); (B and C) Gemini (GE) site. Suspended sediments are found in the massive ice, giving it a translucid color, with several small pockets of sediments and some ~0.5- to 2-cm-thick horizontal fine-grained sediment bands.
The geochemical facies of the massive ice vary according to site elevation (Figs. 4 and 5). At Dump and Station Creek (40–45 m asl), the ice is characterized by an NaCl geochemical facies, with Na/Cl molar ratios similar to Slidre Fjord and seawater, but with total dissolved solutes (TDS) concentrations of one to two orders of magnitude lower (average TDS of 868 ± 1990 and 310 ± 353 mg L−1, respectively). The TDS variation in the 1.3 m core from Dump remains rather constant with depth (60–200 mg L−1), except at 0.4–0.5 m depth, where TDS increase to 8830 mg L−1 due to a higher sediment content. At Station Creek, the TDS remain near 200 mg L−1 in the uppermost 0.75 m, then progressively increase to 1200 mg L−1 at 1.5 m depth. The ice at the higher-elevation site (Gemini, 115 to 120 m asl) is characterized a Ca-SO4 geochemical facies with lower TDS (average of 28 ± 15 and 37 ± 10 mg L−1, respectively). At both coring sites (G1 and G3), the TDS show little variation with depth. The Na/Cl molar ratios of the massive ice at Gemini is higher than at Dump and Station Creek, but more similar to that of snow.
Volumetric ice content (fVi), δ18O, D-excess, and concentrations of major ions in the tabular massive ice, Eureka Sound Lowlands, high Arctic Canada. (A and B) Gemini site, G3 and G1, respectively; (C) Station Creek (SC) site; (D) Dump site (DS).
Geochemical composition of the tabular massive ice compared with Slidre Fjord and snow in Eureka Sound Lowlands, high Arctic Canada. (A) Piper diagram of Cl−, SO42−, and Na/Cl; (B) Piper diagram of Mg2+, Ca2+, and (Na+K); (C) scatter plot of Na/Cl and (Na+K)/Cl molar ratios. G3 and G1 from Gemini site.
Similar to the chemistry, the δD-δ18O composition of the ice also varies by site elevation (Fig. 6). The higher-elevation site (Gemini) has the highest average δ18O composition (−30.9 ± 0.4‰), whereas the sites near the coast (Dump and Station Creek) have lower average δ18O composition (−32.3 ± 0.7‰ and −31.1 ± 1.6‰, respectively). At Dump, the δ18O profiles of the ice in the DS-1 and DS-2 adjacent cores were similar, and the samples were distributed along linear regression slope values of 7.8 and 7.3, respectively. At Station Creek, the δ18O composition of the massive ice ranged from −34.9‰ to −25.5‰ and a linear regression slope value of 7.5 was calculated between δD and δ18O (δD = 7.5 δ18O − 5.4; r2 = 0.98). At Gemini, two cores were collected at ~15 m vertical distance. The δ18O composition of the G3 core (near the surface) is ~3‰ higher (−30.8 ± 0.4‰) than that of the G1 core (−33.9‰ to −32.1‰), and both cores showed a decreasing trend with depth. The δD-δ18O slopes of G3 and G1 are 8.4 (δD = 8.4 δ18O + 22.7; r2 = 0.96) and 7.5 (δD = 7.5 δ18O − 7.8; r2 = 0.85), respectively.
(A) δD-δ18O composition of the tabular massive ice compared with Slidre Fjord and snow in Eureka Sound Lowlands, high Arctic Canada. Also shown is the Eureka local meteoric water line (LMWL: δD = 7.4, δ18O − 9.1; IAEA/WMO, 2015). (B) Na/Cl and δ18O composition of the tabular massive ice compared with Slidre Fjord and snow in Eureka Sound Lowlands. G3 and G1 from Gemini site.
The DOC concentration in the ice at all sites ranges from 0.9 to 2.9 mg L−1, with δ13CDOC averaging −26.4 ± 2.1‰. The DOC/Cl molar ratios average 1.66 at Gemini, 0.31 at Station, and 0.14 at Dump. The 14CDOC ages at Gemini and Dump sites, ranging from 26,150 ± 330 (30,460 cal yr BP) to 8850 ± 45 14C yr BP (9955 cal yr BP) (Table 1), are older than the timing of deglaciation and marine regression (Fig. 7) (Note that the 14CDOC ages reflect that of the DOC in the source(s) of water and not the timing of freezing.) The 14CDOC ages show a general trend of increasing ages toward lower elevation.
Comparison of 14CDOC of the tabular massive ice and ice wedges sampled at different elevation with the regional isostatic curves. Marine regression curves are from Hodgson and Nixon (Reference Hodgson and Nixon1998) and Bell and Hodgson (Reference Bell, Hodgson, Garneau and Alt2000). 14CDOC ages of ice wedges are from Campbell-Heaton et al. (Reference Campbell-Heaton, Lacelle, Fisher and Pollard2021).
Dissolved organic carbon (DOC), δ13CDOC, and 14CDOC results for samples of tabular massive ice collected at three sites in Eureka Sound Lowlands, high Arctic Canada.a
a The 14C ages were calibrated using OxCal and the IntCal20 calibration curve (Bronk Ramsey, Reference Bronk Ramsey2009; Reimer et al., Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey and Butzin2020).
DISCUSSION
Origin of tabular massive ground ice in the raised marine-deltaic sediments
The origin of the bodies of massive ice can be inferred by combining various approaches (i.e., cryostratigraphic observations, chemistry, and δD-δ18O of the ice). At the three sites, and others in the ESL, the massive ice is found stratigraphically below laminated fine-grained marine sediments, and the upper contact between the ice and the overlying marine sediments is conformable and gradational. Marine sediments also occur as thin layers within the upper part of the ice. The internal structure of the massive ice is horizontal to subhorizontal and concordant with the structure of the overlying sediments. These characteristics are more consistent with a segregation origin instead of burial of glacial ice (French, Reference French1998; French and Shur, Reference French and Shur2010; Lacelle and Vasil'chuk, Reference Lacelle and Vasil'chuk2013). The chemical composition of the ice also supports a segregation origin. The concentration of major ions in the tabular massive ice (tens to thousands parts per million range) is orders of magnitude higher than expected for firnified glacial ice (i.e., nearby Agassiz Ice Cap has ionic concentrations in parts per billion range; Fisher et al., Reference Fisher, Koerner and Reeh1995; Koerner et al., Reference Koerner, Fisher and Goto-Azuma1999). The geochemical facies (NaCl) and Na/Cl and DOC/Cl molar ratios in the massive ice suggest a contribution of seawater to the groundwater reservoir, especially for the lower-elevation sites (Dump and Station Creek). However, the low δ18O values (−35‰ to −25‰) of the ice indicate that, despite the addition of seawater, the main source of water for the ice is glacial meltwater. Ice formed by equilibrium freezing of water will typically have δD-δ18O regression slope values between 4 and 6, with the value largely depending on isotopic composition of initial water and freezing rate (Souchez and Jouzel, Reference Souchez and Jouzel1984; Lacelle, Reference Lacelle2011). However, the massive ice at the three sites has co-isotope slope values between 7.3 and 8.4, which at first sight would suggest a firnified ice origin (Fig. 6).
Freezing of water with δD-δ18O slope value near the global meteoric water line
In an open system where there is equilibrium freezing and recharge of the reservoir, Souchez and Jouzel (Reference Souchez and Jouzel1984) suggested that the value of the freezing slope (S) from an initial reservoir (r) that receives input of water (i) could be estimated from:
where β = fractionation factor of 18O ice water, α = fractionation factor of D ice water, δDi = δD composition of input, δDr = δD composition of reservoir, δ18Oi = δ18O composition of input, δ18Or = δ18O composition of reservoir, A = coefficient for input, and S = coefficient for freezing.
This equation shows that the freezing slope depends on the initial δD and δ18O of the reservoir (which influences the value of the fractionation factor) and the difference in δD and δ18O values between the input and the reservoir. In most natural systems, the difference in δD and δ18O between the input and the reservoir is small, and the forming ice would develop a freezing slope (values between ca. 4 and 6, and lower than the meteoric water line). However, Souchez and de Groote (Reference Souchez and de Groote1985) found basal ice at Gruben-gletscher in Valais, western Switzerland, with a δD-δ18O slope of 8.3, much higher than predicted for open-system freezing but similar to the meteoric water line. Souchez and de Groote (Reference Souchez and de Groote1985) simulated the evolution of δD and δ18O in ice if the reservoir was being recharged with water having an isotopically lighter composition than the initial reservoir using Eq. 1. It was found that the δD-δ18O of the ice would plot along a meteoric slope, albeit with slightly lower D-excess values. It was therefore suggested that the glacial meltwater received additional input from an ice-dammed lake near the glacier margin that had an isotopically lighter composition. Another example of meteoric slope produced under an open freezing system that receives an isotopically light input is at the Mammuthöhle cave in the Alps of central Austria (Kern et al., Reference Kern, Fórizs, Pavuza, Molnár and Nagy2011).
The landscape model for the development of the tabular bodies of massive ice in the ESL would fit this condition of open-system freezing and where the groundwater reservoir is being recharged by water with a substantially different δD-δ18O composition. However, there is one difference with Souchez and de Groote (Reference Souchez and de Groote1985). In the ESL, the reservoir initially recharged by glacial meltwater is likely receiving some input of brackish seawater with a higher isotopic composition. In any case, the reservoir would be recharged by water whose δD-δ18O composition is either (1) being progressively depleted as the melting of the IIS progresses to the deeper and older ice with lower isotopic composition (i.e., δ18O transition from the Early Holocene to late Pleistocene age ice at Agassiz Ice Cap is in the order of 7–10‰; Fisher et al., Reference Fisher, Koerner and Reeh1995; Zdanowicz et al., Reference Zdanowicz, Fisher, Clark and Lacelle2002) or (2) enriched over the initial reservoir composed of glacial meltwater as it receives input of seawater, likely on the order of 15–20‰. Support for scenario 1 is provided from the 14CDOC ages obtained from the tabular massive ice, where an increase in age is observed toward the sites at lower elevation (Fig. 7). This would suggest that as the land emerged, the reservoir was being recharged by the melting of the deeper-older Innuitian ice layers (i.e., the late Pleistocene ice layers). However, the lower-elevation sites also show a higher seawater contribution with respect to the Gemini site. Based on the TDS, the Dump and Station Creek sites likely received 1–10% contribution of seawater in the reservoir.
To demonstrate the effect of both scenarios on the value of the δD-δ18O regression slope of the forming ice, we used the Souchez and Jouzel (Reference Souchez and Jouzel1984) model (Eq. 1) and FREEZCH9, an isotope-augmented version of FREZCHEM with a recursive capability to mix water of fixed chemistry and δD-δ18O composition to the remaining reservoir (Faucher et al., Reference Faucher, Lacelle, Fisher, Weisleitner and Andersen2020; Fisher et al., Reference Fisher, Lacelle, Pollard and Faucher2020). The equilibrium fractionation factors for D and 18O during freezing (αi-w) were those of O'Neil (Reference O'Neil1968), and we increased the difference in δ18O values between the input and the reservoir from 0‰ to 25‰ (equivalent of 0‰ to 200‰ for δD). For scenario 1 (recharge by an isotopically depleted source), if |δ18Oi − δ18Or| is small, the regression slope value is near the one expected for freezing; however, if |δ18Oi − δ18Or| increases to >10‰, the regression slope value during freezing approaches that of the meteoric water line (S = 8) (Fig. 8). For scenario 2 (recharge by an isotopically enriched source), if |δ18Oi − δ18Or| is small, the regression slope value is also near the one expected for freezing; however, if |δ18Oi − δ18Or| is 3‰, the regression slope value reaches ~14, and if |δ18Oi − δ18Or| continues to increase, regression slope value approaches that of the meteoric water line. The relative contribution of input to recharge under equilibrium conditions does not affect the findings; for example, the same results are obtained for 99% or 50% of residual water remaining in the reservoir and being recharged to the initial volume. At our three sites in the ESL, the tabular massive ice has co-isotope slope values between 7.3 and 8.4. This would suggest that |δ18Oi − δ18Or| was likely >10‰ and that possibly both scenarios could have occurred for the growth of the tabular massive ground ice bodies.
Numerical modeling of the effect of equilibrium freezing under open-system condition and recharge from an input (δi) with a different δD-δ18O composition relative to that of the initial reservoir (δr) on the value of the δD-δ18O regression slope. To demonstrate the effect of δi > δr and δi < δr on the value of the δD-δ18O regression slope of the forming ice, we used both the Souchez and Jouzel (Reference Souchez and Jouzel1984) model (Eq. 1) and FREEZCH9 (Faucher et al., Reference Faucher, Lacelle, Fisher, Weisleitner and Andersen2020; Fisher et al., Reference Fisher, Lacelle, Pollard and Faucher2020). For the simulations, δ18Or and δDr = −25‰ and −190‰, respectively, and δi was increased or decreased accordingly. The equilibrium fractionation factors for D and 18O during freezing (αi-w) were those of O'Neil (Reference O'Neil1968).
Holocene permafrost evolution in the ESL and the growth of tabular massive ice in the raised marine-deltaic sediments
The Fosheim Peninsula was occupied by the IIS, which reached a thickness of ~1.6 km (Dyke et al., Reference Dyke, Andrews, Clark, England, Miller, Shaw and Veillette2002; England et al., Reference England, Atkinson, Bednarski, Dyke, Hodgson and Ó Cofaigh2006). Deglaciation in the ESL (10,300–8800 14C yr BP [~12,100–9800 cal yr BP]; England et al. Reference England, Atkinson, Bednarski, Dyke, Hodgson and Ó Cofaigh2006) was followed by seawater incursion and the deposition of up to 35 m of marine-deltaic sediments over the weathered sandstone bedrock (8800–7800 14C yr BP [~9800–8600 cal yr BP]; Bell and Hodgson, Reference Bell, Hodgson, Garneau and Alt2000). Marine regression and isostatic rebound later exposed the marine sediments to the atmosphere. Paleotemperature reconstructions from the Agassiz Ice Cap suggest that mean 25-yr air temperature was 1–3°C warmer than today between 7800 and 5000 14C yr BP (~8600–5700 cal yr BP), thus in the −17°C to −15°C range for the ESL (Lecavalier et al., Reference Lecavalier, Fisher, Milne, Vinther, Tarasov, Huybrechts and Lacelle2017). The 14CDOC measurements of ice wedges indicate that the air temperatures were sufficiently cold to allow for the aggradation of permafrost in the marine sediments immediately following marine regression (Campbell-Heaton et al., Reference Campbell-Heaton, Lacelle, Fisher and Pollard2021; Fig. 7). In addition, thermal modeling of permafrost aggradation during the Early to Middle Holocene in marine sediments of a fjord in Svalbard using mean annual air temperature in the −6 to −4°C range showed that permafrost could aggrade tens of meters in a relatively short time period (~200 yr) (Rotem et al., Reference Rotem, Lyakhovsky, Christiansen, Harlavan and Weinstein2022).
The development of tabular massive ice in the ESL raised marine-deltaic sediments involves glacial meltwater recharging an aquifer in the poorly consolidated weathered sandstone beneath the Holocene marine sediments. When ice sheets overlie permeable materials, such as the unconsolidated Tertiary sandstones and shales of the ESL, infiltration of subglacial meltwater into the bed and into the groundwater reservoir is expected (e.g., Boulton and Dobbie, Reference Boulton and Dobbie1993; Lemieux et al., Reference Lemieux, Sudicky, Peltier and Tarasov2008). Although cold-based glaciers are largely frozen to the bed, some may have areas of basal melting (Paterson, Reference Paterson1969). In addition, recent work from the Greenland Ice Sheet showed the importance of moulins for transferring surface meltwater to the subglacial zone (Gulley et al., Reference Gulley, Benn, Müller and Luckman2009; Catania and Neumann, Reference Catania and Neumann2010). For the ESL, the IIS meltwater recharging the aquifer might have been sourced from both surface and basal melt.
The groundwater in the aquifer would then migrate under a hydraulic and thermal gradient to the freezing front as permafrost aggrades in the marine sediments, allowing for the formation of thick tabular bodies of massive ice. However, permafrost aggradation and formation of massive ice in the marine sediments in the ESL would follow a spatiotemporal gradient because of the marine regression and isostatic rebound. Permafrost would first begin to aggrade at the Gemini site (120–125 m asl) between 8500 and 8000 14C yr BP (~9500 and 8900 cal yr BP). At that time, the aquifer would have been recharged by meltwater of Early Holocene–age ice (δ18O: −27‰ to −25‰) and late Pleistocene meltwater (δ18O: −34‰ to −31‰) as the receding IIS was located in proximity to the ESL (Fig. 1). The massive ice at Gemini is very pure (few to no sediment inclusions) and has a thinner sediment overburden (<2 m). The TDS in Gemini ice are the lowest and have nearly the same major ion molar ratios as glacial ice from the Greenland Ice Sheet (Yde et al., Reference Yde, Knudsen, Hasholt and Mikkelsen2014). The fact that marine effects are minimal at this site is likely due to it being close to marine limit and still in proximity to a receding IIS at the time of permafrost aggradation.
Relative sea-level curves place the Station Creek and Dump sites (40–50 m asl) as emerging at approximately 5000 14C yr BP (~5700 cal yr BP) (Fig. 7). The slightly lower δ18O values of the tabular massive ice at these sites would suggest that the groundwater reservoir was decreasing in volume as freezing of water imparts a progressive decrease in δ18O in the residual water. The decreasing reservoir was caused by the receding IIS (by 5000 14C yr BP [~5700 cal yr BP], the IIS was constrained to the mountains along eastern Ellesmere Island) and continued permafrost aggradation at higher elevations. For example, at the elevation of Gemini (120–125 m asl), permafrost has been aggrading for ca. 3000 yr, and it probably reached a few hundred meters thick by 5000 14C yr BP (~5700 cal yr BP). In addition to a decreasing groundwater reservoir, the combination of the IIS being located along the eastern mountains and thick permafrost likely caused a stoppage in glacial meltwater seepage to the reservoir. Instead, the groundwater reservoir was receiving leakage of seawater into the residual aquifer. The tabular massive ice at both Station Creek and Dump sites has Na/Cl molar ratios close to seawater, and TDS concentrations suggest a contribution of seawater on the order of 1–10%. The addition of seawater to the reservoir with higher δ18O values would also explain the near meteoric regression slope values of the massive ice. Overall, the groundwater aquifer in the poorly consolidated sandstone beneath the marine sediments appears to have experienced a switch from seepage of glacial meltwater during the Early Holocene to slow leakage of seawater during the mid-Holocene as a result of permafrost aggradation and isostatic rebound.
CONCLUSIONS
The aim of this study was to assess the origin of tabular massive ground ice in the ESL from a series of sites located between the marine limit and modern coastline. Based on the results, three conclusions can be reached.
First, the cryostratigrahic observations and geochemical composition of the tabular massive ice at the three sites are consistent with an ice segregation origin. However, the massive ice has δD-δ18O regression slope values that are near meteoric. Such co-isotope slope values can be produced in ground ice that formed in an open system with an input that has a significantly different isotopic signature than an initial reservoir.
Second, the development of tabular massive ice in the ESL raised marine-deltaic sediments involves glacial meltwater recharging an aquifer in the poorly consolidated sandstone bedrock beneath the Holocene marine sediments. Based on the low δ18O values, the water recharging the aquifer might have been sourced from both surface and basal melt of the IIS.
Third, massive ground ice in the ESL formed synchronously following marine regression. However, permafrost aggradation and formation of massive ice in the marine sediments in the ESL would follow a spatiotemporal gradient because of the marine regression and isostatic rebound. Near the marine limit, permafrost began to aggrade at 8500–8000 14C yr BP (~9500–8900 cal yr BP), and these sites preserve a geochemical signature similar to glacial meltwater, likely due to their continued proximity to a receding IIS. At modern coastal sites (40 m asl), permafrost began to aggrade near 5000 14C yr BP (~5700 cal yr BP). During that time, the IIS was constrained to the mountains along eastern Ellesmere Island, and permafrost probably reached a few hundred meters thick at the higher-elevation sites. These conditions likely caused a stoppage in glacial meltwater seepage to the reservoir, but instead, the groundwater reservoir was receiving leakage of seawater into the residual aquifer. The tabular massive ice at both the Station Creek and Dump sites has Na/Cl molar ratios close to seawater, and TDS concentrations suggest a contribution of seawater in the order of 1–10%.
This study is relevant to other localities where ice is forming from a reservoir receiving inputs of water of different δD-δ18O composition.
Acknowledgments
We would like to thank the Eureka Weather Station for accommodations. We thank S. Wetterich, an anonymous reviewer, and the associate editor (W. Oswald) for their constructive comments. This project was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grants to DL and WP with logistical support provided by NSERC Northern Supplement, Polar Continental Shelf Program and the Northern Scientific Training Program.
Data Availability Statement
Data are available as supplementary online material.
Supplementary Material
The supplementary material for this article can be found at https://doi.org/10.1017/qua.2023.30
Open access
