Early Holocene ice on the Begguya plateau (Mt. Hunter, Alaska) revealed by ice core ¹⁴C age constraints

Investigating North Pacific climate variability during warm intervals prior to the Common Era can improve our understanding of the behavior of ocean–atmosphere teleconnections between low latitudes and the Arctic under future warming scenarios. However, most of the existing ice core records from the Alaskan and Yukon region only allow access to climate information covering the last few centuries. Here we present a surface-to-bedrock age scale for a 210 m long ice core recovered in 2013 from the summit plateau of Begguya (Mt. Hunter; Denali National Park, Central Alaska). Combining dating by annual layer counting with absolute dates from micro-radiocarbon dating, a continuous chronology for the entire ice core archive was established using an ice flow model. Calibrated ¹⁴C ages from the deepest section (209.1 m, 7.7 to 9.0 ka cal BP) indicate that basal ice on Begguya is at least of early Holocene origin. A series of samples from a shallower depth interval (199.8 to 206.6 m) were dated with near-uniform ¹⁴C ages (3 to 5 ka cal BP). Our results suggest this may be related to an increase in annual net snow accumulation rates over this period following the Northern Hemisphere Holocene Climate Optimum (around 8 to 5 kyr BP). With absolute dates constraining the timescale for the last >8 kyr BP, this paleo-archive will allow future investigations of Holocene climate and the regional evolution of spatial and temporal changes in atmospheric circulation and hydroclimate in the North Pacific.

1 Introduction

Arctic surface temperatures have increased more than twice as fast as global temperature during the early 20th century and since the 1970s (Bengtsson et al., 2004; Tokinaga et al., 2017; Svendsen et al., 2018). Recent modeling results suggest that during the early 20th century, as the Pacific Decadal Oscillation (PDO) transitioned to a positive phase, there was a concomitant deepening of the Aleutian Low that warmed the Arctic through poleward low-level advection of extratropical air (Svendsen et al., 2018). The impact of Pacific multidecadal variability on Arctic warming has considerable implications for sea ice extent (Screen and Francis, 2016) and hence the possible linkage between Arctic amplification, sea ice loss, and enhanced midlatitude winter variability (Cohen et al., 2014, 2018, 2019; Francis et al., 2017; Screen et al., 2018; Blackport et al., 2019). Whether the present positive PDO conditions will persist and contribute to Arctic warming at an even higher rate in the future remains a fundamental question (Svendsen et al., 2018). A longer-term perspective on Pacific decadal variability and the teleconnection between the tropical Pacific, North Pacific, and the Arctic, particularly during warm intervals in the Holocene outside those captured in the instrumental record, would be an important contribution to this problem (e.g., Park et al., 2019). High-mountain ice cores in the North Pacific region have the advantage of sampling atmospheric moisture (e.g., snow) and aerosol deposition, as well as preserving physical characteristics (e.g., melt), all of which can be related to Pacific climate processes (Zdanowicz et al., 2014; Osterberg et al., 2017; Winski et al., 2018) if Holocene (or greater) length records can be recovered.

Figure 1 Location map of North Pacific ice core sites and the age of the oldest ice dated from each location: Begguya (Mt. Hunter; this study), McCall Glacier (Klein et al., 2016), Aurora Peak (Tsushima, 2015), Mt. Wrangell (Yasunari et al., 2007), Bona-Churchill (Porter et al., 2019), Mt. Logan (Fisher et al., 2008), Eclipse Icefield (Yalcin et al., 2007), and Mt. Waddington (Neff et al., 2012).

The general timing of deglaciation in Alaska (Brooks Range, Central Alaska Range, and southern Alaska) was determined based on terrestrial cosmogonic radionuclides, lichenometry, and radiocarbon dating to between 10 and 20 kyr BP (Dortch, 2007). Following the Last Glacial Maximum (LGM), glaciers in the Brooks Range retreated up valley to, or even within, their modern limits by ca. 15 ka (Pendleton et al., 2015). Given the small extent of the Brooks Range glaciers prior to the Holocene thermal maximum, during which some glaciers in southern Alaska disappeared entirely (Barclay et al., 2009), it is possible that the Brooks Range glaciers may have disappeared as well. In the Central Alaska Range, reaching much higher altitudes and considering today's glacier extent, this is rather unlikely. Nevertheless, it is unclear where preserved ice from the early Holocene (or older) can be found in basal layers of these glaciers. Most of the ice cores recovered from the Alaskan and Yukon region did not reach bedrock and are thus limited in the time covered, reaching back a few centuries only (Fig. 1). The Prospector-Russell Col (PRCol) ice core from Mt. Logan is an exception, having an estimated bottom age of ∼20 kyr BP based on the assumption that the significant depletion in the water stable isotope ratios observed in the very bottom section of the core is a signal of the LGM cold conditions (Fisher et al., 2008). The PRCol chronology is further constrained by a large δ¹⁸O minimum and coeval increases in deuterium excess and Ca²⁺, which are assigned to the 4.2 ka BP event (Walker et al., 2019), as well as tephra from the large Alaskan eruption of Aniakchak (3.6 kyr BP; Walker et al., 2019). The PRCol record serves as a global auxiliary stratotype for the Middle and Late Holocene subdivision boundary (Walker et al., 2019). However, there are no chronologic tie points in the PRCol record prior to the 4.2 ka BP event (Walker et al., 2019).

