Glacial erosion and history of Inglefield Land, northwestern Greenland

We used mapping of bedrock lithology, bedrock fractures, and lake density in Inglefield Land, northwestern Greenland, combined with cosmogenic nuclide (¹⁰Be and ²⁶Al) measurements in bedrock surfaces, to investigate glacial erosion and the ice sheet history of the northwestern Greenland Ice Sheet. The pattern of eroded versus weathered bedrock surfaces and other glacial erosion indicators reveal temporally and spatially varying erosion under cold- and warm-based ice. All of the bedrock surfaces that we measured in Inglefield Land contain cosmogenic nuclide inheritance with apparent ¹⁰Be ages ranging from 24.9 ± 0.5 to 215.8 ± 7.4 ka. The ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios require minimum combined surface burial and exposure histories of ∼ 150 to 2000 kyr. Because our sample sites span a relatively small area that experienced a similar ice sheet history, we attribute differences in nuclide concentrations and ratios to varying erosion during the Quaternary. We show that an ice sheet history with ∼ 900 kyr of exposure and ∼ 1800 kyr of ice cover throughout the Quaternary is consistent with the measured nuclide concentrations in most samples when sample-specific subaerial erosion rates are between 0 and 2 × 10⁻² mm yr⁻¹ and subglacial erosion rates are between 0 and 2 × 10⁻³ mm yr⁻¹. These erosion rates help to characterize Arctic landscape evolution in crystalline bedrock terrains in areas away from focused ice flow.

1 Introduction

The Greenland Ice Sheet currently contributes more to sea level rise than any other single ice mass and is predicted to continue to melt at an accelerated rate throughout the next century (e.g., Aschwanden and Brinkerhoff, 2022; Goelzer et al., 2020). Ice sheet models used to predict future ice sheet evolution are aided by knowledge of the long-term ice sheet history and patterns of past ice flow variability, as these play a large role in modulating ice sheet mass balance (e.g., Hubbard et al., 2009). Uncertainty in ice sheet model parameters can be lessened when paleo ice sheet simulations are performed alongside geologic constraints (e.g., Briner et al., 2020; Cuzzone et al., 2016; Patton et al., 2017). Observations of many ice sheet processes are sparse because it is difficult to access the bed of modern ice sheets; thus, investigating the beds of former ice sheets (i.e., glaciated landscapes) provides information on past ice sheet dynamics and basal processes such as glacial erosion.

The distribution of glacial erosional features at the ice sheet scale, including striations, lakes, lateral meltwater channels, and sculpted bedrock across glaciated landscapes, has been used to map past basal thermal conditions and relative ice velocities delineating zones of warm-bedded and/or fast-flowing ice where erosional features are abundant from areas of cold-based ice and ancient (Quaternary or pre-Quaternary) landscape preservation where these features are absent (Andrews et al., 1985; Daly, 1902; Flint, 1943; Margold et al., 2015; Sugden, 1978). Identifying areas of differential erosion in glaciated landscapes is particularly useful for mapping paleo ice streams, as knowing their extent is helpful for understanding the mass balance of former ice sheets through ice sheet models. Ice streams play a large role in modulating the volume of the modern Greenland Ice Sheet, and their stability is directly linked to overall ice sheet mass balance (e.g., Khan et al., 2022; Mouginot et al., 2015). These features are large areas of fast-moving ice, with onset zones that feed into main channels, which eventually calve into the ocean and are constrained by topography or slow-moving ice (Benn and Evans, 2014). Paleo ice stream onset zones (i.e., the heads of ice streams) are found at the transitions between warm- and cold-based ice, which can be delineated via mapping of erosion imprints in a landscape and cosmogenic nuclide analysis (Briner et al., 2008; Margold et al., 2015, 2018).

Cosmogenic nuclides are produced mostly in the upper few meters of Earth's surface when it is exposed to the cosmic ray flux, and they are routinely used in glaciated landscapes to quantify erosion rates and the timing of past ice sheet fluctuations (Gosse and Phillips, 2001). In areas covered by warm-based, highly erosive glaciers, nuclides that have accumulated in the upper 2–3 m of bedrock are often removed through efficient glacial erosion. In areas of minimal to no glacial erosion (i.e., covered by cold-based glaciers or short-lived erosive ice), bedrock surfaces contain inventories of cosmogenic nuclides from multiple periods of exposure (known as cosmogenic nuclide inheritance; Bierman et al., 1999). Patterns of inheritance (or a lack thereof) across a landscape can be used to delineate ice streaming where faster-flowing, erosive ice depleted nuclide inventories (Briner et al., 2006; Corbett et al., 2013; Roberts et al., 2013).

Cosmogenic nuclide inheritance resulting from minimal glacial erosion can pose an issue for single-nuclide surface exposure dating, where the apparent exposure age will be anomalously old (Ivy-Ochs and Briner, 2014). However, some cosmogenic nuclides are radioactive and decay when they are shielded from the cosmic ray flux by ice (e.g., ¹⁰Be and ²⁶Al: no significant production beneath tens of meters of ice; Miller et al., 2006). The production ratio of ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ is 7.3:1 in quartz at Earth's surface across Greenland (Corbett et al., 2017). Because ²⁶Al has a shorter half-life (705 kyr) than ¹⁰Be (1388 kyr), departure from this ratio indicates surface burial (Korschinek et al., 2010; Nishiizumi, 2004). By measuring multiple nuclides in bedrock surfaces near the modern Greenland Ice Sheet margin, researchers have exploited cosmogenic nuclide inheritance to investigate periods of ice sheet minima, glacial erosion, and longer-term ice sheet history (e.g., Corbett et al., 2013; Skov et al., 2020; Young et al., 2021).

Recently, attention has been drawn to retrieving samples from the bed of modern extant ice sheets (i.e., the Greenland and Antarctic ice sheets; Briner et al., 2022; Johnson et al., 2025; Spector et al., 2018). The information contained at the contemporary ice–bed interface provides valuable, and rare, terrestrial constraints on previous ice sheet minima and long-term ice sheet histories. Cosmogenic nuclide and luminescence analysis of sediment and bedrock samples collected from the modern glacier bed have provided insight into the stability of the ice sheet throughout the Quaternary (Balco et al., 2023; Bierman et al., 2024; Christ et al., 2020; Christ et al., 2021; Christ et al., 2023; Schaefer et al., 2016). Studying the landscape at the fringes of the Greenland Ice Sheet allows us to systematically investigate large areas of the former ice sheet bed, without the need to drill through the ice sheet, providing complementary results for efforts focused on obtaining material from under the ice. We mapped bedrock features and used cosmogenic nuclide analysis to determine the extent and magnitude of erosion and the ice sheet history across Inglefield Land, northwestern Greenland, over the Quaternary.

2 Inglefield Land

