Delayed and rapid deglaciation of alpine valleys in the Sawatch Range, southern Rocky Mountains, USA

We quantify retreat rates for three alpine glaciers in the Sawatch Range of the southern Rocky Mountains following the Last Glacial Maximum using ¹⁰Be ages from ice-sculpted, valley-floor bedrock transects and statistical analysis via the BACON program in R. Glacier retreat in the Sawatch Range from at (100 %) or near (∼83 %) Last Glacial Maximum extents initiated between 16.0 and 15.6 ka and was complete by 14.2–13.7 ka at rates ranging between 35.6 and 6.8 m a⁻¹. Deglaciation in the Sawatch Range commenced ∼2–3 kyr later than the onset of rising global CO₂ and prior to rising temperatures observed in the North Atlantic region at the Heinrich Stadial 1–Bølling transition. However, deglaciation in the Sawatch Range approximately aligns with the timing of Great Basin pluvial lake lowering. Recent data–modeling comparison efforts highlight the influence of the large North American ice sheets on climate in the western United States, and we hypothesize that recession of the North American ice sheets may have influenced the timing and rate of deglaciation in the Sawatch Range. While we cannot definitively argue for exclusively North Atlantic forcing or North American ice sheet forcing, our data demonstrate the importance of regional forcing mechanisms for past climate records.

1 Introduction

Alpine glaciers worldwide underwent substantial retreat in response to climate warming during the last deglaciation (Shakun et al., 2015; Palacios et al., 2020). However, the general trend of warming through the last deglaciation was interrupted by internally forced and regionally heterogeneous climate changes such as the cool Heinrich Stadial 1 (17.5–14.7 ka), abrupt warming into the Bølling–Allerød period (14.7–12.9 ka), and the Younger Dryas cold period (12.9–11.7 ka) all centered in the North Atlantic region (NGRIP members, 2004; Rasmussen et al., 2014). To thoroughly characterize the influence of these climatic oscillations, their expression throughout the Northern Hemisphere is often investigated using records of mountain glaciation (Ivy-Ochs et al., 2006; Schaefer et al., 2006; Young et al., 2011; Shakun et al., 2015; Marcott et al., 2019; Young et al., 2019). Mountain glacier deposits serve as suitable archives since mountain glaciers are particularly sensitive to changes in climate (e.g., Oerlemans, 2005; Roe et al., 2017). Furthermore, where deposits are carefully mapped and dated, quantitative retreat or thinning rates of glaciers can be compared to records of climatic forcings. Statistical approaches to quantify retreat and thinning rates have been previously applied to ice sheets (e.g., Johnson et al., 2014; Jones et al., 2015; Koester et al., 2017; Small et al., 2018; Lesnek et al., 2020) but only for a few mountain glaciers (e.g., Hofmann et al., 2019).

In the western United States (US; Fig. 1), mountain glaciers expanded out of the high elevations of the Rocky Mountains, the Sierra Nevada, and many other smaller ranges during the Last Glacial Maximum (LGM; Porter et al., 1983; Pierce, 2003). During the last deglaciation, many glaciers retreated from their extended LGM positions and eventually melted from their cirques during the late glacial to early Holocene (e.g., Munroe and Laabs, 2017; Marcott et al., 2019). Yet, the temporal and spatial patterns of retreat throughout the western US and their relationship to hemispheric and global forcing are still a subject of debate. Glaciers in the western US may have retreated in response to increasing global temperature forced by rising atmospheric CO₂ concentrations, thus broadly synchronous with other mountain glaciers around the world (e.g., Shakun et al., 2015; Marcott et al., 2019). However, some evidence suggests a delay of deglaciation until the Bølling due to either persistent stadial conditions (e.g., Young et al., 2011) or as a response to increased local moisture supply to some glaciers from nearby pluvial lakes (e.g., Laabs et al., 2009).