New surface-to-bedrock ice cores were recovered from the Begguya plateau (Mt. Hunter; Denali National Park, Alaska; 62.93^∘ N, 151.09^∘ W; Fig. 1) in 2013 at 3900 m elevation (Winski et al., 2017). The two surface-to-bedrock cores (DEN-13A, DEN-13B) reached depths of 211.2 and 209.7 m, respectively. Analysis of the upper 190 m of DEN-13B (2013 to 810 CE) revealed that snow accumulation at the drilling site has doubled since ∼1840 CE, coeval with the warming of western tropical Pacific sea surface temperatures (Winski et al., 2017) and intensification of the Aleutian Low system (Osterberg et al., 2014, 2017). The same core also shows a 60-fold increase in water equivalent of total annual melt between 1850 CE and present, which suggests a summer warming rate of 1.92±0.31 ^∘C per century during the last 100 years in the altitude range of 3900 m (Winski et al., 2018). The Begguya melt layer record is significantly correlated with surface temperatures in the central tropical Pacific through a Rossby-wave-like pattern that enhances temperatures over Alaska (Winski et al., 2018). Taken together, these hydroclimate changes are consistent with linkages between Pacific decadal variability and Arctic hydroclimate changes seen in the observational record (Svendsen et al., 2018), and they demonstrate that the North Pacific hydroclimate response since 1850 CE is unprecedented in the past millennium.

The annual layer-counting-based chronology of the Denali core results in an ice age of 1203±41 years at a depth of 190 m (152.8 m w.e., water equivalent; Winski et al., 2017). Below that depth, annual layering was less consistent due to the loss of seasonal resolution caused by the glacier-flow-induced thinning of layers. However, based on previously reported age–depth scales of ice cores from cold, high-elevation glaciers frozen to bedrock, the bottom 20 m of ice may contain most of the record in terms of time, covering the Holocene and potentially even reaching into the Last Glacial period (Uglietti et al., 2016; Licciulli et al., 2020). The Denali ice core therefore provides the possibility of establishing a new Holocene North Pacific hydroclimate record reaching beyond the Common Era if a precise and absolute-dated chronology can be established in the bottom 20 m of the core. The water-insoluble organic carbon (WIOC) and dissolved organic carbon (DOC) ¹⁴C dating method has been validated and applied for multiple midlatitude ice cores (e.g., Jenk et al., 2009; Uglietti et al., 2016; Hou et al., 2018; Fang et al., 2021). The technique makes use of the transport and deposition of carbonaceous aerosols onto the glacier. Before the industrial period, carbonaceous aerosols were mainly emitted from the living biosphere and from biomass burning. Consequently, this carbon reflects the contemporary atmospheric ¹⁴C content (Jenk et al., 2006). After deposition, the WIOC and DOC are incorporated into glacier snow, firn, and ice and undergo radioactive decay with a half-life of 5730 years (Godwin, 1962). Here we report results from ¹⁴C analysis of the bottom 60 m of the Denali ice core. These absolute dates extend the existing late Holocene Begguya chronology (Winski et al., 2017), providing the first high-latitude Northern Hemisphere ice core chronology based on absolute dates from radiometric methods. We discuss our results in relation to Holocene ice extent and climate in the North Pacific region.

2 Methods

2.1 Annual layer counting

Two surface-to-bedrock ice cores (DEN-13A, DEN-13B) were drilled in 2013 at 3900 m elevation above sea level (a.s.l.) from the saddle between the north and middle peaks of Begguya (Mt. Hunter), Alaska (Winski et al., 2017, 2018; Osterberg et al., 2017; Polashenski et al., 2018). The annual layer counting for DEN-13B was previously published (Winski et al., 2017) and is only briefly described here. The timescale from 2013 to 1777 CE was determined by counting annual oscillations in δ¹⁸O (summer peak), melt layers (summer peak), magnesium (spring peak), dust (spring peak), liquid conductivity (summer peak), ammonium (summer peak), and methanesulfonic acid (MSA; late summer–fall peak), consistent with previous North Pacific ice cores (Yasunari et al., 2007; Osterberg et al., 2014; Tsushima et al., 2015). Between 1777 and 1500 CE annual layer counting is based on annual oscillations of δ¹⁸O, δD, dust concentration, and liquid conductivity that were measured at higher resolution than the other analytes, while conductivity and dust concentrations were exclusively used to date the ice core from 1500 back to 810 CE. For this study, the counting based on these two parameters has been extended back in time (see Sect. 3.2).

Table 1 The ¹⁴C results of the Denali ice core samples (DEN-13B), given as F¹⁴C, ¹⁴C ages, and calibrated ¹⁴C ages. For ¹⁴C calibration, chronological layering was assumed (sequential deposition; see main text). Samples were dated using the WIOC fraction, except for section 235 in which the DOC fraction was analyzed. Numbers of the carbon amount available for ¹⁴C AMS analysis as well as the concentration of WIOC (DOC) in the sample are also provided. Additionally shown is the range of the dating based on ALC (range from top to bottom depth of section) and the final age scale (inverse ice flow model; provided is the uncertainty range at mid-depth of the analyzed section).

^a Following recommendations, samples with a carbon mass of significantly less than 10 µg C were not considered (Uglietti et al., 2016). ^b Results from the DOC fraction. ^c After correction for in situ ¹⁴C production (Fang et al., 2021; see main text).

