Terrestrial cosmogenic nuclide bedrock depth profiles used to infer changes in Holocene glacier cover, Vintage Peak, southern Coast Mountains, British Columbia

The majority of glaciers in North America reached their maximum Holocene downvalley positions during the Little Ice Age (1300–1850 CE), and in most cases, this expansion also destroyed earlier evidence of glacier activity. Substantial retreat in the 20th and early 21st centuries exposed bedrock that fronts many glaciers that may record Early to mid-Holocene exposure and later burial by ice, which can be elucidated using multiple-nuclide cosmogenic surface exposure dating. Cores of bedrock allow the measurement of cosmogenic nuclide depth profiles to better constrain potential exposure and burial histories. We collected four bedrock surface samples for ¹⁰Be and ¹⁴C surface exposure dating and shallow (<0.6 m depth) bedrock cores from Vintage Peak, in the southern Coast Mountains of British Columbia, Canada. We apply a Monte Carlo approach to generate combinations of exposure and burial duration that can explain our data. We found that Vintage Peak became uncovered by the Cordilleran Ice Sheet between 14.5 and 9.7 ka, though upper elevations on Vintage Peak retained ice until 10–12 ka before retreating to smaller than modern positions. Glaciers on Vintage Peak advanced within 100 m of Late Holocene maximum positions around 4–6 ka. Poorly constrained subglacial erosion rates, possible inheritance, and variable mass shielding complicate our ability to more robustly interpret bedrock cosmogenic surface exposure histories. Nine ¹⁰Be ages on Late Holocene moraines reveal that glaciers reached their greatest Holocene extents by ca. 1300 CE. Our results agree with other regional glacier records and demonstrate the utility of surface exposure dating applied to deglaciated bedrock as a technique to help construct a record of Holocene glacier activity where organic material associated with glacier expansion may be absent or poorly preserved. Further work to increase exposure and/or burial history modeling complexity may help to better constrain complex exposure histories in glaciated alpine areas.

1 Introduction

Moraines provide the most unequivocal evidence to gauge the magnitude of a glacier advance; however, Holocene moraines deposited prior to the Little Ice Age (LIA, 1300–1850 CE) in North America are rarely preserved (Barclay et al., 2009; Luckman, 2000; Menounos et al., 2009; Solomina et al., 2015). Glacier activity, particularly the timing of glacier advances, prior to the LIA can be elucidated through lateral moraine stratigraphy, overridden wood in glacier forefields, and proglacial lake sediments, among other proxy records (Hawkins et al., 2021; Luckman et al., 2020; Tomkins et al., 2008). Alternative methods are needed to determine the duration of ice cover at alpine sites to place current rates of glacier retreat into a long-term perspective. Cosmogenic nuclide exposure analysis offers a useful method to document past glacier activity (Balco, 2011).

Cosmogenic nuclide exposure analysis measures the concentration of rare isotopes in rock at Earth's surface to calculate the duration the surface has been exposed and/or buried (Balco, 2011; Gosse and Phillips, 2001). For Late Pleistocene and Holocene studies, the combination of ¹⁴C ($t_{\mathrm{1}/\mathrm{2}}=\mathrm{5730}\text{years}$) and ¹⁰Be ($t_{\mathrm{1}/\mathrm{2}}=\mathrm{1.38}\mathrm{Myr}$) provides several advantages over the measurement of a single nuclide alone. Due to the relatively short half-life of ¹⁴C, times of previous exposure, such as during Pleistocene interglacials, are quickly “forgotten” through radioactive decay. This short half-life has been used in many studies to elucidate glacier extent change and/or subglacial erosion rates (e.g., Goehring et al., 2011, 2013; Graham et al., 2019; Schimmelpfennig et al., 2022; Schweinsberg et al., 2018). Hippe (2017) demonstrated that 10 kyr of burial (such as under a thick ice sheet) accompanied by 2.4 m of subglacial erosion would remove any ¹⁴C inventory to below detection limits. Conversely, ¹⁰Be decays slowly enough that it is effectively stable over Holocene timescales. In situations of negligible erosion, ¹⁰Be acts to record total exposure duration at a site, while ¹⁴C can reveal periods of burial. Periods of Holocene burial will cause depth profiles of ¹⁴C to steepen and shift towards lower concentrations as ¹⁴C decays. Hippe (2017) provides a clear overview of the concepts behind using paired ¹⁴C and ¹⁰Be to investigate Holocene surficial processes.

Cosmogenic nuclide production is rapidly attenuated with depth in bedrock and is the combination of spallogenic and muogenic reactions (Gosse and Phillips, 2001). At the surface with no overlying mass, spallation reactions account for nearly 99 % of ¹⁰Be and 77 % of ¹⁴C production, with the remaining nuclide production coming from fast muon interactions and negative muon capture (Dunai, 2010). Assuming a rock density of 2.65 g cm⁻³, 1 m below the bedrock surface, the nuclide production by spallation accounts for ∼97 % of ¹⁰Be and just ∼47 % of ¹⁴C. The ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ ratio should increase with depth, as muons are less affected by mass shielding and contribute to a larger proportion of ¹⁴C production than ¹⁰Be with depth. The ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ ratio at the surface is simply set by the nuclide's relative production rate and the duration of exposure. Erosion of the bedrock surface will serve to lower the concentration of both nuclides and elevate the ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ ratio. A nuclide depth profile in bedrock would shift to the left (less total nuclide) in an erosive scenario compared to one of constant exposure and no erosion (Fig. 1). Nuclide depth profiles in bedrock provide additional nuclide measurements that constrain possible durations of exposure, burial, and erosion more than a single paired nuclide measurement at the surface.

