Early Holocene ice retreat from Isle Royale in the Laurentian Great Lakes constrained with ¹⁰Be exposure-age dating

The timing of the Laurentide Ice Sheet's final retreat from North America's Laurentian Great Lakes is relevant to understanding regional meltwater routing, changing proglacial lake levels, and lake-bottom stratigraphy following the Last Glacial Maximum. Recessional moraines on Isle Royale, the largest island in Lake Superior, have been mapped but not directly dated. Here, we use the mean of 10 new ¹⁰Be exposure ages of glacial erratics from two recessional moraines (10.1 ± 1.1 ka, one standard deviation; excluding one anomalously young sample) to constrain the timing of Isle Royale's final deglaciation. This ¹⁰Be age is consistent with existing minimum-limiting ¹⁴C ages of basal organic sediment from two inland lakes on Isle Royale, a sediment core in Lake Superior southwest of the island, and an estimated deglaciation age of the younger of two subaqueous moraines between Isle Royale and Michigan's Keweenaw Peninsula. Relationships between Isle Royale's landform ages and Lake Superior bottom stratigraphy allow us to delineate the retreat of the Laurentide ice margin across and through Lake Superior in the early Holocene. We suggest that Laurentide ice was in contact with the southern shorelines of Lake Superior later than previously thought.

1 Introduction

Following the Last Glacial Maximum (LGM), the southwest margins of the Labrador Dome of the Laurentide Ice Sheet (LIS) retreated, exposing the Laurentian Great Lakes of North America (hereafter referred to as “the Great Lakes”; Fig. 1). Initial retreat of the LIS from LGM margins in the Great Lakes region is associated with increased summer boreal insolation (Clark et al., 2009; Ullman et al., 2015), and much of meltwater from the melting LIS drained through the Great Lakes throughout the late Pleistocene and early Holocene (Carlson et al., 2007; Fisher, 2020; Teller and Thorleifson, 1983; Teller et al., 1983). Overall, meltwater from the LIS accounts for ∼ 50 %–60 % of the late-Pleistocene sea-level budget (Clark and Mix, 2002). Constraining spatial and temporal retreat patterns of the LIS retreat, including along its southern margins through the Great Lakes, improves our understanding of meltwater routing from proglacial lakes to the oceans and proglacial lake organization (Breckenridge, 2013; Broecker, 2006; Broecker et al., 1989; Carlson et al., 2007; Farrand, 1969; Fisher, 2020; Fisher and Breckenridge, 2022; Leydet et al., 2018; Teller, 1990; Teller and Mahnic, 1988; Teller et al., 2002, 2005).

Figure 1 The Laurentian Great Lakes of North America and the Laurentide Ice Sheet (LIS) margins at the Last Glacial Maximum (LGM; solid blue line; Ehlers et al., 2011). Inset shows full extent of the North American LGM ice cover inclusive of the Cordilleran (CIS), Greenland (GIS), and Innuitian (IIS) ice sheets (Ehlers et al., 2011), as well as the Labrador Dome (LD) where ice that forms the LIS's southern margins originates. Existing sample sites are shown for ¹⁴C data (black dots; Dalton et al., 2020), ¹⁰Be data (yellow triangles; Ceperley et al., 2019; Kelly et al., 2016; Leydet et al., 2018; Lowell et al., 2021; Ullman et al., 2015), and optically simulated luminescence (magenta squares; Fisher et al., 2020; Schaetzl et al., 2017). Only US states and Canadian provinces covered by the LIS are shown. The state of Michigan comprises the Lower Peninsula (LP) and Upper Peninsula (UP). Locations of Isle Royale (IR) and Keweenaw Peninsula (KP) are shown.

Although previous work has investigated the age and rate of ice retreat through the Great Lakes region (Dalton et al., 2020; Dyke, 2004), geochronological constraints on the timing of ice retreat are sparse, especially within lake basins, and Dalton et al. (2020) note that the chronology of Laurentide ice retreat through the Great Lakes has not received as much attention by researchers as other margins of the LIS over the past 20 years. The existing chronology of ice retreat from the Great Lakes is constrained primarily by minimum-limiting ¹⁴C ages, most of which come from the southern LGM margins of the LIS (Dyke, 2004; Dalton et al., 2020; Fig. 1). Additional constraints on the timing of ice retreat from the Great Lakes region include ¹⁰Be exposure-age dating (Fig. 1; Ceperley et al., 2019; Colgan et al., 2002; Leydet et al., 2018; Lowell et al., 2021; Ullman et al., 2015) and in some cases, optically stimulated luminescence burial-age dating (e.g., Fisher et al., 2020; Schaetzl et al., 2017). As a note, all ages in this study with units of cal ka BP are published ¹⁴C calibrated ages and were rounded to the nearest 100 years; none are significantly different than those recalibrated using CALIB rev. 8 (Stuiver and Reimer, 1993).

Lake Superior is the largest and northernmost of the Great Lakes and the largest freshwater lake by surface area in the world (IAGLR, 2012). Despite decades of geochronological efforts mapping ice margin positions across Lake Superior (e.g., Bajc et al., 1997; Black, 1976; Leydet et al., 2018; Lowell et al., 1999, 2005, 2009, 2021; Saarnisto, 1974), the long distances between shorelines and the paucity of dateable organic material retrieved from lake-bottom sediments has hampered development of a detailed history of LIS retreat across the lake (Breckenridge, 2007; Breckenridge et al., 2004; Hyodo and Longstaffe, 2011). Thus, uncertainties about LIS retreat across the lake following the LGM, and hence meltwater routing and proglacial lake organization, remain. Ice retreat from western Lake Superior appears to have been mostly continuous across the late-Pleistocene-to-Holocene transition (Lowell et al., 2021) with some clear exceptions of moderate ice re-advance during the Younger Dryas stadial (Carlson, 2010; Loope, 2006; Lowell et al., 1999, 2009). The timing of this advance, known as the Marquette Readvance, is well constrained at the Lake Gribben forest-bed site in Michigan's Upper Peninsula (UP) to ∼ 11.9 cal ka BP (Fig. 2a; Lowell et al., 1999), and it is known to have produced the Grand Marais I moraine complex across the UP and the Marks moraine in Ontario (Fig. 2; Carlson, 2010; Hobbs and Breckenridge, 2011; Loope, 2006; Lowell et al., 1999, 2009). There is less agreement regarding how far ice advanced in the western Lake Superior basin during the Marquette Readvance (Fig. 2b; Black, 1976; Clayton and Moran, 1982; Colman et al., 2020; Drexler et al., 1983; Farrand and Drexler, 1985; Peterson, 1985). This uncertainty results from little dateable material in Lake Superior sediments, the lack of well-defined ice-marginal deposits associated with the Marquette Readvance around the western Lake Superior basin, and ongoing discussion about how some ¹⁴C ages and the stratigraphy they come from should be interpreted (Colman et al., 2020; Hobbs and Breckenridge, 2011).