2.2 Denali ice core ¹⁴C analysis

A total of 16 samples were selected from the lower portion of DEN-13B (Table 1). Because WIOC concentrations at this site were assumed to be low, ice samples of at least 1 kg of mass were cut, aiming for extracted yields of carbon that allow dating with a reasonable uncertainty of 10 %–20 % (>10 µg C; Uglietti et al., 2016). In order to process such large sample volumes, a splitting of the sample for melting was required and the overall filtration time had to be increased. Using artificial ice produced from ultra-pure water, the adapted procedures were tested to reach low blanks similar to the ones previously achieved for smaller samples (Jenk et al., 2009; Uglietti et al., 2016; Fang et al., 2019). Otherwise, the samples for WIOC ¹⁴C dating were prepared following the protocol described in Uglietti et al. (2016) with a brief summary provided here. In order to remove potential contamination in the outer layer of the ice core, pre-cut samples from the inner part of the core were rinsed with ultra-pure water. After the melting of the sample in a pre-cleaned jar (1 L, PETG, Semadeni), due to the size split in two, the carbonaceous particles contained as impurities in the sample ice were filtered onto a prebaked quartz fiber filter (Pallflex Tissuquartz 2500QAT-UP). Potential particulate carbonates also remaining on the filter were removed by acidifying three times with 0.5 µL of 0.2 M HCl. These initial steps were performed in a class 100 laminar flow box to ensure clean conditions. At the University of Bern (Laboratory for the Analysis of Radiocarbon with AMS – accelerator mass spectrometer – LARA laboratory) the WIOC samples were then combusted in a thermo-optical OC/EC (organic carbon and elemental carbon) analyzer (Sunset Modeldoc4L, Sunset Laboratory Inc, USA) with a non-dispersive infrared sensor for CO₂ quantification, using the established Swiss 4S protocol for OC/EC separation (Zhang et al., 2012). Being coupled to a 200 kV compact accelerator mass spectrometer (AMS, mini radiocarbon dating system MICADAS), equipped with a gas ion source and a gas interface system (GIS; Ruff et al., 2007; Synal et al., 2007; Szidat et al., 2014), the LARA Sunset–GIS–AMS system (Agrios et al., 2015, 2017) allowed for final, online ¹⁴C measurements of the CO₂ produced from the WIOC fraction.

For the deepest sample from ∼209 m depth (Denali235) the available amount of ice was very limited (∼200 g). To ensure a sufficient mass of carbon for final AMS analysis, the ¹⁴C dating was performed on the DOC fraction, for which a higher concentration compared to the WIOC fraction is expected (Legrand et al., 2013). By a catalyzed UV oxidation in a dedicated system, DOC was converted to CO₂ which was then cryogenically trapped and flame sealed in glass ampules for final AMS analysis. Details can be found in Fang et al. (2019).

All ¹⁴C results are expressed as fraction modern (F¹⁴C), which is the ¹⁴C $/$ ¹²C ratio of the sample divided by the same ratio of the modern standard referenced to the year 1950 CE (NIST standard oxalic acid II, SRM 4990C), both being normalized to −25 ‰ in δ¹³C to account for isotopic fractionation. Daily AMS calibration was performed using sets of modern (NIST oxalic acid II, SRM 4990C, F¹⁴C $=\mathrm{1.3407}\pm \mathrm{0.0005}$) and fossil standards (sodium acetate, Sigma-Aldrich, no. 71180, F¹⁴C $=\mathrm{0.0018}\pm \mathrm{0.0005}$). Final values presented in Table 1 are the AMS F¹⁴C raw data after corrections accounting for constant contamination and cross contamination in the Sunset–GIS–AMS system (or GIS–AMS system for DOC, respectively) and the overall procedural blank contribution introduced from the preparation of ice samples to final AMS analysis. F¹⁴C of DOC was corrected for the contribution from ¹⁴C in situ production following Fang et al. (2021). The applied small shift in F¹⁴C of 0.019±0.010 was derived using an in situ production rate of 260.9 ¹⁴C atoms g${}_{\mathrm{ice}}^{-\mathrm{1}}$ yr⁻¹ as the best estimate for the site latitude and elevation (Lal et al., 1987; Lal and Jull, 1990; Lal, 1992), an average accumulation rate of 1.0±0.5 m w.e. (a best initial guess based on the annual values from Winski et al., 2017, ranging from 0.2 to 2.0 m w.e. for the time period 810 to 2013 CE), and assuming an average incorporation into DOC of 18±7 % (Hoffmann, 2016). This correction shifts the calibrated age to be 300±200 years older, with uncertainty being fully propagated as for all other ages. Note that the upper estimate does not exceed the achieved dating precision defined by the analytical uncertainty (see Table S1 in the Supplement). For all samples, calibrated radiocarbon ages were derived by calibrating final F¹⁴C values using OxCal v4.4.4 (Ramsey, 2021) with IntCal20 (the Northern Hemisphere calibration curve; Reimer et al., 2020) and the OxCal in-built sequence model (Bayesian-approach-based deposition model; Ramsey, 2008, 2017). All calibrated ¹⁴C ages are presented as the 1σ range in years before present (cal BP, with BP referring to the year 1950 CE).

3 Results

3.1 Englacial stratigraphy