Figure 1 Key moraine chronologies from the southern Rocky Mountains and locations of glaciation centers and large pluvial lakes in the western US following the Last Glacial Maximum (LGM). LL: Lake Lahontan, LB: Lake Bonneville, A: Colorado Front Range, B: Sangre de Cristo Mountains, C: San Juan Mountains, and D: Winsor Creek valley, New Mexico. The largest star corresponds to our field site in the Sawatch Range. LGM ice limits are from Dalton et al. (2020). The inset is of the western portion of North America. CIS: Cordilleran Ice Sheet, LIS: Laurentide Ice Sheet.

Over a decade of work has resulted in detailed moraine chronologies in three adjacent alpine valleys in the Sawatch Range of central Colorado (Fig. 2; Briner, 2009; Young et al., 2011; Shroba et al., 2014; Leonard et al., 2017b; Schweinsberg et al., 2020). While these studies primarily focused on mapping and dating the range-front moraines and associated outwash terraces, a transect of ages from bedrock samples in Lake Creek valley (Fig. 2) documented rapid retreat between 15.6±0.7 and 13.7±0.2 ka (Leonard et al., 2017b; Schweinsberg et al., 2020). Schweinsberg et al. (2020) suggested a possible link between North Atlantic climate forcing and the rapid deglaciation observed in Lake Creek valley, but similar transects from adjacent valleys are lacking to bolster or refute this hypothesis.

Figure 2 Ice-sculpted bedrock ¹⁰Be ages from Lake Creek (LC; orange), Clear Creek (CC; green), and Pine Creek valleys (PC; blue). New ages are open circles and previously published ages are closed. Bedrock ages with italicized labels are suspected outliers. Included are LGM and recessional moraines (solid colored lines and labels with black text) with reported ages for the LGM moraines in all three valleys, including the younger mode in Pine Creek valley at 16.0±0.9 ka (n=7; Young et al., 2011) and a recessional moraine in Lake Creek valley at 15.6±0.7 ka (n=5; Schweinsberg et al., 2020). There is a similar, undated recessional moraine in Clear Creek valley that we hypothesize is also ∼16 ka. Ice-sculpted bedrock ages reported here include analytical uncertainty, and moraine ages are reported as the mean and 1 standard deviation. Glaciers at their mapped LGM extents are delineated in gray.

Here, we combine 12 new cosmogenic ¹⁰Be exposure ages with 10 previously published ¹⁰Be ages from bedrock samples along transects in three adjacent alpine valleys in the Sawatch Range, southern Rocky Mountains (Fig. 2). By dating bedrock sites along valley transects, we characterize the timing and pace of glacier retreat during the last deglaciation. We calculate rates of deglaciation for each valley with best-fit time–distance plotting using the R program BACON (Fig. 4). Our results suggest that glaciers in the Sawatch Range may have been influenced more heavily by regional forcing than by global CO₂ concentrations.

2 Setting

