Introduction
Responses of the cryosphere to climate fluctuations since the last deglaciation provide important analogs to what may occur under future carbon emission scenarios. Of particular importance are the glacial dynamics of under-constrained regions like Alaska, which accounted for 22% of global glacier mass loss over 2000–2023 CE, the largest contributor of the 19 Randolph Glacier Inventory (RGI) 7.0 regions (Zemp et al., Reference Zemp, Jakob, Dussaillant, Nussbaumer, Gourmelen and Dubber2025). Between the last deglaciation and the onset of the Holocene, global surface temperatures increased by ∼6°C (Osman et al., Reference Osman, Tierney, Zhu, Tardif, Hakim, King and Poulsen2021). The Bølling-Allerød warm period (BA; 14.7–12.8 ka) ended with a return to colder conditions during the Younger Dryas (YD; 12.8–11.7 ka) before the onset of Holocene warmth. While temperatures were relatively steadier during the Holocene interglaciation than the latest Pleistocene, climate fluctuations did occur, with some climate reconstructions suggesting an estimated temperature anomaly of ∼2°C for the Holocene (11.7 ka–present) in the extratropical Northern Hemisphere (Kaufman et al., Reference Kaufman, McKay, Routson, Erb, Dätwyler, Sommer, Heiri and Davis2020). The magnitude of these temperature anomalies following the last deglaciation span the range of projected warming by 2100 CE across various emissions scenarios (IPCC, 2023). Evaluating past glacier response to these temperature fluctuations enables a better understanding of the coupled cryosphere/atmosphere/ocean system and may help to more accurately predict sea-level response to the projected future warming.
There is much existing research constraining glacial history in southeast Alaska and surrounding northwestern North America following the last glacial maximum (LGM; locally ∼20–17 ka; Lesnek et al., Reference Lesnek, Briner, Lindqvist, Baichtal and Heaton2018). During the LGM, the Cordilleran Ice Sheet (CIS) spanned from the southern Yukon down to the northwestern corner of the contiguous United States, and extended laterally from the Canadian Rocky Mountains in the east to a margin somewhere on the Pacific Ocean continental shelf, as evidenced through a combination of marine sediment core proxies and cosmogenic 10Be dating of glacial landforms along the outer coast (Barrie and Conway, Reference Barrie and Conway1999; Lesnek et al., Reference Lesnek, Briner, Lindqvist, Baichtal and Heaton2018; Walcott et al., Reference Walcott, Briner, Baichtal, Lesnek and Licciardi2022; Dalton et al., Reference Dalton, Dulfer, Margold, Heyman, Clague, Froese and Gauthier2023; Fig. 1). LGM ice extents for the CIS were time-transgressive, with the Puget Lobe to the south reaching its maximum extent at ∼17 ka (Porter and Swanson, Reference Porter and Swanson1998), while the CIS reached its maximum between Haida Gwaii and Vancouver Island before 18.1 ka (Darvill et al., Reference Darvill, Menounos, Goehring, Lian and Caffee2018). In the Alexander Archipelago, LGM extents were achieved sometime between ∼20 and 17 ka (Lesnek et al., Reference Lesnek, Briner, Lindqvist, Baichtal and Heaton2018).
Figure 1. Representative radiocarbon (Baichtal et al., Reference Baichtal, Lesnek, Carlson, Schmuck, Smith, Landwehr and Briner2021) and cosmogenic nuclide (Lesnek et al., Reference Lesnek, Briner, Baichtal and Lyles2020; Walcott et al., Reference Walcott, Briner, Baichtal, Lesnek and Licciardi2022) deglaciation ages of this sector of the Cordilleran Ice Sheet (CIS) in southeast Alaska.
As temperatures increased during the latest Pleistocene, marine-terminating outlet glaciers of the CIS began to retreat. Stratigraphy of marine sediment cores from Dixon Entrance, located to the south of the Alexander Archipelago (Fig. 1), suggest that an extensive glacier occupying the entrance during the LGM deglaciated between 15.5 and 15.0 cal ka BP (Barrie and Conway [Reference Barrie and Conway1999] reported 14C ages calibrated using Marine20). Farther to the south, marine sediment cores indicate that the continental shelf just west of Queen Charlotte Sound was completely deglaciated by ∼13 cal ka BP (Blaise et al., Reference Blaise, Clague and Mathewes1990). Within the Alexander Archipelago, cosmogenic 10Be ages reveal the retreat of the CIS from the continental shelf and onto the coast of Baranof and Kruzof Islands by 15.1 ± 0.9 ka (Walcott et al., Reference Walcott, Briner, Baichtal, Lesnek and Licciardi2022). Additional cosmogenic 10Be ages from the Alexander Archipelago interior, along with 14C ages from uplifted marine sediments, suggest that ice recession continued inland between ∼15 and 14 ka, followed by subsequent dispersion into independent centers of alpine glaciers that later deglaciated by 11.3 ± 1.1 ka (Lesnek et al., Reference Lesnek, Briner, Baichtal and Lyles2020; Baichtal et al., Reference Baichtal, Lesnek, Carlson, Schmuck, Smith, Landwehr and Briner2021; Fig. 1). By ∼14 ka, the CIS had retreated inland to British Columbia, where centers of alpine glaciation persisted for some millennia (Menounos et al., Reference Menounos, Goehring, Osborn, Margold, Ward, Bond and Clarke2017). Cosmogenic 10Be dating of moraines documents glacial advances during the BA and YD following CIS retreat (Menounos et al., Reference Menounos, Goehring, Osborn, Margold, Ward, Bond and Clarke2017), and glacial erratics in the Telkwa Range, located just east of the Coast Range in central British Columbia, date to 12.4 ± 0.7 ka, recording a persistent ice cap occupying the Coast Range into the YD (Margold et al., Reference Margold, Stroeven, Clague and Heyman2014).
The deglaciation of the CIS in southern Alaska and western Canada has been of particular interest as a means of refuting or validating posited routes of human migration through the region. It has been theorized, largely based on studies documenting plant and animal distributions, that glacial refugia (non-glaciated areas) existed in southeast Alaska throughout the last glacial period, allowing humans to effectively move down the coast into what is now the contiguous United States (e.g., Heaton et al., Reference Heaton, Talbot and Shields1996; Carrara et al., Reference Carrara, Ager and Baichtal2007). While recent work using cosmogenic 10Be dating across the Alexander Archipelago has refuted many of the purported refugia zones (Lesnek et al., Reference Lesnek, Briner, Lindqvist, Baichtal and Heaton2018, Reference Lesnek, Briner, Baichtal and Lyles2020; Walcott et al., Reference Walcott, Briner, Baichtal, Lesnek and Licciardi2022), a route following the coastline was perhaps more feasible after 17 ka (when the CIS would have retreated from its LGM position), and humans may have then migrated along this route (Lesnek et al., Reference Lesnek, Briner, Lindqvist, Baichtal and Heaton2018).
The Holocene extent of Alexander Archipelago glaciers, following dispersal into local centers of glaciation after 11.3 ± 1.1 ka (Lesnek et al., Reference Lesnek, Briner, Lindqvist, Baichtal and Heaton2018), is not well constrained. While there is disagreement in the precise timing, records from surrounding regions suggest the period of maximum Holocene warmth and minimum glacial extents occurred sometime in the Early–Middle Holocene (∼11–5 ka), followed by glacial readvance that culminated in Holocene maximum extents during the Little Ice Age (LIA; 1770–1850 CE) (Mann and Streveler, Reference Mann and Streveler2008; Barclay et al., Reference Barclay, Wiles and Calkin2009; Menounos et al., Reference Menounos, Osborn, Clague and Luckman2009; Kaufman et al., Reference Kaufman, Axford, Henderson, McKay, Oswald, Saenger and Anderson2016). Whether or not Alexander Archipelago glaciers survived the Early-Middle Holocene climate before Late Holocene cooling remains unclear.
Although the Coast Range hosts extensive upland icefields and heavy valley glaciation, glaciers within the archipelago itself today are largely constrained within a collection of small (<1 km2) alpine cirques, the majority of which are located on Baranof Island in the northern part of the archipelago (Molnia, Reference Molnia2007). Remote sensing-based reconstructions of glacial change in southeast Alaska using aerial and satellite imagery reveal substantial and increasing retreat since the twentieth century. In a comprehensive reconstruction of glacial change on the Juneau Ice Field since the LIA, Davies et al. (Reference Davies, McNabb, Bendle, Carrivick, Ely, Holt, Markle, McNeil, Nicholson and Pelto2024) document continuous glacial retreat that accelerated fivefold from 2015 to 2019 CE compared with 1979 to 1990 CE. Bevington and Menounos (Reference Bevington and Menounos2022) used Landsat imagery to quantify glacier area change from 1984 to 2020 CE across western Canada, documenting a nearly sevenfold acceleration in area loss from 1984 to 2010 CE compared with 2011 to 2020 CE. These records demonstrate the acute sensitivity of maritime glaciers in the North Pacific to ongoing climate warming. To date, no local-scale record exists of how the small alpine glaciers currently occupying Baranof Island have responded to climate change in the recent historical past and whether that response has been congruent with the larger ice fields inland of the Alexander Archipelago.
Here, we seek to elucidate the glacial history of a local center of alpine glaciation within the Alexander Archipelago, using Baranof Island as a study site. Specifically, we aim to reconstruct: (1) the timing of Late Pleistocene deglaciation of Baranof Island, (2) glacier changes spanning climate events of the latest Pleistocene, and (3) glacier evolution spanning the Holocene up to the present. We create these reconstructions using cosmogenic 10Be dating of erratic boulders on a transect toward modern glaciers, a continuous record of Holocene glacial fluctuations from a proglacial lake sediment archive, and historical imagery spanning the last century.
Study area
Southeast Alaska is a collection of glacially sculpted and incised islands rising steeply from the coast. We focus here on the glacial history of Baranof Island, a large island in the Alexander Archipelago that contains the heaviest contemporary glaciation within the archipelago. The bedrock of Baranof Island is a heterogeneous accreted-terrane complex, with abundant quartz-rich granitic lithologies (Wilson et al., Reference Wilson, Hults, Mull and Karl2015). The climate is maritime, characterized by wet, mild winters and cool summers, with an average annual temperature of 7.4°C and average precipitation of ∼2200 mm/yr (1981–2010 CE average from https://wrcc.dri.edu/cgi-bin/cliMAIN.pl?ak8494, last accessed 6 January 2025).
In the summer of 2024, we collected samples for cosmogenic 10Be dating from granitic erratic boulders on local topographic highs next to Necker Bay fjord (west-central Baranof Island) and Baranof Lake (east-central Baranof Island) (Fig. 2). Necker Bay cuts into the western flank of Baranof Island, approximately 45 km south of Sitka. This region is dominated by Sitka Graywacke and Paleogene granite, and all 10Be samples were taken from quartz-rich granites/granodiorites. We collected samples from three erratic boulders from each of the two sites within the fjord—one at the mouth where it is abutted by Jamboree Bay and one at the head near Benzeman Lake—and three samples from erratics were also collected from the site on the eastern side of the island near the outflow of Baranof Lake.
Figure 2. Beryllium−10 sample locations on Baranof Island. JB, Jamboree Bay; BL, Benzeman Lake. Three boulders were sampled at each site. One outlier (24SEAK-01; 18 ± 1.2 ka) was not included in the mean site age calculation (shown in bold italics). Elevation data are from the Interferometric Synthetic Aperture Radar (IFSAR) digital terrain model and glacier outlines come from the Randolph Glacier Inventory (7.0) (RGI Consortium, 2023).
Baranof Lake (57.08731°N, 134.87033°W, 44 m above sea level [m asl]) is a glacially fed lake drained by a 0.5-km-long outflow stream into Chatham Strait (Fig. 3). The lake’s 83 km2 catchment currently (2023 CE) has ∼7% glacial coverage, consisting largely of cirque glaciers mantling upper ridges and basins. The lake is 2.9 km2 and has three subbasins (mapped in the field via an echo-sounder), the deepest of which reaches 100 m and is situated on the inflow side. We collected sediment cores from depths ranging from 12 to 25 m from a bench between two subbasins (∼60 and ∼70 m depths) at the outflow side of the lake (Fig. 3).
Figure 3. Baranof Lake study area. (a) Interferometric synthetic aperture radar (IFSAR) digital terrain model of the Baranof Lake watershed (50 m contour interval), with glacier outlines from the Randolph Glacier Inventory (7.0) (RGI Consortium, 2023). (b) Bathymetry of Baranof Lake with inset map showing locations where sediment cores were extracted.
Methods
Cosmogenic 10Be exposure dating
All samples had minimal apparent postglacial weathering, and bedrock surfaces were highly sculpted with sporadic glacial grooving still preserved. Sample elevations in Necker Bay ranged from 77 to 162 m asl, which is above the local marine limit of ∼10 m in western Baranof Island (Baichtal et al., Reference Baichtal, Lesnek, Carlson, Schmuck, Smith, Landwehr and Briner2021). Near the outflow of Baranof Lake, sample elevations ranged from 81 to 105 m asl, which we believe is above the marine limit here based on a highest dated marine limit of ∼36 m on adjacent Admiralty Island (Baichtal et al., Reference Baichtal, Lesnek, Carlson, Schmuck, Smith, Landwehr and Briner2021). We used a hammer and chisel to sample the top few centimeters of each boulder surface. Sample locations were recorded in GAIA GPS (vertical uncertainty ± 5 m).
We processed samples at the Cosmogenic Isotope Laboratory at the University at Buffalo following established quartz cleaning procedures (Corbett et al., Reference Corbett, Bierman and Rood2016). After being purified, samples were dissolved with a known quantity of 9Be carrier (PRIME Lab 2023.10.11 #7/#3; concentration 1033.4 ppm). Isolated beryllium was packed into target cathodes and sent to Lawrence Livermore National Laboratory, where samples were measured via accelerator mass spectrometry against the 07NKSTD standard (ratio 10Be/9Be = 2.85 × 10−12). Measured 10Be/9Be ratios were finalized after subtracting a process blank (1.495 ± 0.28 × 10−15). The AMS error is between 2% and 3% for all samples except 24SEAK-01, which has an error of 6.5%.
We calculated exposure ages using the online CRONUS-Earth exposure age calculator version 3 (https://hess.ess.washington.edu/, last accessed 13 February 2025). We used the Arctic 10Be production rate (Young et al., Reference Young, Schaefer, Briner and Goehring2013) and the time-dependent (Lm) scaling scheme of Lal (Reference Lal1991). Lm scaling was used to maintain consistency with other cosmogenic nuclide studies in the region (e.g., Lesnek et al., Reference Lesnek, Briner, Lindqvist, Baichtal and Heaton2018, Reference Lesnek, Briner, Baichtal and Lyles2020; Young et al., Reference Young, Briner, Schaefer, Zimmerman and Finkel2019; Walcott et al., Reference Walcott, Briner, Baichtal, Lesnek and Licciardi2022).
Lake sediment coring
We used a self-made Nesje-style piston coring system to extract two 5-m-long, 7.62-cm (3-inch)-diameter sediment cores and a Universal Corer from Aquatic Research Instruments to collect three shorter surface cores (Fig. 3). Coring locations were chosen based on a bathymetric map of the lake bottom created from a grid of transects made using a Humminbird Fishfinder echo-sounder and associated software package. The Fishfinder uses software called AutoChart Live to create a bathymetric map in real time while driving repeated transects, including over multiple days. Our survey lines covered the entirety of the lake, except for some gaps shown in the areas in white on Figure 3b. Cores were packaged in the field for transportation and were split at the University at Buffalo. We selected the paired piston (24BNF-3A) and surface (24BNF-3B) cores for detailed analysis and radiocarbon dating, because they contained the most complete sediment sequence among our cores. Split core segments were processed on the Geotek MSCL-XYZ Split-Core Logger at the Continental Scientific Drilling Facility at the University of Minnesota for high-resolution core imaging, magnetic susceptibility, and grayscale reflectance. Loss-on-ignition and dry-bulk density were measured at the University at Buffalo at variable increments, ranging from 2.5 cm in the upper 40 cm of 24BNF-3A and most of 24BNF-3B (5 to 39 cm), 3 cm through the organic-rich portion of 24BNF-3A (241.5 to 296 cm), and 10 cm throughout the rest of 24BNF-3A (296 to 407 cm). The upper 5 cm of 24BNF-3B was sampled for loss-on-ignition at 0.5-cm increments to capture the transition from silt to gyttja, which was not captured in the 24BNF-3A piston core. Wet samples were weighed, dried overnight at 105°C, reweighed to calculate water content, then dried at 550°C for 4 hours, and reweighed again to measure % loss-on-ignition (%LOI). Values of dry-bulk density were calculated by dividing the sediment weight after drying at 105°C by the sediment volume.
Eleven samples of terrestrial macrofossils (wood, unidentified species) were isolated from 24BNF-3A (n = 10) and 24BNF-3B (n = 1) for radiocarbon dating. Each sample was wet-sieved from 1 cm3 of sediment at 850 μm and 125 μm and sent to the ACE Isotope Lab at Northern Arizona University for AMS measurement. We used the Bayesian age–depth modeling R software Bacon (v. 3.2.0) to construct an age–depth model using the radiocarbon ages. A turbidite was manually prescribed from 3.07 to 3.12 m where a sand plug with sharp contacts with adjacent sediment is found near the base of 24BNF-3A. Well-defined, coupled laminations from 3.45 m to the base of 24BNF-3A (4.07 m) were assumed to be annual varves (see “Results” for justification of this assumption) and used to prescribe accumulation rates for the basal portion of the core. Counting and measuring the laminations for accumulation rate calculations were done using the high-resolution images from the Geotek MSCL-XYZ Split-Core Logger at the Continental Scientific Drilling Facility at the University of Minnesota.
Remote sensing of glacial change
We derived historical glacier outlines from 1986 to 2023 CE using images from the Landsat Thematic Mapper (Landsat 4-5; 1984–2013 CE) and Operational Land Imager (Landsat 8-9; 2013 CE–present). Images used had less than 10% land cloud cover and were from the end of the melt season (August/September) to minimize misclassification of clouds or perennial snow cover as glacial ice. The 11 study years chosen reflect the availability of images fitting these criteria. All spatial analyses were completed in QGIS 3.34.13. We calculated the glacier band ratio (GBR) for each Landsat image to emphasize glacier delineations (Hall et al., Reference Hall, Ormsby, Bindschadler and Siddalingaiah1987):
\begin{equation}GBR{\text{ }} = {\text{ }}NIR/SWIR1\end{equation}where both the NIR and SWIR1 bands have a spatial resolution of 30 m. A series of 1:63,360-scale USGS topographic maps covering central Baranof Island based on 1948 CE aerial imagery were used to manually digitize 1948 CE glacial extents. These polygons were cross-referenced with those of the RGI 7.0 (RGI Consortium, 2023), a collection of glacier outlines for approximately the year 2000 CE; any glaciers present in the RGI vector but missing from the 1948 CE vector were appended to the 1948 CE vector under the assumption that glaciers did not grow on Baranof Island between 1948 CE and ca. 2000 CE. We deem this assumption valid under the finding by Molnia (Reference Molnia2007) that all glaciers in the Alexander Archipelago were thinning and retreating throughout the first decade of Landsat collection, as well as the finding by Davies et al. (Reference Davies, McNabb, Bendle, Carrivick, Ely, Holt, Markle, McNeil, Nicholson and Pelto2024) that Juneau Ice Field glaciers have been monotonically retreating since 1675 CE.
We manually digitized glaciers from their GBR rasters. Nunataks were also digitized and subtracted from resulting glacier polygons. We removed any glaciers smaller than 0.05 km2 and calculated glacier area and fragmentation (defined here as the number of distinct glacier polygons) for each sample year.
We used State of Alaska High-Resolution Imagery (2020 CE; https://geoportal.alaska.gov, last accessed 6 May 2025) to manually identify a moraine at the head of the Baranof Lake valley that we infer to have formed during the LIA based on it being sharp crested and mostly unvegetated. The glacial extent demarcated by this moraine is compared with our twentieth- and twenty-first-century glacial extents, temporally extending this remote sensing–based record back to the timing of moraine stabilization.
Results
Cosmogenic 10Be exposure ages
The site adjacent to westernmost Necker Bay fjord (Necker Outer site) is a bedrock knoll with erratic boulders between Necker and Jamboree Bays (Figs. 2 and 4). The three boulders from this knoll date to 18.2 ± 1.2 ka (24SEAK-01), 14.3 ± 0.3 ka (24SEAK-02), and 14.2 ± 0.3 (24SEAK-03) (all ages reported to 1σ; Table 1). The mean age of these three samples is 15.6 ± 2.3 ka (n = 3, 1 SD). However, the high uncertainty in the 10Be/9Be ratio of 24SEAK-01 suggests it is less reliable (Table 1), and in our discussion, we take the mean landform age of this site to be the average of 24SEAK-02 and 24SEAK-03, which is 14.3 ± 0.1 ka (n = 2, 1 SD).
Figure 4. Photos of cosmogenic 10Be sampling sites. Ages reported with 1σ internal uncertainty.
Table 1. Beryllium−10 exposure age data.
a AMS results standardized with 07KNSTD (Ratio 10Be/9Be = 2.85 × 10−12; Nishiizumi et al., Reference Nishiizumi, Imamura, Caffee, Southon, Finkel and McAninch2007) and corrected with process blank (1.495 ± 0.28 × 10−15)
b Beryllium−10 ages reported with 1σ internal uncertainty; external uncertainty in parentheses. Ages are calculated using CRONUS-Earth v. 3 with Arctic production rate (Young et al., Reference Young, Schaefer, Briner and Goehring2013) and Lal (Reference Lal1991) scaling scheme.
Farther inland in Necker Bay (Necker Inner site), we collected three boulder samples from a bedrock ridge between Necker Bay and Benzeman Lake (Figs. 2 and 4). These samples date to 11.1 ± 0.3 ka (24SEAK-05), 10.8 ± 0.3 (24SEAK-06), and 11.7 ± 0.3 ka (24SEAK-07) (Table 1). The mean age of these samples is 11.2 ± 0.5 (n = 3, 1 SD).
At Baranof Lake, we collected three boulder samples from the bedrock ridge separating the outflow end of Baranof Lake from Warm Springs Bay (Figs. 2 and 4). These samples date to 10.4 ± 0.2 (24SEAK-09), 10.0 ± 0.3 (24SEAK−11), and 12.1 ± 0.4 (24SEAK−12). These samples average 10.8 ± 1.1 ka (n = 3, 1 SD).
Lake sediment record
While we collected a total of five lake sediment cores from Baranof Lake (Fig. 3), we focused our investigations on the piston core 24BNF-3A and its paired surface core 24BNF-3B, because they hold the most complete stratigraphic record of the retrieved cores. We first report the results from our sedimentological analysis of the cores, including sediment descriptions and values for magnetic susceptibility, dry-bulk density, and %LOI, and then report radiocarbon ages from extracted macrofossils and the resulting age–depth model produced for these two paired cores.
Sedimentological analysis
We captured 4.07 m of sediment in the piston core 24BNF-3A and 0.39 m in the paired surface core 24BNF-3B (Fig. 5). We did not hit resistance when extracting 24BNF-3A but were forced to stop sampling when the core tube became full. For this reason, we cannot be sure that our basal sediment age marks the age of deglaciation and lake initiation. The upper 2 cm of 24BNF-3B is a brown gyttja, while the rest of the core is light-brown clayey silt. Visual inspection of the core following extraction from the lake showed a clear water column above the sediment, suggesting that the sediment–water interface was captured intact. Additional surface cores that we extruded in the field confirm the uppermost few centimeters are composed of brown gyttja. The clayey silt that dominates the rest of 24BNF-3B also dominates the top ∼2.40 m of 24BNF-3A, below which the sediment gradually transitions to a darker brown gyttja. There is no upper layer of gyttja preserved at the surface of 24BNF-3A, suggesting that the surface sediment–water interface was not captured in the piston core. At 2.95 m, the dark brown gyttja in 24BNF-3A abruptly transitions to a blue-gray silt that continues downward to ∼3.45 m. A 5-cm plug of gray sand is found at 3.07–3.12 m. At ∼3.45 m, the silt begins to gradually transition to visibly distinct coupled tan silt/gray sand laminations. These laminations are millimeter-scale from 3.45 to 3.78 m before becoming centimeter-scale down to the base of the core at 4.07 m.
Figure 5. Photos and simplified schematics of cores 24BNF-3A and 24BNF-3B. Red circles represent terrestrial macrofossil samples that were radiocarbon dated (ka). Ages were calibrated using CALIB 8.2 and are reported in cal ka BP with 2σ uncertainty.
We find the sedimentological variation observed visually in the sediment cores reflected in changes in magnetic susceptibility, dry-bulk density, and %LOI (Fig. 6). The base of 24BNF-3A shows magnetic susceptibility fluctuating from high values below 3.90 m (24–28), to dipping slightly between 3.90 and 3.44 m (averaging 18), before increasing again to its maximum values between 3.44 and 2.95 m. This pattern is mimicked by dry-bulk density values, which also reach their maximum (1.4 g/cm3) between 3.44 and 2.95 m. The lowest %LOI values are found at the base of the core and increase gradually toward the 2.95 m transition from blue-gray silt to gyttja. At this transition, %LOI abruptly increases to its highest value (7.5%), while magnetic susceptibility and dry-bulk density decrease sharply. Both magnetic susceptibility and dry-bulk density reach their lowest values (averaging 13 and 0.8 g/cm3, respectively) in this brown gyttja section. At the transition from gyttja to lighter brown silt around 2.40 m, magnetic susceptibility increases slightly and remains relatively stable throughout the rest of the core, while dry-bulk density continues to increase and %LOI decreases. In the upper ∼2 cm of core 24BNF-3B, where the light brown silt transitions to a less dense, darker brown sediment, %LOI increases abruptly, while magnetic susceptibility decreases.
Figure 6. Magnetic susceptibility, dry-bulk density, and organic matter percentage of core 24BNF-3A. Magnetic susceptibility and organic matter of the corresponding surface core (24BNF-3B) are shown in red dotted lines. Age scale based on age–depth model shown in Figure 7.
Age–depth model
We correlate 24BNF-3B and 24BNF-3A via a distinct peak in the magnetic susceptibility curves of the two cores. This correlation indicates that ∼3 cm of sediment from the top of 24BNF-3A was not captured, and depths in the age–depth model were accordingly modified to create a composite core depth. Our 11 radiocarbon ages range from 0.6 ± 0.1 to 10.6 ± 0.5 cal ka BP (Table 2). There were two stratigraphic reversals in the 14C ages. We constructed an age–depth model with all 11 samples and removed the two samples that forced an abrupt shift in accumulation rate (with no explanatory change in sedimentological characteristics); the age–depth model reported here does not include these two samples (Fig. 7). These anomalously old samples may be from organic material that resided on the landscape or submerged in shallow edge areas of the lake for some centuries before being deposited into the lake basin. We assume that the sediment–water interface in the surface core was intact and assign it an age of −74 cal yr BP. We manually assign a hiatus for the sand section at 3.07–3.12 m, interpreted to be a turbidite with rapid deposition. Because the lowest macrofossil was sampled from 2.95 cm, we have no direct age constraint on the accumulation rate of the lower laminated portion of the core. However, we use the context of our cosmogenic 10Be ages at the distal edge of the lake as a maximum limit to the core bottom age, inferring based on the sediment morphology that this portion of the core (below 2.95 cm) was not deposited subglacially, to extend our age–depth model to the core bottom. Instead of extrapolating the accumulation rate from the adjacent gyttja section into this minerogenic section, which would produce a basal age of >14 ka, we manually prescribe accumulation rates below 2.95 m based on counting laminations and assuming they are varves. While we cannot definitively classify these layers as varves, we justify this assumption through a few rationales: First, the depositional environment is conducive to varve production, having a steep-walled catchment (reducing wind-fetch) and a deeper basin with a relatively smaller surface area (Zolitschka et al., Reference Zolitschka, Francus, Ojala and Schimmelmann2015). Second, the morphology and sedimentology of the layers are consistent with typical varve structure, having a coarser lamina indicative of high-energy deposition capped by a clay lamina indicative of lower-energy deposition (Zolitschka et al., Reference Zolitschka, Francus, Ojala and Schimmelmann2015). Finally, when we assume these layers are varves and count them as such, we produce an extrapolated basal age for the core that agrees within uncertainty with our 10Be age for the outflow ridge. Therefore, while we cannot definitively characterize them as such, we maintain our interpretation of these layers as varves. The total length of the centimeter-scale varve section from 3.78 to 4.07 m was measured and divided by the number of varves, yielding an average accumulation rate of 2.9 mm/yr. We used this same method on millimeter-scale varves present from 3.45 to 3.78 m to calculate an accumulation rate of 1.3 mm/yr for this section. This accumulation rate was also assigned to the blue-gray silt section from 2.95 to 3.78 m, where millimeter-scale laminations are visible but too indistinct to count. This extrapolation is done under the assumption that the accumulation rate of this highly minerogenic blue-gray section is likely closer to the varved portion than to the organic-rich portion just above 2.95 m.
Table 2. Radiocarbon ages from Baranof Lake sediment cores.
a Not included in age–depth model
Figure 7. Age–depth model based on nine terrestrial macrofossil samples from cores 24BNF-3A (n = 8) and 24BNF-3B (n = 1) output from Bacon (v. 3.2.0). Two anomalously old samples were not included in the age–depth model and are plotted as their calibrated ages from CALIB 8.2.
Combining our calculated accumulation rates from the basal sediments with the radiocarbon ages yields our age–depth model, which has a basal age of 10.6 ± 0.3 cal ka BP. We date the transition from silt to gyttja at 2.95 m to 10.4 ± 0.2 cal ka BP, indicating a rapid ∼200 year deposition of the basal minerogenic section. The more gradual transition from the dark brown gyttja to light brown silt at ∼2.40 m dates to 7.9 ± 0.4 cal ka BP, bracketing the organic-rich section of the record to between ∼10.4 and ∼7.9 cal ka BP. A pair of 14C-dated macrofossil samples just above the transition to organic-rich sediments document the significant decrease in sedimentation rate with the transition from glacial sediment to brown gyttja, decreasing from 10 mm/yr in the blue-gray silt portion to 0.4 mm/yr in the gyttja. This rate then decreases to its lowest value (0.2 mm/yr) for the majority of the brown gyttja section. Sedimentation rates across the remainder of the core above the organic-rich sediment section increase and are relatively stable, fluctuating between 0.2 and 0.4 mm/yr.
Remote sensing
After removing images with more than 10% land cloud cover or excess perennial snow cover, we used 11 years of Landsat images between 1986 and 2023 CE to delineate glacial extents in the catchment of Baranof Lake (Fig. 8, Table 3). Minor noise in the trends in both glacier area and fragmentation could arise from difficulty in identifying glacier margins due to cloud-cover and perennial snow with a small study area and limited sample size of available images. Even when manually digitizing glacier outlines, applying the 30 m spatial resolution of Landsat imagery to a small spatial domain leads to some uncertainty. Despite these complications, overall trends in these delineations are valuable to understand historical glacial change within the Baranof Lake watershed. We find that net glacier area loss between 1948 and 2023 CE was 10.6 km2, a 34% reduction, with 9.5 km2 (89%) of this loss occurring between 1986 CE and 2023 CE (Table 3). The rate of glacier area loss between 1948 CE and 1986 CE was slow (−0.03 km2/yr) before increasing by an order of magnitude between 1986 CE and 2023 CE (−0.29 km2/yr). Glacier fragmentation (i.e., the number of separate ice masses) followed the same trend; while fragmentation increased fourfold across the total sampling period (increasing from 4 in 1948 CE to 16 in 2023 CE), no fragmentation occurred between 1948 and 1986 CE (Fig. 9b, Table 3). The moraine mapped at the head of the valley that we attribute to the LIA sits 1.7 km laterally from the 1948 CE extent (along the strike of the valley) and is ∼1000 m lower in elevation (Fig. 8).
Figure 8. Glacial extents in the Baranof Lake catchment since 1948 CE. The 1948 CE extent is based on USGS historical topographic maps (1:63,360) based on aerial photographs; all other years are derived from Landsat imagery. Representative years are displayed for clarity; all sample years are shown in Table 3.
Figure 9. (a) Mean summer (JJA) temperatures for Sitka (red) from the Sitka Rocky Gutierrez Airport Station (https://www.wunderground.com/history/monthly/us/ak/sitka/PASI/date/2014-6, last accessed 28 February 2025). (b) Glacier fragmentation in the Baranof Lake watershed (blue) and glacier disconnections (separation of ablation and accumulation areas) for the Juneau Ice Field (gray). (c) Glacier area for Baranof Lake (black) and the Juneau Ice Field (gray). All Juneau Ice Field data come from Davies et al. (Reference Davies, McNabb, Bendle, Carrivick, Ely, Holt, Markle, McNeil, Nicholson and Pelto2024).
Table 3. Glacier area and fragmentation from 1948 to 2023 CE.a
a The 1948 CE extent is derived from USGS historical topographic maps (1:63,360) based on aerial photographs, all other extents are derived from Landsat imagery.
Discussion
Latest Pleistocene and Early Holocene
At each of our three cosmogenic 10Be dating sites, we interpret the mean age of the erratics to represent the timing of ice retreat from that area of the landscape. The deglaciation age of the mouth of Necker Bay fjord at 14.3 ± 0.1 ka suggests that ice retreated from the mouth of Necker Bay fjord at the same time or just after retreat from nearby Biorka Island to the north (14.8 ± 0.8 ka) and shortly after retreat from sites farther southwest on Baranof Island (15.4 ± 1.1 ka) (Walcott et al., Reference Walcott, Briner, Baichtal, Lesnek and Licciardi2022; Fig. 1). The terminus of the glacier filling Necker Bay continued to occupy the fjord (possibly retreating and advancing) until deglaciating the Benzeman Lake area by 11.2 ± 0.5 ka, roughly coincident with the onset of the Holocene. We find that these pulses of retreat, first across the western coastline of Baranof Island between ∼15 and ∼14 ka and subsequently from the head of Necker Bay at the onset of the Holocene, coincide with pulses of warming atmospheric and sea-surface temperatures (SST) at these times (Fig. 10). Proxy-based global temperature records show an increase of ∼1°C between 15 and 14 ka and an increase of ∼1.5°C from the YD to 11 ka (Shakun et al., Reference Shakun, Clark, He, Marcott, Mix, Liu, Otto-Bliesner, Schmittner and Bard2012) and this pattern is echoed in a regional (53–60°N, 126–144°W) surface air temperature (SAT) record from the global paleotemperature product of Osman et al. (Reference Osman, Tierney, Zhu, Tardif, Hakim, King and Poulsen2021; Fig. 10a and b). Alkenone-inferred mean annual SSTs document warming on the order of 4–5°C in the adjacent Gulf of Alaska at both the onset of the BA and the Holocene (Praetorius et al., Reference Praetorius, Mix, Walczak, Wolhowe, Addison and Prahl2015; Fig. 10d); these are regional temperature swings that could control glacial extent on Baranof Island. These warming events are also evident in a planktonic foraminifera δ18O record from the Gulf of Alaska (Fig. 10c) and have been linked to contemporaneous large-scale warming events in the North Atlantic (recorded by NGRIP δ18O), with signals between these two basins being most highly correlated at the onset of the Bølling (14.7–14.1 ka) and Holocene (11.7–10 ka) (Praetorius and Mix, Reference Praetorius and Mix2014).
Figure 10. Records from lake core 24BNF-3A and 10Be ages plotted against paleoclimate records of the last 15,000 years. All shaded ribbons indicate 1σ uncertainty. Timing of the Younger Dryas (YD; 12.8–11.7 ka) is shown with a gray rectangle. (a) Global surface air temperature (SAT) anomaly from Shakun et al. (Reference Shakun, Clark, He, Marcott, Mix, Liu, Otto-Bliesner, Schmittner and Bard2012) (orange) (relative to 11–6.5 ka) and Marcott et al. (Reference Marcott, Shakun, Clark and Mix2013) (blue) (relative to 1000–1850 CE). (b) Average North Pacific (53–60°N, 126–144°W) SAT derived from the Last Glacial Maximum Reanalysis (Osman et al., Reference Osman, Tierney, Zhu, Tardif, Hakim, King and Poulsen2021). (c) Neogloboquadrina pachyderma sinistral δ18O composite for the Gulf of Alaska from cores EW0408-26JC, EW0408-66JC, and EW0408-85JC from Praetorius and Mix (Reference Praetorius and Mix2014) (pink, y-axis reversed). (d) Alkenone (UK37’)-inferred sea-surface temperature (SST) of the Gulf of Alaska from Praetorius et al. (Reference Praetorius, Condron, Mix, Walczak, McKay and Du2020) (purple). (e and f) Loss-on-ignition percentage (black) and magnetic susceptibility (aqua, y-axis reversed) of piston core 24BNF-3A and surface core 24BNF-3B (red dotted). (g) Constructed conceptual time–distance record for the Alexander Archipelago based on all other records in this study. Blue circles represent 10Be sample sites, with 1σ error bars.
At the time of these SST increases, sea-level was rising concurrently with the collapse of the forebulge that propped up the outer Alexander Archipelago (Baichtal et al., Reference Baichtal, Lesnek, Carlson, Schmuck, Smith, Landwehr and Briner2021). According to the relative sea-level compilation from Baichtal et al. (Reference Baichtal, Lesnek, Carlson, Schmuck, Smith, Landwehr and Briner2021), Necker Bay fjord would have been isolated above sea level until flooding by the ocean at ∼10 ka. That high SST coinciding with timing of ice retreat on Baranof Island suggests that regional atmospheric temperature may have echoed the SST pattern. This suggested covariance between ocean and air temperatures is supported in the pollen records of this region, with a record from Hummingbird Lake in southwestern Baranof Island (Fig. 1) showing some evidence for vegetation response to cooling during the YD via an increase in the cold-tolerant Alnus sp. and a decrease in Pinus contorta, which prefers the warmer lowlands (Ager, Reference Ager2019). Additionally, a lake sediment record from Pleasant Island (just north of Baranof Island; Fig. 1) supports an ecosystem response to the YD, with increased sedimentation and higher counts of tundra-plant pollen grains from 12.6 to 11.3 cal ka BP (Hansen and Engstrom, Reference Hansen and Engstrom1996). During the YD, decreased air and ocean temperatures may have caused local glaciers to come to a standstill (or readvance) in central Baranof Island, supporting an ice cap with an outlet glacier that extended down Necker Bay fjord before the Holocene warming initiated its retreat inland. Our dates from Necker Bay fjord do not preclude or affirm a YD advance; we interpret them to confirm a glacier terminated in Necker Bay fjord between 14.3 ± 0.1 ka and 11.2 ± 0.5 ka, with the possibility of this glacier advancing to terminate farther down-fjord in response to the YD cooling event.
The post-YD local glacial retreat found at the head of Necker Bay fjord is also documented on the eastern side of Baranof Island by the deglaciation of the Baranof Lake outflow at 10.8 ± 1.1 ka. These ages from both the head of Necker Bay fjord and Baranof Lake together suggest that post-YD glacial retreat toward the central spine of the island was generally coeval across mid-Baranof Island. Collectively, the ages presented here demonstrate that while the CIS had retreated into British Columbia ∼14 ka (Barrie and Conway, Reference Barrie and Conway1999; Menounos et al., Reference Menounos, Goehring, Osborn, Margold, Ward, Bond and Clarke2017; Darvill et al., Reference Darvill, Menounos, Goehring, Lian and Caffee2018)—by which time the inland sounds of the Alexander Archipelago had opened—local centers of glaciation persisted on Baranof Island (and possibly in the mountainous areas of other Alexander Archipelago islands; cf. Lesnek et al., Reference Lesnek, Briner, Baichtal and Lyles2020) for several millennia until the Holocene. This inference is supported by the relative presence and absence of tephra from the latest Pleistocene Mt. Edgecumbe eruptive sequence recorded by Riehle et al. (Reference Riehle, Mann, Peteet, Engstrom, Brew and Meyer1992) across Baranof Island and the surrounding area; the timing of this sequence was later dated to between 14.6 and 13.1 ka (Praetorius et al., Reference Praetorius, Mix, Jensen, Froese, Milne, Wolhowe, Addison and Prahl2016), effectively providing a minimum age for ice retreat from these areas where tephra is present. We use these tephra locations, available 10Be ages (Walcott et al., Reference Walcott, Briner, Baichtal, Lesnek and Licciardi2022; this study), and topographic constraints to create a conceptual depiction of an ice cap during the YD in central Baranof Island (Fig. 11).
Figure 11. Conceptual depiction of an extent of a local Baranof Island ice cap during the Younger Dryas (YD). Presence/absence locations of tephra are based on the map from Riehle et al. (Reference Riehle, Mann, Peteet, Engstrom, Brew and Meyer1992). Kruzof Island, Biorka Island, and southwestern Baranof Island 10Be ages come from Walcott et al. (Reference Walcott, Briner, Baichtal, Lesnek and Licciardi2022). Minimum tephra age of 13.1 ka is constrained by Praetorius et al. (Reference Praetorius, Mix, Jensen, Froese, Milne, Wolhowe, Addison and Prahl2016).
The Baranof Island ice cap may be representative of other islands in the Alexander Archipelago during the YD chronozone following retreat of the CIS. Beryllium−10 ages across the southern portion of the archipelago provide support for persistent local ice, with evidence for Early Holocene glacial retreat on southern Kuiu, northern Prince of Wales, and Zarembo Islands (Lesnek et al., Reference Lesnek, Briner, Baichtal and Lyles2020; Fig. 1). The ages reported by Lesnek et al. (Reference Lesnek, Briner, Baichtal and Lyles2020) include a subset with a mean of 11.5 ± 0.9 ka, corresponding to the onset of the Holocene and perhaps reflecting the timing of final deglaciation of local ice. However, each of these sites also had samples with significantly older ages (∼17–15 ka), and multiple ages from single sites are complicated to interpret.
Holocene glaciation
Our sediment cores from Baranof Lake can be used to continue our reconstruction of Baranof Island glaciation into the Holocene. Although core 24BNF-3A did not reach bedrock or a till deposit at its base, and we therefore cannot assign our lowest age of 10.6 ± 0.3 cal ka BP as the timing of lake initiation, the age does closely align with the forefield deglaciation of 10.8 ± 1.1 ka based on 10Be ages of erratic boulders. Throughout this discussion, we interpret our proxies of magnetic susceptibility, dry-bulk density, and %LOI to be indicators of glacial extent in the catchment. While it is important to note that these are proxies, not directly measuring glacial extent and susceptible to non-glacial catchment changes that could influence their values through our record, we feel that using these three proxies in tandem and evaluating the record they produce within the context of paleoclimate and glacial records from the surrounding region bolsters our interpretation. We interpret the highest magnetic susceptibility, lowest %LOI section composed of distinctly laminated couplets and blue-gray silt (indicative of a high sedimentation rate environment) at the bottom of 24BNF-3A (2.95–4.07 m; Fig. 5), all of which are signals correlated with increased glacial activity, to represent the timing of farthest downvalley glacial extent represented in the sediment core. We infer that a glacier may have terminated within Baranof Lake at 10.6–10.4 ka. The laminations that we interpret to be varves in this section of the core decrease in thickness from centimeter-scale at 4.07–3.75 m to millimeter-scale at 3.75–3.46 m, before becoming indistinct. This decrease in thickness through time supports that a nearby glacier margin (terminating within the upper reaches of the lake) was retreating during this portion of the record. The transition at 10.4 ka from this ice-proximal sediment to more organic-rich sediments suggests that glaciers in the catchment retreated out of the lake and farther upvalley rapidly at this time, leading to a sedimentation rate decrease of an order of magnitude. This transition may be the response of the ice cap to elevated global air temperatures in the earliest Holocene (Praetorius and Mix, Reference Praetorius and Mix2014; Ager, Reference Ager2019), with the central valley glacier quickly disintegrating into upland cirques and maybe even disappearing altogether at this time.
The Early Holocene period of reduced (or absent) glacial activity spans from ∼10.4 to ∼7.9 ka, marked by a period with relatively high %LOI and low magnetic susceptibility and dry-bulk density (Fig. 6). If we do, in fact, take the moraine mapped at the head of lake valley to be the LIA moraine, the much lower (higher) magnetic susceptibility and dry-bulk density (%LOI) from ∼10.4 to ∼7.9 ka compared with the LIA suggests that ice extents during this time were retracted relative to this moraine (Fig. 3). This presumed LIA extent therefore provides a maximum constraint on glacial extents during the Early Holocene and documents a significant glacial retreat at 10.6 ka (Fig. 3). Indeed, that %LOI values are higher from ∼10.4 to ∼8 ka than the postindustrial record (shown in 24BNF-3B) suggests that glaciers may have been smaller during the Early Holocene than they are today (Fig. 10e). It is important to again note, however, that this moraine has not been dated, and we are simply speculating that it is the LIA moraine, with this speculation supported by our lake sediment proxies and records from the surrounding region showing Holocene maxima being reached during the LIA (e.g., Barclay et al., Reference Barclay, Wiles and Calkin2009; Menounos et al., Reference Menounos, Osborn, Clague and Luckman2009; Hawkins et al., Reference Hawkins, Menounos, Goehring, Osborn, Pelto, Darvill and Schaefer2023, Reference Hawkins, Goehring and Menounos2025).
Records from western Canada and elsewhere in southern Alaska document warm conditions and restricted glacial extents during the Early Holocene. Remnants of subfossil trees situated above the modern timberline in the Canadian Rockies suggest that temperatures were higher than today between 9.1 and 6.7 ka (Luckman and Kearney, Reference Luckman and Kearney1986; Menounos et al., Reference Menounos, Osborn, Clague and Luckman2009), largely coincident with the global Holocene Thermal Maximum (Fig. 10a), and Early Holocene wood emerging from currently retreating glaciers in the southern Coast Mountains and Athabasca glacier in the Canadian Rockies suggests relatively restricted Early Holocene glacial extents (Menounos et al., Reference Menounos, Osborn, Clague and Luckman2009). Sediment cores from Kodiak Island in southwest Alaska and Pleasant Island just north of Baranof Island additionally support warmer conditions during the Early Holocene (Hansen and Engstrom, Reference Hansen and Engstrom1996; Peteet et al., Reference Peteet, Nichols and Mann2019). The timing of warmth and reduced glacial extents these records document during the Early Holocene largely coincides with our record from Baranof Lake, suggesting a regional glacial minimum ∼10–8 ka.
We attribute the transition from the organic-rich sediments into a lighter brown silt at 2.40 m accompanied by increasing magnetic susceptibility and decreasing %LOI at ∼8 ka (Figs. 5 and 6) to glaciers growing within the catchment beginning around this time. This timing of glacial growth is in alignment with evidence for readvance across western Canada from 8.6 to 8.1 ka, which has been suggested to be in response to the 8.2 ka cooling event (Menounos et al., Reference Menounos, Koch, Osborn, Clague and Mazzucchi2004; Rasmussen et al., Reference Rasmussen, Andersen, Svensson, Steffensen, Vinther, Clausen and Siggaard-Andersen2006). Paleoclimatic and paleoenvironmental records of the 8.2 ka event in Beringia are scarce, however, and Kaufman et al. (Reference Kaufman, Axford, Henderson, McKay, Oswald, Saenger and Anderson2016) suggest that more finely resolved records from this region are needed to determine whether this North Atlantic climate event was expressed in regions of the North Pacific.
We infer continued glacial expansion throughout the remainder of the Holocene up until the last century from decreasing %LOI and increasing dry-bulk density over this period. In a compilation of Holocene climate proxy records across eastern Beringia, Kaufman et al. (Reference Kaufman, Axford, Henderson, McKay, Oswald, Saenger and Anderson2016) found expression of maximum Holocene temperatures between 7 and 5 ka across all proxy types, with glaciers inferred to be at their Holocene minima during this period. In contrast, our data suggest that Baranof Island glaciers began to expand by 7 ka and continued expanding through the Middle Holocene. An earlier Holocene thermal maximum on Baranof Island, as indicated by high %LOI and low magnetic susceptibility, which in turn reflect glacier length, is echoed, however, in the timing of peatland initiation and lake formation in eastern Beringia, with the highest rates of initiation occurring from 11 to 10 ka before decreasing and remaining low through the remainder of the Holocene (Kaufman et al., Reference Kaufman, Axford, Henderson, McKay, Oswald, Saenger and Anderson2016). Additionally, 14C-dated glacially sheared stumps from the southern Coast Mountains, BC, document a period of glacial expansion between 7.0 and 5.6 cal ka BP (Menounos et al., Reference Menounos, Osborn, Clague and Luckman2009), overlapping with our record of advancing glaciers through the Middle Holocene. A lake sediment record from Hallet and Greyling Lakes in the Chugach Mountains of southern Alaska documents an onset of Holocene warmth at 10 ka, largely in agreement with our record; however, their archive documents this warm period continuing until 6 ka based on %LOI and dry-bulk density values (McKay and Kaufman, Reference McKay and Kaufman2009). After 6 ka, the %LOI trends from these lakes closely align with ours, with %LOI steadily and gradually decreasing throughout the remainder of the Holocene.
Records across southern Alaska document widespread glacial growth by ∼4 ka and continuing until the LIA (Mann and Streveler, Reference Mann and Streveler2008; Barclay et al., Reference Barclay, Wiles and Calkin2009). Glacial advances are documented by detrital wood and glacially sheared stumps across parts of southern Alaska and the Canadian Rockies at 3.3–2.9 cal ka BP, 2.2–2.1 cal ka BP, and 2.0 cal ka BP, followed by so-called first millennium advances ∼1.2 ka (Barclay et al., Reference Barclay, Wiles and Calkin2009). Three main phases of LIA advance are documented by dendrochronology in the Kenai, Chugach, and Coast Mountains during the 1180s–1320s, 1540s–1710s, and 1810s–1880s CE, with nearly all records documenting Holocene maxima during the latest LIA phase (Barclay et al., Reference Barclay, Wiles and Calkin2009). Cosmogenic nuclide dating of a moraine in the forefield of Vintage Peak glacier in the southern Coast Mountains, BC, dates the Late Holocene maximum extent to 1090–1530 CE (Hawkins et al., Reference Hawkins, Goehring and Menounos2025). The sequential advance of glaciers throughout the Late Holocene culminating in a late LIA maximum agrees with our record showing increasing magnetic susceptibility and dry-bulk density and decreasing %LOI culminating during the LIA (Fig. 6).
Historical glacial change
Most of Alaska’s glaciers are temperate alpine glaciers, a class of glaciers that responds quickly and acutely to climate fluctuations (Molnia, Reference Molnia2007). We interpret the abrupt increase in %LOI paired with a decrease in magnetic susceptibility in the upper 1.5 cm of core 24BNF-3B, which corresponds to the last ∼30 years according to our age–depth model, as the retreat of glaciers in the Baranof Lake catchment in response to anthropogenic climate warming (Fig. 6). This retreat inferred from lake sediments is echoed in remotely sensed imagery, with a 1.7 km lateral retreat of the central valley glacier from its purported LIA maximum to the 1948 CE ice margin and a 34% decline in the area of glaciers in and around the Baranof Lake watershed between 1948 and 2023 CE. It is important to note that the change in trend occurring coincident to the change in image source type (from aerial imagery in 1948 CE to Landsat imagery in 1986 CE) raises the possibility of this change being a methods artifact rather than a real signal, with the lower spatial resolution of the Landsat images potentially introducing greater misattribution of glacial pixels. To quantify the error in the Landsat extent for 1986 CE (the first year that Landsat imagery was used), we assume a systematic overcounting/undercounting; that is, we assume that all partially glaciated pixels were counted as fully glaciated (upper limit) or non-glaciated (lower limit). This calculation yields an area estimate of 29.838 ± 1.597 km2 for 1986 CE. The calculated area for 1948 CE (30.991 km2) is within these error margins, making it possible that no change in glacial extent occurred in the catchment from 1948 to 1986 CE. Additionally, these error margins suggest the possibility of 0.4 km2 glacial growth during this interval; however, we deem this possibility unlikely, given the climatic context and the trend of monotonic retreat observed in the nearby Juneau Ice Field over this interval (Davies et al., Reference Davies, McNabb, Bendle, Carrivick, Ely, Holt, Markle, McNeil, Nicholson and Pelto2024).
The increase in the rate of Baranof glacier area loss and fragmentation starting between 1980 and 1990 CE as recorded in Landsat imagery is coincident with increasing mean summer (JJA) temperatures for Sitka (measured by the Sitka Rocky Gutierrez Weather Station) during the same decade (Fig. 9a). Mean summer temperature in Sitka from 1949 to 1991 CE was 11.7 ± 1.2°C, increasing to 12.9 ± 1.0°C from 1991 to 2024 CE. From 1991 to 2023 CE, glacier area and summer temperature are negatively correlated (ρ= –0.68, P = 0.02, n = 8) while fragmentation and summer temperature are positively correlated (ρ = 0.83, P < 0.001, n = 8).
These trends in Baranof Island glacier loss are reflected elsewhere in southeast Alaska and western Canada. Extensive retreat from the Late Holocene maximum ice margin is echoed in the findings of Hawkins et al. (Reference Hawkins, Goehring and Menounos2025) for Vintage Peak Glacier in the southern Coast Mountains, BC, where glacier area decreased by 95% between 1090–1530 CE and 2020 CE. Additionally, Hawkins et al. (Reference Hawkins, Menounos, Goehring, Osborn, Pelto, Darvill and Schaefer2023) find that in the Mackenzie and Selwyn mountain ranges to the northwest of our site, along the border of the Yukon and Northwest Territories, Late Holocene moraines sat 0.7 to 2 km downvalley of 2020 CE ice extents. Glacier area loss since the LIA for the nearby Juneau Ice Field (see Fig. 1) was quantified by Davies et al. (Reference Davies, McNabb, Bendle, Carrivick, Ely, Holt, Markle, McNeil, Nicholson and Pelto2024) using a combination of geomorphological mapping, aerial photographs, and satellite imagery. This record documents trends similar to our findings in the Baranof Lake watershed, with glacier area loss increasing from 0.18%/yr in 1979–1990 CE to 0.39%/yr in 1990–2005 CE. They also observe a marked acceleration in area loss in 2015 CE, reaching its highest values of the record at 0.96%/yr. Davies et al. (Reference Davies, McNabb, Bendle, Carrivick, Ely, Holt, Markle, McNeil, Nicholson and Pelto2024) additionally measured “disconnections”—defined as regions where a glacier’s ablation and accumulation zones become separated—which align closely with our glacier fragmentation record on Baranof Island (Fig. 9b). While we are limited in delineating these more acute trends in the Baranof Lake watershed due to the small total area of glacial coverage and the limited availability of high-resolution, cloudless imagery, we still observe strong correlation between our record and that of the Juneau Ice Field in both area loss (ρ = 0.92, P < 0.001, n = 11) and disconnection/fragmentation (ρ = 0.89, P < 0.001, n = 11) after 1986 CE. In western Canada, Bevington and Menounos (Reference Bevington and Menounos2022) used Landsat imagery to derive glacial changes from 1984 to 2020 CE and found that rate of area loss increased from −78.7 km2/yr between 1984 and 2011 CE to −481 km2/yr from 2011 to 2020 CE, showing largely the same trend in retreat rate documented in southeast Alaska.
Glaciers elsewhere across Alaska echoed this rapid twentieth- and twenty-first-century retreat. In a systematic analysis of glacial change since the LIA across all the Alaskan mountain ranges, Molnia (Reference Molnia2007) found that 98% of Alaskan glaciers are thinning or retreating, with the small portion of glaciers not currently undergoing retreat predominantly being current or former tidewater glaciers, which are well-known to follow the climatically decoupled “tidewater glacier cycle” (Brinkerhoff et al., Reference Brinkerhoff, Truffer and Aschwanden2017). Global glacier mass loss inventories since 2000 CE show that Alaska is the largest contributor out of the 19 RGI regions, estimated to account for 22% of the global glacier mass loss of 273 ± 16 Gt/yr from 2000 to 2023 CE (Zemp et al., Reference Zemp, Jakob, Dussaillant, Nussbaumer, Gourmelen and Dubber2025).
Continued global glacial recession through the next century is projected under all emissions scenarios (Rounce et al., Reference Rounce, Hock, Maussion, Hugonnet, Kochtitzky, Huss and Berthier2023). Extrapolation into the future of the 1990–2023 CE rate of glacier area loss in our field area calculated here projects a complete disappearance of glaciers from the Baranof Lake watershed (and, thus, likely the island as a whole) within this century. This projection is echoed in the Mackenzie and Selwyn mountain region, where the Open Global Glacier Model predicts 85–97% glacier volume loss by the end of the century across various emissions scenarios (Hawkins et al., Reference Hawkins, Menounos, Goehring, Osborn, Pelto, Darvill and Schaefer2023). The modern historical retreat of Beringian glaciers has already contributed significantly to eustatic sea-level rise, with southeast Alaska and British Columbia alone contributing 0.04 ± 0.01 mm/yr to sea-level rise from 1948 to 2000 CE (Larsen et al., Reference Larsen, Motyka, Arendt, Echelmeyer and Geissler2007), and Alaska as a whole contributing 0.12 ± 0.02 mm/yr from 1962 to 2006 CE (Berthier et al., Reference Berthier, Schiefer, Clarke, Menounos and Rémy2010). With anticipated continued retreat, 17 ± 4 mm of sea-level rise contribution is expected from Alaskan glaciers under emission scenario RCP 4.5 (Hock et al., Reference Hock, Rasul, Adler, Cáceres, Gruber, Hirabayashi, Jackson, Pörtner, Roberts, Masson-Delmotte, Zhai, Tignor, Poloczanska and Mintenbeck2019).
Conclusion
Here, we provide for the first time a post-deglaciation record of local ice cover in southeast Alaska, helping to minimize the spatial data gap in glacial history records for this region. We conclude that:
• A local ice cap persisted on Baranof Island for millennia following the retreat of the CIS, with pulses of retreat (∼15−14 ka and ∼11 ka) occurring concurrently with increasing sea-surface and atmospheric temperatures.
• Baranof Island glaciers were likely smaller than their modern extents during the Early Holocene (from ∼10.4 to ∼8 ka).
• In the historical record, the rate of glacier mass loss increased by an order of magnitude after 1986 CE, echoing a pattern documented elsewhere in Alaska and western Canada.
• Extrapolation of the rate of glacier area loss since the 1980s suggests that glaciers could be entirely absent from Baranof Island within this century, returning to conditions last seen on Baranof Island ∼8 ka.
Future work is needed to more precisely delineate the proposed Baranof Island ice cap boundary and its history between ∼15 ka and the Early Holocene. For example, it remains unclear whether glaciers retreated when the CIS collapsed during the Bølling and then readvanced during the colder Allerød and/or YD, or if they persisted throughout this interval without large changes in extent. Further work using glacial extent indicators such as %LOI, dry-bulk density, and grain-size analysis on sediment cores from proglacial lakes that became uncovered during the Bølling could help clarify this portion of the record. The 10Be ages and topography of Baranof Island suggest that the southern tip or outer western coast of the island may hold such archives. Deglaciation from the southwestern coast of Baranof Island by ∼15 ka as well as the lower elevations in the southern portions of the island support the existence for more temporally extensive lake sediment records in this region (Fig. 11).
Acknowledgments
We thank Katherine Prussian, Kuldeep More, and Michaela Richter for field assistance; Maya Lasker for lab assistance; and the United States Forest Service for logistical assistance. This work was done on the ancestral homelands of the Lingít tribe. Gunalchéesh to the Sitka Tribe of Alaska and the broader community of Sitka. This work was supported by NSF-2221988. Finally, we thank the associate editor and two anonymous reviewers for help improving this article.
Competing Interests
The authors declare that they have no conflicts of interest.
Open access