Around the Begguya drill site, no folding was observed in ground-penetrating radar (GPR) data, and the bedrock geometry appears to be uncomplicated (Campbell et al., 2013). New radar data were collected in 2022. Ice thickness, bed topography, and internal stratigraphy of the core site were mapped using GPR (10 MHz center frequency radar system, Blue Systems Integration). Standard processing techniques were applied to the data: clipping stationary periods, applying horizontal stacking, bandpass filtering, and correction for antenna separation (Lilien et al., 2020). Data were interpolated for standard trace spacing and then migrated using the SeisUnix sumigtk routine. Clear, visible layering is evident in the majority of the ice column; however, the interpretation of the stratigraphy at depth is complicated by sidewall reflections produced from the trough beneath the ice core site. There is no conclusive evidence from these data of either stratigraphic continuity or discontinuity in the bottom-most 10 m of ice (Fig. 2). Future measurements using the millimeter-precision capabilities of autonomous phase-sensitive radar (Brennan et al., 2014) would be beneficial to resolve englacial stratigraphy close to the bedrock.

Figure 2 Ground-penetrating radar profile collected with 10 MHz BSI radar across the Begguya plateau in 2022. Standard processing techniques were applied to the data using ImpDAR (Lilien et al., 2020). The two-way travel time (TWTT) is plotted on the y axis on the right side. The Denali ice core drilling (DEN-13B) is indicated by the vertical blue line.

3.2 Annual layer counting

Annual layer counting (ALC), previously published in Winski et al. (2017) back to 810 CE (Sect. 2.1.), was extended back to 339 CE, i.e., for the top 197 m. The uncertainty in the ALC chronology back to 810 CE was estimated through statistical comparisons among individual layer positions indicated by three individuals (see Winski et al., 2017, for details). By 1900 CE, uncertainty estimates are ±4 years, increasing to ±10 years at 1500 CE and ±30 years by 810 CE (190.05 m). Only one individual (DW) performed ALC below 190 m, prohibiting a similar approach to estimate uncertainties, but we estimate an uncertainty of around ±60 years at 339 CE. These estimates are for ALC only and do not consider additional, constraining information from time horizons. There is no offset between the timescale and inferred volcanic eruptions as indicated by peaks in sulfate, chloride, and conductivity during the 19th and 20th centuries, indicating that an accuracy within ±2 years throughout the last 200 years is likely. The sulfate and chloride peaks in the 18th century used for chronology validation (inferred as Laki, 1784 CE, and Pavlof, 1763 CE) were offset by 1 year from the ALC chronology. Additionally, ¹³⁷Cs concentrations in the Denali core strongly peak in the layer assigned to the year 1963 CE, 1 year after the most extensive atmospheric testing of nuclear weapons, which matches the ¹³⁷Cs residence time in the atmosphere (Tian et al., 2007; Winski et al., 2017).

3.3 Denali ice core ¹⁴C data

Air masses leading to precipitation on Begguya (∼3900 m a.s.l.) originate predominantly from the Pacific and contain relatively low organic aerosol concentrations (Haque et al., 2016; Choi et al., 2017). The WIOC concentration in the Denali core is thus significantly lower than in ice cores from the Alps. The WIOC concentrations range from 6 to 31 µg C kg⁻¹ ice with an average of 13±7 µg C kg⁻¹ (Table 1). This is slightly higher than in snow at Summit, Greenland (4.6 µg C kg⁻¹; Hagler et al., 2007), but only about half of the pre-industrial WIOC concentrations in European Alpine ice cores, with 24±9 µg C kg⁻¹ (Legrand et al., 2007) and 32±18 µg C kg⁻¹ (Jenk et al., 2009) from Colle Gnifetti, Monte Rosa, Switzerland, and 24±7 µg C kg⁻¹ from Fiescherhorn Glacier (Jenk et al., 2006). In agreement with findings from previous studies (Legrand et al., 2007), the concentration of DOC (80 µg C kg⁻¹), measured in the deepest sample, was significantly higher than the concentration of WIOC.

The ¹⁴C calibration was performed using the OxCal in-built sequence model (Ramsey, 2008, 2017; see Methods). The assumption that samples are in chronological order allows statistical constraints for the most likely age distribution of the individual samples in the sequence. This assumption of chronological ordering will be discussed below. Samples containing less than 10 µg C are generally characterized by a wide range of age probability. A reduction in the dating precision for those samples is expected due to the small carbon amount available for analysis. Small amounts on the one hand cause reduced AMS measurement precision (lower ¹²C current and fewer ¹⁴C counts) and a lower, unfavorable signal-to-noise ratio (i.e., the ratio between size of sample and procedural blank) on the other hand. Combined, this leads to a larger overall analytical uncertainty, finally translating into a wider range of possible ages. Although we used a considerable amount of ice for each sample (∼1 kg), the total carbon amount in five samples was significantly below this 10 µg C threshold recommended to obtain a reliable dating with a final uncertainty <20 % for samples older than around 1000 years (Uglietti et al., 2016). These samples will thus not be discussed in the following (but can be found in the Supplement, together with calibration results without sequence constraints).