Figure 1 Theoretical ¹⁴C concentrations after 10 kyr of continuous exposure with depth below bedrock (2.65 g cm⁻³). Panel (a) shows the expected nuclide concentration profile under three scenarios: (1) continuous exposure with no mass shielding (snow or till cover), (2) 0.002 cm yr⁻¹ subaerial erosion, and (3) no erosion but with 220 cm of 0.38 g cm⁻³ snow cover for half of the year. Note that both mass shielding and erosion reduce the total nuclide concentration. Erosion also steepens the nuclide concentration curve. Panel (b) shows the relative contribution of spallogenic vs. muogenic ¹⁴C production with depth. Muons are a significant contributor to total ¹⁴C production but only contribute 1 %–3 % of total ¹⁰Be production at mass depths less than 250 g cm⁻².

At locations where a well-constrained glacier chronology exists, the cosmogenic nuclide inventory within previously glaciated bedrock surfaces can constrain subglacial erosion rates (Balter-Kennedy et al., 2021; Goehring et al., 2013). Subglacial erosion rates, however, can vary significantly across previously glaciated surfaces (Graham et al., 2023; Koppes, 2022; Magrani et al., 2022). Past studies of subglacial erosion rates have sampled multiple sites across latitudinal and longitudinal transects of a glacier's flowline to assess the variability in subglacial erosion across zones a differing slope, ice thickness, and resulting ice velocity (Graham et al., 2023; Rand and Goehring, 2019; Wirsig et al., 2016). The rate of subglacial erosion must be independently estimated in cases where there are few bedrock samples across the formerly glaciated landscape and past glacier activity is unknown or poorly constrained.

In formerly glaciated alpine sites of western Canada, bedrock surfaces that lie beyond limits of Late Holocene glacier expansion would have been exposed to cosmogenic nuclide flux since the retreat of the Cordilleran Ice Sheet (CIS) following the termination of the Last Glacial Maximum (LGM; Menounos et al., 2009, 2017; Seguinot et al., 2016). The CIS, large decaying masses of the CIS, or even cirque glaciers may have readvanced during the Younger Dryas cold interval, after which the ice sheet retreated rapidly (Menounos et al., 2017). Alpine glaciers in western Canada were likely smaller than their 2000–2010 CE extents for much of the Early Holocene (Koch et al., 2007, 2014; Menounos et al., 2004) and may have begun periodically advancing as early as 8 ka, with more evidence for glacier advances of generally increasing magnitude by ca. 6 ka. At most sites in western Canada, glaciers reached their greatest Holocene extent during the Little Ice Age, between 700 and 150 years before 2000 CE (Menounos et al., 2009; Mood and Smith, 2015a). What remains unclear is the size of alpine glaciers prior to neoglacial expansion in the mid-Holocene and the amount of time glaciers were near their Late Holocene maximum extents or retreated to within their modern (ca. 2024 CE) extents.

We seek to investigate the record of Late Pleistocene through Holocene glacier change at an alpine site in the southern Coast Mountains of British Columbia, Canada, through several different applications of cosmogenic nuclide surface exposure dating. To determine the duration of ice cover at this site, we make necessary assumptions of the impact of snow cover, subglacial and subaerial erosion, and inheritance on the calculated exposure and burial duration at each bedrock site. This study highlights the utility of paired nuclide dating and bedrock cores to more closely constrain complex exposure histories and the limitations of bedrock exposure dating; it also provides new numerical ages for end moraines in the southern Coast Mountains.

2 Study area

Vintage Peak (1876 m above sea level – m a.s.l.) is located at the head of Powell Lake on the traditional territory of the Tla'amin Nation in the southern Coast Mountains of British Columbia, Canada (50.25° N, 124.30° W; Fig. 2). The mountain is composed of a mid-Cretaceous granodiorite pluton that forms broad sloping surfaces and steep rock walls that extend down to the valley bottom (Bellefontaine et al., 1994). A small (<0.05 km²) north-facing glacierette herein referred to as “Lockie's Glacier” is located at 1780 m a.s.l. in a cirque west of Vintage Peak, below a sub-peak called Lockie's Table. Another north-facing glacierette resides on the north side of Vintage Peak. The cirque headwalls are variably steep, though large portions of the headwall give way to rolling, convex ridgelines (Fig. 2).

Figure 2 Location of bedrock and boulder sampling locations on Vintage Peak. ¹⁰Be surface exposure ages are listed in text boxes with analytical errors. The green star in the inset map is the study location. The base hillshade is from a 2017 aerial lidar survey. Contour interval is 100 m.

The study site is 50 km east of the Strait of Georgia, with the coastal climate producing high precipitation and moderate temperatures. Between 1998 and 2022, the mean annual air temperature and mean annual precipitation at a BC Ministry of Forests weather station at Toba Camp (52 m a.s.l.), 38 km northeast of Vintage Peak, were 10.1 °C and 1401 mm, respectively (in 2022). Assuming a lapse rate of 6.5 °C km⁻¹, the average annual air temperature at 1600 m on Vintage Peak is close to 0 °C. There are two province of British Columbia snow survey sites with near-annual observations since 1938 CE located just 2 km away from Lockie's Glacier at 882–910 and 1002–1040 m a.s.l. (Fig. 2) that together have a mean annual snow water equivalent (SWE) at the end of March or early April of 847 ± 399 mm water equivalent.