Much of the recent geochronological work around Lake Superior focuses on retreat of the LIS's Superior Lobe along the northwest shore (the North Shore; Fig. 2a) leading up to and following the Marquette Readvance, rather than establishing cross-lake relationships between the North Shore and Wisconsin or Michigan's UP (Fig. 2). For example, near-constant retreat of the Superior Lobe along the Lake Superior North Shore uncovered eastward-draining outlets from glacial Lake Agassiz to the Lake Superior basin via Lake Nipigon (Fig. 2; Fisher, 2020; Kelly et al., 2016; Leydet et al., 2018; Lowell et al., 2009, 2021; Teller and Mahnic, 1988). Organic ¹⁴C and ¹⁰Be exposure ages of landforms along the North Shore constrain deglaciation of these outlets from Lake Nipigon to Lake Superior at Black Bay at ∼ 10.7 ka (mean of ¹⁰Be ages; Fig. 2a; Leydet et al., 2018; Lowell et al., 2021; Teller and Mahnic, 1988). Meltwater from proglacial lakes flowed through western Lake Superior to eastern Lake Superior and drained first to Lake Michigan through the Au Train–Whitefish outlets across Michigan's UP and later through the St. Mary's River (Fig. 2; Breckenridge, 2013). However, the path meltwater took through Lake Superior during this drainage is unclear because ice margin positions remain less well constrained since there is less recent research and few numerical deglaciation ages across Michigan's UP and none along the Keweenaw Peninsula (Fig. 2a; Drexler, 1981; Huber, 1973; Hughes, 1963).

In the absence of quantitative dating across much of the Lake Superior region, much of our knowledge of the chronology of LIS retreat here comes from analyses of the Lake Superior lake-bottom stratigraphy, which primarily consists of pre-Marquette Readvance outwash or till, then red glaciolacustrine sediment (i.e., varves; Breckenridge, 2007), gray varves, and lastly post-glacial Holocene sediments (Farrand, 1969). The older red varves are sourced from red clay, till, and/or bedrock in the southwest Lake Superior region, compared with the younger gray varves, which are sourced from the bedrock of the Canadian Shield (Breckenridge et al., 2021; Farrand, 1969). Lake Agassiz drainage into the Lake Superior basin is associated with the onset of red varve deposition in areas that were ice-free at the time (Fig. 2b; Breckenridge and Johnson, 2009), most likely when lake levels were between Duluth and Minong stages (Breckenridge et al., 2004; Colman et al., 2020; Teller and Mahnic, 1988). In the western Lake Superior basin, red varves are identified north of Isle Royale, in some inland lakes on Isle Royale's northern margin (Figs. 2b, 3a), and down-ice of the IR2 moraine (Fig. 2b; Breckenridge, 2007; Breckenridge et al., 2004; Colman et al., 2020; Maher, 1977; Raymond et al., 1975; Teller and Mahnic, 1988). Red varves are also observed in parts of the eastern Lake Superior basin (Fig. 2b; Breckenridge et al., 2004; Fisher and Whitman, 1999; Mothersill, 1985), indicating the western and eastern parts of the Lake Superior basin were hydrologically linked. Red varves have not been clearly or consistently reported for areas of the northern Lake Superior basin. Varves are ubiquitous in the northern Lake Superior basin, however, but are gray in color, if any color is reported at all, and are usually associated with ice retreat from the lake basin altogether (Fig. 2b; Breckenridge, 2007; Breckenridge et al., 2004; Colman et al., 2020; Dell, 1973, 1976; Farrand, 1969; Halfman and Johnson, 1984; Hyodo and Longstaffe, 2011; Johnson, 1980; Johnson and Fields, 1984; Kemp et al., 1978; Landmesser et al., 1982; Maher, 1977; Mothersill, 1979, 1985, 1988; Mothersill and Fung, 1972; O'Beirne, 2013; Raymond et al., 1975; Teller and Mahnic, 1988; Thomas and Dell, 1978; Yu et al., 2010). Where red varves are absent in Lake Superior, gray varves directly overlie pre-Marquette Readvance till or outwash, including at a shallow rise just northeast of Isle Royale (gray square in Fig. 2b; Mothersill and Fung, 1982).

Isle Royale is the largest island in Lake Superior (Fig. 2), and it is proximal to, and strikes parallel with, the well-constrained North Shore deglaciation chronology (Fig. 2; Kelly et al., 2016; Leydet et al., 2018; Lowell et al., 2009, 2021). Isle Royale's position within Lake Superior offers the opportunity to extend the limits of the existing deglacial chronology and draw chronological relationships between the north and south shores of the Lake Superior basin and the lake bottom stratigraphy therein. In this study, we present 11 new ¹⁰Be exposure ages of glacial erratics from the crests of recessional moraines on Isle Royale (Huber, 1973), which are broadly constrained by one existing ¹⁴C age from each of two inland lakes on Isle Royale northeast and southwest of these moraines (Flakne, 2003). We place our interpretation of this new ¹⁰Be exposure-age dataset in the context of the existing Lake Superior deglaciation chronology, various Lake Superior bottom stratigraphy records, relationships drawn between regional subaerial and subaqueous glacial landforms, and the organization of proglacial lakes that occupied the western Lake Superior basin following the Marquette Readvance.