Calibrated ¹⁴C ages range from 0.3±0.3 ka cal BP at 166.2 m (131.4 m w.e.) depth to 8.4±0.6 ka cal BP for the deepest sample (Denali235; 209.1 m or 169.8 m w.e.), the last sample above bedrock (0.6 m). These results show the characteristic exponential increase in age with depth, expected for a cold glacier archive due to the associated ice flow dynamics (e.g., Dansgaard and Johnsen, 1969; see also Sect. 4.1.), and most importantly, reveal ice of early Holocene origin in the Denali ice core (Table 1 and Fig. 3). The absolute ages from radiocarbon dating are in agreement with the independently derived ages from the ALC reported in Winski et al. (2017), extended back to 339 CE in this study (see Sect. 3.2 “Annual layer counting”). For the youngest sample, Denali183 from a depth of 166.2 m (131.4 m w.e.), and for Denali214 from 192.6 m (155.0 m w.e.), the 1σ age range is 4–679 and 958–1410 yr cal BP, respectively – similar to the corresponding ALC-derived ages of 340–380 and 1200–1410 yr BP. The 1σ¹⁴C age range for Denali210-211 at 189.5 m (152.3 m w.e.) is 527–930 yr cal BP and with a possible age of 930 yr cal BP only slightly younger than the ALC-derived age range of 1020–1200 yr BP (in agreement within the 2σ range of 317–1174 yr cal BP).

Figure 3 Calibrated ¹⁴C age probability distributions for samples from the Denali ice core (DEN-13B), as derived in OxCal v4.4.4 using the IntCal 20 radiocarbon calibration curve (Ramsey, 2021; Reimer et al., 2020). Light green areas indicate the a priori age probabilities and the dark green areas the a posteriori probabilities when sequential ordering of samples is assumed (see main text). The agreement index (A) indicates overlap between these two probability functions. A value <60 indicates poor agreement (see main text). The 1σ and 2σ range is indicated by the lines below the probability distribution areas.

Samples of indistinguishable ages, with regard to the dating uncertainty achieved (i.e., analytical precision), were observed in the depth interval from 200.3 to 206.2 m (161.9 to 167.2 m w.e.). This interval corresponds to a time period from around 3.2 to 4.3 kyr BP. For the respective samples (Denali223, Denali224-225, Denali229-230, Denali231), a low agreement index (denoted as A in OxCal) resulted for the applied ¹⁴C calibration approach. A indicates the level of agreement between the probability function derived by the ordinary calibration approach (a priori distribution) and the calibration with an additional constraint (a posterior distribution; see OxCal and Ramsey, 2008, 2017, for more details). Distributions are shown in Fig. 3. A value of 100 indicates no alteration in the distribution (100 % or unity), while a value lower than 60 indicates a warning to check for the validity of the underlying assumption, i.e., (i) a non-sequential layering of samples or (ii) the presence of analytical outliers. It is apparent from Fig. 3 that the two samples with lowest A (<10), Denali223 and Denali231, are also characterized by an exceptionally large uncertainty. For the batch of samples with AMS lab IDs BE-10013.1.1 to BE-10022.1.1 (Table 1; see also Supplement Fig. S1 and Table S1), the contribution to the final overall uncertainty from AMS analysis alone was around twice as much as what typically can be achieved for samples of that carbon mass. For that measurement day, we also observed above-average uncertainties for the measured sets of AMS calibration standards, with a slight elevation in the fossil standard value (+0.02 in F¹⁴C; see Methods). This is an indication of non-ideal AMS conditions due to sub-optimal instrument tuning on the one hand and an elevated, potentially non-stable background that day on the other hand. Thus, neither ¹⁴C ages nor the englacial stratigraphy gives sufficient evidence to conclude a non-sequential ordering of samples (i.e., an age reversal in the Denali ice core). Additionally, there is evidence from other studies from the region suggesting hydrological changes between around 4 to 2 kyr BP (e.g., increased lake levels and precipitation; see Discussion), which coincides with the time period in question here. Because increased accumulation rates would lead to a reduced increase in age per unit depth, an unambiguous resolving of the sequence then depends on the achievable analytical uncertainty. Having pushed the limit of the analytical method with the small amounts of carbon available for ¹⁴C analysis and considering all the above, we thus exclude assumption (i) and are confident that the applied ¹⁴C sequence calibration approach does provide us with the most accurate dates.

4 Discussion

4.1 Denali ice core chronology

Modeling the age scale in high-elevation mountain ice cores can be attempted either by applying rather simple glaciological one-dimensional (1D) flow models (e.g., Nye, 1963; Dansgaard and Johnsen, 1969; Bolzan, 1985) or by much more complex three-dimensional (3D) models based on a suite of observational data from the glaciological survey (e.g., Campbell et al., 2013; Licciulli et al., 2020). Independent of model complexity, age scale modeling, particularly of mountain glaciers, is strongly challenged to provide accurate or even conclusive ages along the profile at a specific point on the glacier (e.g., the ice core drill site; Campbell et al., 2013; Licciulli et al., 2020). This is especially the case close to bedrock, where ice flow can become highly complex and because past annual net accumulation rates with potential variations over time are unknown. Layers of known age along a glacier depth profile, e.g., from ice core dating, can provide crucial model constraints, allowing free model parameters to be tuned for a best fit between observations and model output. For a defined point, moving along a single axis (bed to surface), 1D models benefit from their simplicity to do so (less parameters). 1D models have been applied for decades to obtain continuous age–depth relationships at sites on polar ice sheets (e.g., Dansgaard and Johnsen, 1969), thereby also accounting for past changes in accumulation rates by inverse modeling approaches (e.g., Buiron et al., 2011; Buchardt and Dahl-Jensen, 2008). However, applications to sites from high-mountain glaciers are more recent (e.g., Jenk et al., 2009; Uglietti et al., 2016).

