Latest Pleistocene glacial chronology and paleoclimate reconstruction for the East River watershed, Colorado, USA

Abstract

Reconstructing Pleistocene glaciation timing and extent is vital for understanding paleoclimate. Whereas late Pleistocene glaciation has been studied extensively in western North American mountain ranges, the glacial history of the western Elk Range in Colorado remains understudied, particularly in the East River watershed, a site of intense scientific focus. Here we use cosmogenic nuclide exposure and depth–profile dating methods to determine the timing of glaciation in the East River watershed. We use glacier modeling to reconstruct paleoglacier extents and quantify past climate conditions. Our findings indicate that the East River glacier retreated from its maximum position approximately 17–18 ka, moving to recessional positions between 13 and 15 ka, before experiencing more substantial retreat to high-elevation cirques around 13 ka. Glacier modeling suggests that the maximum ice extents at 17–18 ka could have been sustained by temperature depressions of approximately −6.5°C compared to modern conditions, assuming consistent precipitation. Additionally, the ice position at 13–15 ka could have been supported by temperature depressions of around −4.0°C. These results offer insights into the deglaciation timeline in the East River watershed and broader western Elk Range as well as paleoclimate conditions during the late Pleistocene, which may aid future research on critical zone evolution in the East River watershed.

Introduction

Glaciers are sensitive to changes in climate due to the influences of temperature and snow accumulation on ice mass balance. The current widespread retreat of mountain and alpine glaciers highlights the sensitivity of glaciers to modern climate change (Dyurgerov and Meier, 2000; Oerlemans, 2005; Zhang et al., 2015; Roe et al., 2017). Glaciers often leave behind evidence of changes in their extent and ice volume in geomorphic and stratigraphic records. Thus, the geologic record of mountain and alpine glacier fluctuations has been used as a robust indicator of past terrestrial climate change (Agassiz and Bettannier, 1840; Gilbert, 1890; Putnam et al., 2013; Shakun et al., 2015; Marcott et al., 2019), provided the timing of glaciation can be established.

Records of Pleistocene glaciations and past climates in the Rocky Mountains of western North America have been the focus of a diverse array of studies for over a century (e.g. Gilbert, 1890; Ray, 1940; Leonard, 1989; Murray and Locke, 1989; Locke, 1990, 1995; Pierce, 2003; Quirk et al., 2022). In particular, numerous records detailing the chronology and associated climates of Pleistocene mountain glaciation have been developed across the Colorado Rocky Mountains (Benson et al., 2004; Brugger, 2007; Guido et al., 2007; Briner, 2009; Ward et al., 2009; Dühnforth and Anderson, 2011; Young et al., 2011; Leonard et al., 2017; Brugger et al., 2019a, b; Schweinsberg et al., 2020; Tulenko et al., 2020). These studies broadly indicate that mountain glaciers in Colorado reached their maximum extents near the end of or shortly after the global last glacial maximum (LGM; 19.0–26.5 ka; Clark et al., 2009) followed by several thousand years of punctuated deglaciation, and, in some cases, minor re-advances, which is similar to results from other parts of the Rocky Mountains (e.g., Licciardi and Pierce, 2008, 2018; Quirk et al., 2020, 2022) and western North America (Laabs et al., 2020, and references therein). The glacier climate reconstructions available for the Rocky Mountains of Colorado (Brugger, 2010, 2019a, b; Leonard et al., 2017; Brugger et al., 2021) generally suggest temperature depressions of 7–9°C colder than present and only modest changes in precipitation. However, compared to other ranges in the Colorado Rocky Mountains, the Elk Range (Fig. 1) has received less attention and currently lacks any absolute chronologies of Pleistocene glaciations, despite a well-preserved geomorphic record of glaciations (Gaskill et al., 1967, 1986, 1991). Thus, the geologic record of recent Pleistocene glaciations preserved in the Elk Range of the Colorado Rockies presents a unique opportunity to further our understanding of regional glacier and Pleistocene climate change.

Figure 1. (A) Regional elevation map of Colorado indicating the locations of mountain ranges discussed in the text. The location of the East River watershed is denoted by the red circle. Inset map showing of North and Central America with Colorado outlined in red. (B) Hill shade map overlying high-resolution elevation data for the East River watershed study area. The reconstructed LGM extent for the East River glacier is shown in white and double arrows indicate locations of potential ice-flow over drainage divides. Figure S1 shows the sites targeted for cosmogenic exposure and depth profile dating in detail. Colored lines show the locations of mapped moraine crests. CB = Mount Crested Butte, SG = Snodgrass Mountain, GM = Gothic Mountain, WG = Washington Gulch. The longitude and latitude indicate the locations of Castle Peak (northern crosshair) and Mount Emmons (southern crosshair).