Inglefield Land is an ice-free area in northwestern Greenland situated between the Greenland Ice Sheet and the coastline 30 km to the northwest. It is bounded by Smith Sound to the west, Kane Basin to the north, Prudhoe Dome and the main body of the northern Greenland Ice Sheet to the south, and the Hiawatha Glacier sector of the Greenland Ice Sheet to the east (Fig. 1). Today, much of the ice bordering Inglefield Land is cold-based, and ice velocities are low, ranging from 10 to ∼ 120 m yr⁻¹ (Fig. 1; Joughin et al., 2018; MacGregor et al., 2022). Inglefield Land is characterized by relatively low-relief uplands that reach 700 m a.s.l. near the ice margin with valleys incised along the northern coast. Ice sheet meltwater drains into these valleys and into four embayments (from west to east): Force Bay, Rensselaer Bay, Marshall Bay, and Dallas Bay (henceforth the unnamed valleys crossing Inglefield Land will be referred to by the bays into which they drain).

Figure 1 Overview map of Inglefield Land and key geographic locations. Ice velocity from Joughin et al. (2018); 25 m DEM from Korsgaard et al. (2016).

During the Last Glacial Maximum (LGM; 26–19 ka), the northern Greenland Ice Sheet covered Inglefield Land completely as it advanced into Kane Basin, where it coalesced with the Innuitian Ice Sheet and flowed southward, eventually terminating in an ice shelf spanning northern Baffin Bay (England, 1999; Couette et al., 2022; Batchelor et al., 2024). Although the timing and duration of the LGM ice advance remain unknown, retreat onto the modern coast of Inglefield Land is constrained to ∼ 8.6–7.9 ka based on radiocarbon dating of organic material in raised beach deposits and in situ ¹⁴C ages from erratic boulders (Blake et al., 1992; Mason, 2010; Nichols, 1969; Søndergaard et al., 2020). The ice sheet continued to decay, arriving at the modern margin by ∼ 7 ka and maintaining a smaller-than-modern position between ∼ 5.8 and 0.3 ka, an estimate based on radiocarbon ages from reworked wood fragments at the modern ice margin (Søndergaard et al., 2020).

The pre-Holocene ice sheet history and dynamics of Inglefield Land are not known, though cosmogenic nuclide inheritance (¹⁰Be) found in boulders across Inglefield Land indicates that much of the landscape was covered by cold-based ice during the last glacial cycle (Søndergaard et al., 2020). Limited ice sheet terrestrial records suggest that there were major deglaciation events during the Quaternary as captured by cosmogenic nuclide and luminescence analysis of sub-ice sediment and bedrock at Camp Century and Summit Camp, during which Inglefield Land was almost certainly ice-free (Christ et al., 2021; Christ et al., 2023; Schaefer et al., 2016). Additionally, offshore marine sediment records imply that the ice sheet underwent several cycles of advance and retreat, though the inland extent of ice sheet recession in northern Greenland during interglacials is poorly constrained (Bierman et al., 2016; Colville et al., 2011; Hatfield et al., 2016; Knutz et al., 2019; Reyes et al., 2014).

3 Methods

We used mapping, field observations, and ¹⁰Be and ²⁶Al measurements from bedrock surfaces at 11 sites to assess the ice sheet history and erosive conditions of the northwestern Greenland Ice Sheet across Inglefield Land. Mapping and field observations allowed us to broadly identify variable subglacial and subaerial erosion across the landscape, while cosmogenic nuclide measurements helped us quantify the ice sheet history and erosion.

3.1 Geologic mapping – GIS and field observations

3.1.1 Bedrock lithology

The relatively simple bedrock geology of Inglefield Land allowed us to investigate the long-term pattern of erosion and landscape evolution of the area using pre-existing geologic maps. Inglefield Land is underlain by crystalline paragneiss and capped with near-horizontally bedded sedimentary rocks; therefore, generally speaking, outcrops of basement rock indicate areas where cap rocks have been removed by erosion (Fig. 2). Some of this cap rock removal likely took place prior to Quaternary glaciation (e.g., Krabbendam and Bradwell, 2014), but the patterns seen in the geologic map – with crystalline lithologies in glacial troughs – hint at the role of past glacial erosion in shaping the landscape. To identify the removal of sedimentary cap rocks, we created a simplified geologic map of Inglefield Land using the 1:500 000 Greenland-wide geologic map in ArcGIS Pro (Kokfelt et al., 2023) and classified units as crystalline, sedimentary, Quaternary sediments, or lakes. In the field, we noted the lithology at each of our sample sites and whether bedrock outcrops exhibited evidence of ice sculpting (e.g., sculpted bedrock, striations, and glacial polish) or weathering (e.g., grussification and weathering pits; Figs. 2 and 3).

Figure 2 Simplified geologic map of Inglefield Land (Kokfelt et al., 2023).

Figure 3 Representative images of bedrock surfaces at sample sites in Inglefield Land. The photo locations are shown in Fig. 1.

3.1.2 Mapping bedrock fractures

Bedrock fractures in glaciated terrains are exposed through the removal of regolith by glacial erosion (Gordon, 1981; Skyttä et al., 2023; Sugden, 1974, 1978). Areas with a high density of exposed bedrock fractures within crystalline bedrock areas have been used to identify intense glacial erosion (e.g., Sugden, 1978). This is in contrast to areas where bedrock fractures are obscured by sediment cover, a regolith mantle, or, in some cases, subhorizontal sedimentary bedrock units. Our field area has a heterogenous pattern of mappable fractures that are easily identifiable at the landscape level (>25 m), which is indicative of areas where regolith has been stripped via glacial erosion. We outlined zones of bedrock fractures at the landscape scale in ArcGIS Pro using a 25 m resolution digital elevation model (Korsgaard et al., 2016). We created hillshade images with a 3 times vertical exaggeration for our mapping (Fig. 4). Our digital elevation model resolution of 25 m was fine enough to capture areas of obvious bedrock fractures on a landscape scale but coarse enough to avoid mistaking smaller fractures for other landscape features.

Figure 4 Mapped fracture zones in Inglefield Land. The inset shows a zoom-in of the fracture zone in the Rensselaer Valley.

3.1.3 Mapping lake density

Lake density has also been used as a proxy for glacial erosion (e.g., Andrews et al., 1985; Briner et al., 2008; Sugden, 1978). Areas of regolith cover are generally topographically smooth at the landscape scale and contain few lakes. Meanwhile, areas of high erosion rates beneath ice sheets can remove this overlying regolith and expose the underlying bedrock, after which the bedrock is susceptible to ice sheet scouring and erosion. These bedrock basins then fill with water following glacial retreat, forming lakes across glaciated areas, and thus the density of lakes can be used as an indicator of past ice sheet erosion. We created an inventory of lakes in Inglefield Land to calculate lake density. While the 1:500 000 geologic map of Greenland shows larger lakes, smaller lakes are not included due to the relatively coarse map resolution (Fig. 2; Kokfelt et al., 2023). Therefore, we used a semi-automated process in ArcGIS Pro to map all of the lakes in the study area. First, we used cloud-free LANDSAT8 images to visualize the surface water (Band 5; visible blue–green). We explored a range of threshold values to extract cells with surface water from this image (Fig. 5; cells with a value of less than the threshold of 8000 were water, though we note that this value may not be suitable for any location) and evaluated these against the original LANDSAT8 images to determine a suitable threshold value that adequately captured lakes and streams. We converted this raster to polygons of lake and river extents. Finally, we conducted manual quality control by removing rivers and any lakes that were by ice sheets and then merged our new lake polygons with those from the geologic map. To calculate lake density across Inglefield Land, we determined the percentage of each cell in a 1 × 1 km grid covered by lake polygons.