In the case of the Denali ice core, accurate dating by ALC supported with independent time horizons for the upper two-thirds of the core and absolute-dated horizons for the deep section of the core (¹⁴C dates) is available. Winski et al. (2017) developed a well-defined age scale for the upper part of the core based on ALC supported by distinct time horizons. Since age–depth relationships are less challenging to model in the upper 90 % of the ice core, because of relatively moderate layer thinning and little if any influence from bedrock, Winski et al. (2017) used a combination of 1D modeling and a 3D glacier flow model developed for this site (Campbell et al., 2013) to determine a significant increase in accumulation rates since around 1850 CE. Therefore, significant changes in net accumulation rates at the Denali ice core drill site should be expected to a have also occurred in the more distant past.

Figure 4 Denali ice core (DEN-13B): annual layer counting (ALC), dating horizons (¹⁴C, volcanos, ¹³⁷Cs peak), and modeled, continuous age–depth relationship (1D ice flow model; see main text). (a) Model output for constant accumulation rates (b, in m w.e. yr⁻¹). (b) Modeled age–depth relationship for variable b (inverse model). (c) Zoom of the deepest part. All error bars indicate the 1σ uncertainty.

Due to its simplicity, we used the 1D two-parameter model (2p-model; Bolzan, 1985) to provide a first, best estimate for a continuous age–depth relationship from surface to bedrock, building on the available data points presented. The 2p-model is based on a simple analytical expression for the decrease in the annual layer thickness with depth and has 2 degrees of freedom, the mean annual net accumulation rate b and the thinning parameter p, characterizing the strain rate function; both are assumed to be constant over time. Knowing the glacier thickness of 209.7 m from the ice core length (supported by ground-penetrating radar data; 170.4 m w.e.) and with all depths converted from meters to meters water equivalent based on the ice core density profile allowed us to find the best solutions for b and p by fitting the model (least squares approach, as described in Fahnestock et al., 2001) through the time horizons in the Denali ice core (volcanos, ¹³⁷Cs, ¹⁴C). The derived value for p was 0.79±0.01. The resulting value of b of 1.5±0.1 m w.e. yr⁻¹, representing the mean annual net accumulation rate for the entire period covered by the ice core, is similar to the recently observed 21st century values. It is however significantly higher than the average value of around 0.5 m w.e. yr⁻¹ previously determined for the last 810 years (ranging from around 0.3 to 1.5 m w.e. yr⁻¹; Winski et al., 2017). This is no surprise, considering the likelihood that similar variations may also have occurred further back in time. As a consequence of being constrained by the age of dated layers, the model results are in agreement with the observational data for the total time period covered within the ice column. However, at various depths along the depth profile, a significant offset between model output and data can be observed (Fig. 4a). Again, this is not unexpected considering the fact that the accumulation rate was kept constant in the model, while significant changes over time are known to have occurred (Winski et al., 2017). In Fig. 3a, the effect on model results for variations in b is illustrated (runs with b equal to 0.5, 1.5, and 2 m w.e. yr⁻¹, with p determined as before).

To achieve our final goal, obtaining a continuous age–depth relationship based on the absolute dating presented, we next applied a simple inverse modeling approach. We tightly fit the model to the experimental data by numerically solving for the exact value of b for each depth with a determined age (p and H as before). To reduce and account for potential noise in the data, an uncertainty-weighted three-point running mean to obtain the non-steady-state values for b was calculated (starting from top, then reversed from bottom, thereby propagating the values for continuity). These values, interpolated for depths between the dated layers, were finally used for model input, yielding a continuous age–depth relationship (Fig. 4b and c). All uncertainties have been fully propagated throughout the calculations (from analysis to modeling). We derived annual net accumulation rates of 0.5±0.1 m w.e. yr⁻¹ at around 1000 CE, eventually increasing to a 20th century average value of 1.1±0.2 m w.e. yr⁻¹ (Fig. S2). This is in good agreement with what was determined previously by Winski et al. (2017) for the corresponding periods, based on results from different models investigated (for the 3D model considered best, 0.25 m w.e. yr⁻¹ around 1000 CE, with models ranging from 0.05–0.7 m w.e. yr⁻¹, and 1.1±0.3 m w.e. yr⁻¹ for the 20th century average). During the Holocene Climate Optimum (around 8 to 5 kyr BP, Kaufman et al., 2016) we obtained net accumulation rates of 1.2±0.3 m w.e. yr⁻¹, similar to the average rate observed since 1950 CE, followed by higher rates of 1.7±0.4 m w.e. yr⁻¹ from around 4.3 to 3.2 kyr BP. Then, the rates decrease over the next 500 to 1000 years to around 0.4±0.2 m w.e. yr⁻¹. See Sect. 4.3 for further discussion. Our derived age–depth scale results in ages of 9–14 kyr BP at 0.5 m above bedrock, strongly suggesting the presence of at least early Holocene ice at the Denali ice core drill site.

Table 2 Overview of existing North Pacific ice cores.

^a McCall Glacier (Klein et al., 2016). ^b Aurora Peak (Tsushima, 2015). ^c Begguya (this study). ^d Mt. Wrangell (Yasunari et al., 2007; Sasaki et al., 2016). ^e Bona-Churchill (Porter et al., 2019). ^f Mt. Logan (Fisher et al., 2008). ^g Eclipse Icefield (Yalcin et al., 2007). ^h Mt. Waddington (Neff et al., 2012).