The East River watershed, a headwater tributary of the Colorado River within the Elk Mountains, is a United States Department of Energy Watershed Function Scientific Focus Area (SFA) that serves as a natural laboratory for furthering our understanding of watershed function and biogeochemical interactions (Hubbard et al., 2018). Research in the SFA is focused on improving our understanding of how mountainous watersheds store, partition, and release solutes into the environment under variable climate regimes (Hubbard et al., 2018). Because the development of soils, subsurface weathering fronts, and other aspects of the critical zone that influence watershed function are directly tied to the timing of deglaciation, it is necessary to establish a detailed chronology of post-LGM ice retreat within the East River watershed.

We used exposure and depth profile dating using in situ-produced cosmogenic ¹⁰Be to establish the first glacial chronology for the East River watershed. We then used the glacial chronology to reconstruct local late Pleistocene climate using a coupled energy–mass balance and ice-flow modeling approach. The glacier chronologies and climate modeling provide quantitative estimates of latest Pleistocene climate conditions for the central Rocky Mountains and can aid investigations of subsurface biogeochemical process rates during the time following the last glaciation.

Study area

The East River watershed is located in the Elk Range in the central Colorado Rocky Mountains (~38.9°N, 106.9°W). The watershed is west of the continental divide in the headwaters of the Gunnison River, a major tributary of the Colorado River. Bedrock in the watershed is primarily Paleozoic and Mesozoic sedimentary formations that have been intruded by Paleocene- to Miocene-age laccoliths and stocks with quartz monzonite porphyry and granodiorite composition (Gaskill et al., 1967, 1991). Peaks in the Elk Range are associated with igneous intrusions and often exceed 4000 m in elevation, whereas the elevation of the East River valley bottom is near 2700 m. The Elk Range is characterized by extensive glacial and periglacial features, including cirques, U-shaped valleys, moraine deposits, polished and striated bedrock, and active rock glaciers (Gaskill et al., 1991). The modern continental, subarctic climate of the East River watershed is characterized by short, cool summers and long, cold winters (Hubbard et al., 2018). Mean annual temperature within the watershed is around 0°C and the higher elevation areas receive approximately 0.7–1.2 m of precipitation per year on average, the majority of which falls as snow (Hubbard et al., 2018; Fang et al., 2019).

Quaternary sediments in the watershed are dominated by moraine, till, glacial outwash, and landslide deposits. Gaskill et al. (1967, 1986, 1991) characterized and mapped the moraine deposits targeted for exposure dating in this study. We collected samples for cosmogenic exposure and depth–profile dating at six sites from across the study area, including the East River terminal moraine complex, moraines at Washington Gulch south and north, moraines at the mouth and headwaters of Copper Creek, and striated bedrock at Schofield Pass (Fig. 1; LIDAR-derived hill-shade maps with exposure age sample locations are shown in Fig. S1).

The East River terminal moraine complex is located to the southeast of and across the river valley from Mount Crested Butte and is comprised of at least three curvilinear ridges (Fig. 1; Fig. S1). We collected samples for exposure dating from boulders deposited on the two outermost moraine crests preserved on the south side of the valley. Moraine deposits on the north side of the valley likely have been eroded by the East River (Fig. 1; Fig. S1). Till deposits immediately outboard of the terminal moraine complex that we sampled have been mapped as likely corresponding to the penultimate Bull Lake glaciation (Gaskill et al., 1986).

Moraine deposits indicate that ice flowed laterally to the southwest out of the East River Valley and into Washington Gulch at two locations where there are topographic lows on the drainage divide (Gaskill et al., 1967, 1991). The Washington Gulch south site is associated with an ice lobe that exited the East River valley through the topographic gap between Mount Crested Butte and Snodgrass Mountain (Fig. 1). The moraine is characterized by at least two nearly continuous curvilinear, latero-end moraine ridges (Fig. 1; Fig. S1) that bound broad flat topographic surfaces. The Washington Gulch north site was similarly deposited by ice flowing through the pass between Snodgrass and Gothic mountains (Fig. 1). The moraine complex is exceptionally well preserved, featuring two continuous latero-end moraine ridges with a broad flat area separating the two ridges (Fig. S1).

The Gothic moraine site is located in the Gothic townsite near the mouth of Copper Creek (Fig. 1) and is composed of two curvilinear ridges that can be traced to lateral moraine counterparts along the northern walls of the Copper Creek valley (Fig. S1). Portions of the moraine have been eroded by the south-flowing East River. Farther up the Copper Creek valley, we sampled a set of boulder-rich linear till ridges (Fig. 1; Fig. S1), which we interpreted as moraines (Copper Lake moraine site). Near the headwaters of the East River in Schofield Pass, we observed a broad, flat outcrop of Dakota Sandstone quartzite with pervasive glacially polished surfaces and areas with clear striations parallel to the primary south-southeast ice-flow direction (Figs. 1, S1).