Figure 2 Maps showing the Lake Superior region of North America and the position of the Laurentide Ice Sheet at ∼ 10.5 ka (Breckenridge, 2007), the location of new ¹⁰Be data presented in this study (magenta star). Abbreviated locations of the Lake Superior basin mentioned in this study include Thunder Bay (TB), Black Bay (BB), and Nipigon Bay (NB); regional cities include Duluth (D), Houghton (H), Marquette (M), and Thunder Bay (T). Background imagery is 3 arcsec Bathymetry of Lake Superior dataset, which shows elevations in meters above sea level and depths below present-day lake level (NOAA Great Lakes Environmental Research Lab, 1999). Panel (a) shows the Lake Superior basin with locations of moraines known to be associated with the Marquette Readvance (solid black lines). Moraine names are shown with unformatted black text (Isle Royale moraines 1 and 2: IR1 and IR2, respectively; Colman et al., 2020), and ice margin positions thought to be associated with the maximum extent of ice during the Marquette Readvance are shown: dashed purple line (Black, 1976; Clayton and Moran, 1982; Drexler et al., 1983; Hughes and Merry, 1978); stippled yellow line (Peterson, 1985); and cyan dashed–dotted line (Colman et al., 2020). Black dots and yellow triangles show existing ¹⁴C and ¹⁰Be data locations, respectively (Dalton et al., 2020; Leydet et al., 2018; Lowell et al., 2021; Ullman et al., 2015). The Lake Gribben forest-bed site is specifically shown by a white circle with a black X (Lowell et al., 1999). Panel (b) shows the Lake Superior basin with locations of cores where varves (red, gray, or no color reported) have been identified and the LS00-3P sediment core (Breckenridge, 2007). The gray square northeast of Isle Royale is a shallow rise where gray varves sit directly over pre-Marquette Readvance till (Mothersill and Fung, 1972). Red and gray shaded areas south of Isle Royale show the extent of red and gray varve deposition relative to the IR1 and IR2 moraines (Colman et al., 2020).

2 Study location: Isle Royale and the Lake Superior basin

The landscape of Isle Royale is characterized by glacial sediments that blanket the ∼ 1.1 Ga Portage Lake Volcanics and Copper Harbor Conglomerate on the southwest third of the island (Fig. 3a; Davis et al., 2022; Elling et al., 2022; Huber, 1973). Glacial striae are oriented parallel to bedrock ridges and indicate northeast-to-southwest Laurentide ice flow where bedrock is exposed, but striae and crag-and-tail structures (mapped as drumlins) indicate east-to-west ice flow where bedrock is covered by glacial till. Striae and crag-and-tail structures point to, and terminate at, arcuate ice marginal deposits that mark recessional positions of Laurentide ice following the Marquette Readvance (Fig. 3; Huber, 1973).

The topographical expression of ice-marginal landforms (Fig. 4) is an indication they are recessional moraines that were formed subaerially and not reworked by wave processes (Huber, 1973). These moraines have never been directly dated, but they are bounded by minimum-limiting ¹⁴C ages on bulk organic sediments extracted from the bottoms of cores at Lily Lake and Lake Ojibway on Isle Royale's southwest and northeast ends, respectively (Fig. 3a; Flakne, 2003), which record ice-free conditions on Isle Royale by ∼ 10.7 cal ka BP (Lily Lake; median age; 2σ age range of 11.1–10.5 cal ka BP) and ice retreat from Isle Royale by ∼ 9.8 cal ka BP (Lake Ojibway; median age; 2σ age range of 10.1–9.7 cal ka BP). Paleomagnetic reconstructions and varve counting for a sediment core southwest of Isle Royale (LS00-3P; Fig. 2b) indicate this site was ice-free by at least ∼ 10.7 ka (Fig. 2; Breckenridge, 2007). Dated strandlines that were correlated across the western Lake Superior basin, including Isle Royale, constrain the drawdown of glacial Lake Duluth to post-Minong levels to a two-century period between ∼ 10.8 and 10.6 cal ka BP, draining first to the Lake Michigan basin via the Au Train–Whitefish outlets and ultimately to the St. Mary's River (Figs. 2, 3b; median ages as presented in Breckenridge, 2013).

Two prominent subaqueous moraines (IR1 and IR2) span the western Lake Superior basin from Isle Royale to Michigan's Keweenaw Peninsula but have never been associated with any land-based moraines nor have they been dated directly (Figs. 2, 3; Breckenridge, 2013; Colman et al., 2020; Landmesser et al., 1982). Formation of the IR2 moraine is estimated at ∼ 10.5–10.2 ka (Colman et al., 2020) based on the presence of red varves found only on the down-ice side of IR2, the transition from red to gray varves at the LS00-3P core (Fig. 2; Breckenridge, 2007), and correlations to the dated red-to-gray varve transition in cores from the eastern Lake Superior basin (Breckenridge et al., 2004; Hyodo and Longstaffe, 2011).

Figure 3 Panel (a) shows bedrock and surficial geological map of Isle Royale (NPS, 2008), locations of subaqueous moraines IR1 and IR2 (Colman et al., 2020), and inland lakes (Flakne, 2003; Raymond et al., 1975). ¹⁴C ages are median ages of the 2σ calibrated age range. Panel (b) shows Isle Royale's recessional moraines (Huber, 1973) and named strandlines (Breckenridge, 2013). Locations of samples from this study are shown by circles; IR-18-11 is shaded gray because it is anomalously younger than other samples and not included in the mean age of deglaciation we present for Isle Royale. Panel (c) shows the locations of samples relative to positions of Isle Royale's highest moraines and strandlines (see Figs. 4, 5). Bold numbers are sample IDs following the format IR-18-XX. Solid gray lines (W–Z) are topographic profiles in Fig. 4. LiDAR imagery courtesy of Isle Royale National Park.

Figure 4 Elevation along topographic profile lines in Fig. 3c: (a) profile W, (b) profile X, (c) profile Y, and (d) profile Z. Elevation data comes from LiDAR imagery courtesy of Isle Royale National Park. Elevations of mapped strandlines are from Breckenridge (2013) and coloration of strandlines is the same as in Fig. 3; minor strandlines are indicated by gray dashed lines. Vertical exaggeration (VE) on each profile is indicated.

3 Methods