4.2 Ice core chronologies in Eastern Beringia

So far, existing ice cores from Eastern Beringia (Table 2) were dated with ages covering less than the last millennium except for the Denali core discussed in this study and the 188 m long PRCol core (Fig. 1) drilled to the bed surface on the summit plateau of Mt. Logan in 2001 and 2002. The older part of the PRCol core was dated based on a signal interpreted as the Younger Dryas to Holocene transition (sudden reduction in electrical conductivity coinciding with a drop in δ¹⁸O and an increase in various chemical species) and a bottom age estimate from an ice flow model of about 20 ka (Fisher et al., 2008). Another ice core from Mt. Logan (King Col; 60.59^∘ N, 140.60^∘ W; 4135 m a.s.l.) was drilled in 2002 reaching a depth of 220.5 m. This core was not dated, but a potential age range of 0.5 to 1.3 ka was estimated based on modeling results (Shiraiwa et al., 2003). The 152 m ice core drilled in 2008 on the McCall Glacier was dated by using a combination of ALC and specific horizons. The upper 37 m of ice date back to 65 years, and the full 152 m core was estimated to cover more than 200 years, but no actual dating of the lower section was performed (Klein et al., 2016). The Aurora Peak site is located southeast of Mt. Hayes, and the ice core was also drilled in 2008. The total ice thickness at the drilling site is 252±10 m, but this core (180.17 m) did not approach the bed surface. By annual layer counting, the estimated bottom age of the Aurora Peak core is about 274 years (Tsushima et al., 2015). Two cores were collected at Eclipse Icefield in 2002. The chronology of these cores is based on multi-parameter ALC of seasonal oscillations in the stable isotope (δ¹⁸O) and major ion records (Na⁺) supported by the identification of volcanic horizons. The longest core (345 m) covers the period 1000 to 2002 CE (Yalcin et al., 2007) but did not reach bedrock. In 2004, a 212 m ice core was drilled from Mt. Wrangell. The ice depth in the summit caldera is probably over 900 m, but the definite bottom has not yet been detected (Benson et al., 2007). For this core, a short 12-year record of dust and δD was reported in Yasunari et al. (2007), and dating was later extended back to 1981 (23 years) at the depth of 100.1 m (Sasaki et al., 2016). The record from Mt. Waddington only covers a period of 1973–2010 CE (Neff et al., 2012). The total length of the Mt. Waddington core is 141 m, but the total ice thickness at the drilling site is about 250 m. The ice core from Bona-Churchill reached bedrock at a depth of 460 m, but the age–depth scale has only been established for the last ∼800 years (depth of 399 m); the deepest ice is estimated to exceed 1500 years in age (Porter et al., 2019).

Because none of the cores from the Eastern Beringia region was drilled to the bed surface and no ice close to the bed was ever dated by an absolute dating method, no conclusive evidence about the age of the oldest glacier ice preserved in this region exists so far. In this study, we achieved a first complete and absolute (radiometric) dating through the first application of ¹⁴C analysis on a high-latitude Northern Hemisphere ice core from Begguya, which reached the bed surface. Our results, with calibrated ¹⁴C ages of 7.7 to 9.0 kyr BP close to the bottom (0.61 m above bedrock) and model-based indication for potentially even older ice further below (>12 kyr BP), clearly indicate that glaciers in this region can be of early Holocene or even Pleistocene origin.

Table 3 Regional paleoclimate events.

Figure 5 Regional paleoclimate changes as reported in previous studies (Anderson et al., 2005a, b, 2006, 2016, 2019; Broadman et al., 2020; Clegg and Hu, 2010; Denton and Karlén, 1977; Heusser et al., 1985; Jones et al., 2019; Kelly et al., 2013; Schiff et al., 2009; Solomina et al., 2015) and the period of elevated annual net accumulation rates indicated in the Denali ice core DEN-13B (this study; see main text).

4.3 Possible implications for Holocene hydroclimate in Eastern Beringia

In recent decades, extensive work has been done in the North Pacific region to characterize Holocene hydroclimate (Table 3, Fig. 5). Following a modest Holocene thermal maximum that was 0.2–0.5 ^∘C warmer than the last millennium average (Kaufman et al., 2016), although 1.7 ^∘C cooler than present (Porter et al., 2019), glaciers across the region advanced synchronously at about 4.5 kyr BP (Solomina et al., 2015). This Neoglaciation continued through 3.5 to 2.5 kyr BP in the Yukon Territory based on past tree line variations, lake levels, and carbonate oxygen isotopes (Denton and Karlén, 1977; Anderson et al., 2005a). While a mid-Holocene temperature decrease may have played a role, Denton and Karlen (1977) hypothesized that an increase in regional precipitation contributed to the regional Neoglaciation, a conclusion also reached by later studies (e.g., Anderson et al., 2011).