3 Methods

3.1 Field methods

To constrain the age of the inferred Late Holocene moraines in the study area, we sampled nine large boulders on or near the crest of moraines or those resting directly on exposed bedrock. We sampled six boulders from the Lockie's Glacier moraine and three boulders from the Vintage Peak moraine (Fig. 2). In addition, we sampled one erratic (17-VP-07), approximately 60 m beyond the Lockie's Glacier end moraine. At each site, we recorded boulder dimensions, GPS position and elevation, local horizon, and boulder self-shielding using a Brunton compass and inclinometer. We used a gas-powered concrete saw and hammer and chisel to collect ∼1 kg of rock from each boulder.

A common practice in cosmogenic surface exposure dating on moraines is to sample large boulders from on or near the crest of the moraine, with the aim to reduce the chance of sampling a boulder that experienced significant snow cover, exhumation, or post-depositional movement (Gosse and Phillips, 2001; Heyman et al., 2016). The Hakai Airborne Coastal Observatory completed an aerial lidar and orthoimagery survey of Vintage Peak on 15 September 2017 and again on 16 May 2023 (Fig. 6). Snow is retained in May on Vintage Peak, while the September survey was largely snow-free. The elevation difference between the co-registered, 1 m spatial resolution digital elevation models (DEMs) gives an approximation of the snow depth in May 2023 across the survey area, notably around the sampled moraines and erratic. Further background on the lidar acquisition and data processing can be found in Donahue et al. (2023).

We collected four bedrock core samples within or just beyond the Late Holocene extent of Lockie's Glacier (Fig. 3). Core sample locations were chosen to capture a variety of potential bedrock exposure histories. At each location, we used a Shaw Tool gasoline-powered portable backpack drill to collect a pair of 0.35–0.55 m deep, 41 mm outer-diameter bedrock cores. Cores were not taken deeper than 0.55 m due to the slow penetration rate of the drill (supplied drill bits for the rental drill were for sedimentary rather than intrusive igneous rock). A pair of cores within 0.10–0.15 m of each other were taken at each site to have enough material at each depth for nuclide measurements. Similar to boulder samples, we recorded sample coordinates, elevation, surface orientation, and horizon shielding at each bedrock site. Samples were transported in plastic tubes, with core orientation labeled and maintained throughout transport and laboratory sampling (Fig. 3). Below, we refer to the pairs of cores at each site as a single core.

Figure 3 Bedrock core from site 17-VPC4. Photo scale converted to digital representation.

One core (17-VPC4) was collected from a bedrock knob, approximately 150 m beyond the presumed Late Holocene extent of Lockie's Glacier. Cores 17-VPC3 and 17-VPC2 are within Late Holocene glacier limits, with 17-VPC3 ∼50 m inside the Late Holocene moraine and 17-VPC2 ∼40 m beyond the modern (2017 CE) ice limit (Fig. 2). Core 17-VPC1 is located just off a broad ridgeline on the northeast side of the Lockie's Table cirque and is presumed to have been covered by minimally erosive ice when Lockie's Glacier was near its maximum Holocene extent. Our interpretation that the site of core 17-VPC1 experienced negligible erosion is due to a lack of headwall at this point in the ridge, which would have allowed accumulating ice to drape over the ridgeline, forming a minimally erosive flow divide near the site of the core. At all bedrock core sites, we attempted to sample bedrock with minimal glacier erosion by sampling from areas with oxidized surfaces and preserved striations. We avoided surfaces on the lee of bedrock steps that may be the result of subglacial plucking (Graham et al., 2023; Rand and Goehring, 2019).

3.2 Laboratory processing

All samples were processed at the Tulane University Cosmogenic Nuclide Laboratory. We followed the procedure of Nichols and Goehring (2019) to physically and chemically isolate quartz from each sample. Extraction of ¹⁰Be from the sampled boulders and the top 5 cm of each bedrock core followed standard chemical isolation procedures (Ditchburn and Whitehead, 1994). Each batch of approximately eight samples included a sample blank (Table S1 in the Supplement). Each bedrock core was divided into five, 0.05 m long sections. While both ¹⁰Be and ¹⁴C were measured at the top 0.05 m of each core, we only measured ¹⁴C along the remaining length of each core. We measured both ¹⁰Be and ¹⁴C on the one erratic boulder, 17-VP-07. Extraction of ¹⁴C followed the procedures of Goehring et al. (2019).

Concentrations of ¹⁰Be were measured by accelerator mass spectrometry (AMS) at Purdue University's PRIME Laboratory, with all ¹⁰Be measurements normalized to the 07KNSTD standard (Nishiizumi et al., 2007). We measured ¹⁴C concentrations via AMS at the Woods Hole Oceanographic Institution's National Ocean Sciences Accelerator Mass Spectrometry (NOSAMS) facility. Sample blank corrections for the ¹⁴C samples followed the procedure of Balco et al. (2023), where blank corrections are treated as a lognormal distribution. Exposure ages were calculated using version 3 of the online exposure age calculator formerly known as CRONUS-Earth and LSDn (Lifton–Sato–Dunai) scaling (Balco et al., 2008; Borchers et al., 2016; Lifton et al., 2014). Our production rate calibration dataset is based on the measurement of saturated CRONUS-A data presented in Goehring et al. (2019). Topographic and self-shielding corrections were applied to all moraine and bedrock sample surfaces. Corrections for snow cover were applied to all bedrock and moraine boulder samples separately and are discussed further below. Moraine ages are presented as a median ± interquartile range. This choice of summary statistics does not rely on the assumption of a normal distribution in our dataset. Median ages and individual sample ages are presented in years (a) before the year of sampling (2017 CE).