Methods

Cosmogenic ¹⁰Be exposure dating

Moraine boulders and glacially striated bedrock in the western Elk Range were targeted for in-situ-produced cosmogenic ¹⁰Be exposure dating (Fig. S1). We primarily sampled quartz monzonite and quartzite boulders with broad horizontal surfaces that were deposited on moraine crests. Boulders sampled for cosmogenic exposure dating were, on average, 53 ± 23 cm (1σ; Table S1.A) tall (as measured above the ground surface). At the Schofield Pass site, we collected two bedrock samples with glacial striae and polish. Samples were collected using a portable rotary hammer drill, sledgehammer, and chisel. Sample thickness and mass ranged from 2–4 cm and 1–3 kg, respectively. Topographic shielding data were collected in the field by measuring clear-sky obstructions (e.g., valley walls, mountain peaks, etc.) as the angle above the horizon at 20° azimuthal increments.

Twenty-seven samples were prepared for in-situ cosmogenic ¹⁰Be measurement at the University of Massachusetts Amherst Cosmogenic Nuclide Lab. Samples were crushed and sieved to a target grain size of 250–850 μm. Quartz grains were isolated using a rare earth hand magnet, Outotec magnetic roll separator, froth flotation, and dilute acid treatment. The quartz purification process was accomplished by repeated etching in dilute hydrofluoric acid (Kohl and Nishiizumi, 1992). We tested aliquots of each sample to determine quartz purity via inductively coupled plasma–optical emission spectroscopy (ICP–OES) analysis, including measurements of Al, Be, Ca, K, Mg, Na, Fe, and Ti. We used a ⁹Be carrier solution to add ~250 μg of Be to each sample and procedural blank. The Be fraction of each sample was chemically isolated using ion chromatography following Ditchburn and Whitehead (1994) and Stone (2001), calcined by heating in a quartz crucible with a gas torch until BeO glowed orange for 1 minute, and packed with Nb powder into targets for ¹⁰Be/⁹Be measurement by accelerator mass spectrometry (AMS) at the Lawrence Livermore National Laboratory. The ratio of ¹⁰Be/⁹Be was measured against AMS beryllium standard 07KNSTD with a known ratio of 2.85 × 10⁻¹² (Nishiizumi et al., 2007).

Sample ¹⁰Be/⁹Be ratios ranged from 1.41 × 10⁻¹³ to 7.92 × 10⁻¹³ (Table 1). The ¹⁰Be/⁹Be ratios in procedural blanks ranged from 4.15 × 10⁻¹⁶ to 5.28 × 10⁻¹⁶ (Table 1). Therefore, the blank corrections were minimal because less than 1% of the measured ¹⁰Be is attributed to the blank. The range of AMS measurement uncertainties (1σ) for all samples was approximately 1–2%.

Table 1. In situ-produced ¹⁰Be sample data

*Latitude, longittude, and elevation were measured with a decimeter-scale precision GPS (Bad Elf Flex with Satellite-based Augmentation System (SBAS) differential correction).

All samples use AMS standard 07KNSTD