Samples were collected in 2018 from the tops of quartz-bearing glacial erratics situated on the highest-elevation recessional moraines on Isle Royale (Fig. 3; Table 1; Appendix A). Nine samples were collected from the Mt. Desor moraine (IR-18-01 through IR-18-09), one sample was positioned on the small Sugar Mountain moraine (IR-18-10), and another (IR-18-11) was positioned on an upland till plain between the Sugar Mountain and Mt. Desor moraines. (Since there are no known oral traditions or archival records of Indigenous Ojibwemowin names for these moraines, we refer to them by their proximity to local geographic names.) Sampled moraines are relatively higher than the surrounding landscape (Fig. 4) and targeting the highest-elevation moraines on Isle Royale minimizes the potential of prolonged subaqueous erratic emplacement under lowering proglacial lake levels (Breckenridge, 2013). To this end, all samples collected are higher in elevation than Lily Lake, which was isolated from proglacial lakes in the Lake Superior basin by ∼ 10.7 cal ka BP (median age; Flakne, 2003) while lake levels dropped from glacial Lake Duluth to post-Minong levels ∼ 10.8–10.6 cal ka BP (Figs. 3, 4; median ages presented in Breckenridge, 2013). If any of the samples were initially exposed from beneath retreating ice while underwater, they experienced less than ∼ 100 years of subaqueous exposure, which is not long enough to significantly affect ¹⁰Be production and thus exposure ages. The tops of all erratics were >25 cm above the forest floor and 8 of the 11 were >40 cm above the forest floor (Table 1).

Figure 5 Average elevation profile and traces of major strandlines identified on Isle Royale (adapted from Breckenridge, 2013). Drawdown of glacial Lake Duluth to post-Minong stages, isolation of Lily Lake from proglacial lakes (blue square; median age; Flakne, 2003), and exposure of erratics used in this study (circles; gray circle is IR-18-11, which is not included in the mean age of deglaciation presented in this study) occurred within the ∼ 200-year period from ∼ 10.8 to 10.6 cal ka BP (median ages as presented in Breckenridge, 2013).

Table 1 Sample locations and characteristics.

^a Elevations are present day, meters above sea level, and are not corrected for glacial isostatic adjustment. ^b IR-18-11 was collected along the sloped top side of the erratic, which was adjacent to an even larger mafic erratic. In addition to topographic shielding, shielding for IR-18-11 accounts for self-shielding of the erratic's sloped surface (strike 308^∘, dip 12^∘) and local horizon shielding from the adjacent erratic ([170^∘, 0.64], [180^∘, 0.86], [205^∘, 0.64]) using Balco et al.'s (2008) shielding calculator (wrapper 2.0, skyline 2.0).

Quartz was isolated from samples using froth-flotation and acid etching at the Purdue Rare Isotope Measurement (PRIME) Laboratory. Be was isolated and purified from quartz at the National Science Foundation/University of Vermont Community Cosmogenic Facility following established procedures (Corbett et al., 2016). We dissolved between 6.1 and 21.5 g of quartz and isotopically diluted each sample with ∼ 250 µg ⁹Be using an in-house carrier, termed UVM-SPEX, created from dilution of SPEX 10 000 ppm Be standard, with a resulting concentration of 304 µg mL⁻¹ (Table 2). The 11 samples were processed alongside one process blank. ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratios were measured at PRIME in March 2022 and normalized to primary standard 07KNSTD3110, with an assumed ratio of 2.850 × 10⁻¹² (Nishiizumi et al., 2007). We corrected the measured ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratios with the ratio and uncertainty of the one process blank and propagated uncertainties in quadrature (Table 2).

Exposure ages of erratics were determined in February 2023 using Balco et al.'s (2008) exposure-age calculator (version 3, default production rate) based on present-day elevations, and Lifton et al.'s (2014) LSDn ¹⁰Be production rate scaling scheme. Topographic shielding for each sample was calculated in ArcGIS Pro using Li's (2018) point-based shielding tool and LiDAR imagery; the sloped surface of IR-18-11 required self-shielding calculations (see footnote in Table 1). ¹⁰Be exposure ages presented in this study assume no post-glacial erosion. Although a 70-year monthly snowfall record exists for Isle Royale (see Appendix B), its location is not near sampled sites, and it comprises <1 % of the duration of the Holocene, so we do not apply a snow shielding correction on reported ages. Erratic exposure via exhumation through degrading moraine crests is not considered to be significant on Isle Royale because the moraines are small (<30 m in height), younger than the Younger Dryas, and are relatively prominent and geomorphically well defined above the surrounding landscape, all of which greatly reduce the likelihood of long-term erratic exhumation (Putkonen and Swanson, 2003). The normality of the dataset was assessed using a Shapiro–Wilk W test in JMP software (version 16), and we subsequently use Chauvenet's criterion to identify potential age outliers from the dataset (Clark et al., 2009; Rinterknecht et al., 2006).

4 Results

Measured concentrations of ¹⁰Be from the 11 erratics range from (3.98 to 7.35) × 10⁴ atoms g⁻¹ (Table 2), which correspond to exposure ages of 7.0 ± 0.8 to 12.4 ± 1.0 ka (1σ internal uncertainties; Fig. 6; Table 2). Samples from the Mt. Desor moraine (n=9; IR-18-01 through IR-18-09) yield exposure ages ranging from 9.0 ± 0.7 to 12.4 ± 1.0 ka (1σ, internal). There is no relationship between erratic height above the surrounding landscape and exposure age (R²=0.22, p=0.15), and the fact that the lowest-lying erratic (IR-18-02, 25 cm) returns an exposure age that is older than the tallest two erratics (IR-18-07 and IR-18-11, 75 cm) suggests that erratics were not exposed by exhumation through degrading moraine crests. The exposure-age dataset passes the Shapiro–Wilk test for normality (W=0.92; p=0.35), but one sample, IR-18-11, was identified as an outlier using Chauvenet's criterion. IR-18-11 was collected between the Mt. Desor and Sugar Mountain moraines, and we did not observe any evidence in the field of IR-18-11 having been rolled, tipped, buried, exhumed, or exposed due to spalling in the past. The age of IR-18-11 places it out of chronological and geographic sequence between the otherwise similar exposure ages of erratics on the Mt. Desor and Sugar Mountain moraines. As a result, we do not include IR-18-11 in subsequent statistics. After omitting the age of IR-18-11, the mean of 10 samples from the Mt. Desor and Sugar Mountain moraines is 10.1 ± 1.1 ka (1 standard deviation (SD); Fig. 6).