3.3 Bedrock exposure history modeling

The measured nuclide concentration in each of the bedrock cores represents an integrated signal of exposure, burial by ice and snow, and erosion. While this integrated signal precludes us from being able to constrain the precise timing of Holocene glacier advance and retreat, we can apply reasonable estimates of the subglacial erosion rate and seasonal snow cover at each site and then determine the most likely cumulative duration of burial and exposure for each site that can be explained by our data.

We assume that Cores 1 and 4 experienced negligible erosion since they were exposed following the termination of the Last Glacial Maximum. Core 1 lies on a broad bedrock shoulder and would have been covered by thin, minimally erosive ice when Lockie's Glacier was at or near its Holocene maximum extent. Core 4 lies outside of the Lockie's Glacier's Late Holocene moraines and would have experienced no erosion by ice over the course of the Holocene. Subaerial erosion rates in alpine environments are typically low (Elkadi et al., 2022; Lehmann et al., 2020), rarely more than 0.001 cm yr⁻¹. As such, we assume that during periods of exposure, there was negligible subaerial erosion of bedrock surfaces.

Cores 2 and 3 would have been covered by Holocene ice and are expected to have experienced subglacial erosion during periods of burial. Subglacial erosion rates have been shown to be highly variable across regions and within single glacier forefields (i.e., Koppes et al., 2015; Magrani et al., 2022; Rand and Goehring, 2019). While others have imposed a “known” glacier history to model likely subglacial erosion rates (Schimmelpfennig et al., 2022; Steinemann et al., 2020, 2021; Wirsig et al., 2017), we chose to fix a likely subglacial erosion rate in order to estimate Holocene glacier activity. Given our sampling strategy that avoided evidence of Late Holocene plucking, subglacial erosion would have primarily been through abrasion. Abrasion rates have been estimated by others to be 0.009–0.035 cm yr⁻¹ in the Puget Lowland of Washington State under the Cordilleran Ice Sheet (Briner and Swanson, 1998), 0.002–0.033 cm yr⁻¹ at the Rhone Glacier in Switzerland (Goehring et al., 2011), and $\mathrm{0.013}\pm \mathrm{0.008}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}$ at Jakobshavn Isbrae in Greenland (Graham et al., 2023). As striations on the bedrock at Cores 2 and 3 clearly indicate there was a non-zero amount of erosion of the bedrock surface, we present likely cumulative exposure histories for these sites given a subglacial erosion rate of 0.005, 0.01, and 0.02 cm yr⁻¹.

We used the exposure age of an erratic boulder (17-VP-07) outside of Late Holocene glacier limits to constrain the maximum possible exposure duration for all sites, assuming the exposure age of the erratic represents the timing of Late Pleistocene deglaciation on Vintage Peak. The minimum exposure duration at each sample location was constrained by the ¹⁴C exposure age of bedrock surface samples. We plotted the measured ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ ratios on a two-isotope diagram with burial isochrons to determine the apparent burial history of each sample (Fig. 7). We note that the apparent burial history is dependent on model parameters; any amount of erosion will result in longer apparent periods of burial and exposure. We later compared the apparent burial history to the model exposure and burial histories discussed below.

All sites are assumed to have been seasonally covered by snow throughout their exposure history. We used the average snow depth and density from nearby snow monitoring sites to apply a snow cover correction for 6 months of the year. The average (±1 standard deviation) snow density is $\mathrm{0.38}\pm \mathrm{0.07}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$ and average depth is 220 ± 102 cm (Table S2). During periods of burial by ice, we assume the ice was sufficiently thick to cease nuclide production.

For each site, we run a Monte Carlo simulation 10 000 times per discrete erosion rate (ϵ) of 0.0, 0.005, 0.01, or 0.02 cm yr⁻¹, with each run picking a uniform, randomly distributed value of exposure duration (t_e) and burial duration (t_b) in years. This simulation results in 40 000 runs per site. Our method is similar to the method presented in Jones et al. (2023), except our model looks at a cumulative burial and exposure history rather than developing the exposure history in a stepwise fashion over the Holocene. The model starts by calculating the nuclide production rate at each depth below the surface, beginning at the air–snow interface through to 100 cm depth within the bedrock. The production rate of each nuclide is calculated from the sum of spallation, fast muon, and negative muon contributions (Balco et al., 2008). We chose to make no assumptions about the likelihood of exposure or burial over the course of the Holocene, simply limiting to maximum and minimum possible exposure and burial durations at each site. The exposure duration is a value between the maximum possible exposure duration of 15 kyr and the minimum exposure duration, which is the ¹⁴C age of the bedrock surface for each site. The burial duration is a value between zero and the difference between the maximum exposure duration and the minimum exposure duration at each site. The theoretical nuclide concentration at each depth in the model space was calculated using Eq. (1):

$$\begin{array}{}\text{(1)}& N_{i,z}={\displaystyle \frac{P_{i,z}}{{\mathit{\lambda }}_i}}(\mathrm{1}-e^{(-t_{\mathrm{e}}{\mathit{\lambda }}_i)})\left(e^{-{\mathit{\lambda }}_i{t}_{\mathrm{b}}}\right)\left(e^{\frac{-\mathit{\rho }\mathit{ϵ}{t}_{\mathrm{b}}}{\mathrm{\Lambda }}}\right),\end{array}$$