The high peaks of south-central Colorado and northern New Mexico compose the southern end of the Rocky Mountain Range in North America and were home to many alpine glaciers during multiple glaciations throughout the Pleistocene (Fig. 1; Pierce, 2003; Leonard et al., 2017b; Marcott et al., 2019; Laabs et al., 2020; ages discussed below are recalculated using the Promontory Point production rate calibration of Lifton et al., 2015, and the LSDn scaling model of Lifton et al., 2014). Transects of ¹⁰Be ages from bedrock along valley axes exist for a few valleys in the upper Boulder Creek drainage in the Front Range, Colorado (Benson et al., 2004; Ward et al., 2009; Dühnforth and Anderson, 2011). While some evidence from the Boulder Creek drainages may suggest delayed deglaciation, chronologic scatter in the ages makes it difficult to determine the exact timing and how quickly glaciers retreated to their cirques. Existing ages from one valley, the Sangre de Cristo Range in south-central Colorado, suggest that a glacier there remained at or readvanced to near its LGM terminus at ∼16 ka but then retreated to its cirque in a period of ∼2 kyr (Leonard et al., 2017a). In the Animas River valley of the San Juan Mountains, southwest Colorado, existing ¹⁰Be ages indicate that glacier retreat began as early as ∼19 ka, with complete retreat of nearly 70 % of the total valley length beginning ∼16 ka and finishing by ∼12.7 ka (Guido et al., 2007). Relatively early initial retreat of the glacier in the Animas River valley is contingent on dating at a single site. Near Baldy Peak in northern New Mexico, LGM moraines and what appear to be cirque moraines have been surveyed in the Winsor Creek valley (Armour et al., 2002; Marcott et al., 2019). ¹⁰Be ages from the cirque, ∼4 km up-valley from the LGM moraines, range from 15.8 to 14.3 ka, suggesting that the glacier retreated to near its cirque within that interval. The recessional and LGM moraines remain undated so it is difficult to know when the glacier began retreating. In summary, while there is some chronologic scatter in ages from these sites, there is evidence to suggest that some glaciers in the southern Rocky Mountains remained relatively expanded through the beginning of the last deglaciation and were delayed in their retreat. However, once retreat was underway, all sites observed thus far reveal that glaciers completely retreated at least up to their cirques prior to the Younger Dryas cold period with no evidence for subsequent moraine deposition.

Prominent moraines originally mapped as part of the surficial geologic map of the Granite 7.5^′ quadrangle (updated by Shroba et al., 2014) exist at the mouths of multiple glacially sculpted valleys within the Sawatch Range (e.g., Briner et al., 2009; Young et al., 2011; Brugger et al., 2019a; Schweinsberg et al., 2020). Of these, moraines deposited at the mouths of three adjacent valleys, Lake Creek, Clear Creek, and Pine Creek, have been thoroughly surveyed and dated (Fig. 2; Briner, 2009; Young et al., 2011; Schweinsberg et al., 2020). The moraine chronologies reported thus far reveal that following the LGM (which culminated between ∼22 and 19 ka), a recessional moraine at 82 % of the LGM position sampled in the Lake Creek system was deposited at 15.6±0.7 ka (Schweinsberg et al., 2020). There is a similar-appearing moraine at 83 % of the LGM position in Clear Creek valley. Although it is undated, we tentatively correlate this moraine in Clear Creek valley with the moraine dated to 15.6±0.7 ka in Lake Creek valley. Finally, there is no recessional moraine in Pine Creek valley, but a cluster of ages at 16.0±0.9 ka from the LGM moraine suggest that the glacier readvanced to or remained at its LGM extent until nearly the same time when glaciers in the other two valleys deposited recessional moraines (Briner, 2009; Young et al., 2011). Young et al. (2011) argued that since all three glaciers are east-facing and in close proximity – yet show differences in the timing of LGM culmination between the valleys – it is possible that non-climatic factors, such as glacier hypsometry, may have influenced the timing and extent of LGM culminations. While there are preexisting ¹⁰Be ages measured in a transect along Lake Creek Valley that track the retreat of the glacier through the last deglaciation, the other two valleys have not yet been surveyed. As such, it remains unclear if glacier hypsometry also influenced the timing and pace of deglaciation between all three adjacent valleys.

3 Methods and materials

Sample collection for ¹⁰Be dating from Clear Creek and Pine Creek valleys was conducted in the summers of 2017 and 2018. A total of 12 samples were collected from exposed, glacially sculpted bedrock surfaces along the Clear Creek (n=8) and Pine Creek (n=4) valley floors, spanning from just within range-front moraines up to each respective cirque (Figs. 2 and 3). Bedrock surfaces located in the bottoms of valley floors – where glacial erosion is maximized – were specifically targeted since the potential for incomplete scouring of these surfaces can lead to inherited nuclides and ages that are older than expected.

