Massive ground ice of glacial meltwater origin in raised marine-deltaic sediments, Fosheim Peninsula, high Arctic Canada

Abstract

In the Canadian high Arctic, tabular massive ground ice is found extensively throughout the Eureka Sound Lowlands (ESL). This study evaluates the development of tabular massive ice in raised marine-deltaic sediments of the ESL based on new cryostratigraphic data from sites found between the coastline and the Holocene marine limit. At all sites, massive ice is found below laminated fine-grained marine sediments, and the upper contact between the ice and the overlying marine sediments is conformable and gradational. The concentration of major ions in the massive ice is orders of magnitude higher than expected for glacial ice, but Na/Cl molar ratios vary following elevation: the higher-elevation site has ratios similar to glacial ice, but sites at lower elevations have ratios closer to seawater. The δ¹⁸O values of the ice indicate that the main source of water is glacial meltwater but the δD-δ¹⁸O regression slope values suggest that the ice formed in an open system while receiving an influx that had a substantially different isotopic signature than the initial reservoir. The development of massive ice in the marine-deltaic sediments involves glacial meltwater recharging an aquifer beneath the Holocene marine sediments with a contribution of 1–10% of seawater.

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, 1990; Dyke and Savelle, 2000; Ingólfsson and Lokrantz, 2003; Lacelle et al., 2018; Coulombe et al., 2019). These permafrost terrains are, as a result, vulnerable to thaw due to their high ground ice content (Rudy et al., 2017). 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., 2017). 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, 2000; Lakeman and England, 2012) 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, 1988; Harry et al., 1988; Mackay and Dallimore, 1992; Lacelle et al., 2004). Based on evidence from sites in the Mackenzie delta and the Yukon Coastal Plain, Rampton (1988) 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, 1998). 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, 2000a; Ward Jones et al., 2019). 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, 1991, 2000b; Pollard et al., 2015). Pollard (2000b) 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 (1988) 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-δ¹⁸O regression slope values that were near the meteoric water line, thus leaving a buried glacial ice origin a possibility (Pollard, 1991).

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-δ¹⁸O, and ¹⁴C_DOC 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-δ¹⁸O 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-δ¹⁸O composition.

Figure 1. 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 ¹⁴C yr BP, respectively. Marine regression curves are from Hodgson and Nixon (1998) and Bell and Hodgson (2000). Recessional glacial limits are from England et al. (2006). 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, 1996; Hodgson and Nixon, 1998).

The Late Quaternary glacial history of the northernmost cluster of islands in Canada's Arctic archipelago has long been debated with contrasting views (Hodgson, 1985; Bell, 1996; England et al., 2006). The most recent glacial history reconstruction suggests that the IIS covered most of this region and reached its maximum extent at about 15,800 ¹⁴C yr BP (~19,100 cal yr BP), with the ice being ca. 1.6 km thick in the ESL (Dyke et al., 2002; England et al., 2006). 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., 2006). Deglaciation started between 10,300–8800 ¹⁴C yr BP (~12,100–9800 cal yr BP) (England et al., 2006), and remnants of the ice sheet are still found at the nearby Agassiz Ice Cap. Marine transgression followed until 8800–7800 ¹⁴C yr BP (~9800–8600 cal yr BP) with the Holocene marine limit found between 139 and 150 m asl (Bell, 1996). 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, 1998; Bell and Hodgson, 2000. 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, 2000).

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⁻¹, 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., 2016). Recently disturbed areas are largely vegetation free and covered by efflorescence salts (Kokelj and Lewkowicz, 1999). Holocene air temperatures were reconstructed from δ¹⁸O 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., 2017). Presently, air temperatures are rapidly warming and are now at their warmest in the past 6800–7800 yr (Lecavalier et al., 2017; Copland et al., 2020).

The climate and vegetation conditions ensure that permafrost is cold, continuous, and ~500 m thick (Pollard, 2000a). 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., 2015). 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, 1998). The landscape is dominated by ice wedge polygons, and they occupy ~50% of the ESL (Couture and Pollard, 1998; Bernard-Grand'Maison and Pollard, 2018). Based on ¹⁴C_DOC dating, the ice wedges started to develop immediately after marine regression (Campbell-Heaton et al., 2021). 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, 1998; Pollard, 2000a, 2000b). 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, 2000a, 2000b). 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 ¹⁴C yr BP (~10,000 cal yr BP), while the two lower-elevation sites became exposed at ca. 5000 ¹⁴C yr BP (~5700 cal yr BP) (Fig. 1).

Figure 2. 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-δ¹⁸O, major ions, dissolved organic carbon (DOC), and ¹⁴C_DOC 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, 2005). The volumetric ice content and excess ice contents were calculated following Pollard and French (1980). The concentrations of major cations (Ca²⁺, K⁺, Mg²⁺, Na⁺) and anions (SO₄²⁻, Cl⁻, NO₃⁻) 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 HNO₃ acid. Analytical reproducibility for solute analysis is ±5%.