where N_i,z is the nuclide concentration (atoms g⁻¹) of the isotope in question at a specific mass depth (¹⁴C or ¹⁰Be), P_i,z is the depth-specific total production rate ($\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$), λ_i is the nuclide-specific decay constant (s⁻¹), ρ is the rock density (assumed to be 2.65 g cm⁻³), and Λ is the spallation attenuation length (160 g cm⁻²). While the production rate calculation does account for spallation and muon-produced nuclides, we acknowledge that this equation does not account for muon-produced nuclide concentrations during periods of burial and subglacial erosion. However, this effect is small relative to the total expected nuclide concentration.

We ran the model setup described above assuming no inherited ¹⁴C or ¹⁰Be. Since inherited ¹⁰Be is more likely, we ran our model again for each site assuming a duration of possible inheritance as indicated by the excess ¹⁰Be, which can be explained by ∼3 kyr of inheritance. This produced another 40 000 total simulations per site.

We tested the influence of a 0.001 cm yr⁻¹ subaerial erosion rate on our simulations and found that for Cores 1–3 there was no significant difference in the modeled burial or exposure ages, and at Core 4, inclusion of subaerial erosion yielded no simulations with exposure and burial histories that matched our observations. Thus, we present our results with no subaerial erosion.

We retained simulations that produced modeled ¹⁴C profiles within the 2σ error bounds of the measured ¹⁴C profiles and a modeled surface ¹⁰Be concentration within the 2σ errors in the measured ¹⁰Be concentration. We then calculated the χ² statistic to evaluate goodness of fit between the modeled and measured nuclide concentrations with depth. The solution with the lowest χ² value was selected as the most likely burial and exposure history for each site.

4 Results and discussion

4.1 Late Holocene moraine ages

At their Late Holocene maximum, the glacier on the north aspect of Vintage Peak covered a maximum area of nearly 1.9 km², while Lockie's Glacier covered approximately 1.5 km². We estimate that glacier ice would have been at least 35–40 m thick when Lockie's Glacier was near its Late Holocene maximum extent. This estimate comes from assuming a linear surface slope between the headwall and toe of the glacier of 16°and a yield stress of 100 kPa (Cuffey and Paterson, 2010). By 2020 CE, both glaciers were less than 0.1 km² and are now likely too thin to flow.

Figure 4 Lockie's Holocene glacier extent change. The boundary of Lockie's Glacier when Core 1 was exposed is uncertain but must have at least covered the site of Core 3. Four time periods are shown in this figure: (1) the Little Ice Age maximum extent of Lockie's Glacier, (2) the intermediate extent of Lockie's Glacier that must have covered the site of Core 3 but left Core 1 exposed, (3) the 2006 CE extent, and (4) the 2017 CE extent of Lockie's glacier taken from satellite imagery. Figure basemap is an oblique hillshaded digital elevation model from a 2017 aerial lidar survey.

Moraines sampled in this study were bouldery and draped over bedrock, with crests that were no more than 5 m above the bedrock surface (Sect. S1 “Data” in the Supplement). Collectively, the nine moraine boulders sampled on Vintage Peak returned a median ¹⁰Be age of 705 ± 219 a (1090–1530 CE; Fig. 5). Moraine boulders on Lockie's Glacier and Vintage Glacier moraines date between 934 and 538 a, except a single boulder on the Lockie's Glacier moraine that dates to 2849 ± 88 a.

Figure 5 (a) Probability density function of ¹⁰Be moraine ages on Vintage Peak (blue line) with kernel density plot showing individual sample ages. Black markers are the median ± interquartile range ¹⁰Be ages from moraines fronting Gilbert Glacier, 65 km north of Vintage Peak (Hawkins et al., 2021). (b) Normalized probability from organic ¹⁴C samples from overridden wood in glacier forefields and composite moraines in western Canada, interpreted to closely constrain periods of glacier advance (modified from Hawkins et al., 2021). The black line is the mean reconstructed temperature anomaly with respect to the 1800–1900 CE average for the 30–60° N latitude band from Kaufman et al. (2020). Blue shaded region is the Little Ice Age period from 1300–1850 CE.

Figure 6 Snow depth on Vintage Peak moraines. Elevation change between aerial-lidar-derived DEMs collected on 16 May 2023 (snow-covered) and 15 September 2017 (snow-free). Note the moraine crest within panel A is diffuse, while the moraine is more sharply crested in panels B and C.

Repeat lidar surveys completed in the late spring and autumn reveal the depth and distribution of snow across Vintage Peak (Fig. 6). In our survey data, moraine crests were covered by up to several meters less snow than adjacent bedrock areas. While annual snow depth will vary, the distribution of snow over the landscape will be broadly consistent each year. Manual snow depth measured at the nearby snow survey sites was ∼25 % of normal on 30 March 2023, while an automated snow survey site (3A25P, 1340 m a.s.l.) 80 km southeast of Vintage Peak similarly recorded the snowpack to be 25 % of normal in mid-May. While the lidar surveys support the assumption that sampling boulders for surface exposure dating on moraine crests will reduce the influence of snow cover on moraine exposure ages, we do still apply a snow cover correction to moraine boulders and bedrock as discussed in the Methods section.