Figure 3 Field photos of ice-sculpted bedrock surfaces from selected locations. Clockwise from top left: 18CC-04, 18PC-01, 18PC-05, 17CC-08. The color scheme for ages matches Figs. 2 and 4: Clear Creek valley samples are in green; Pine Creek valley samples are in blue.

Samples were processed at the University at Buffalo Cosmogenic Isotope Laboratory following the versions of quartz purification and beryllium extraction procedures refined at the University of Vermont (Corbett et al., 2016). After quartz purification, samples were dissolved in acid along with a ⁹Be carrier spike in two batches each with a process blank. Beryllium was then purified and extracted, oxidized, and packed into targets for measurement at the Center for Accelerated Mass Spectrometry at Lawrence Livermore National Laboratory. ¹⁰Be∕⁹Be ratios were measured and standardized to the reported 07KNSTD3110 ratio of $\mathrm{2.85}\times {\mathrm{10}}^{-\mathrm{12}}$ (Nishiizumi et al., 2007). For samples collected in 2018, the process blank ¹⁰Be∕⁹Be ratio was $\mathrm{2.96}\times {\mathrm{10}}^{-\mathrm{15}}$, and for samples collected in 2017 the process blank ¹⁰Be∕⁹Be ratio was $\mathrm{9.56}\times {\mathrm{10}}^{-\mathrm{16}}$ (see Table 1 for details on sample collection dates). Our 12 ages and 10 previously published ages were calculated using the Cronus Earth online calculator (developmental version 3; https://hess.ess.washington.edu/math/index_dev.html, last access: 12 April 2020; Balco et al., 2008). We calculate ages using the Promontory Point production rate (Lifton et al., 2015) and the LSDn scaling model (Lifton et al., 2014) – a combination used extensively throughout the western US (e.g., Licciardi and Pierce, 2018; Quirk et al., 2018; Brugger et al., 2019b; Schweinsberg et al., 2020). Below, we discuss in more detail how different production rate calibrations and scaling schemes impact our results. We do not attempt to make any corrections for snow cover or post-depositional bedrock surface erosion.

Table 1 Sample data and ¹⁰Be ages.

Notes: rock density for all samples is 2.65 g cm³; zero surface erosion rate applied to all samples. Sample names of previously published ages are italicized, and sample names for new ages are in bold. ^a Ages calculated using the Promontory Point production rate calibration (Lifton et al., 2015) and LSDn scaling (Lifton et al., 2014). ^b Ages calculated using the northeast North America production rate calibration (Balco et al., 2009) and Lm scaling (Lal, 1991; Stone, 2000). ^c Process blank for the following samples was $\mathrm{9.56}\times {\mathrm{10}}^{-\mathrm{16}}$ (atoms g⁻¹). ^d Process blank for the following samples was $\mathrm{2.96}\times {\mathrm{10}}^{-\mathrm{15}}$ (atoms g⁻¹).

To calculate retreat rates, we used the BACON program in R (Blaauw and Christen, 2011). This program generates age–depth models for stratigraphic records based on chronologic constraints at various depths. Here, we use the ¹⁰Be ages and their 1-sigma internal uncertainties measured in each valley as the age input and the geographic coordinates of each age as the depth inputs. The position along the valley floor is scaled such that the toe of the glacier at the LGM is the starting point (e.g., 100 % or maximum length), and the base of each valley's cirque wall is the end point (e.g., 0 % or minimum length). The model then interpolates between each point using Bayesian analysis and the geologic principle of superposition to build an age–length model with an unweighted statistical treatment of uncertainty. The interpolation between points is smoothed (i.e., nonlinear) based on retreat rates at previous positions. The retreat rates presented here are net retreat rates because it is possible there may have been short-lived readvances that did not lead to significant moraine deposition. BACON outputs a time series of age–length points and non-Gaussian 95 % confidence intervals. Calculated retreat rates are assumed to be linear, and we report the 95 % uncertainty range.

4 Results

The 12 new sculpted-bedrock ¹⁰Be ages reported here range 15.8±0.3–13.7±0.3 ka (Fig. 2; Table 1). Combined with the previously published samples in our study area, all 22 sculpted-bedrock ¹⁰Be ages, which span from immediately inboard of the innermost moraine to the cirque floors, range between 16.0±0.4 and 13.5±0.3 ka (Fig. 2, Table 1). In Lake Creek valley, seven ages span 67 %–1 % of the distance of the valley floor, ranging between 15.2±0.4 and 13.5±0.3 ka. Nine ¹⁰Be ages spanning 68 %–1 % in Clear Creek valley range between 15.3±0.2 and 13.7±0.2 ka. In Pine Creek valley, six ¹⁰Be ages span 78 %–3 % and range between 16.0±0.4 and 14.2±0.3 ka.

Most ages in each valley are in stratigraphic order and fall within the 95 % confidence interval calculated in BACON, except for four ages (Fig. 4). Ages from Lake Creek valley suggest the glacier retreated from its recessional moraine position (82 %) at 15.6±0.7 ka and reached its cirque (∼1 %) by 13.7±0.2 ka. Clear Creek valley ages suggest the glacier retreated from its recessional moraine position (83 %) at 15.6±0.7 ka and reached its cirque (∼1 %) by 13.7±0.3 ka. Finally, Pine Creek valley ages suggest the glacier was at its LGM extent (100 %) until 16.0±0.9 ka and then retreated to its cirque (∼3 %) by 14.2±0.2 ka.

Figure 4 Summary plots of ¹⁰Be ages and BACON statistical analysis results. (a) Lake Creek valley (orange), (b) Clear Creek valley (green), and (c) Pine Creek valley (blue). Ages in solid black fill at the bottom of each transect are from recessional moraine ages (Young et al., 2011; Schweinsberg et al., 2020). BACON results are mean (color lines) and 95 % confidence intervals (gray shading). The left y axes are total valley-floor distances from the LGM moraine to the base of each respective cirque headwall (note that scales are different because valley lengths are different). The right y axes are the same but normalized values; 1 indicates the LGM moraine position and 0 the base of the cirque headwall. (d) Distribution of all ages using the PPT (Lifton et al., 2015) and LSDn (Lifton et al., 2014) as well as the NENA (Balco et al., 2009) and Lm (Lal, 1991; Stone, 2000) production rate calibration and scaling scheme combinations. All other reasonable combinations mentioned in the text produce ages that fall somewhere between the ages calculated using the PPT–LSDn and NENA–Lm combinations.

Results from BACON analysis suggest that the net retreat rate for the glacier in Lake Creek valley between 15.6±0.7 ka (Schweinsberg et al., 2020) and 13.7±0.2 ka ranges 35.6–13.8 m a⁻¹ at 95 % confidence (Fig. 4). The net retreat rate calculated from BACON for the glacier in Clear Creek valley between 15.6±0.7 and 13.7±0.3 ka ranges 15.5–8.2 m a⁻¹ at 95 % confidence. Finally, the net retreat rate for the glacier in Pine Creek valley from the LGM position at 16.0±0.9 ka (Young et al., 2011) to 14.2±0.3 ranges 18.3–6.8 m a⁻¹ at 95 % confidence. Removal of potential outliers reduces retreat rates by 1.7 %, 2.7 %, and 6 % for Lake Creek, Clear Creek, and Pine Creek valleys, respectively. The calculated average valley gradients for each valley – measured as the elevation change divided by the horizontal length of each valley bottom transect from the LGM moraine up to the base of each respective cirque – are 29 m km⁻¹ for Lake Creek valley, 37 m km⁻¹ for Clear Creek valley, and 65 m km⁻¹ for Pine Creek valley.

5 Discussion

5.1 Reliability of bedrock ages

While most bedrock ages along each valley transect are in stratigraphic order, we find four ages that do not comply with stratigraphic order and fall outside the 95 % confidence interval of the retreat rates calculated from BACON. For example, in Lake Creek valley, one age at 13.5±0.3 ka is younger than all up-valley ages, which average 13.8±0.2 ka (excluding one possible outlier outside the BACON 95 % confidence interval). In Clear Creek valley, the age from the farthest down-valley site of 14.8±0.3 ka may be a possible outlier because the next three ages up-valley are all older and in stratigraphic order, the oldest of which is 15.3±0.2 ka. Finally, one age from the Pine Creek cirque of 15.3±0.3 ka may be an outlier because it is older than the next age down-valley (14.6±0.3 ka) and a second sample from the cirque of 14.2±0.3 ka. Two suspected outliers are older than expected, which may have been caused by insufficient glacial erosion. The two remaining potential outliers are younger than expected, which could have resulted from excessive soil and snow cover, enhanced post-depositional bedrock surface erosion, or a combination of these factors.

Although we interpret our results using the Promontory Point production rate calibration site (Lifton et al., 2015) and the LSDn scaling scheme (Lifton et al. 2014), we calculate exposure ages with other commonly used calibration sites (e.g., northeastern North America NENA from Balco et al., 2009, and the “global” production rate from Borchers et al., 2016) and another commonly used scaling scheme (Lal/Stone Lm scaling; Lal, 1991; Stone, 2000). Samples used for the NENA production rate calibration range in elevation between ∼50 and 400 m a.s.l. and are located ∼3000 km northeast of the Sawatch Range. This combination produces the oldest ages, given the previously mentioned reasonable production rate calibrations and scaling schemes, that are between 9 % and 12 % older than when using PPT–LSDn (all other combinations fall somewhere between; Fig. 4; Table 1). We do not feel confident in calculating exposure ages using other production rate calibration sites since the sites in closest proximity likely shared the most similar exposure histories. Ultimately, we favor the Promontory Point production rate calibration site (Lifton et al., 2015) because the site is closest in both location (the site is ∼600 km from the Sawatch Range) and elevation (sample elevations are ∼1600 m a.s.l.) to our study area.

5.2 The last deglaciation of the Sawatch Range and the southern Rocky Mountains

The pattern of deglaciation in both Clear Creek valley and Pine Creek valley appears to follow the pattern previously observed in Lake Creek valley (Young et al., 2011; Leonard et al., 2017b; Schweinsberg et al., 2020). All three glaciers remained at – or readvanced to – (100 %) or near (82 %–83 %) their LGM extents between 16.0 and 15.6 ka, after which all three glaciers rapidly retreated to their cirques within the next ∼2 kyr at rates ranging between 35.6 and 6.8 m a⁻¹. It is possible that the glacier in Pine Creek valley began retreating ∼500 years earlier than the other two glaciers and likewise completely deglaciated ∼500 years earlier. Pine Creek valley is shorter and steeper than the other two valleys. Thus, it is possible that variations in valley hypsometry between Pine Creek valley and the other two valleys may have caused the slight difference in their deglaciation histories (Young et al., 2011). We also observe that Pine Creek valley has the steepest average valley gradient and generally the slowest net retreat rate, which is predictably a direct result of valley hypsometry since glacier lengths in steeper valleys generally adjust less to equivalent changes in equilibrium line altitude (ELA). On the other hand, glaciers occupying the lower-gradient Lake Creek and Clear Creek valleys experienced generally higher reconstructed rates of retreat. Regardless, we find that while Pine Creek may have initiated ∼500 years sooner than the other two, all three valleys were in a period of ∼1–1.5 kyr synchronous retreat once the other two glaciers began retreating. We conclude that while there may have been some hypsometric influences on the timing of deglaciation across our study site, evidence suggests that these influences did not keep the glaciers from synchronously retreating during a majority of their deglaciation. We find that all three valley glaciers did not begin significantly retreating until ∼5–6 kyr after the culmination of the LGM in the Sawatch Range (since we assume that boulder ages on both LGM and recessional moraines represent the timing of moraine abandonment). However, once glacier retreat initiated, deglaciation was completed within ∼2 kyr.

From the existing records in the southern Rocky Mountains synthesized above, we find that the pattern of deglaciation observed in the Sawatch Range was consistent in a few but not all sites across the region. Collecting more records of alpine deglaciation in the southern Rocky Mountains may be necessary to further test which pattern, if any, is the dominant pattern of deglaciation in the region.

5.3 Drivers of southern Rocky Mountain deglaciation

Records of global climate change over the last deglaciation suggest a link between rising CO₂ concentrations and global temperature (Denton et al., 2010; Shakun et al., 2012; Putnam et al., 2013; Fig. 5). However, there is noticeable spatial heterogeneity in both the timing and magnitude of warming through the last deglaciation that cannot be attributed to global CO₂ forcing alone (e.g., Clark et al., 2012). We find that the initiation of significant deglaciation in some locations across the southern Rocky Mountains lagged rising CO₂ concentrations by as much as ∼2–3 kyr (Fig. 5), which suggests these glaciers were more likely influenced by regional forcings rather than global CO₂.

Figure 5 Deglaciation of the Sawatch Range compared to other climate proxies. From top to bottom: atmospheric CO₂ concentrations (Bereiter et al., 2015); global and Southern Hemisphere temperature stacks (Shakun et al., 2012); synthesized Greenland temperature from ice cores (Buizert et al., 2014); lake level reconstructions of Lake Bonneville (LB; black dashed line) and Lake Lahontan (LL; gray dashed line) from Reheis et al. (2014); normalized BACON plots from Lake Creek (LC; orange), Clear Creek (CC; green), and Pine Creek valleys (PC; blue). Vertical lines correspond to the onset of CO₂ rise beginning ∼18 ka and the Heinrich Stadial 1–Bølling transition at 14.7 ka.

Ice core records – among other records – reveal a complex pattern of abrupt warming and cooling events that occurred in the North Atlantic region during the last deglaciation (Fig. 5; Buizert et al., 2014). Despite rising CO₂ concentrations beginning ∼18 ka, North Atlantic records reveal that cold conditions persisted until 14.7 ka, known as Heinrich Stadial 1 (HS-1). Following these sustained cold conditions, an abrupt transition to warmer conditions is marked by the HS-1–Bølling boundary at 14.7 ka (Buizert et al., 2014). We find that deglaciation at some locations in the southern Rocky Mountains encompasses the HS-1–Bølling transition. Furthermore, the relatively rapid and short-lived nature of retreat for glaciers in the Sawatch Range – and some others across the southern Rocky Mountains – appears to be more consistent with the abrupt manner of warming observed in the North Atlantic. However, glaciers were already retreating prior to the abrupt HS-1–Bølling transition at ∼14.7 ka. Therefore, it is difficult to argue that North Atlantic warming alone forced glacier retreat in the southern Rocky Mountains.

In addition to the alpine glaciers that existed in the mountainous regions of the western US during the late Pleistocene, large pluvial lakes such as Lake Lahontan and Lake Bonneville existed across the Great Basin (Fig. 1; Gilbert, 1890; Russell, 1885; Orme, 2008). These lakes could have been sustained by increased precipitation delivery to the southwestern US (e.g., Munroe and Laabs, 2013; Oster et al., 2015; Lora and Ibarra, 2019) or were maintained simply by colder temperatures persisting throughout the region (e.g., Benson et al., 2013). Recent syntheses of past Great Basin lake levels reveal that Lahontan and Bonneville lakes resided at relative high stands between 15.5 and 14.5 ka (Benson et al., 2013; Reheis et al., 2014; Oviatt, 2015). After this time, each lake experienced notable declines in lake level (Fig. 5), which could have been the result of reduced precipitation due to rearranging storm tracks, warming temperature, or a combination of both (Benson et al., 2013; Oster et al., 2015; Lora and Ibarra, 2019).

Recent modeling efforts have highlighted how North American ice sheets likely influenced atmospheric circulation and regional climate throughout the Pleistocene (COHMAP members, 1988; Löfverström et al., 2014; Liakka and Lofverstrom, 2018; Lora and Ibarra, 2019). Specifically, there appears to have been a drastic shift in climatologies over the western US when the Cordilleran (CIS) and Laurentide (LIS) ice sheets separated (Lofverstrom et al., 2014; Lora et al., 2016; Tulenko et al., 2020). For example, during the last deglaciation, once the CIS and LIS separated, some model results suggest the western US became warmer and drier (Lora et al., 2016). The latest synthesis of the last deglaciation of the major North American ice sheets suggests that the separation occurred between 16 and 15 ka (Dalton et al., 2020). Thus, it is possible that the saddle collapse and separation of the CIS and LIS and resulting atmospheric reorganization may have led to both drastic pluvial lake level reductions and the rapid deglaciation of some glaciers in the southern Rocky Mountains.

Between North Atlantic forcing and North American ice sheet forcing, it is difficult to conclude what the primary driver of deglaciation in the Sawatch Range was; it may be a combination of both forcings. We find that the approximate timing and rate of deglaciation observed in the Sawatch Range point to abrupt warming and/or drying; this is supported by pluvial lake level records in the western US, which have also been tied to both North Atlantic forcing and North American ice sheet forcing (Munroe and Laabs, 2013; Benson et al., 2013; Lora and Ibarra, 2019). Regardless, the data synthesized here underscore the dominance of regional forcing mechanisms over global forcing mechanisms in some climate records in the western US.

6 Conclusions

We constrain the timing and rate of deglaciation in three alpine valleys in the Sawatch Range, southern Rocky Mountains. Beryllium-10 ages from ice-sculpted bedrock in each valley reveal that the significant retreat of glaciers from their LGM extents (100 %) or near (82 %–83 %) their LGM extents was initiated shortly after 16.0–15.6 ka, despite ∼2–3 kyr of increasing global temperature forced by rising atmospheric CO₂. Glaciers in three adjacent valleys retreated rapidly to their cirques within ∼2 kyr, culminating at ∼ 14.2–13.7 ka, at rates ranging between 35.6 and 6.8 m a⁻¹. We recognize that using the NENA production rate and Lm scaling produces ages 9 %–12 % older than the ages reported herein, which might change the interpretation of the dataset. However, we favor the PPT–LSDn combination because the PPT calibration site is closest in proximity and elevation to the Sawatch Range.

We hypothesize that one of two – or a combination of both – possible regional climatic mechanisms was responsible for driving the pattern of deglaciation for some glaciers in the southern Rocky Mountains. First, we find that for some alpine glaciers in the region, the relatively rapid, short-lived, and synchronous nature of retreat – including those in the Sawatch Range – across the southern Rocky Mountains is more consistent with the abrupt manner of warming observed in the North Atlantic than with steadily increasing global temperature forced by CO₂ rise. However, evidence suggests that glaciers were already retreating prior to the HS-1–Bølling transition. Alternatively, lake level records reveal that both Bonneville and Lahontan lakes lowered nearly in step with some retreating alpine glaciers across the southern Rocky Mountains. Previous studies have linked Great Basin pluvial lake regression to warming and the migration of prevailing storm tracks due to atmospheric reorganization that may have been forced by the separation of North American ice sheets. Thus, warming and drying induced by abrupt atmospheric reorganization at the time of LIS and CIS separation may have driven both Great Basin lake level lowering and rapid alpine glacier retreat in some valleys in the southern Rocky Mountains. While we cannot conclude that either one of the aforementioned forcing mechanisms was solely responsible for deglaciation of the Sawatch Range, we suggest that one or both were stronger controls than increasing global temperature forced by CO₂ rise.