Although all sampled erratics were chosen for the suitability for ¹⁰Be exposure-age dating, and even though exposure ages are normally distributed, we present the mean age with 1 SD uncertainties (∼ 11 % SD) because it is a conservative estimate that reflects the geological uncertainty of our samples and field area. The standard deviation of the mean age of our dataset is larger than the average shift in exposure ages that results from using any of the following: (i) a different exposure-age calculator (i.e., the Ice-TEA exposure-age calculator; 6.7 % shift; Jones et al., 2019), (ii) the Northeast North America (NENA) ¹⁰Be production rate calibration (0.5 % shift; Balco et al., 2009), (iii) a ¹⁰Be production rate scaling scheme other than LSDn (0.4 %–2.2 % shift), (iv) glacial isostatic adjustments (3.8 % shift; Jones et al., 2019), or (v) snow shielding corrections (2.9 % shift). (See Appendix B for details on age comparisons.)

Table 2 ¹⁰Be isotopic data and exposure ages.

^a 9Be was added using an in-house carrier, termed UVM-SPEX, created from dilution of SPEX 10 000 ppm Be standard, with a resulting concentration of 304 µg mL⁻¹. ^b Isotopic analysis was conducted at PRIME Laboratory on 8 March 2022; ratios were normalized against standard 07KNSTD3110 with an assumed ratio of 2.850 × 10⁻¹² (Nishiizumi et al., 2007). ^c Background corrections were made using one measurement blank sample. UVM batch number: 692; AMS cathode number: 162873; AMS ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ blank sample ratio: $\mathrm{4.73}\times {\mathrm{10}}^{-\mathrm{15}}$; ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ blank sample ratio uncertainty: $\mathrm{8.03}\times {\mathrm{10}}^{-\mathrm{16}}$. ^d Exposure ages calculated using Balco et al.'s (2008) exposure age calculator (v. 3, default production rate), Lifton et al.'s (2014) scaling scheme (LSDn), a rock density of 2.7 g cm⁻³, present-day elevations, topographic shielding (Li, 2018), and self-shielding for one sample (IR-18-11; see Table 1). Other production rate calibrations, calculators, production rate scaling schemes, snow shielding, and glacial isostatic adjustment corrections were considered but ultimately not incorporated into the ages presented here (see Appendix B).

Figure 6 Distribution of ¹⁰Be exposure ages from this study (exclusive of IR-18-11; see Results). Ice-free conditions at Lily Lake and Lake Ojibway (Flakne, 2003) are shown by median ¹⁴C ages and their 1σ and 2σ distributions (dark and light blue, respectively; ¹⁴C age 2σ ranges were recalibrated using CALIB rev. 8; Stuiver and Reimer, 1993; median ages are unchanged from original publication). Minimum ice-free conditions southwest of Isle Royale at the LS00-3P core site (magenta triangle) and the age estimates for the subaqueous IR2 moraine (orange rectangle) south of Isle Royale are also shown (Breckenridge, 2007; Colman et al., 2020).

5 Discussion

New ¹⁰Be ages of erratics from Isle Royale's Mt. Desor and Sugar Mountain moraines indicate they share a similar history of exposure at ∼ 10.1 ± 1.1 ka (1 SD). This date provides new insights to the timing of LIS retreat from Isle Royale and the western Lake Superior basin following the Marquette Readvance (Fig. 7a). We interpret this age as the mean ¹⁰Be exposure age of the Mt. Desor and Sugar Mountain moraines, recognizing that exposure ages of individual erratics, or the whole dataset, could be shifted older if factors like snow shielding or exhumation could be reasonably accounted for. Exposure ages could likewise be shifted younger if our measurements include inherited ¹⁰Be, which is possible (Briner et al., 2016).

The ¹⁰Be-constrained timing of LIS retreat from the Mt. Desor and Sugar Mountain moraines is chronologically and geographically consistent with the timeline of presently known deglaciation events on and around Isle Royale. For instance, the stratigraphy of the LS00-3P core southwest of Isle Royale has been interpreted to suggest ice-free conditions by ∼ 10.7 ka (Breckenridge, 2007), which is coeval with the timing of Lily Lake's isolation from proglacial lakes at ∼ 10.7 cal ka BP (minimum-limiting median age; 11.1–10.5 cal ka BP, 2σ age range; Fig. 6; Flakne, 2003). The mean ¹⁰Be age of ice retreat from the Isle Royale's moraines (∼ 10.1 ka) is consistent with the estimated timing of ice retreat from the subaqueous IR2 moraine south of Isle Royale (∼ 10.5–10.2 ka; Colman et al., 2020), which links the island to the Keweenaw Peninsula on the Michigan UP mainland. The earlier timing of ice-free conditions at the LS00-3P core and Lily Lake sites is reasonable given that they are positioned down-ice of the Mt. Desor, Sugar Mountain, and subaqueous IR2 moraines. We take this regional consistency among datasets to be an indication that, despite the 1.1 kyr uncertainty we report, the mean ¹⁰Be age of ice retreat of ∼ 10.1 ka is meaningful within the context of existing datasets and the deglaciation history of the western Lake Superior basin around Isle Royale.