Glaciers on Vintage Peak reached their Holocene maxima at the beginning of the classic LIA. The timing of glaciers reaching their Late Holocene maxima on Vintage Peak accords with other nearby glacier records (Fig. 5). Nearby glaciers advanced multiple times from 1.4 to 0.3 ka (Hawkins et al., 2021; Koehler and Smith, 2011; Mood and Smith, 2015b; Reyes and Clague, 2004; Ryder and Thomson, 1986). The regional radiocarbon record provides evidence of glacier advance but not the timing of glaciers reaching their maximum positions and depositing moraines. While there is good correspondence among records of glacier advance preceding the age of moraines on Vintage Peak, most glaciers in the region continued to advance after 0.65 ka. Though many glaciers in western Canada reached their greatest Holocene extents between 0.4 and 0.1 ka (Menounos et al., 2009), there is increasing evidence that glaciers were near their Late Holocene maximum positions several times as early as 2 ka (Hawkins et al., 2021; Maurer et al., 2012). Stochastic variability in glacier response to climate forcing or differing response times to climate perturbations may explain why glaciers on Vintage Peak reached their greatest Holocene extents slightly earlier than many other glaciers in western Canada (e.g., Roe, 2011).

4.2 Paired ¹⁴C and ¹⁰Be on erratic and bedrock

The erratic boulder (17-VP-07) lies approximately 65 m beyond the Late Holocene moraine that fronts Lockie's Glacier. Due to the boulder's distance from adjacent rock walls and lack of additional boulders in the area, we interpret that the boulder is most likely an erratic deposited during deglaciation of the Cordilleran Ice Sheet rather than a rockfall deposit. We do note, however, that 17-VP-07 was subangular and lacked any obvious striations that would have been more indicative of glacial transport. While we cannot definitively conclude the boulder is an erratic, we proceed with the interpretation of glacial origin. The boulder 17-VP-07 has a calculated ¹⁰Be age of 11.65 ± 0.2 ka and a ¹⁴C age of 9.74 ± 0.9 ka. When plotted on a paired isotope diagram, the erratic falls within the field of continuous exposure (Fig. 7). The discrepancy between the apparent ¹⁰Be age and ¹⁴C age could be the result of inheritance. The erratic may contain the equivalent of ∼2 kyr of inherited ¹⁰Be, which is more likely than the boulder having any significant quantity of inherited ¹⁴C. Conversely, another way to produce the observed difference in ¹⁰Be and ¹⁴C ages on the erratic is by mass shielding under seasonal snow cover. Mass shielding, in this case by snow, increases the calculated age of a surface. ¹⁴C is more sensitive to thin mass shielding, which means that the apparent ¹⁴C age increases faster than the apparent ¹⁰Be age with increasing shielding correction (Hippe, 2017). With a shielding factor that would be produced by 190 cm of 0.38 g cm⁻³ snow for 6 months of the year, both nuclides yield an equivalent exposure age of ∼14.5 ka. Given the average snow depth of 220 ± 100 cm and average density of $\mathrm{0.38}\pm \mathrm{0.67}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$ (n=123, mean ± standard deviation) at nearby snow monitoring sites which are nearly 600 m lower than 17-VP-07, this snow cover value is reasonable. Work by Menounos et al. (2017) and Seguinot et al. (2016) provides cosmogenic ages and ice sheet modeling that suggest the area around Vintage Peak was likely deglaciated around 14 ka, giving further confidence that the snow-corrected age of the erratic is an accurate deglaciation age for Vintage Peak. However, as we cannot rule out the possibility of inheritance, we constrain the age of deglaciation on Vintage Peak to 14.5–9.7 ka.

Figure 7 Paired ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ nuclide diagram with y axis normalized to the production ratio of ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ and x axis (atoms g^{−1 10}Be) normalized to the site-specific production ratio of ¹⁰Be. Dashed lines represent burial isochrons, the uppermost indicating 2 kyr burial, with each descending line representing an additional 2 kyr of burial. The uppermost solid black line is the line of continuous exposure with zero erosion, and the lower black line is the line of steady-state equilibrium; samples plotting between these two lines were continuously exposed and eroding at a given rate.

The erratic boulder and lack of apparent Younger Dryas moraine(s) indicates that glaciers on Vintage Peak have not extended beyond their Late Holocene positions since Late Pleistocene deglaciation that may have commenced as early as 14.5 ka. Some cirque glaciers in western Canada that were above the decaying CIS advanced during the Younger Dryas and deposited moraines on the landscape (Menounos et al., 2017). However, we did not observe any moraines beyond the Late Holocene moraines at Vintage Peak. Our findings on Vintage Peak further support that, on average, glacier advances during the Little Ice Age were the most extensive since the termination of the LGM.

Both ¹⁰Be and ¹⁴C bedrock surface concentrations at each core site follow expected trends based on the geomorphic position relative to the modern glacierette (Fig. 8). The greatest nuclide concentrations are found at Core 4 (17-VPC4), outside of the Late Holocene moraines, with a surface ¹⁰Be concentration of $\mathrm{1.669}\pm \mathrm{0.0167}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}$ and ¹⁴C concentration of $\mathrm{1.886}\pm \mathrm{0.445}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}$. Nuclide concentrations at Core 3 (17-VPC3, just inside of the Late Holocene moraine) and Core 2 (17-VPC2, just distal from the 2017 CE glacier terminus) are equivalent within errors (Table 1). This suggests that, within our measurement uncertainty, the bedrock surface at Cores 2 and 3 experienced nearly equivalent exposure histories over the Holocene. Cores 2 and 3 have 63 %–66 % less ¹⁰Be and 55 %–57 % less ¹⁴C than Core 4. We interpret equivalent nuclide concentrations at Cores 2 and 3 to indicate that Lockie's Glacier responded to climate perturbations in a roughly binary manner: existing near the glacier's Late Holocene maximum or retracted close to modern extents (Fig. 4).