Figure 5 LANDSAT Band-6 image (a) of the portion of Inglefield Land in Fig. 1, showing the workflow of automatically selecting surface water with a threshold value of 8000 (b) and manually removing rivers (c).

3.2 Cosmogenic nuclide measurements

3.2.1 Sampling approach

We collected samples for cosmogenic nuclide measurements from bedrock surfaces along two SSE–NNW transects from the ice margin to the coast in Inglefield Land (Fig. 1; Table S1 in the Supplement). We sampled roughly 1 kg of crystalline rock from 15 bedrock surfaces and one boulder using a handheld angle grinder, hammer, and chisel during the summer of 2022. We collected nine samples along our western transect (beginning at the mouth of the Rensselaer Valley): two from bedrock surfaces and one from a boulder (2.5 m long × 1.6 m wide × 1.1 m tall) close to the ice margin, one bedrock sample from each of three separate inland sites, one bedrock sample from a higher-elevation coastal site (255 m a.s.l.), and two bedrock samples from a lower-elevation coastal site (∼ 100 m a.s.l.). Along our eastern transect (beginning at the mouth of the western Marshall Valley), we collected bedrock samples: three from near the ice margin, one from each of three different inland sites, and one from a low-elevation coastal site (∼ 100 m a.s.l.).

3.2.2 Laboratory procedures