The ¹⁸O/¹⁶O 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 ¹⁸O/¹⁶O and D/H ratios measurements of H₂O and verified for spectral interference contamination. The results are presented using the δ-notation (δ¹⁸O and δD), where δ represents the parts per thousand differences for ¹⁸O/¹⁶O or D/H in a sample with respect to Vienna Standard Mean Ocean Water (VSMOW). Analytical reproducibility for δ¹⁸O and δD is ±0.3‰ and ±1‰, respectively. Based on the δD and δ¹⁸O values, the D-excess was calculated according to the equation provided in Dansgaard (1964).

For DOC and ¹⁴C_DOC 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 (2003) 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 ¹⁴C analysis, is described in Grinter et al. (2019) and Murseli et al. (2019). Graphitized samples were analyzed on a 3 MV tandem mass spectrometer (Kieser et al., 2015). The ¹⁴C/¹²C ratios are expressed as fraction of modern carbon (F¹⁴C) and corrected for spectrometer and preparation fractionation using the AMS measured ¹³C/¹²C ratio (Crann et al., 2017). Radiocarbon ages are calculated as −8033ln(F¹⁴C) and reported in ¹⁴C yr BP (BP = AD 1950) (Stuiver and Polach, 1977). In the text, radiocarbon ages with a standard error between 50 and 1000 yr are rounded to the nearest 10 yr. Our ¹⁴C_DOC ages were calibrated using OxCal and the IntCal20 calibration curve (Bronk Ramsey, 2009; Reimer et al., 2020). 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).

Figure 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⁻¹, respectively). The TDS variation in the 1.3 m core from Dump remains rather constant with depth (60–200 mg L⁻¹), except at 0.4–0.5 m depth, where TDS increase to 8830 mg L⁻¹ due to a higher sediment content. At Station Creek, the TDS remain near 200 mg L⁻¹ in the uppermost 0.75 m, then progressively increase to 1200 mg L⁻¹ at 1.5 m depth. The ice at the higher-elevation site (Gemini, 115 to 120 m asl) is characterized a Ca-SO₄ geochemical facies with lower TDS (average of 28 ± 15 and 37 ± 10 mg L⁻¹, 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.

Figure 4. Volumetric ice content (fVi), δ¹⁸O, 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).

Figure 5. 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⁻, SO₄²⁻, and Na/Cl; (B) Piper diagram of Mg²⁺, Ca²⁺, 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-δ¹⁸O composition of the ice also varies by site elevation (Fig. 6). The higher-elevation site (Gemini) has the highest average δ¹⁸O composition (−30.9 ± 0.4‰), whereas the sites near the coast (Dump and Station Creek) have lower average δ¹⁸O composition (−32.3 ± 0.7‰ and −31.1 ± 1.6‰, respectively). At Dump, the δ¹⁸O 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 δ¹⁸O 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 δ¹⁸O (δD = 7.5 δ¹⁸O − 5.4; r² = 0.98). At Gemini, two cores were collected at ~15 m vertical distance. The δ¹⁸O 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-δ¹⁸O slopes of G3 and G1 are 8.4 (δD = 8.4 δ¹⁸O + 22.7; r² = 0.96) and 7.5 (δD = 7.5 δ¹⁸O − 7.8; r² = 0.85), respectively.

Figure 6. (A) δD-δ¹⁸O 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, δ¹⁸O − 9.1; IAEA/WMO, 2015). (B) Na/Cl and δ¹⁸O 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⁻¹, with δ¹³C_DOC 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 ¹⁴C_DOC ages at Gemini and Dump sites, ranging from 26,150 ± 330 (30,460 cal yr BP) to 8850 ± 45 ¹⁴C yr BP (9955 cal yr BP) (Table 1), are older than the timing of deglaciation and marine regression (Fig. 7) (Note that the ¹⁴C_DOC ages reflect that of the DOC in the source(s) of water and not the timing of freezing.) The ¹⁴C_DOC ages show a general trend of increasing ages toward lower elevation.

Figure 7. Comparison of ¹⁴C_DOC 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 (1998) and Bell and Hodgson (2000). ¹⁴C_DOC ages of ice wedges are from Campbell-Heaton et al. (2021).

Table 1. Dissolved organic carbon (DOC), δ¹³C_DOC, and ¹⁴C_DOC results for samples of tabular massive ice collected at three sites in Eureka Sound Lowlands, high Arctic Canada.^a

^a The ¹⁴C ages were calibrated using OxCal and the IntCal20 calibration curve (Bronk Ramsey, 2009; Reimer et al., 2020).

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-δ¹⁸O 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, 1998; French and Shur, 2010; Lacelle and Vasil'chuk, 2013). 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., 1995; Koerner et al., 1999). 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 δ¹⁸O 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-δ¹⁸O regression slope values between 4 and 6, with the value largely depending on isotopic composition of initial water and freezing rate (Souchez and Jouzel, 1984; Lacelle, 2011). 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-δ¹⁸O 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 (1984) suggested that the value of the freezing slope (S) from an initial reservoir (r) that receives input of water (i) could be estimated from:

(1)$$S = \displaystyle{{{\rm \alpha }\left[{( {{\rm \alpha \;}-{\rm \;}1} ) ( {1000{\rm \;} + {\rm \;\delta Dr}} ) {\rm \;}-{\rm \;}\left({\displaystyle{{\rm A} \over S}} \right)( {{\rm \delta Di\;}-{\rm \;\delta Dr}} ) } \right]} \over {{\rm \beta }[ {{\rm \beta \;}-{\rm \;}1} ) ( {1000{\rm \;} + {\rm \;}{\rm \delta }^{18}{\rm Or}} ) -( {{\rm A}/S} ) ( {{\rm \delta }^{18}{\rm Oi\;}-{\rm \;}{\rm \delta }^{18}{\rm Or}} ) ] }}$$

where β = fractionation factor of ¹⁸O ice water, α = fractionation factor of D ice water, δDi = δD composition of input, δDr = δD composition of reservoir, δ¹⁸Oi = δ¹⁸O composition of input, δ¹⁸Or = δ¹⁸O 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 δ¹⁸O of the reservoir (which influences the value of the fractionation factor) and the difference in δD and δ¹⁸O values between the input and the reservoir. In most natural systems, the difference in δD and δ¹⁸O 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 (1985) found basal ice at Gruben-gletscher in Valais, western Switzerland, with a δD-δ¹⁸O slope of 8.3, much higher than predicted for open-system freezing but similar to the meteoric water line. Souchez and de Groote (1985) simulated the evolution of δD and δ¹⁸O 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-δ¹⁸O 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., 2011).

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-δ¹⁸O composition. However, there is one difference with Souchez and de Groote (1985). 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-δ¹⁸O 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., δ¹⁸O transition from the Early Holocene to late Pleistocene age ice at Agassiz Ice Cap is in the order of 7–10‰; Fisher et al., 1995; Zdanowicz et al., 2002) 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 ¹⁴C_DOC 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-δ¹⁸O regression slope of the forming ice, we used the Souchez and Jouzel (1984) model (Eq. 1) and FREEZCH9, an isotope-augmented version of FREZCHEM with a recursive capability to mix water of fixed chemistry and δD-δ¹⁸O composition to the remaining reservoir (Faucher et al., 2020; Fisher et al., 2020). The equilibrium fractionation factors for D and ¹⁸O during freezing (α_i-w) were those of O'Neil (1968), and we increased the difference in δ¹⁸O 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 |δ¹⁸Oi − δ¹⁸Or| is small, the regression slope value is near the one expected for freezing; however, if |δ¹⁸Oi − δ¹⁸Or| 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 |δ¹⁸Oi − δ¹⁸Or| is small, the regression slope value is also near the one expected for freezing; however, if |δ¹⁸Oi − δ¹⁸Or| is 3‰, the regression slope value reaches ~14, and if |δ¹⁸Oi − δ¹⁸Or| 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 |δ¹⁸Oi − δ¹⁸Or| was likely >10‰ and that possibly both scenarios could have occurred for the growth of the tabular massive ground ice bodies.

Figure 8. Numerical modeling of the effect of equilibrium freezing under open-system condition and recharge from an input (δi) with a different δD-δ¹⁸O composition relative to that of the initial reservoir (δr) on the value of the δD-δ¹⁸O regression slope. To demonstrate the effect of δi > δr and δi < δr on the value of the δD-δ¹⁸O regression slope of the forming ice, we used both the Souchez and Jouzel (1984) model (Eq. 1) and FREEZCH9 (Faucher et al., 2020; Fisher et al., 2020). For the simulations, δ¹⁸Or and δDr = −25‰ and −190‰, respectively, and δi was increased or decreased accordingly. The equilibrium fractionation factors for D and ¹⁸O during freezing (α_i-w) were those of O'Neil (1968).

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., 2002; England et al., 2006). Deglaciation in the ESL (10,300–8800 ¹⁴C yr BP [~12,100–9800 cal yr BP]; England et al. 2006) was followed by seawater incursion and the deposition of up to 35 m of marine-deltaic sediments over the weathered sandstone bedrock (8800–7800 ¹⁴C yr BP [~9800–8600 cal yr BP]; Bell and Hodgson, 2000). 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 ¹⁴C yr BP (~8600–5700 cal yr BP), thus in the −17°C to −15°C range for the ESL (Lecavalier et al., 2017). The ¹⁴C_DOC 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., 2021; 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., 2022).

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, 1993; Lemieux et al., 2008). Although cold-based glaciers are largely frozen to the bed, some may have areas of basal melting (Paterson, 1969). 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., 2009; Catania and Neumann, 2010). 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 ¹⁴C yr BP (~9500 and 8900 cal yr BP). At that time, the aquifer would have been recharged by meltwater of Early Holocene–age ice (δ¹⁸O: −27‰ to −25‰) and late Pleistocene meltwater (δ¹⁸O: −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., 2014). 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 ¹⁴C yr BP (~5700 cal yr BP) (Fig. 7). The slightly lower δ¹⁸O 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 δ¹⁸O in the residual water. The decreasing reservoir was caused by the receding IIS (by 5000 ¹⁴C 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 ¹⁴C 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 δ¹⁸O 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-δ¹⁸O 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 δ¹⁸O 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 ¹⁴C 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 ¹⁴C 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-δ¹⁸O composition.