Figure 8 ¹⁰Be and ¹⁴C concentrations in bedrock cores with elevation. Nuclide concentration is normalized relative to Core 4, which had the lowest production rate of the sample locations.

Table 1 Vintage Peak ¹⁰Be sample information.

^a Be carrier 2014.10.20 used for all samples and contains a Be concentration of 1014 ppm. ^b All Be concentrations were measured at the Purdue Rare Isotope Measurement Laboratory (PRIME). ¹⁰Be / ⁹Be ratios are not corrected for ¹⁰Be detected in procedural blanks. ^c Ages are calculated using version 3 of the online exposure age calculator formerly known as the CRONUS-Earth online exposure age calculator found at https://hess.ess.washington.edu/ (last access: 26 August 2020) (wrapper 3.0.2; muons: 1A; constants as of 26 August 2020). All ages are calculated using the Lifton–Sato–Dunai “LSDn” scaling and the default production rate. ^d Shielding by snow cover applied to bedrock core surfaces assuming 220 cm of 0.38 g cm⁻³ snow cover for 6 months of the year.

The ¹⁰Be and ¹⁴C surface concentrations of Core 1 (17-VPC1), which is assumed to have experienced variable cover by minimally erosive ice but less total ice cover than Core 2 or Core 3, are $\mathrm{1.179}\pm \mathrm{0.0351}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}$ and $\mathrm{1.401}\pm \mathrm{0.1795}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}$, respectively. Core 1 has 29 % less ¹⁰Be and 26 % less ¹⁴C than Core 4.

4.3 ¹⁴C bedrock depth profiles

Five measurements of ¹⁴C concentration along the depth of each core show decreasing or consistent, within uncertainties, nuclide concentrations. At Core 4, the ¹⁴C concentration is $\mathrm{1.886}\pm \mathrm{0.445}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}$ in the upper 0.05 m of the core and decreases to $\mathrm{1.464}\pm \mathrm{0.225}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}$ at the bottom 0.40–0.45 m of core (Fig. 9, Table S3). There is no clear pattern to measurement uncertainty with depth in each of the four cores, with measurement uncertainty ranging from 10 % to 47 % (average 25 %) of the mean concentration. Cores 2 and 3 have quasi-vertical ¹⁴C concentration profiles of equivalent magnitude, within errors (Fig. 9). Core 1 has a slight decrease in ¹⁴C concentration in the upper section of the core, but within uncertainties it is quasi-vertical as well.

Figure 9 Bedrock ¹⁴C profiles with depth. The height of each box represents the amalgamated depth of each 5 cm long sample taken every 5–10 cm along each core. The width of the boxes is the 1σ ¹⁴C measurement error. Lines connect the mean ¹⁴C concentration measurement along each core.

4.4 Monte Carlo simulations

Monte Carlo simulations of possible burial and exposure histories provide a range of solutions that adequately explain our data. Subglacial erosion and possible inheritance from previous bedrock exposure significantly influences the apparent exposure age of the bedrock. Though we are unable to independently determine the magnitude of inheritance or subglacial erosion at each site, we make reasonable estimates of these values given each site's geomorphic position.

Core 1, assuming no subglacial erosion, experienced 9.3 kyr of exposure and 4.9 kyr of burial. Assuming complete resetting of the cosmogenic nuclide inventory during the LGM, this means that Core 1 became ice-free following the LGM around 14.3 ka.

Cores 2 and 3 likely experienced some amount of subglacial erosion during the Holocene, as evidenced by striations on the sampled bedrock surface. However, if we assume negligible erosion, Cores 2 and 3 experienced 4.5 and 5.4 kyr of exposure and 5.5 and 5.9 kyr of burial, respectively. If we invoke a subglacial erosion rate of 0.005 cm yr⁻¹, then Core 2 would have experienced 6.4 kyr of exposure and 4.5 kyr of burial, while Core 3 experienced 7.6 kyr of exposure and 4.6 kyr of burial. At both sites, increasing the subglacial erosion rate notably increases the modeled exposure duration and decreases the burial duration. With a subglacial erosion rate of 0.01 cm yr⁻¹, the measured nuclide concentration at Core 2 is best explained by 8.5 kyr exposure and 3.8 kyr of burial, while Core 3 potentially experienced 9.9 kyr of exposure and 3.8 kyr of burial. Core 4 would have experienced no subglacial erosion since deglaciation following the LGM. Our nuclide measurements at Core 4 are best explained by 14.9 kyr of exposure and 2.9 kyr of burial.

While the apparent exposure age of Core 4 is coincident with the timing of regional deglaciation, the nearly 3 kyr of burial is at odds with our assumption that Core 4 has experienced continuous exposure following deglaciation. The apparent burial signal may be the result of inheritance due to insufficient resetting of the nuclide inventory in the bedrock during the LGM, mass shielding from transient till cover following deglaciation and/or more snow cover than we have accounted for, or a combination of these factors. If the apparent burial was caused by till cover, assuming a till density of 2.4 g cm⁻³, there would need to be ∼1.9 m of till to effectively reduce nuclide production below detection limits. Testing the effect of thin mass (from till or shallow ice) shielding on the bedrock nuclide inventory is beyond the scope of this study. Samples from bedrock sites with a comparable geomorphic context to Core 4 would help indicate if bedrock in this area was fully reset during the LGM or if areas that did not experience ongoing Holocene glaciation contain inherited nuclides. Monte Carlo simulations given 3 kyr of apparent ¹⁰Be inheritance estimate nearly 13 kyr of exposure and ∼2.5 kyr of burial at Core 4 under a zero subglacial erosion scenario (Fig. S7 in the Supplement).