Concurrent with this Neoglaciation, effective moisture rose across much of the region. Based on pollen reconstructions, Heusser et al. (1985) inferred a doubling of southern Alaskan mean annual precipitation from around 3.9 to 3.5 kyr BP (Fig. 3). Clegg and Hu (2010) found that effective moisture, particularly during winters, increased markedly between 4.0 and 2.5 kyr BP. Hansen and Engstrom (1996) suggested cooler and wetter conditions in Glacier Bay at around 3.4 kyr BP. At Jellybean Lake and Marcella Lake, lake levels were high between 2.0–4.0 kyr BP (Anderson et al., 2005a, b), which was attributed to changes in the strength and positions of the Aleutian Low (Anderson et al., 2005b), consistent with the more recent interpretation of hydroclimate changes from the Denali ice core (Winski et al., 2017; Osterberg et al., 2017). Records from Mica Lake (Schiff et al., 2009) and Sunken Island Lake (Broadman et al., 2020) showed wetter conditions associated with a stronger Aleutian Low at 4 kyr and 4.5 kyr BP, respectively. Greenpepper Lake experienced high lake levels from 2–5 kyr BP (Anderson et al., 2019), and a major shift from moss to sedge occurred at Horse Trail Fen, concurrent with a large isotopic anomaly at 3–4 kyr BP (Jones et al., 2019). At the same time, paleoenvironmental records showed a decrease in wildfire (Anderson et al., 2006; Kelly et al., 2013).

Together, previous work indicates an enhanced flux of moisture into the region, likely associated with a strengthened Aleutian Low, sometime near 4 kyr BP (Anderson et al., 2016). The Denali ice core may provide corroborating evidence for this mid-Holocene shift in hydroclimate (Table 3, Fig. 5). As presented before, samples of indistinguishable ages, at least for the achieved analytical precision, were observed in the depth interval from 200.3 to 206.2 m (161.9–167.2 m w.e.), corresponding to the modeled time period from 4.3±0.5 to 3.2±0.5 kyr BP (see Sects. 3.3 “Denali ice core ¹⁴C data” and 4.1 “Denali ice core chronology”). Elevated snow accumulation provides a possible explanation for this clustering of dates and would support many previous studies. While our model results based on ¹⁴C ages are consistent with existing interpretations of mid-Holocene changes in regional precipitation, applying other independent dating methods using the remaining parallel ice sections from this depth interval (e.g., from DEN-13B) could be used, and additional geophysical and modeling approaches are needed to rigorously test this hypothesis.

Importantly, some hydroclimate studies do not show a shift toward wetter conditions at 4 kyr BP. On Adak Island, conditions grew cooler and drier at 4.5 kyr BP, which is consistent with the prevailing interpretation of a stronger Aleutian Low advecting warmer, moister air into the Gulf of Alaska and cooler, drier air to the western Aleutians (Bailey et al., 2018). Certain sites located to the north of the Alaska or St. Elias ranges also show a drying trend or no major features around 4 kyr BP (Lasher et al., 2021; Finney et al., 2012; King et al., 2022; Chakraborty et al., 2010), emphasizing the idea that orography and rain shadows are critical for controlling the relationship between site precipitation and circulation (Anderson et al., 2016). In fact, Winski et al. (2017) showed that during the instrumental era, Begguya snowfall is highly correlated with precipitation along the Gulf of Alaska but bears little resemblance to nearby precipitation recorded in interior Alaska – a pattern that seems to hold through the mid-Holocene. We note that the Aleutian Low is a wintertime phenomenon such that the role of summertime precipitation may be an important contributor to some of the observed variability among regional paleorecords. Comparing records with different seasonality or with seasonal resolution will be critical in the future given that most of the isotope-based records listed above are dominated by wintertime Aleutian Low dynamics.

5 Conclusions

Although ¹⁴C analysis of ice-incorporated carbonaceous aerosols has allowed radiocarbon dating of various high-elevation ice cores from low and midlatitudes, this technique has not been applied before for high-latitude ice cores because of the generally lower carbon content. The ¹⁴C results from the Denali ice core are the first from a high-latitude ice core. These were achieved by small adaptations in the ice sample preparation procedures for the WIOC ¹⁴C dating method which allowed the processing of larger ice samples up to >1 kg of ice and the application of a new technique for ¹⁴C dating of the DOC fraction, which benefits from higher concentration levels in ice compared to the WIOC fraction (by around a factor of 3). Combining dating by ALC to a depth of 197.2 m (159.2 m w.e.; ∼1674 years BP or 339 CE); volcanic tie points from sulfate, chloride, and conductivity; and the new ¹⁴C-dated horizons, a complete continuous chronology over the entire core was established using a simple inverse ice flow modeling approach. For the overlapping sections, ages based on ALC are confirmed by the agreement with the absolute, radiometric ¹⁴C dates.

The ¹⁴C dating of a sample from just 0.61 m above bedrock at around 209 m depth yielded the first absolute date for near-bedrock ice in the region. Dated to be 7700 to 9000 years old, our result clearly indicates this ice at the very bottom is of early Holocene age. The additional model results indicate a high likelihood of even older ice below (>12 ka). The old ice at the bottom of the Denali core confirms that at least some glacier ice in the central Alaska Range survived the Holocene thermal maximum. Future, independent dating methods would be beneficial to further constrain and improve the timescale presented here. Our results show the applicability and great potential of ¹⁴C dating on low-carbon-content samples from North Pacific and Arctic ice cores. While they indicate the Denali ice core is currently one of the few existing archives in the North Pacific region providing an opportunity to reconstruct Holocene hydroclimate variability, we do expect that similar or even longer paleo-ice-core records can be recovered from North Pacific glaciers if bedrock can be reached.