We calculated exposure ages using the Balco et al. (2008) online exposure age calculator, version 3.0 (http://hess.ess.washington.edu/math/). The online calculator implements several algorithms for scaling cosmogenic nuclide production with respect to geography and time, and we report results from the nuclide-specific, time-dependent scaling model of Lifton et al. (2014). Additionally, we used the ¹⁰Be global production rate calibration dataset presented by Borchers et al. (2016) to facilitate comparison against published glacial chronologies calculated in a similar manner. Input data used for the exposure age calculations are included in Tables S1.A and S1.B. We calculated exposure ages assuming no erosion of the boulder surfaces and demonstrated that this assumption had limited effect on calculated exposure ages and interpretations by recalculating ages using a modest erosion rate of 1 mm kyr⁻¹. The difference in calculated exposure ages using zero and 1 mm kyr⁻¹ erosion rates is less than 2% (Table S1.C). Exposure ages from previous work that are presented here have been recalculated using identical methods.

Landform ages (i.e., moraine ages) were calculated as the arithmetic mean of all samples, excluding outliers, from a given landform, and we report landform age uncertainty as the standard error of the mean. Uncertainty associated with scaling schemes and production-rate calibration, which is likely on the order of 3–4%, is suggested to be systematic with little geographic bias (Balco, 2020). We report external uncertainty in Table 2 to facilitate comparison of the exposure ages we report against those from different chronometers. Further, to assess how well analytical uncertainty alone explains the observed variance, we calculated a reduced chi-squared statistic (χ²R) from the distribution of ages for each landform (Table 2). We interpret χ²R values greater than 1 to indicate that the observed variance cannot be explained solely by analytical errors, and instead must include other sources of error (e.g., Schaefer et al., 2009; Rood et al., 2011).

Table 2. In situ-produced ¹⁰Be exposure age results

†Landform age uncertainty calculated as the standard error of the mean. The standard error in ages for the Washington Gulch South site is zero because the two ages are the same; however we report and uncertainty of 0.1 ka for this site.

*All ages calculated using version 3.0 of the Balco et al. (2008) calculator with the Lifton-Sato-Dunai time dependent scaling model (Lifton et al., 2014) and the primary ¹⁰Be calibration data set (Borchers et al., 2016), and are rounded to the nearest 100 years.

Cosmogenic ¹⁰Be depth profile

We collected five samples from a soil pit on the Gothic moraine for in situ cosmogenic ¹⁰Be measurements (Fig. S2). We opted to use an in-situ-produced ¹⁰Be depth profile due to the lack of suitable boulders on the Gothic moraine crest. Samples and corresponding bulk density measurements were collected from 25–145 cm below the moraine crest surface. We used the Monte Carlo algorithm developed by Hidy et al. (2010) to simultaneously estimate denudation rates, cosmogenic ¹⁰Be inheritance, and landform exposure age. The model also allows for incorporation of independent geologic information, such as previously established age and erosion constraints, into the calculation of these parameters.

We used a site-specific, time-averaged (10 and 20 ka) spallation surface production rate of 32.7 atoms g⁻¹ yr⁻¹, which was calculated using the same nuclide-specific, time-dependent scaling model as for exposure age calculations (Lifton et al., 2014) and includes the effects of topographic shielding. In all simulations, we used a depth-averaged bulk density of 1.9 g cm⁻³ based on measurements of the till matrix and the abundance of clasts in the profile. Details of the input parameters are provided in Table S2.

¹⁰Be/⁹Be ratios for the Gothic moraine depth profile samples ranged from 7.53 × 10⁻¹⁴ at 145–150 cm depth to 2.65 × 10⁻¹³ near the surface at 25–30 cm depth. The samples were prepared alongside a procedural blank with a ratio of 3.74 × 10⁻¹⁶ (Table S2). The measured ¹⁰Be/⁹Be ratios for all five samples preserve stratigraphic-¹⁰Be concentration order within analytical uncertainty. However, we excluded the 115–120 cm depth sample from our analysis because the measured concentration (1.06 × 10⁵ ± 2.2 × 10³ atoms ¹⁰Be g⁻¹) had a very similar concentration to the next lowest sample from 145–150 cm depth (9.86 × 10⁴ ± 2.4 × 10³ atoms ¹⁰Be g⁻¹), which resulted in a departure from the theoretical (and otherwise observed) exponential decay of concentration with depth and made it difficult to fit modeled profiles to the measurements.

The Monte Carlo model inputs included a range of likely limiting estimates of moraine exposure age and total permissible erosion. Preliminary results from a simulation where landform age and total erosion were allowed to vary without constraints indicated that model solutions did not include any erosion rates > 20 cm ka⁻¹. Thus, we limited maximum erosion to 20 cm ka⁻¹ for all experiments. All experiments were run until 10,000 simulation solutions were found within a 3σ uncertainty envelope. We ran multiple models with setups that included combinations of unconstrained exposure age (0–50 ka), unconstrained maximum erosion (10 meters), constrained exposure age (12.9–18.0 ka), and constrained erosion threshold (0–40 cm). Limits on the moraine exposure age are based on ages for other landforms that we exposure dated. The estimate of 40 cm for maximum erosion threshold is a conservative estimate based on the observation of maximum possible soil erosion around boulders, even though there are no clear signs of boulder exhumation in the exposure age data (i.e., Putkonen and Swanson, 2003; Applegate et al., 2012).

Glacier and climate reconstructions

The coupled energy–mass balance and ice-flow models used in this study were originally developed by Plummer and Phillips (2003) and have been applied to estimate paleoclimate conditions for paleoglaciers in a variety of settings (Harrison et al., 2014; Rowan et al., 2014; Leonard et al., 2017; Laabs et al., 2006; Quirk et al., 2018, 2020). The energy–mass balance model calculates monthly snow accumulation and ablation across the modeling domain using a digital elevation model (DEM) and monthly climate parameters (e.g., precipitation, temperature, wind speed, etc.) as inputs. The monthly mass balance data are used to calculate an annual mass balance distribution, effectively defining the annual average zones of accumulation and ablation, which are used as inputs to the ice-flow model along with the DEM and associated ice-flow parameters. The ice-flow model used here is similar to the finite-element ice-sheet model of Fastook and Chapman (1989) and uses the common shallow-ice approximation. The iterative modeling approach involves calculating a mass-balance distribution across the topographic surface for a given climate scenario, calculating the ice-flow and resulting ice-thickness distributions, and comparing to known, mapped ice extents based on moraines. In subsequent simulations, climate parameters such as precipitation and temperature are varied to determine the best fit with mapped ice extents. Importantly, because we matched modeled ice extents to moraine deposits in this study, the ice-flow model was run until the simulated glacier reached steady-state.

We primarily focused on matching two paleoglacier extents: the East River terminal moraine complex and the Gothic recessional moraine. The model configuration used here is nearly identical to that presented by Quirk et al. (2018, 2020), including the use of PRISM (Parameter-elevation Regression on Independent Slopes Model) (Daly et al., 2008) mean monthly precipitation and temperature (1971–2000 CE) values as inputs to the energy–mass balance model. A detailed description of additional climate parameters used as model input and data sources are in Table S3. We used the model uncertainty as determined by Quirk et al. (2020) via sensitivity analyses of ± 1.0°C and 30% for precipitation. We applied uniform temperature depressions and precipitation factors to the entire monthly input data.

Results

Cosmogenic ¹⁰Be exposure ages

The exposure age of the East River terminal moraine complex is 17.8 ± 0.1 ka (n = 10), with no statistical difference in age between the two outermost moraine crests (Fig. 2). The Washington Gulch south and north moraines exhibit landform exposure ages of 17.3 ± 0.1 (n = 2) and 17.9 ± 0.2 ka (n = 10), respectively. Glacially polished and striated quartzite bedrock samples from near Schofield Pass have an exposure age of 14.9 ± 0.8 ka. The moraine ridge in the headwaters of Copper Creek yielded the youngest age found in our study of 13.0 ± 0.1 ka (n = 3).

Figure 2. Kernel density estimates of cosmogenic ¹⁰Be exposure ages for each landform in the study. The thin, gray curves show individual exposure ages, whereas the bold curves in color show the cumulative distribution for each landform.

Cosmogenic ¹⁰Be depth profile

We conducted depth–profile modeling to limit the likely range of values for erosion rate, exposure age, and inheritance (Fig. 3). The completely unconstrained simulation, where age and total-erosion threshold are allowed to vary over a large range of values, generally indicates low erosion rates and a latest Pleistocene age for the Gothic moraine with median values of 7.0 cm ka⁻¹ and 29.6 ka, respectively. The stratigraphic position of the Gothic moraine indicates that it must be younger than the East River terminal moraine and Washington Gulch end moraine and that it must be older than the moraines in the headwaters of Copper Creek. Thus, we used the chronological data to constrain age limits used in the depth profile calculations. Experiments constraining either age (12.9–18.0 ka; using bracketing exposure ages presented here) or the conservatively estimated erosion threshold (0–40 cm) yielded median erosion rates and exposure ages ranging from 1.6–3.6 cm ka⁻¹ and 12.7–15.4 ka, respectively (Fig. 3B). Constraining both age and the erosion threshold yielded a median erosion rate and age of 2.3 cm ka⁻¹ and 13.5 ka, respectively. Predicted inheritance values are relatively consistent across all simulations with median values ranging from 2.5–2.8 × 10⁴ atoms per gram of quartz. These results suggest ice retreated from the Gothic moraine around 13–15 ka. For the purposes of glacier modeling, the differences in insolation across this narrow time interval have only modest effects on our results and are within model uncertainty. Model results for all four simulations are included in Table S1.C.

Figure 3. (A) Modeled ¹⁰Be concentration depth profiles results (red curves) fit to the Gothic moraine soil pit sample in-situ-produced ¹⁰Be concentrations (blue circles). Model depth profiles do not substantially vary across our simulation parameters. The single outlier sample collected from 115 cm depth (indicated by light-gray square) was not included in the model fit. All error bars indicate 3σ uncertainty; x-axis concentration ×10⁵. (B) Modeled exposure ages (red points) plotted against erosion rates for simulation with unconstrained age limits (0–50 ka) and 40 cm maximum erosion threshold. (C) Same as in (B) but for simulation constraining both age (12.8–18 ka; dashed horizontal lines) and 40 cm maximum erosion threshold. The green circle highlights the median solution with exposure age and erosion rate of 13.5 ka and 2.3 cm ka⁻¹, respectively. In both (B) and (C) the best 100 chi-square fits are highlighted as blue points, and the 40 cm maximum erosion threshold is indicated by a black curve.

Glacier modeling

Model results indicate temperature depressions and precipitation change combinations for the stadia represented by the East River terminal moraine complex range from −8.6 to −2.2°C and 50–300%, respectively (Fig. 4). Similarly, for the Gothic moraine stadia these values range from −6.3 to −1.5°C and 50–200%, respectively (Fig. 4). Model solutions for each of the stadia are presented in full in Table S4.

Figure 4. Temperature depression and precipitation factor combinations that satisfactorily reproduce the East River maximum ice extent (black dots) and recessional ice position at the Gothic moraine (red dots). Modeled exponential fits, corresponding equations, and R ² values for the East River terminal and Gothic moraine positions are shown in black and red, respectively.

Modeled glacier extents match the geomorphic constraints for both the maximum terminal and recessional Gothic moraine positions (Fig. 5). Simulations matching the East River terminal moraine positions generally agree with the mapping and chronology presented here and predict the presence of the outlet glaciers that deposited the Washington Gulch north and south moraines. However, the model appears to simulate a slight excess of ice flowing down the main Washington Gulch tributary, which coalesces with the Washington Gulch north moraine lobe (Fig. 5).

Figure 5. Predicted ice thickness for the East River maximum (top) and Gothic recessional (bottom) glacier terminus positions at 17–18 and 13–15 ka, respectively.

Discussion

Chronology of Pinedale Glaciation in the East River watershed

Our findings provide new insight on the timing of glaciation and ice-retreat during the last glacial episode in the Elk Range. In the context of the Rocky Mountain glacial chronology, all the landforms investigated here belong to the most recent Pinedale glaciation (Blackwelder, 1915) which is broadly correlative with Marine Oxygen Isotope Stage 2 (MIS 2). Unless stated otherwise, landform ages presented here are interpreted to represent the time of moraine abandonment by retreating ice. The assumption that the ages represent abandonment, in lieu of data suggesting otherwise, is reasonable considering that the boulders we sampled were at the tops of moraine crests and are likely the youngest glacial deposits within the moraine deposit and thus closest to deglaciation in age.

The terminal moraine complex in the East River valley and the moraines in Washington Gulch exhibit similar exposure ages and suggest that these three positions were abandoned by ice around the same time at ca.17–18 ka. Similar to mountain ranges across western North America, ice retreat from the maximum extent in the East River watershed appears to have occurred shortly after the end of the global LGM. Based on exposure age chronology alone, it is difficult to determine the relative timing of ice retreat up-valley to Schofield Pass (14.9 ± 0.8 ka) and abandonment of the Gothic recessional moraine (13–15 ka). However, the modeled glacier reconstructions (Fig. 5) indicate that the Schofield Pass site was ice covered when the glacier that occupied Copper Creek terminated at the Gothic moraine. However, given that the East River glacier was likely not in steady-state as it retreated above the striated bedrock at Schofield Pass, this result must be interpreted with caution. In any case, the available data suggest that much of the East River valley was deglaciated and that the main glacier had retreated and separated into several discrete ice bodies by around 13–15 ka. Lastly, the exposure age for the Copper Lake moraine indicates that ice in the catchment had retreated above 3400 m in elevation to a shaded cirque position by ca. 13.0 ka.

Comparison to regional glacier chronologies

Numerous cosmogenic exposure chronologies for the most recent glaciation of the central and southern Rockies have been developed over the last two decades. These include records from the Sawatch and Mosquito ranges located less than 100 km northeast of the East River watershed (Brugger, 2007; Briner, 2009; Young et al., 2011; Brugger et al., 2019a, b; Schweinsberg et al., 2020). Terminal moraine ages vary from 16.8–22.9 ka in the Sawatch Range. Similarly, in the Mosquito Range, terminal moraines vary in age from 17.6–22.3 ka. In the East River, maximum glacier positions appear to have been occupied until ca. 17–18 ka. However, we did not observe terminal moraine ages that coincided with the end of the global LGM as observed in the nearby Mosquito and Sawatch ranges. One hypothesis that could explain this observation is that the East River glacier occupied the same or similar ice-extent throughout much of last glaciation, and only the most recent recession around 17–18 ka is recorded. Alternatively, ice could have been positioned farther up the East River valley prior to advancing to the terminal position shortly before 17–18 ka, thus overriding any previous deposits. Yet another explanation for the apparent lack of LGM-aged moraines in the East River is that they lie beyond the margin of the East River terminal moraine. For example, the moraines immediately outboard of the East River terminal moraine complex previously were mapped by Gaskill et al. (1986) as penultimate or Bull Lake in age. However, numerical age limits for these moraines have not yet been established. Our cursory investigation of these moraines did not yield conclusive indications, such as relative weathering indicators, of the landform age, although a dramatic reduction in boulder size and frequency was noted. Thus, it remains possible that glacial deposits coinciding with the end of the LGM are present within the East River watershed.

In the Front and Park ranges to the north of the East River (Benson et al., 2004; Ward et al., 2009; Dühnforth and Anderson, 2011), terminal moraines exhibit a slightly older age distribution, from 18.6–22.8 ka, than observed in the Elk, Sawatch, and Mosquito ranges. South of the Elk Range, in the San Juan Mountains (Benson et al., 2004) and Sangre de Cristo Mountains (Leonard et al., 2017), terminal moraines vary in age from 19.1–20.1 and 18.2–19.2 ka, respectively. Thus, the data suggest that retreat from ice-distal positions in the Elk, Sawatch, and Mosquito ranges spanned a greater range of time, and tend to be younger, than in other ranges in the Colorado Rocky Mountains to the south and north. These observed differences could reflect a spatially limited paleoclimatic change (e.g., localized decreases in summer temperatures or increases in winter precipitation). Alternatively, the ice dynamics of individual glaciers or glacier hypsometry could influence the observed age distributions.

To understand the chronology of ice retreat in the Elk Mountains, we calculated normalized length and normalized elevation change of the glacier toe along the flow path from the Copper Creek headwaters (0) to the terminal moraine complex (1). Similarly, we calculated normalized, headwater to terminus, glacier elevation and length changes for several ranges in the Colorado Rocky Mountains with existing chronologies (Fig. 6). Assessing both changes in glacier length and toe elevation can help to better understand the effects of hypsometry on glacier response at different sites (Shakun et al., 2015). For example, our data indicate a substantial length decrease of approximately 60% during the early phases of deglaciation while the elevation data indicate only an approximately 15% reduction from the East River terminal to the Gothic moraine sites. This observation is consistent with documented glacier responses globally and is likely controlled by the relationship between catchment slope and temperature change (Oerlemans, 2005; Shakun et al., 2015) where glaciers such as the East River flow from areas of high slopes in the headwater to lower slopes near the terminus.

Figure 6. Colorado Rocky Mountain glacier changes over time shown as normalized glacier elevation (A) and normalized glacier length (B), relative to glacier headwater (0) and terminus (1). The data suggest variable initial deglaciation ages from maximum positions (i.e., normalized elevation and length = 1), followed by relatively consistent regional retreat to less than 20% maximum lengths and to high elevations by ca. 13–15 ka.

Our data suggest glacier length and elevation retreat rates of approximately 725 and 46 m kyr⁻¹, respectively, in the west Elk Range from ca. 18–14 ka. Glacier elevation retreat rates for the ca. 14–13 ka time interval accelerated compared to the ca. 18–14 ka period by an order of magnitude to more than 500 m kyr⁻¹ while length changes slowed to just under 500 m kyr⁻¹. The temporal pattern of ice retreat in the Elk Range is similar to that observed in the nearby Sawatch Range (Young et al., 2011; Tulenko et al., 2020) and in the San Juan (Guido et al., 2007) and Sangre de Cristo mountains (Leonard et al., 2017) to the south as well as in the Front Range (Ward et al., 2009; Dühnforth and Anderson, 2011) to the northeast (Fig. 6). Generally, we observe the same deglacial pattern as Quirk et al. (2022) identified, with variable initial deglaciation ages (ca. 20–16 ka) and nearly uniform subsequent deglaciation after ca. 16 ka.

Glacier model paleoclimate interpretations

The glacier reconstructions provide new insight on the spatiotemporal pattern of glaciation and ice-extents in the Elk Range as well as the climate that produced them. Glaciers predominantly respond to changes in temperature during the melt or ablation season (i.e., summer) and accumulation season (i.e., winter) precipitation. Thus, while our results for the East River glacier and the results of others are reported as mean annual temperature and precipitation changes, it is likely that the results more closely reflect seasonal changes in these parameters.

Brugger (2010) reconstructed the maximum extent of the East River Pinedale-age glacier using detailed field mapping and estimated the necessary climate conditions to sustain the mapped ice extent using a degree-day model. The modeled reconstructions presented here match the field-based mapping and reconstructions of Brugger (2010) well. For example, our model predicts the concurrent ice margins at the terminal moraine and the moraines in Washington Gulch and is consistent with the chronology developed here. Further, both the energy–mass balance model used here and the Brugger (2010) degree-day model indicate that, with little to no change in precipitation, glaciers in the Elk Range could have been sustained by temperature depression around −7°C (our result: −6.5 ± 1.0°C). General circulation model ensemble averages, interpolated for the Sawatch and Elk Range (Oster et al., 2015), substantiate the assumption of only modest change in LGM precipitation (112 ± 47% of modern). The degree-day and coupled ice flow–mass balance models have not been directly compared before, and although limited to a single site, our results suggest the two methods yield comparable results.

Reconstructed glacial systems from across the central Rocky Mountains suggest similar climatic changes during the last glaciation. To the north in the Sawatch Range, slightly greater temperature depressions of around −8 to −9°C were required to maintain maximum Pinedale glacier extents (Brugger et al., 2019b). In the Mosquito Range, temperature depressions of approximately −7.5 to −8°C were required to maintain maximum ice extents during the last glacial (Brugger et al., 2019a). In the Sangre de Cristo to the south, Leonard et al. (2017) and Brugger et al. (2021) found Pinedale maximum glaciers required approximately −5 and −8.5°C temperature changes relative to present, respectively—possibly underscoring the importance of local precipitation variability (e.g., micro-climate) or ice-dynamic (e.g., ice rheology, transient flow) effects. Comparisons are further complicated by the fact that reconstructions enumerated here are not necessarily coeval with one another, because the timing of maximum ice extent varies by several thousand years across the Rocky Mountains in Colorado. However, in aggregate, the glacial reconstructions based on our work and prior studies provide a broadly consistent record of temperature change ca. 21–17 ka during local Pinedale maxima in the central Rocky Mountains and indicate temperature depressions were about −7 to −9°C.

Recessional ice configurations are reconstructed less often than maximum Pinedale glacier extents. As such, there are fewer data and proxy records with which we can compare our modeling results for the Gothic moraine stadial between 13 and 15 ka. However, if we assume that precipitation did not vary drastically between the Pinedale maxima (ca. 17–18 and 13–15 ka), our results suggest a temperature increase of approximately 2.5°C (0.5–4.5°C, accounting for model uncertainty) over this time interval. Alternatively, if precipitation increased, for example from Heinrich Event 1 (Munroe and Laabs, 2013), and drove the culmination of this stadial, temperature increases could have been greater. For example, a temperature increase of 5°C (3–7°C, accounting for model uncertainty) would be required had precipitation increased two-fold at this time.

Implications for weathering zone development in the East River valley

Predicting the depth, rates, and extent of bedrock weathering is an ongoing focus of research in the East River watershed due to the role that subsurface chemical weathering plays in solute generation and nutrient export (e.g., Winnick et al., 2017; Wan et al., 2019, 2021). At a hillslope with Mancos Shale bedrock located near the valley floor and downvalley from the Gothic moraine, weathering fronts for soluble minerals (pyrite, calcite, dolomite) are 3–4 m below the ground surface (Wan et al., 2019, 2021). Assuming deep glacial erosion of the valley floor, these weathering fronts have propagated to their current depth since deglaciation 13–15 ka at rates on the order of 0.2–0.3 mm yr⁻¹. Such rates are minimum estimates because they do not account for post-glacial erosion of the hillslope. Further, the modeled glacier reconstructions predict that many areas remained ice-free, such as the low-elevation summit of Snodgrass Mountain (3400 m), which may have experienced limited glacier erosion during the Pleistocene. Hence, in addition to bedrock properties and topography (e.g., Uhlemann et al., 2022), factors such as the depth of glacial erosion and exposure to periglacial weathering regimes may influence the character of the regolith and weathering zone within the watershed.

Conclusions

¹⁰Be exposure dating indicates the East River glacier abandoned maxima positions at the East River terminal and Washington Gulch moraines around 17–18 ka. Exposure dating and a ¹⁰Be depth profile indicate ice retreat from a recessional moraine near Gothic, Colorado, occurred ca. 13–15 ka and that ice retreated farther to high-elevation cirques by ca. 13 ka. Exposure ages from terminal moraines in the East River are somewhat younger than those from neighboring mountain ranges in Colorado. However, our results are consistent with the relatively rapid deglaciation after 16 ka observed throughout the Rocky Mountains. Our glacier reconstructions place plausible limits on climate conditions at ca. 17–18 and 13–14 ka, with estimated annual or summer temperature depressions of −6.5°C and −4.0°C, respectively, assuming no change in precipitation relative to present. These results are consistent with other glacier reconstructions from the central Rocky Mountains and general circulation model simulations for the LGM. The deglacial chronology also provides constraints on the rates of weathering-front propagation and development of the regolith within the East River watershed.