Figure 10 Reduced χ² values from Monte Carlo modeling results, assuming zero inheritance. Each plot shows the exposure and burial duration with the lowest χ² value (red cross and values printed in plot). The colored region in each plot indicates the χ² values within the range of possible exposure and burial histories for each core and subglacial erosion rate. The total number of simulations out of 10 000 combinations of burial and exposure histories that fit within measurement errors is shown by “n sims”.

In summary, Monte Carlo simulations provide insight to the effect of varying subglacial erosion on the cumulative burial and exposure histories that can be explained by our data. At Cores 2 and 3, a small increase in subglacial erosion would have to be explained by greater total exposure duration and less burial by ice. Preserved striations on the bedrock surfaces of Cores 2 and 3 demonstrate that some amount of subglacial erosion did occur. A subglacial erosion rate of 0.005 cm yr⁻¹ indicates that ice may have covered Cores 2 and 3 until 10.9–12.2 ka. Core 4 records 14.9 kyr of exposure but contains an apparent burial signal of nearly 3 kyr that can be explained by inheritance, transient mass shielding, or a combination of the two.

4.5 Results synthesis and regional context

Based on our results and comparison to published glacier records from the region, we conclude that Vintage Peak was fully covered by the erosive Cordilleran Ice Sheet during the Last Glacial Maximum. As early as ∼14.5 ka, the Cordilleran Ice Sheet had thinned to expose most of the bedrock surface around Vintage Peak, though alpine glaciation kept bedrock within the Late Holocene glacier maxima covered. Following deglaciation, the site of Core 4 may have been covered by till for ∼3 kyr before the till was eroded away. While we cannot rule out the possibility of former till cover, we do not observe preserved Late Pleistocene till in nearby proglacial environments, which suggests this explanation is unlikely. Over the course of the Holocene, Lockie's Glacier was larger than present for 4.5–4.6 kyr and smaller than present for 6.4 kyr. Glaciers on Vintage Peak reached their greatest Holocene extent 700 ± 220 a (1090–1530 CE) and have since retreated to minimal extents.

Decreasing insolation in the Northern Hemisphere strongly influenced the growth of glaciers in western Canada from the mid-Holocene until the culmination of the Little Ice Age (Menounos et al., 2009; Solomina et al., 2015). Holocene temperature reconstruction by Kaufman et al. (2020) shows a peak in Holocene temperatures at around 6.8 ka, followed by a consistent decrease in temperature until around 1850 CE (Fig. 5). While our bedrock exposure data do not allow us to determine distinct times of glacier advance or retreat, regional temperature reconstructions and nearby glacier records suggest Lockie's Glacier likely advanced to a position near its Late Holocene maximum extent during the last 6 kyr.

Our results are at odds with the interpretations of other surface exposure dating studies, namely Menounos et al. (2017) and Darvill et al. (2022). Without correcting for snow cover, the ¹⁰Be ages for the bedrock surface at Core 4 and the erratic agree with the median age of cirque moraines across western Canada formed in response to Younger Dryas cooling (Menounos et al., 2017). However, our paired ¹⁴C–¹⁰Be dating of Core 4 and the erratic reveals a more complex history of exposure that is best explained by the significant influence of snow and in the case of Core 4 inheritance and/or transient mass shielding following deglaciation. Measurements of ¹⁴C on boulders presented in Menounos et al. (2017) would show whether the boulders have experienced a complex exposure history.

5 Conclusions

Our study used paired ¹⁴C–¹⁰Be dating on alpine bedrock surfaces and bedrock cores to elucidate glacier activity prior to the Late Holocene maximum extent of an alpine glacier in the southern Coast Mountains of British Columbia. We estimate Late Pleistocene deglaciation at Vintage Peak to have occurred around 14.5 ka. Bedrock beyond Holocene ice limits suggests limited burial that cannot be explained by snow cover alone and may be the result of previous cover by till or inherited ¹⁰Be. A broad ridgeline that would have been covered by minimally erosive ice during maximum ice extents indicates around 9.3 kyr of exposure and ∼4.9 kyr of burial since Late Pleistocene glaciation. Bedrock exposure histories close to the Late Holocene maximum extent and the much less extensive 2017 CE ice margin are nearly equivalent and indicate that the bedrock at this site may not have become exposed following Late Pleistocene deglaciation until 10–12 ka. Ice was more extensive than the 2017 CE glacier extent for around 4.5–5.9 kyr. Glaciers on Vintage Peak were likely episodically larger than present from the neoglaciation until the culmination of the LIA, in response to cooling beginning in the mid-Holocene. Nine ¹⁰Be ages from boulders on two moraines that front Vintage Peak Glacier yield a median exposure age of 700 ± 220 a (1090–1530 CE), dating the Late Holocene maximum extent of glaciers on Vintage Peak to the earlier portion of the Little Ice Age. Future work to test additional simulated exposure histories to investigate the impact of transient till cover and variable pre-LGM inheritance could further constrain glacier change throughout the Holocene. Bedrock exposure ages are influenced by several variables (exposure, burial, inheritance, mass shielding, and erosion); confident determinations of actual exposure and burial durations require independent constraints on several of these variables.