Figure 7 Maps showing inferred Laurentide Ice Sheet margin positions at selected times. Panel (a) indicates ∼ 12.9–11.6 ka, during the Marquette Readvance (Dalton et al., 2020; Leydet et al., 2018; Lowell et al., 1999, 2005, 2009, 2021). Panel (b) indicates ∼ 10.7 ka (Breckenridge, 2007; Flakne, 2003; Leydet et al., 2018; Lowell et al., 2021), when deposition of red varves are first observed across the western Lake Superior basin (Colman et al., 2020; Maher, 1977; Teller and Mahnic, 1988). Panel (c) indicates ∼ 10.1 ka (this study), by which point ice retreats north across the eastern Lake Superior basin (Breckenridge, 2007; Fisher and Whitman, 1999; Saarnisto, 1974), the transition from red to gray varve deposition occurs (Colman et al., 2020), and ice last traverses the western Lake Superior basin linking the North Shore with the Keweenaw Peninsula via Isle Royale (this study; Colman et al., 2020). Ice continues to cover northern Lake Superior basin as indicated by gray varves directly over till (Mothersill and Fung, 1982). Panel (d) indicates ∼ 9 ka, by which point gray varves are ubiquitous across Lake Superior, terminating at ∼ 9 ka (Dell, 1976; Breckenridge 2007; Breckenridge et al., 2004; Colman et al., 2020; Fisher and Whitman, 1999; Halfman and Johnson, 1984; Hyodo and Longstaffe, 2011; Johnson and Fields, 1984; Maher, 1977; O'Beirne, 2013; Raymond et al., 1975; Teller and Mahnic, 1988; Yu et al., 2010). Minimum-limiting ¹⁴C ages (cal ka BP) of Rous Lake and Ring Lake are median ages of organic sediment (Bajc et al., 1997; Saarnisto, 1974). Background imagery is 3 arcsec Bathymetry of Lake Superior dataset, which shows elevations in meters above sea level and depths below present-day lake level (NOAA Great Lakes Environmental Research Lab, 1999). Color shading for the background imagery of all panels is the same as in Fig. 2.

Based on the similarity between the mean ¹⁰Be age of Isle Royale's moraines reported here (∼ 10.1 ka) and the age of formation of the subaqueous IR2 moraine south of the island (∼ 10.5–10.2 ka; Colman et al., 2020), we suggest that the LIS retreat from the western Lake Superior basin south of the island occurred later than it did along Lake Superior's North Shore (∼ 10.7 ka; mean of ¹⁰Be ages from the Nipigon moraine and nearest sites west of Lake Nipigon; Fig. 7; Leydet et al., 2018; Lowell et al., 2021) as these sample sites and ours are likely to share similar snow cover and inheritance histories. Near-constant modeled ice retreat rates of ∼ 40–60 m yr⁻¹ are reported along Lake Superior's North Shore (Lowell et al., 2021); a similar retreat rate of ∼ 45 m yr⁻¹ can be inferred between Lily Lake and Lake Ojibway on Isle Royale (∼ 42 km) from minimum-limiting ¹⁴C bulk-sediment ages (Fig. 3; Flakne, 2003). However, ice retreat south of Isle Royale from the LS00-3P core site (∼ 10.7 ka; Fig. 2; Breckenridge, 2007) to the Mt. Desor and IR2 moraines (∼ 10.1 ka; ∼ 60 km; this study; Colman et al., 2020) indicates a mean ice retreat rate of ∼ 100 m yr⁻¹. Although the ice retreat rate south of Isle Royale appears to outpace the North Shore retreat rate and that between inland lakes on Isle Royale, ice retreat south of the island must have paused long enough to construct the ≤ 20 m high subaerial Mt. Desor moraines (profile X in Fig. 4) and the ∼ 60–100 m high subaqueous IR1 and IR2 moraines during a time when proglacial lake levels were similar to those of glacial Lake Duluth (Colman et al., 2020).

Subaqueous moraines were not identified in Lake Superior north of Isle Royale during multiple seismic surveys (Colman et al., 2020; Johnson, 1980; Landmesser et al., 1982). Neither the ∼ 45 m yr⁻¹ retreat rate between Isle Royale's inland lakes nor the ∼ 100 m yr⁻¹ retreat rate south of Isle Royale we infer based on new ¹⁰Be exposure ages account for periods of standstill during which time these large moraines were constructed. Therefore, the estimates of ice retreat south of Isle Royale should be considered minimum average retreat rates, with rapid retreat before and after periods of standstill, perhaps being facilitated by high rates of iceberg calving over deeper water (Colman et al., 2020; Hanson and Hooke, 2000). This interpretation of ages and retreat rates in the western Lake Superior basin supports the previously proposed hypothesis that the topography of Isle Royale and the Keweenaw Peninsula provided lateral stability to the retreating ice margin in relatively shallow waters, which delayed its overall retreat, relative to ice retreat long the North Shore (Colman et al., 2020).

Early work, based on the relative positions of subglacial, ice marginal, and proglacial landforms on Isle Royale, suggests that moraines on Isle Royale formed in a subaerial environment, likely associated with the Beaver Bay stage of glacial Lake Duluth (Huber, 1973). The mean ¹⁰Be exposure age (∼ 10.1 ka), however, implies that moraine exposure may have been later, postdating the drawdown of glacial Lake Duluth beyond the Beaver Bay stage, at least to post-Minong stages (∼ 10.8–10.6 cal ka BP, median ages; Breckenridge, 2013). In support of this reinterpretation is the observation that strandlines of early proglacial lakes (e.g., Highbridge to Shelter Bay; Figs. 3b, 4) are present at high elevations on the island, but only on northwest-facing slopes, whereas strandlines of younger proglacial lakes (e.g., Huron Mountain and younger; Figs. 3b, 4) are found on northwest- and southeast-facing slopes at lower elevations (Breckenridge, 2013; Farrand, 1969). All striae, crag-and-tail structures, and ice-marginal landforms are found at elevations above the younger strandlines and were thus formed by ice that covered southeast-facing slopes late into the drawdown history of proglacial lakes (Huber, 1973). Similarly, strandlines older than post-Minong stages appear to be absent from northwest-facing slopes on the Keweenaw Peninsula across from Isle Royale (Breckenridge, 2013), though they do exist down-ice of the IR1 moraine (Hughes, 1963).

The apparent absence of strandlines older than post-Minong on southeast-facing slopes on Isle Royale and northwest-facing slopes along the Keweenaw Peninsula further points to long-lasting ice cover between the two landforms until ice retreat in the early Holocene. Post-Minong and younger strandlines are found along Isle Royale, including northeast (up-ice) of mapped moraines. The emergence of Lake Ojibway on Isle Royale's northeast end at ∼ 9.8 cal ka BP (minimum-limiting median age; 10.1–9.7 cal ka BP, 2σ age range, Flakne, 2003) places a minimum age limit on ice retreat from the island's moraines and the subaqueous IR2 moraine. Together, these ages imply that the timing of LIS retreat from Isle Royale's moraines, abandonment of the IR2 moraine, and the presence of post-Minong stages of glacial Lake Duluth are roughly coeval sometime between ∼ 10.6 and 10.1 ka (this study; Breckenridge, 2013; Colman et al., 2020).

Division of the retreating LIS margin into two separate retreating fronts, northwest and southeast of Isle Royale in the early Holocene, has implications for meltwater routing and the Lake Superior lake-bottom stratigraphy. Red varves associated with influx of meltwater from glacial Lake Agassiz (Breckenridge et al., 2004) were identified down-ice of the subaqueous IR2 moraine south of Isle Royale, but not up-ice (Fig. 2b), and Colman et al. (2020) interpreted this to mean gray varve deposition across Lake Superior was associated specifically with ice retreat from the IR2 moraine, not from the Lake Superior basin itself (Farrand, 1969). Following Colman et al.'s (2020) interpretation, the similar ages of Isle Royale's recessional moraines and IR2 implies that the presence and absence of red varves across the Lake Superior basin delineates a regional ice margin when ice was last present on Isle Royale. Previously drawn maps place the LIS margin across the northern Lake Superior basin by ∼ 10.5 ka (Fig. 2; Breckenridge, 2007), but red varves do not appear to be present in Nipigon Bay nor the northern Lake Superior basin, and gray varves are identified directly overlying till (i.e., no red varves) at a shallow rise northeast of Isle Royale (Fig. 2; Mothersill and Fung, 1972).

The absence of red varves in the northern Lake Superior basin is most likely a result of earlier ice retreat to the Nipigon moraine by ∼ 10.7 ka north of Isle Royale (Leydet et al., 2018; Lowell et al., 2021) where the ice margin paused in its retreat along the divide between Black Bay and Nipigon Bay (Fig. 7b; Teller and Mahnic, 1988). We suggest that water in the western Lake Superior basin was blocked from draining north around both Isle Royale and the Keweenaw Peninsula by the presence of ice spanning from the Keweenaw Peninsula to the Nipigon moraine via the outer subaqueous IR1 moraine and Isle Royale (Fig. 7b).

In this case, the initial influxes of meltwater from glacial Lake Agassiz, which carried the reddish varve-forming sediment, would have first drained into the western Lake Superior basin at Black Bay and then drained across the Keweenaw Peninsula to the eastern Lake Superior basin via the Portage and Allouez outlets (Fig. 7b, c). The Portage and Allouez gaps are known to have previously carried water, but neither has been dated (Hughes, 1963), and we propose that these outlets were abandoned following the rapid ice retreat from Isle Royale and IR2 after ∼ 10.1 ka (Fig. 7d). Ice retreat from Michigan's UP in the eastern Lake Superior basin occurred by ∼ 10.5–10.2 ka (Breckenridge, 2007; Fisher and Whitman, 1999) and quickly retreated north, perhaps to the vicinity of the Chapleau moraine ∼ 10.2 ka (Saarnisto, 1974).

Following retreat from Isle Royale and the IR2 moraine, the LIS ice retreated to the Nakina moraine, exposing the entirety of the northern Lake Superior basin and small lakes between Lake Superior's North Shore and the Nakina moraine in less than 1 kyr (Fig. 7d; Bajc et al., 1997; Breckenridge, 2007, 2013; Breckenridge et al., 2004; Colman et al., 2020; Hyodo and Longstaffe, 2011; Kelly et al., 2016; Lowell et al., 2021; Saarnisto, 1974; Teller and Mahnic, 1988). During this retreat, gray varves were deposited across all of the Lake Superior basin, terminating in a series of 36 anomalously thick varves that are associated with the final pulses of meltwater from the LIS before it retreated from the Nakina moraine at ∼ 9 ka (Breckenridge et al., 2004; Fig. 7). This inferred retreat timing, extrapolated over the >300 km distance from Isle Royale and the IR2 moraine to the Nakina moraine, means ice retreat rates south of Isle Royale nearly tripled to a mean retreat rate of ∼ 270 m yr⁻¹ as the ice margin retreated across deep waters of northern Lake Superior (Hanson and Hooke, 2000; Colman et al., 2020). This rapid retreat accounts for the lack of glacial deposits on the northeasternmost two-thirds of Isle Royale where Precambrian bedrock is exposed (Fig. 3; Huber, 1973).

6 Implications

The interpretations we draw from ¹⁰Be exposure ages from recessional moraines on Isle Royale provide new, quantitative constraints on the timing of ice retreat from Isle Royale. These new ¹⁰Be exposure ages from within Lake Superior expand the geographical coverage of quantitative age constraints on the timing of LIS retreat from the Great Lakes. By inferring that subaerial moraine construction on Isle Royale is contemporaneous with subaqueous moraine construction between Isle Royale and Michigan's Keweenaw Peninsula, we show that the spatial distribution of red varves throughout the Lake Superior basin can be used to delineate the position of the LIS margin in the early Holocene. In doing so, we provide support for the notion that southern margins of the LIS remained in contact with the southern shores of Lake Superior into the early Holocene prior to its final and rapid retreat. Interpreting our new ¹⁰Be-based chronology of ice retreat from Isle Royale within the spatial and geomorphological context of existing proglacial and subglacial landforms and stratigraphy bridges existing ice retreat chronologies along Lake Superior's North Shore with the Lake Superior south shore along the Keweenaw Peninsula. Interpretations drawn from the dataset presented here demonstrate that dating glacial features on islands can be integral to reconstructing regional patterns of deglaciation across bodies of water and drawing relationships between subaerial and subaqueous glacial landforms and regional stratigraphy.

Appendix A: Photographs of sample sites

Glacial erratics are common across Isle Royale. Those erratics suitable for ¹⁰Be exposure-age dating are easily identifiable because they are often granite or metamorphic rocks from the Canadian Shield, which stand in stark contrast to the mafic igneous and volcaniclastic rocks comprising the bedrock of Isle Royale. Most erratics are smaller than those typically found in alpine environments because of the long transport distances from their sources. During the 3 days of fieldwork on Isle Royale, many erratics were identified but only those most suitable for dating were sampled. Sampled erratics were those positioned on moraine crests or broad uplands with minimal topographic shielding and least likely to have rolled down hills or have complex exposure and/or burial histories (Fig. A1).

Figure A1 Photographs of the 11 samples used in this study. Rock hammers, plastic zip bags, GPS units, and people are present for scale.

Appendix B: Sensitivity analyses of exposure-age calculation decisions

In this study, we present exposure ages that were calculated using the online exposure-age calculator presented by Balco et al. (2008; version 3) using the default ¹⁰Be production rate calibration and the LSDn ¹⁰Be production rate scaling scheme (Lifton et al., 2014). Standard ¹⁰Be production shielding corrections were made (Table 1), but we do not incorporate any other production rate scaling options nor use a regional ¹⁰Be calibration dataset. In each section of Appendix B, we explore how ¹⁰Be exposure ages from Isle Royale would have shifted had we made different decisions. Importantly, calculating ages with any of these adjustments only results in minor changes in exposure ages and would not alter the results of this study.

B1 Choice of exposure-age calculator

We considered the use of two online calculators, the first presented by Balco et al. (2008), which is regularly updated, and the newer Ice-TEA calculator (Jones et al., 2019). Both calculators follow the guiding principles of cosmogenic nuclide production (Gosse and Phillips, 2001), but with recognized differences, including how atmospheric pressure is incorporated into age calculations or how easily users can choose regional production rate calibration datasets. Differences in calibration datasets have led to differences in exposure ages from the same input data for at least one previous study in the Great Lakes region (e.g., Lowell et al., 2021). In this study, exposure ages calculated with the Ice-TEA calculator are, on average, 6.7 % older (range: 5.8 %–7.0 %) than ages calculated using Balco et al.'s (2008) calculator (Table B1). We chose to present ages from Balco et al.'s (2008) calculator, primarily because of its longevity and because it is used by other ¹⁰Be exposure-age studies in the Great Lakes region.

Table B1 Exposure ages calculated using Balco et al. (2008) versus Ice-TEA.

B2 Choice of ¹⁰Be production rate calibration

Previous ¹⁰Be exposure-age studies in the Great Lakes region (e.g., Leydet et al., 2018; Lowell et al., 2021) used the Northeast North America (NENA) calibration for ¹⁰Be production rates (Balco et al., 2009), which results in exposure ages that are each 0.5 % younger than those we use in determining the mean exposure age of deglaciation from Isle Royale (Table B2). Although the samples we collected from Isle Royale are ice-marginal landforms found at relative high latitudes (∼ 48^∘ N) and at elevations <1000 m a.s.l., we use the default ¹⁰Be calibration production rate in our age calculations using Balco et al.'s (2008) calculator, primarily because of our distance from NENA calibration sites.

Table B2 Exposure ages calculated in Balco et al. (2008) using the default versus Northeast North America calibration datasets.

B3 Choice of ¹⁰Be production rate scaling scheme

Previous studies present exposure ages using different production rate scaling schemes including non-time-dependent scaling based on atmospheric measurements (St; Lal, 1991; Stone, 2000) and time-dependent scaling based on atmospheric measurements and paleomagnetic reconstructions (Lm; Lal, 1991; Nishiizumi et al., 1989; Stone, 2000). Our preference of ¹⁰Be production rate scaling is to use Lifton et al.'s (2014) LSDn scaling scheme because it accounts for changes in the strength of earth's magnetosphere, changes in solar cosmic ray output, and numerous cosmogenic isotope production pathways (LSDn; Lifton et al., 2014). Exposure ages of erratics on Isle Royale based on non-time-dependent spallogenic scaling (St) are on average 0.4 % older (range: 0.2 %–0.7 %) and exposure ages based on time-dependent spallogenic scaling (Lm) are on average 2.2 % younger (range: 2.0 %–2.8 %) than ages derived using Lifton et al.'s (2014)

LSDn scaling scheme (Table B3).

Table B3 Exposure ages calculated in Balco et al. (2008) using different ¹⁰Be production rate scaling schemes.

B4 Effects of snow shielding

Previous ¹⁰Be studies in the region used sample collection strategies to minimize the potential of snow shielding (e.g., unforested, windswept, high-elevation locations; Leydet et al., 2018; Lowell et al., 2021). This was not an option on Isle Royale. Given its position in the middle of Lake Superior, heavy lake effect snow is probable on Isle Royale, and the forest at the Mt. Desor and Sugar Mountain moraines may prevent fallen snow from blowing away easily. We considered the effects of snow shielding at our sample sites using average monthly snow-depth totals that were measured at Isle Royale's Mott Island Station:

• Station: #205637 (National Climatic Data Center, NCDC)

• Location: 48^∘06^′ N, 88^∘33^′ W

• Elevation: 610 ft. above sea level (∼186 m)

• Period of record: 1 October 1940 to 29 April 2016 (41.8 % of possible observations reported; Table B4)

• Data access: Western Regional Climate Center, https://wrcc.dri.edu/ (last access: 22 February 2023).

Snow shielding calculations are based on Eq. (3.76) in Gosse and Phillips (2001):

$$S_{\mathrm{snow}}={\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}\sum _i^{\mathrm{12}}{e}^{-\left(z_{\mathrm{snow}},i{\mathit{\rho }}_{\mathrm{snow}},i/\mathrm{\Lambda }\right)},$$

where z_snow,i is the monthly average snow thickness above the ground (cm), ρ_snow,i is the monthly average snow density (g cm⁻³), and Λ is attenuation of cosmic rays (160 g cm⁻²). We use an average old-age snow density of 0.27 g cm⁻³ measured from Isle Royale (Peterson, 1977).

Table B4 Monthly average snow depth at Mott Island (October 1940–April 2016).

Including snow shielding led to exposure ages that were, on average, 2.9 % older (range: 2.6 % to 3.3 %) than without snow shielding (Table B5). We chose not to include these ages in our analysis because measured snowfall from 1940 to 2016 may not accurately reflect Holocene averages.

Table B5 Exposure ages calculated in Balco et al. (2008) without versus with snow shielding.

Table B6 Exposure ages calculated in Ice-TEA without versus with glacial isostatic adjustment (GIA) corrections.

B5 Effects of glacial isostatic adjustment

Sample sites in this study are ∼ 300 km up-ice of the Laurentide Ice Sheet's Last Glacial Maximum terminal moraines (Ullman et al., 2015), and strandlines of proglacial lakes that once occupied the western Lake Superior basin have experienced considerable uplift (Breckenridge, 2013; Farrand, 1969). In considering the effects of glacial isostatic adjustment (GIA) on the exposure ages presented in this study, we compare exposure ages from all samples using only the Ice-TEA software package's “Correct for Elevation Change” tool and the ICE-6G GIA model (Jones et al., 2019; Peltier et al., 2015). GIA-corrected ages are compared with uncorrected ages. Exposure ages with GIA corrections were on average 3.8 % older (range: 1.1 %–8.1 %) than exposure ages uncorrected for GIA (Table B6). Following Lowell et al. (2021), we chose to use exposure ages uncorrected for GIA to maintain regional comparability of exposure-age datasets.