We measured ¹⁰Be and ²⁶Al in all 16 of our samples. We isolated quartz and extracted ¹⁰Be and ²⁶Al at the University at Buffalo cosmogenic isotope laboratory (n=12) and the Lamont-Doherty Earth Observatory (LDEO) cosmogenic dating laboratory (n=4) using well-established procedures (Corbett et al., 2016b; Kohl and Nishiizumi, 1992). We spiked our dissolved samples with precisely weighed ⁹Be carrier (PRIME Lab 2017.11.17-Be #3/#4 – ⁹Be concentration of 1074 ± 8 ppm – at the University at Buffalo and LDEO carrier – ⁹Be concentration of 1038.8 ppm – at LDEO). We measured the amount of native ²⁷Al in our dissolved quartz and added varying amounts of ²⁷Al carrier to ensure that each sample had ∼ 2000 mg of ²⁷Al. We measured the total amount of ²⁷Al in aliquots removed after sample digestion with inductively coupled plasma optical emission spectrometry. We sent ¹⁰Be samples processed at the University at Buffalo and all ²⁶Al samples to PRIME Lab and ¹⁰Be samples processed at Lamont-Doherty to Lawrence Livermore National Laboratory, where the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ and ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{27}}\mathrm{Al}$ ratios were measured with accelerator mass spectrometry.

${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratios were measured relative to the 07KNSTD standard (${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratio: 2.85 × 10⁻¹²; Nishiizumi et al., 2007) at both facilities and ²⁶Al samples relative to the KNSTD standard (${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{27}}\mathrm{Al}$ ratio: 1.82 × 10⁻¹²; Nishiizumi, 2004). Analytical uncertainties (1σ) of ¹⁰Be measurements at Lawrence Livermore were between 1.8 % and 1.9 % and ranged from 1.5 % to 6.0 %, with an average of 2.6 ± 1.3 % at PRIME Lab. ²⁶Al measurement uncertainties ranged from 3.6 % to 8.4 % with an average of 4.8 ± 1.5 %. We applied background corrections to both ¹⁰Be and ²⁶Al sample ratios using batch-specific process blank values (Table S1).

We calculated apparent exposure ages using version 3 of the online exposure age calculator of Balco et al. (2008). We used the Arctic ¹⁰Be production rate with a time-independent “St” production rate scaling method to calculate apparent exposure ages (Lal, 1991; Stone, 2000; Young et al., 2013). We did not include production rate uncertainties in our exposure age calculations because we do not compare our age results to other independent dating methods. We used a Greenland-specific ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ surface production ratio in quartz of 7.3 to interpret the exposure and burial histories of both nuclides (Corbett et al., 2017).

3.3 Modeling cosmogenic nuclide accumulation in rock surfaces

In bedrock samples with cosmogenic nuclide inheritance, measurements of two isotopes with different half-lives have been used to calculate Greenland Ice Sheet exposure and burial histories (e.g., Andersen et al., 2020; Beel et al., 2016; Corbett et al., 2013; Knudsen and Egholm, 2018; Knudsen et al., 2015; Skov et al., 2020; Strunk et al., 2017). Although in situ ¹⁴C ages from boulders constrain ice retreat across Inglefield Land to between ∼ 9 and 7 ka, apparent ¹⁰Be ages from boulders date to between 8.3 ± 1.2 and 92.7 ± 1.5 ka, indicating the presence of nuclide inheritance (Søndergaard et al., 2020). Inglefield Land bedrock is also likely to contain nuclide inheritance and therefore yield information about ice sheet history and erosion prior to the last glacial cycle. Given the relatively short distance from the modern ice margin to the coast (<50 km) and the speed at which Inglefield Land deglaciated following the LGM (<2 kyr), we hypothesize that, within ¹⁰Be and ²⁶Al measurement uncertainties, the 15 bedrock samples should have similar ice cover histories on glacial–interglacial timescales (Søndergaard et al., 2020). Therefore, differences in ¹⁰Be and ²⁶Al concentrations across the landscape likely relate to varying sample-to-sample subglacial and subaerial erosion rates and not differences in ice sheet history at each sample location.

To explore this hypothesis, we simulate complex Pleistocene exposure histories using a forward model that calculates cosmogenic ¹⁰Be and ²⁶Al accumulation in rock brought to Earth's surface by subaerial and subglacial erosion. Because little is known about the ice margin history in Inglefield Land prior to the Holocene, we define the pre-Holocene exposure and burial history by applying a threshold value to the benthic δ¹⁸O LR04 stack (Lisiecki and Raymo, 2005) to define a range of plausible exposure and burial scenarios for the last 2.7 Myr, following the approach adopted in prior studies (Balter-Kennedy et al., 2021; Knudsen et al., 2015). Exposure and burial take place at δ¹⁸O values below and above the threshold, respectively. We use a 30 kyr running mean to smooth the δ¹⁸O curve (Knudsen et al., 2015). Given the prevalence of nuclide inheritance across Inglefield Land (Søndergaard et al., 2020), we do not expect ¹⁰Be to give post-LGM deglaciation ages at each sample location, but deglaciation ages from in situ ¹⁴C in boulders at the coast and ice margin are provided by Søndergaard et al. (2020). We therefore estimate site-specific Holocene exposure durations by scaling the deglaciation ages from Søndergaard et al. (2020) to each of our sites based on their relative distance between the coast and ice margin, as calculated along our sample transects. This Holocene exposure scaling does not account for variability in Holocene lateral retreat rates. Ultimately, however, the goal of this exercise is to determine whether a long-term exposure history can explain the measured cosmogenic nuclide concentrations, which is not drastically affected by slight variations in the Holocene exposure history.

The cosmogenic nuclide concentration (atoms per gram), N, per nuclide, i, at the end of each time step, j, is the sum of nuclides inherited from the previous time step (adjusted for radioactive decay) and the new accumulation of nuclides (assumed to be 0 when ice-covered and limited by the surface erosion rate when ice-free):

$$N_{i,j}=N_{i,j-\mathrm{1}}{e}^{-{\mathit{\lambda }}_i{t}_j}+{\int }_{\mathrm{0}}^{t_j}{P}_i(z_{\mathrm{end},j}+{\mathit{\epsilon }}_j\mathit{\tau })e^{-{\mathit{\lambda }}_i\mathit{\tau }}\mathrm{d}\mathit{\tau },$$

where λ is the decay constant (5.00 × 10⁻⁷ for ¹⁰Be; 9.83 × 10⁻⁷ for ²⁶Al) (Korschinek et al., 2010; Nishiizumi, 2004) and t_j is the duration of the time step (years). P_i is the sum of nuclide production (atoms per gram per year) by spallation and muon interactions at a given mass depth (g cm⁻²). Here, mass depth is time-varying and controlled by subaerial and subglacial erosion during the ice-free and ice-covered time steps, respectively. Therefore, the mass depth is defined by the sample depth at the end of each time step, z_end,j, and the erosion rate (g cm⁻² yr⁻¹), ε_j, during that time step. We calculate spallation production rates at Earth's surface using the Arctic ¹⁰Be calibration dataset of Young et al. (2013), a ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratio of 7.3 (Corbett et al., 2017; Young et al., 2021), and the scaling method of Stone (2000). We assume that spallation production decreases exponentially with mass depth at an attenuation length of 160 g cm⁻². We calculate production by muon interactions in MATLAB using the cross sections for ¹⁰Be and ²⁶Al determined by Balco (2017) and implemented in Model 1A in the same reference.

We determine the misfit between the modeled and measured nuclide concentrations for each sample, d, using the error-weighted sum of squares (EWSS):

$$\mathrm{EWSS}={\left({\displaystyle \frac{N_{\mathrm{10},p,d}-N_{\mathrm{10},m,d}}{{\mathit{\sigma }}_{\mathrm{10},m,d}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{N_{\mathrm{26},p,d}-N_{\mathrm{26},m,d}}{{\mathit{\sigma }}_{\mathrm{26},m,d}}}\right)}^{\mathrm{2}},$$

where $N_{i,p,d}$ is the predicted nuclide concentration, $N_{i,m,d}$ is the measured nuclide concentration, and ${\mathit{\sigma }}_{i,m,d}$ is the 1σ measurement uncertainty. An error-weighted sum of squares close to 2 indicates that the difference between the modeled and measured concentrations can be explained by the measurement uncertainty. We consider model runs with error-weighted sums of squares of less than 2.5 to be acceptable fits.

We created 2.7 Myr exposure histories using δ¹⁸O threshold values ranging from 3.60 ‰ to 4.00 ‰ (corresponding to 0.7–1.9 Myr cumulative exposure and 0.8–2 Myr cumulative burial over the last 2.7 Myr) at 0.02 ‰ spacing. We decouple this δ¹⁸O threshold-based exposure timing during Holocene deglaciation (Knudsen et al., 2015), as the timing of Holocene exposure is relatively well-known from the constraints of Søndergaard et al. (2020). For each exposure history, we ran the forward model with subaerial and subglacial erosion rates ranging from 0 to 2.5 × 10⁻¹ mm yr⁻¹ on a log scale (coarse spacing from 0 to 1 × 10⁻⁵ mm yr⁻¹ and finer spacing from 1 × 10⁻⁵ to 2.5 × 10⁻¹ mm yr⁻¹) to capture the potential range in subglacial erosion rates for cold-bedded glaciers and subaerial erosion rates for polar environments (e.g., Cook et al., 2020; Koppes et al., 2015; Portenga and Bierman, 2011). Subglacial and subaerial erosion rates are held constant for all periods of ice cover and exposure, respectively, and therefore we do not consider changes in erosion within glacial cycles or from one glacial cycle to the next. If the total (subaerial + subglacial) erosion through a model run is high, a higher proportion of modeled ¹⁰Be and ²⁶Al accumulates deep in the rock column where production is low, while less erosion results in a larger fraction of the nuclide production near the surface where production is higher. The two erosion rates may therefore trade off. For example, an exposure history may yield a good fit to the data, with a higher subaerial erosion rate and a lower subglacial erosion rate for a certain sample as well as a lower subglacial erosion rate and a higher subaerial erosion rate. Thus, this model allows us to investigate potential Quaternary erosion and ice cover scenarios across Inglefield Land to test the hypothesis that the variability in our cosmogenic nuclide measurements can be explained with a common 2.7 Myr exposure–burial history and sample-specific subaerial and subglacial erosion rates.

4 Results

4.1 Mapping

Paleoproterozoic crystalline paragneiss basement rocks are overlain by subhorizontally bedded Mesoproterozoic to Ordovician sedimentary rocks across Inglefield Land (Fig. 2). Near the ice margin, crystalline basement rocks are exposed, with some areas covered by Quaternary sediments. Along much of the coastline of Inglefield Land, sedimentary rock units extend up to 20 km inland, except where erosion incised the cap rock to expose the underlying basement rocks in the four main valleys. North of Hiawatha Glacier, crystalline basement outcrops from the ice margin to the coastline. In total, 48.2 % of the surface of Inglefield Land is crystalline basement, 29.4 % is sedimentary cap rocks, and 17.4 % is Quaternary sediments, with the remaining 5 % made up of lakes and minor geologic units.

The low-elevation coastal (∼ 100 m a.s.l.) sites exhibit ice-sculpted bedrock with primary glacial erosional features, including striations (Fig. 3). Smoothed bedrock is also present at the higher-elevation coastal site (255 m a.s.l.), though it lacked primary glacial erosional features such as striations or glacial polish. Inland, highly weathered bedrock outcrops, often with abundant gruss, rose only a few meters above the surrounding landscape, which is covered by Quaternary sediments composed of weathered, angular boulders. Near the ice margin sites, we also found highly weathered bedrock outcrops exhibiting lots of gruss and weathering pits on rock surfaces.

Exposed bedrock fractures are sparse across Inglefield Land, except for the areas west and north of Hiawatha Glacier where bedrock fractures are clearly visible at the kilometer scale, as mapped in Fig. 4. There are sizable fracture zones along the Rensselaer Valley and both Marshall valleys; in particular, the fracture zone of the western Marshall Valley extends ∼ 25 km inland. Finally, there are some limited fracture areas in the central inland sector towards the ice sheet.

Lakes are abundant across much of Inglefield Land, covering 161 km² (2 % of the surface area), but they are concentrated in certain sectors (Fig. 6). Regions to the west and north of Hiawatha Glacier have the highest lake densities, with some 1 km² cells 87 % covered by lakes. There are also areas of high lake densities towards the coast, particularly near the coastal valleys. Inland, the highest density of lakes is found in central Inglefield Land. Many grid cells do not contain lakes (grid cells not filled in Fig. 6) or host very small lakes (<1 % of grid cells lake-covered).

Figure 6 Map of lake densities across Inglefield Land. The grid cells are 1 km².

Taken together, there are clear areas of overlap between exposed crystalline basement rock, bedrock fractures, and high lake densities north of Hiawatha Glacier and in the valleys draining into Rensselaer and Marshall bays (Fig. 7). West of Hiawatha Glacier and in western Inglefield Land, there are large areas of exposed basement rocks with lakes but no large-scale fractures. There are also a few areas of sedimentary cap rocks with high lake densities, most notably between Rensselaer, Marshall, and Dallas bays.

Figure 7 Combined map of bedrock lithology, fracture zones, and lake densities with cosmogenic nuclide results. The semitransparent arrows show zones of increased erosion as indicated by mapping proxies. Arrow size does not correspond to amount of erosion; rather, it just marks the areas of erosion. ^b Boulder sample.

4.2 Cosmogenic nuclide analysis

Apparent ¹⁰Be exposure ages are generally youngest at our coastal sites and oldest at the inland sites. ¹⁰Be concentrations are lowest at our ∼ 100 m a.s.l. coastal sites, corresponding to apparent exposure ages between 24.9 ± 0.5 and 36.5 ± 1.7 ka (Figs. 7 and 8; Table S1). These sites exhibit evidence of glacial erosion, including sculpted bedrock and striations, suggesting the presence of warm-based ice. At the higher-elevation coastal site on the western transect (255 m a.s.l.), the apparent ¹⁰Be age is slightly older – 62.2 ± 1.7 ka. Our oldest apparent exposure ages come from the inland sites on both transects, dating from 99.1 ± 2.4 to 213.7 ± 3.5 ka on the western transect and from 87.9 ± 3.8 to 215.8 ± 7.4 ka on the eastern transect. We do not observe any relationship between the apparent ¹⁰Be exposure age of our inland sites and their elevation or distance along the transect. Ages from the ice margin sites generally fall between those from the coastal and inland sites. Apparent ¹⁰Be exposure ages from the western ice margin site range from 50.2 ± 0.9 to 133.7 ± 2.6 ka. At the eastern ice margin site, ¹⁰Be ages are between 49.2 ± 0.9 ka and 62.6 ± 1.2 ka. Our single boulder sample at the western ice margin site has an apparent ¹⁰Be exposure age of 50.2 ± 0.9 ka.

Figure 8 Subset of bedrock sample photos and the resultant cosmogenic nuclide results. Note the smoothness of the coastal outcrops compared to the rough, weathered textures of the inland and ice margin sites.

${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios follow a general pattern that is inverse to the apparent ¹⁰Be ages, with the highest ratios (closest to the production ratio) at the coast, the lowest ratios inland, and ratios from the ice margin in the middle. At the coastal sites, the three samples taken from ∼ 100 m a.s.l. have ratios of 6.38 ± 0.44, 7.58 ± 0.42, and 8.41 ± 0.81 and overlap with the constant exposure isochrone with 95 % confidence (Fig. 9). The 255 m a.s.l. coastal site has a ratio of 5.68 ± 0.30, requiring a minimum burial duration of 500 kyr and at least 25–100 kyr of exposure. In contrast, ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios at our inland sites range from 2.19 ± 0.16 to 5.28 ± 0.30, with the majority falling between the 500 and 1500 kyr burial isochrones and between the 100 and 500 kyr exposure isochrones for total minimum exposure and burial durations of between ∼ 600 and ∼ 2000 kyr. There is one apparent outlier below the 2000 kyr burial isochrone, with a ratio of 2.19 ± 0.16 (22GRO-03). At the western ice margin site, ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios in our three samples range between 4.00 ± 0.18 and 5.14 ± 0.21, corresponding to a minimum glacial history with ∼ 100 and ∼ 250 kyr of cumulative exposure and ∼ 250 and ∼ 1250 kyr of cumulative burial. The boulder has a ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratio of 5.15 ± 0.21, which is the highest ratio at the western ice margin site. At the eastern ice margin site, ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios are between 5.11 ± 0.21 and 5.53 ± 0.26 and represent minimum glacial histories with cumulative exposure durations of ∼ 50–125 kyr and cumulative burial durations of ∼ 500–750 kyr. From both transects, ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios from samples at the ice margin (5.14 ± 0.26; mean ± 1σ) are slightly higher than those from the inland sites (4.47 ± 0.86, excluding the apparent outlier of 2.19 ± 0.16; 22GRO-03), though these overlap at 1σ uncertainty.

Figure 9 ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ two-nuclide diagram. Concentrations are normalized to an Arctic high-latitude sea level ¹⁰Be production rate of 3.96 atoms per gram per year using a Greenland-specific ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ surface production ratio of 7.3 (Corbett et al., 2017; Young et al., 2013).

4.3 Constraints on Quaternary ice cover history and erosion

We modeled cosmogenic nuclide accumulation through the Quaternary to determine whether the measured ¹⁰Be and ²⁶Al concentrations across Inglefield Land can be explained by a common exposure history and differential erosion. We tested exposure histories derived from δ¹⁸O thresholds between 3.6 ‰ and 4.0 ‰, corresponding to 0.7–1.9 Myr of cumulative exposure and 0.8–2 Myr of cumulative burial over the last 2.7 Myr (Fig. 10). For each of these exposure histories, modeled nuclide concentrations yielded a good fit to those we measured for at least one sample. However, only the exposure history constructed with a δ¹⁸O threshold value of 3.74 ‰ yields ¹⁰Be and ²⁶Al concentrations with a good fit to the data in all of the bedrock samples, except for those with ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios above the production ratio (22GRO-01 and 22GRD-CR04-SURF), the sample previously identified as an outlier based on its low ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratio (22GRO-03), and our boulder sample (22GRO-32). However, other scenarios with similar δ¹⁸O threshold values of 3.72 ‰ and 3.76 ‰ fit most of the non-outlier bedrock samples, though they are less than our best-fit δ¹⁸O threshold value of 3.74 ‰. This best-fit exposure history corresponds to ∼ 0.9 Myr of total exposure and ∼ 1.8 Myr of total burial over the last 2.7 Myr (Fig. 11). Given this ice cover scenario, we find model–data agreement, with subglacial and subaerial erosion rates ranging between 0 and ∼ 2 × 10⁻³ mm yr⁻¹ and between 0 and ∼ 2 × 10⁻² mm yr⁻¹, respectively (Fig. 12). For our best-fit exposure scenario using a δ¹⁸O threshold value of 3.74, we were not able to simulate erosion rates for bedrock samples that did not have a good fit to this δ¹⁸O threshold value (22GRO-01, 22GRO-03, 22GRD-CR04-SURF, and 22GRO-32). While there were some exposure and erosion combinations that were compatible with measurements from our single boulder sample (22GRO-32), we prioritized bedrock samples when assessing whether all of the sites could have a common Quaternary ice cover history, as the boulder likely has a different exposure and erosion history than the collocated bedrock. We do not observe clear trends in the modeled subaerial erosion rates, though we note that subglacial erosion rates are modeled to be greater than about 1 × 10⁻⁴ mm yr⁻¹ for ice margin areas, while these exhibit larger variability for inland and coastal samples.

Figure 10 Model–data fits for exposure histories constructed with δ¹⁸O thresholds of 3.60 ‰–4.00 ‰. For each δ¹⁸O threshold, we modeled cosmogenic nuclide concentrations for each sample using subaerial and subglacial erosion rates ranging from 0 to 2.5 × 10⁻¹ mm yr⁻¹ mm yr⁻¹. Colored tiles show the best or lowest error-weighted sum of squares (EWSS) for each sample and the δ¹⁸O threshold across all of the tested erosion rate combinations. We consider an EWSS <2.5 to be an acceptable model–data fit. White tiles indicate that no combination of erosion rates yielded an EWSS <2.5. An exposure history constructed with a δ¹⁸O threshold of 3.74 ‰, outlined in the black box, was the only exposure history we tested that gave an acceptable fit for all non-outlier bedrock samples. ^a Outlier samples identified with ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios. ^b Boulder sample.

Figure 11 Exposure history constructed with a 3.74 ‰ δ¹⁸O threshold (teal line) on the LR04 stack (Lisiecki and Raymo, 2005), below which we considered the site ice-free and above which we considered the site ice-covered. The resulting exposure history is shown in the top panel, where periods of exposure are red and periods of burial are blue. This is the only exposure history we tested that, even when considering site-specific subaerial and subglacial erosion rates, yielded an acceptable fit for all non-outlier bedrock samples. Note that Holocene exposure is decoupled from the δ¹⁸O threshold-based exposure duration calculations and is thus not represented in the plot.

Figure 12 Model–data fits for the tested subaerial (top) and subglacial (bottom) erosion rates for our preferred exposure history constructed with a δ¹⁸O threshold of 3.74 ‰. Colored tiles show the best or lowest error-weighted sums of squares for the tested erosion rate for each sample. White tiles indicate that, for this exposure history, there is no combination of subaerial and subglacial erosion that yields a good fit to the data. For example, in combination with subaerial erosion, only one tested subglacial erosion rate (5 × 10⁻⁴ mm yr⁻¹) yielded an acceptable fit for sample 22GRD-CR02-SURF. However, when applying that subglacial erosion rate, subaerial erosion rates ranging from 0.2 to 1.2 × 10⁻⁴ mm yr⁻¹ yielded an acceptable model–data fit. We tested subaerial and subglacial erosion rates from 0 to 2.5 × 10⁻¹ mm yr⁻¹. For no sample was there an acceptable model–data fit for the preferred exposure history when either erosion rate was >2.5 × 10⁻² mm yr⁻¹. In addition, for any sample that yielded a good model–data fit for an erosion rate of 1 × 10⁻⁵ mm yr⁻¹, we found good fits down to 0 mm yr⁻¹ (no erosion of that type). For clarity, the y axis for both panels is restricted to 1 × 10⁻⁵–2.5 × 10⁻² mm yr⁻¹. ^a Outliers identified from ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios. ^b Boulder sample.

5 Discussion

5.1 Differential erosion, ice streams, and ancient landscapes across Inglefield Land

We identified an ancient landscape, with our oldest sites retaining at least a ∼ 1.5 Myr history. However, the nuclide concentrations vary throughout the landscape despite the fact that the long-term ice margin history should be similar across sites, as evidenced by recent retreat and advance (Søndergaard et al., 2020). After the LGM, the ice sheet retreated over Inglefield Land from its maximum extent in Kane Basin. The modern coast deglaciated ∼ 8.5 ka, and ice retreated behind the present margin by 6.7 ka (Søndergaard et al., 2020), corresponding to our entire <50 km transects becoming ice-free within 2 kyr. We speculate that the LGM advance across Inglefield Land was likely similarly swift given that our field area represents a short distance with respect to the entire ice flowline out to the Greenland Ice Sheet LGM terminus in northern Baffin Bay. These findings suggest that – within the error of our ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$-derived exposure scenarios – all of our sites have experienced similar periods of ice cover and ice-free conditions on glacial–interglacial timescales. This thus begs the question: why do amounts of nuclide inheritance and ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios at our coastal, inland, and ice margin sites differ?

Our combined mappings of bedrock lithology, landscape-scale fractures, and lake density as evidence of past glacial erosion, along with field observations of bedrock surface weathering, provide a first hint that erosion is responsible for differing ²⁶Al and ¹⁰Be concentrations. Though some of the cap rock removal likely occurred before the Quaternary (e.g., Krabbendam and Bradwell, 2014), our mapping of fractures and lakes reveals a clear imprint of differential erosion by the ice sheet. Regions of crystalline bedrock and high lake density inland (with the exception of the areas in front of Hiawatha Glacier) and near the ice margin exhibit weathered bedrock surfaces diagnostic of cold-based ice cover during the last glacial cycle. Overlapping areas of crystalline bedrock, fractures, and high lake densities in the four valleys leading into Force, Rensselaer, and Dallas bays and in front of Hiawatha Glacier, however, contain fresh, unweathered bedrock outcrops with their primary surface features intact. The prevalence of erosion indicators in these areas suggests the presence of erosive ice here with higher velocities than experienced by the surrounding landscape for at least a portion of the last glacial cycle. We posit that these areas were ice stream onset zones during the last glacial cycle and likely earlier glacial cycles as well.

The ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios and lower levels of cosmogenic nuclide inheritance support more erosion at the coastal sites versus the inland and ice margin areas. The samples at 100 m a.s.l. in the Rensselaer and Marshall valleys have the youngest apparent exposure ages and ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios close to the production ratio, suggesting removal of the longer burial signal preserved at the upland sites via erosion. Furthermore, the 255 m a.s.l. coastal site in the Rensselaer Valley has more inheritance and lower ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios compared to the 100 m a.s.l. elevation site, which may reflect differential erosion by an ice stream in the Rensselaer Valley as it increased velocity and erosion downflow. Previous studies on Baffin Island and in Scandinavia used mapping and cosmogenic nuclide measurements to identify areas as ice stream onset zones (Andersen et al., 2018; Briner et al., 2006, 2008; Brook et al., 1996). Ice in these onset zones transitioned from cold-based to warm-based ice and did not erode >3 m during the last glacial cycle, though it is thought that velocities increased downstream from these areas. We find similar cosmogenic nuclide evidence of ice stream onset zones at the mouths of Rensselaer and Marshall bays, suggesting that ice streams may have begun in these areas during glacial maxima.

Cosmogenic nuclide concentrations in weathered bedrock surfaces at our inland and ice margin sites retain combined minimum exposure and burial signals of 500 to 1500 kyr, demonstrating the antiquity of the landscape across Inglefield Land. The preservation of these inland and ice margin landscapes implies past cold-based ice conditions and, in turn, low erosion across Inglefield Land through multiple glacial cycles. This also suggests that the erosional landscapes seen in mapping across the interior of Inglefield Land (i.e., the lack of cap rocks and an exposed crystalline basement, high lake densities, and bedrock fractures) were created prior to the last glaciation.

Such old, low-erosion landscapes have been preserved in other glaciated parts of the Arctic. Our results are similar to studies from southern, western, northwestern, and northeastern Greenland and Baffin Island that identify long-preserved (often >1000 kyr), ancient landscapes with low amounts of erosion and high amounts of cosmogenic nuclide inheritance (Andersen et al., 2020; Beel et al., 2016; Bierman et al., 1999; Briner et al., 2006; Ceperley et al., 2020; Corbett et al., 2016a; Miller et al., 2006; Roberts et al., 2013; Sbarra et al., 2022; Skov et al., 2020; Søndergaard et al., 2020). Similarly, the pattern of significantly lower cosmogenic nuclide inventories in lower, coastal areas in less preserved landscapes is seen elsewhere in Greenland and the Canadian Arctic. In eastern Greenland, for example, Skov et al. (2020) found higher degrees of cosmogenic nuclide inheritance in weathered bedrock uplands and lower concentrations in sculpted bedrock closer to the coast around Dove Bugt, attributing this to a polythermal ice sheet with efficient erosion at lower elevations. Studies on Baffin Island using ¹⁰Be and ²⁶Al measurements found transitions from cold-based and non-erosive ice to warm-based and erosive ice at lower elevations near the coast (e.g., Briner et al., 2006, 2008), supporting our conclusions of transitions from cold-based to warm-based ice near the coast of Inglefield Land.

5.2 Can measured ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratio differences be explained by a common ice cover history and variable erosion rates?

We tested our hypothesis that the disparate nuclide concentrations across Inglefield Land could have resulted from a shared Quaternary exposure history and differential erosion. Indeed, we found that the ¹⁰Be and ²⁶Al concentrations in most of our bedrock samples can be explained by a Pleistocene ice cover history with 0.9 Myr cumulative exposure and 1.8 Myr cumulative burial, together with sample-varying subglacial and subaerial erosion rates.

The range of subaerial erosion rates (0–2 × 10⁻² mm yr⁻¹) that provide good fits to the measured nuclide concentrations for our preferred exposure history is broadly consistent with subaerial erosion rates measured in other cold, arid, and polar environments. There are no clear patterns of spatial variability across the landscape, though estimates from 22GRO-48 at the coast exhibit the largest range of possible subaerial erosion rates. A global synthesis found that the mean subaerial erosion rates at the rock outcrops are 13 × 10⁻³, 8 × 10⁻³, and 4 × 10⁻³ mm yr⁻¹ in cold, arid, and polar environments, respectively (Portenga and Bierman, 2011). A cosmogenic nuclide-based study from the southern Ellsworth Mountains in Antarctica found longer-term subaerial erosion rates of 5.52 ± 0.26 × 10⁻³ mm yr⁻¹ for gneiss (Marrero et al., 2018), similar to those derived from Baffin Island gneiss of 2 × 10⁻³ mm yr⁻¹ (Margreth et al., 2016). Furthermore, subaerial erosion of gneiss in the polar regions is thought to occur primarily through mineral- and granular-scale processes, including abrasion and disintegration, and this is accelerated by wind-blown sediments (Margreth et al., 2016). We noted abundant gruss around our inland and ice margin sites and near-constant winds with wind-transported sediments, suggesting that these processes may be responsible for subaerial erosion in Inglefield Land. Additionally, the presence of weathering features, including gruss and weathering pits across Inglefield Land, indicates that subaerial erosion likely varies not only from outcrop to outcrop as a result of varying lithology, but also across individual bedrock outcrops (e.g., sampled from near an old weathering pit).

Subglacial erosion rates that yield good model–data fits are between 0 and ∼ 2 × 10⁻³ mm yr⁻¹ across our sites, and all of the samples have possible subglacial erosion rates >∼ 1 × 10⁻⁴ mm yr⁻¹. Samples at the ice margin do not have modeled erosion rates below this value, and their maximum modeled erosion rates (∼ 1 × 10⁻³ mm yr⁻¹) are higher modeled maximum erosion rates at the inland sites (∼ 5 × 10⁻⁴ mm yr⁻¹). This indicates that the bedrock surfaces at the modern ice margin are consistently eroded during ice burial and may experience higher subglacial erosion rates. Some inland sites have subglacial erosion rates of 0 mm yr⁻¹, thus highlighting the potential for spatial variability of subglacial erosion rates under a cold-based ice sheet. Furthermore, subglacial erosion rate estimates from our sites with the lowest levels of cosmogenic nuclide inheritance (at the coast) exhibit the largest spread in modeled subglacial erosion rates. These coastal samples also have wide ranges of subaerial erosion rates, and while lower nuclide concentrations indicate that either subglacial or subaerial erosion must be high to limit nuclide accumulation, we cannot determine which kind of erosion is more significant.

Overall, subglacial erosion rate estimates across Inglefield Land are generally lower than those measured elsewhere in Greenland. Many of these previous estimates, however, come from warm-based glaciers, while our results suggest that Inglefield Land was covered primarily by cold-based ice. The upper range of our modeled subglacial erosion rates agrees with that from eastern Greenland derived from sediment flux data of 1–4 × 10⁻² mm yr⁻¹ (Andrews et al., 1994), though this was later updated by Cowton et al. (2012) to 3 × 10⁻¹ mm yr⁻¹ to account for sediment entrained in icebergs (Syvitski et al., 1996). In central western Greenland, suspended sediment measurements from the Watson River yielded average subglacial erosion rates of 5 × 10⁻¹ mm yr⁻¹, with annual erosion rates as high as 4.5 and 4.8 ± 2.6 mm yr⁻¹, all higher than our modeled rates (Cowton et al., 2012; Hasholt et al., 2018; Hogan et al., 2020). Glaciomarine deposits near Petermann Glacier in northwestern Greenland and at the mouth of Jakobshavn Isbræ constrain erosion rates during the last deglaciation to 0.29–0.34 and 0.52 mm yr⁻¹ (Hogan et al., 2012; Hogan et al., 2020). Finally, Balter-Kennedy et al. (2021) determined centennial-scale subglacial erosion rates of 3–8 × 10⁻¹ mm yr⁻¹ and orbital-scale erosion rates of 1–3 × 10⁻¹ mm yr⁻¹ near Jakobshavn Isbræ. These previously derived subglacial erosion rates may be higher than what we report because they come from places with different ice sheet thermal states. Thus, these estimates of subglacial erosion could reflect periods when temperatures were generally warmer or periods of warming, which can lead to higher subglacial erosion rates because of increased basal temperatures and subglacial sliding, perhaps in part due to increased meltwater delivery to the bed in warmer climates (Alley et al., 2019). Additionally, these studies of glacial erosion focus on relatively short timescales (< multi-millennial), while our modeled erosion rates are integrated over the Quaternary, perhaps highlighting the effects of temporal integration as previously noted by Koppes and Montgomery (2009) and Spotila et al. (2004).

Cold-bedded glaciers are commonly considered non-erosive, yet glacial erosional features (striae, scrapes, grooves, and isolated blocks) in landscapes otherwise protected by cold-based ice demonstrate their erosional capacity (Atkins et al., 2002; Cuffey et al., 2000; Sugden, 1978; Sugden et al., 2005; Ugelvig and Egholm, 2018). Our modeled subglacial erosion rates between 0 and ∼ 2 × 10⁻² mm yr⁻¹ in mapped areas of cold-based ice cover point to the erosive nature of cold-based glaciers and are consistent with previous estimates of erosion. Syntheses found modern subglacial erosion rates under frozen-bedded glaciers in the Antarctic Peninsula to be between 1 × 10⁻² and 1 × 10⁻¹ mm yr⁻¹ (Koppes et al., 2015), while on the Meserve Glacier, a cold-based alpine glacier in Victoria Land, Antarctica, subglacial erosion rates have been estimated to be 2 × 10⁻³ mm yr⁻¹ (Cuffey et al., 2000). Studies of tors in the Cairngorms, Scotland, yielded subglacial erosion estimates of ∼ 4.4 × 10⁻³ mm yr⁻¹ during cover by the cold-based Celtic Ice Sheet (Phillips et al., 2006). Various mechanisms have been proposed for erosion by cold-bedded glaciers. Frozen-bedded glaciers entrain debris at their bed and can minorly abrade and pluck bedrock surfaces (Atkins et al., 2002; Cuffey et al., 2000; Sugden, 1978; Sugden et al., 2005; Ugelvig and Egholm, 2018). Extensive studies of the Cairngorms in Scotland using cosmogenic nuclides and geomorphic models of tor formation and erosion found that cold-based ice erosion was potentially capable of significantly modifying pre-existing landforms (Goodfellow et al., 2014; Hall and Glasser, 2003; Hall and Phillips, 2006; Phillips et al., 2006).

In summary, we find that the disparate nuclide concentrations across Inglefield Land can be explained by a common Quaternary exposure–burial history if erosion rates are allowed to vary sample by sample. Erosion rates consistent with the cosmogenic data are similar to previously published estimates in polar areas covered by cold-based glaciers. Variability in erosion rates across the landscape likely reflects differences in lithology as well as subglacial conditions. By allowing our erosion rates to vary at each sampling location, in our model, we are able to account for potential erosional variability resulting from local differences in lithology (e.g., mineralogy and crystal size), landscape position, and sampling location (e.g., within an old weathering pit or a zone of high grussification). However, our model does not account for temporal variability in subglacial erosion. Temporal variability (abrasion versus quarrying) should be diminished when averaged across many glacial cycles, but the imprint of variable glacial erosion during the last glacial cycle (e.g., a site of a cobble-sized block of bedrock removed beneath mostly frozen ice; Atkins et al., 2002; Hall and Phillips, 2006) may lead to differences in the resulting subglacial erosion rates. Furthermore, the basal zone across Inglefield Land may have transitioned from less erosive to more erosive and back to less erosive during switches between warm- and cold-bedded conditions through the thickening and thinning of the ice sheet during a glacial cycle. This could lead to time-varying subglacial erosion rates and shorter periods of increased glacial erosion, as evidenced by the glacial sculpting at our coastal sites.

5.3 Long-term ice sheet fluctuations across Inglefield Land

${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios from our oldest bedrock surfaces suggest that, at a minimum, the Greenland Ice Sheet covered Inglefield Land for 1200 kyr and was smaller than today for 400 kyr over the last 1600 kyr. We determined that a Quaternary exposure history with ∼ 0.9 Myr of cumulative exposure and ∼ 1.8 Myr of cumulative burial (constructed using our best-fit δ¹⁸O threshold of 3.74 ‰) is consistent with the measured cosmogenic nuclide concentrations in nearly all of our bedrock surfaces when erosion rates are allowed to vary for each sample. In this preferred exposure scenario, much of the exposure takes place between 2.7 and 1.2 Ma, before the mid-Pleistocene transition, but this still requires ice-free conditions (and exposure) within the last ∼ 1.2 Myr during major interglacials. These periods of exposure are of particular interest as they indicate times when the Greenland Ice Sheet was at its present extent or smaller.

Our results complement other limited terrestrial studies of Greenland Ice Sheet stability throughout the Quaternary. ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ measurements from the GISP2 bedrock core from under the center of the Greenland Ice Sheet revealed that the ice sheet was present for most of the Quaternary but was nearly completely absent at least once in the last 1.1 Myr (Schaefer et al., 2016). Studies of basal material from Camp Century show similar results, indicating that the northwestern Greenland Ice Sheet was present throughout most of the Pleistocene (Christ et al., 2021). Additionally, when compared to studies of material from under the modern Greenland Ice Sheet, our results suggest that Inglefield Land was ice-free throughout much more of the Quaternary than the interior sectors (i.e., currently covered by the modern ice sheet), which is logical given that Inglefield Land would only be completely covered during glacial maxima. When taken at face value, our preferred exposure history corresponding to a δ¹⁸O threshold of 3.74 ‰ indicates that the northwestern Greenland Ice Sheet persisted at an extent larger than today throughout some Pleistocene interglacials, which is further supported by similar findings from the Laurentide Ice Sheet (LeBlanc et al., 2023).

6 Conclusions

The conclusions of this study are as follows:

• Geologic mapping partitions Inglefield Land into zones of glacial erosion and protection, revealing restricted zones of erosion in the lower Force, Rensselaer, and Marshall valleys and in front of Hiawatha Glacier as ice velocity increased coastward into ice streams; elsewhere, landscapes were likely covered by low-erosion (e.g., frozen-bedded) regimes.

• Patterns of glacial erosion are confirmed by cosmogenic nuclide measurements, revealing an ancient landscape across Inglefield Land. This implies widespread cold-based ice cover across much of Inglefield Land during Quaternary glacial cycles, even during transitions between interglacial and glacial periods. These measurements also indicate that, despite cosmogenic nuclide inheritance, there was perhaps temporary warm-based ice at coastal sites during the last glacial cycle, when the Greenland Ice Sheet extended well beyond Inglefield Land and ice streams originated near the mouths of the modern Force, Rensselaer, and Marshall bays. Mapping and cosmogenic nuclide measurements allow us to differentiate between areas of cold- and warm-based ice, offering a clear look at the polythermal nature of former ice sheets.

• Differential subglacial and subaerial erosion explains the disparate cosmogenic nuclide concentrations found at our inland and ice margin sites. Modeled subaerial erosion rates match those found in other polar regions for similar lithologies. Our modeled subglacial erosion rates are lower than those previously calculated for Greenland, as other studies focused on warm-based parts of the Greenland Ice Sheet. We provide estimates of subglacial erosion under cold-bedded conditions.

• Ice cover durations derived from ${}^{\mathrm{26}}\mathrm{Al}{/}^{\mathrm{10}}\mathrm{Be}$ ratios and cosmogenic nuclide modeling reveal a common ice sheet history, and variable erosion rates indicate that the ice sheet persistently covered Inglefield Land for approximately one-third of the last 2.7 Myr. This exposure occurred mostly early in the Quaternary before the mid-Pleistocene transition. Like other findings (e.g., LeBlanc et al., 2023), when taken at face value, our δ¹⁸O threshold value when applied to the global δ¹⁸O stack implies ice cover during some middle and Late Pleistocene interglacials. Our results also suggest that the Greenland Ice Sheet was at its current extent or smaller during several interglacial periods, adding to the body of literature that indicates a dynamic Greenland Ice Sheet throughout the Quaternary.

• Bedrock surfaces across much of Inglefield Land were not greatly eroded during the LGM and contain a long-term cosmogenic nuclide memory of ice sheet fluctuations through the mid-Pleistocene – at a minimum. Our bedrock sampling strategy and measurements provide an analog for information contained in sub-ice material under the Greenland Ice Sheet, indicating that this material should contain valuable cosmogenic nuclide archives of ice sheet change during minimum phases.