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, Reference Dyurgerov and Meier2000; Oerlemans, Reference Oerlemans2005; Zhang et al., Reference Zhang, Kang, Gabrielli, Loewen and Schwikowski2015; Roe et al., Reference Roe, Baker and Herla2017). 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, Reference Agassiz and Bettannier1840; Gilbert, Reference Gilbert1890; Putnam et al., Reference Putnam, Schaefer, Denton, Barrell, Andersen, Koffman and Rowan2013; Shakun et al., Reference Shakun, Clark, He, Lifton, Liu and Otto-Bliesner2015; Marcott et al., Reference Marcott, Clark, Shakun, Brook, Davis and Caffee2019), 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, Reference Gilbert1890; Ray, Reference Ray1940; Leonard, Reference Leonard1989; Murray and Locke, Reference Murray and Locke1989; Locke, Reference Locke1990, Reference Locke1995; Pierce, Reference Pierce2003; Quirk et al., Reference Quirk, Huss, Leonard, Licciardi, Plummer and Caffee2022). 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., Reference Benson, Madole, Phillips, Landis, Thomas and Kubik2004; Brugger, Reference Brugger2007; Guido et al., Reference Guido, Ward and Anderson2007; Briner, Reference Briner2009; Ward et al., Reference Ward, Anderson, Guido and Briner2009; Dühnforth and Anderson, Reference Dühnforth and Anderson2011; Young et al., Reference Young, Briner, Leonard, Licciardi and Lee2011; Leonard et al., Reference Leonard, Laabs, Plummer, Kroner, Brugger, Spiess, Refsnider, Xia and Caffee2017; Brugger et al., Reference Brugger, Laabs, Reimers and Bensen2019a, Reference Brugger, Ruleman, Caffee and Masonb; Schweinsberg et al., Reference Schweinsberg, Briner, Licciardi, Shroba and Leonard2020; Tulenko et al., Reference Tulenko, Caffee, Schweinsberg, Briner and Leonard2020). 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., Reference Clark, Dyke, Shakun, Carlson, Clark, Wohlfarth, Mitrovica, Hostetler and McCabe2009) 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, Reference Licciardi and Pierce2008, Reference Licciardi and Pierce2018; Quirk et al., Reference Quirk, Moore, Laabs, Plummer and Caffee2020, Reference Quirk, Huss, Leonard, Licciardi, Plummer and Caffee2022) and western North America (Laabs et al., Reference Laabs, Licciardi, Leonard, Munroe and Marchetti2020, and references therein). The glacier climate reconstructions available for the Rocky Mountains of Colorado (Brugger, Reference Brugger2010, Reference Brugger, Laabs, Reimers and Bensen2019a, Reference Brugger, Ruleman, Caffee and Masonb; Leonard et al., Reference Leonard, Laabs, Plummer, Kroner, Brugger, Spiess, Refsnider, Xia and Caffee2017; Brugger et al., Reference Brugger, Leonard, Refsnider and Dolan2021) 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., Reference Gaskill, Godwin and Mutschler1967, Reference Gaskill, Colman, DeLong and Robinson1986, Reference Gaskill, Mutschler, Kramer, Thomas and Zahony1991). 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.
(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., Reference Hubbard, Williams, Agarwal, Banfield, Beller, Bouskill and Brodie2018). 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., Reference Hubbard, Williams, Agarwal, Banfield, Beller, Bouskill and Brodie2018). 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 10Be 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., Reference Gaskill, Godwin and Mutschler1967, Reference Gaskill, Mutschler, Kramer, Thomas and Zahony1991). 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., Reference Gaskill, Mutschler, Kramer, Thomas and Zahony1991). The modern continental, subarctic climate of the East River watershed is characterized by short, cool summers and long, cold winters (Hubbard et al., Reference Hubbard, Williams, Agarwal, Banfield, Beller, Bouskill and Brodie2018). 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., Reference Hubbard, Williams, Agarwal, Banfield, Beller, Bouskill and Brodie2018; Fang et al., Reference Fang, Carroll, Schumer, Harman, Wilusz and Williams2019).
Quaternary sediments in the watershed are dominated by moraine, till, glacial outwash, and landslide deposits. Gaskill et al. (Reference Gaskill, Godwin and Mutschler1967, Reference Gaskill, Colman, DeLong and Robinson1986, Reference Gaskill, Mutschler, Kramer, Thomas and Zahony1991) 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., Reference Gaskill, Colman, DeLong and Robinson1986).
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., Reference Gaskill, Godwin and Mutschler1967, Reference Gaskill, Mutschler, Kramer, Thomas and Zahony1991). 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 10Be exposure dating
Moraine boulders and glacially striated bedrock in the western Elk Range were targeted for in-situ-produced cosmogenic 10Be 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 10Be 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, Reference Kohl and Nishiizumi1992). 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 9Be 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 (Reference Ditchburn and Whitehead1994) and Stone (Reference Stone2001), 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 10Be/9Be measurement by accelerator mass spectrometry (AMS) at the Lawrence Livermore National Laboratory. The ratio of 10Be/9Be was measured against AMS beryllium standard 07KNSTD with a known ratio of 2.85 × 10−12 (Nishiizumi et al., Reference Nishiizumi, Imamura, Caffee, Southon, Finkel and McAninch2007).
Sample 10Be/9Be ratios ranged from 1.41 × 10−13 to 7.92 × 10−13 (Table 1). The 10Be/9Be ratios in procedural blanks ranged from 4.15 × 10−16 to 5.28 × 10−16 (Table 1). Therefore, the blank corrections were minimal because less than 1% of the measured 10Be is attributed to the blank. The range of AMS measurement uncertainties (1σ) for all samples was approximately 1–2%.
In situ-produced 10Be 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. (Reference Balco, Stone, Lifton and Dunai2008) 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. (Reference Lifton, Sato and Dunai2014). Additionally, we used the 10Be global production rate calibration dataset presented by Borchers et al. (Reference Borchers, Marrero, Balco, Caffee, Goehring, Lifton, Nishiizumi, Phillips, Schaefer and Stone2016) 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−1. The difference in calculated exposure ages using zero and 1 mm kyr−1 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, Reference Balco2020). 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 (χ2R) from the distribution of ages for each landform (Table 2). We interpret χ2R 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., Reference Schaefer, Denton, Kaplan, Putnam, Finkel, Barrell and Andersen2009; Rood et al., Reference Rood, Burbank and Finkel2011).
In situ-produced 10Be 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 10Be calibration data set (Borchers et al., 2016), and are rounded to the nearest 100 years.
Cosmogenic 10Be depth profile
We collected five samples from a soil pit on the Gothic moraine for in situ cosmogenic 10Be measurements (Fig. S2). We opted to use an in-situ-produced 10Be 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. (Reference Hidy, Gosse, Pederson, Mattern and Finkel2010) to simultaneously estimate denudation rates, cosmogenic 10Be 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−1 yr−1, which was calculated using the same nuclide-specific, time-dependent scaling model as for exposure age calculations (Lifton et al., Reference Lifton, Sato and Dunai2014) and includes the effects of topographic shielding. In all simulations, we used a depth-averaged bulk density of 1.9 g cm−3 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.
10Be/9Be ratios for the Gothic moraine depth profile samples ranged from 7.53 × 10−14 at 145–150 cm depth to 2.65 × 10−13 near the surface at 25–30 cm depth. The samples were prepared alongside a procedural blank with a ratio of 3.74 × 10−16 (Table S2). The measured 10Be/9Be ratios for all five samples preserve stratigraphic-10Be concentration order within analytical uncertainty. However, we excluded the 115–120 cm depth sample from our analysis because the measured concentration (1.06 × 105 ± 2.2 × 103 atoms 10Be g−1) had a very similar concentration to the next lowest sample from 145–150 cm depth (9.86 × 104 ± 2.4 × 103 atoms 10Be g−1), 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−1. Thus, we limited maximum erosion to 20 cm ka−1 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, Reference Putkonen and Swanson2003; Applegate et al., Reference Applegate, Urban, Keller, Lowell, Laabs, Kelly and Alley2012).
Glacier and climate reconstructions
The coupled energy–mass balance and ice-flow models used in this study were originally developed by Plummer and Phillips (Reference Plummer and Phillips2003) and have been applied to estimate paleoclimate conditions for paleoglaciers in a variety of settings (Harrison et al., Reference Harrison, Rowan, Glasser, Knight, Plummer and Mills2014; Rowan et al., Reference Rowan, Brocklehurst, Schultz, Plummer, Anderson and Glasser2014; Leonard et al., Reference Leonard, Laabs, Plummer, Kroner, Brugger, Spiess, Refsnider, Xia and Caffee2017; Laabs et al., Reference Laabs, Plummer and Mickelson2006; Quirk et al., Reference Quirk, Moore, Laabs, Caffee and Plummer2018, 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 (Reference Fastook and Chapman1989) 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. (Reference Quirk, Moore, Laabs, Caffee and Plummer2018, Reference Quirk, Moore, Laabs, Plummer and Caffee2020), including the use of PRISM (Parameter-elevation Regression on Independent Slopes Model) (Daly et al., Reference Daly, Halbleib, Smith, Gibson, Doggett, Taylor, Curtis and Pasteris2008) 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. (Reference Quirk, Moore, Laabs, Plummer and Caffee2020) 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 10Be 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).
Kernel density estimates of cosmogenic 10Be 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 10Be 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−1 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−1 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−1 and 13.5 ka, respectively. Predicted inheritance values are relatively consistent across all simulations with median values ranging from 2.5–2.8 × 104 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.
(A) Modeled 10Be concentration depth profiles results (red curves) fit to the Gothic moraine soil pit sample in-situ-produced 10Be 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 ×105. (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−1, 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.
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 2 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).
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, Reference Blackwelder1915) 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, Reference Brugger2007; Briner, Reference Briner2009; Young et al., Reference Young, Briner, Leonard, Licciardi and Lee2011; Brugger et al., Reference Brugger, Laabs, Reimers and Bensen2019a, Reference Brugger, Ruleman, Caffee and Masonb; Schweinsberg et al., Reference Schweinsberg, Briner, Licciardi, Shroba and Leonard2020). 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. (Reference Gaskill, Colman, DeLong and Robinson1986) 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., Reference Benson, Madole, Phillips, Landis, Thomas and Kubik2004; Ward et al., Reference Ward, Anderson, Guido and Briner2009; Dühnforth and Anderson, Reference Dühnforth and Anderson2011), 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., Reference Benson, Madole, Phillips, Landis, Thomas and Kubik2004) and Sangre de Cristo Mountains (Leonard et al., Reference Leonard, Laabs, Plummer, Kroner, Brugger, Spiess, Refsnider, Xia and Caffee2017), 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., Reference Shakun, Clark, He, Lifton, Liu and Otto-Bliesner2015). 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, Reference Oerlemans2005; Shakun et al., Reference Shakun, Clark, He, Lifton, Liu and Otto-Bliesner2015) where glaciers such as the East River flow from areas of high slopes in the headwater to lower slopes near the terminus.
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−1, 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−1 while length changes slowed to just under 500 m kyr−1. The temporal pattern of ice retreat in the Elk Range is similar to that observed in the nearby Sawatch Range (Young et al., Reference Young, Briner, Leonard, Licciardi and Lee2011; Tulenko et al., Reference Tulenko, Caffee, Schweinsberg, Briner and Leonard2020) and in the San Juan (Guido et al., Reference Guido, Ward and Anderson2007) and Sangre de Cristo mountains (Leonard et al., Reference Leonard, Laabs, Plummer, Kroner, Brugger, Spiess, Refsnider, Xia and Caffee2017) to the south as well as in the Front Range (Ward et al., Reference Ward, Anderson, Guido and Briner2009; Dühnforth and Anderson, Reference Dühnforth and Anderson2011) to the northeast (Fig. 6). Generally, we observe the same deglacial pattern as Quirk et al. (Reference Quirk, Huss, Leonard, Licciardi, Plummer and Caffee2022) 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 (Reference Brugger2010) 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 (Reference Brugger2010) 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 (Reference Brugger2010) 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., Reference Oster, Ibarra, Winnick and Maher2015), 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., Reference Brugger, Ruleman, Caffee and Mason2019b). 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., Reference Brugger, Laabs, Reimers and Bensen2019a). In the Sangre de Cristo to the south, Leonard et al. (Reference Leonard, Laabs, Plummer, Kroner, Brugger, Spiess, Refsnider, Xia and Caffee2017) and Brugger et al. (Reference Brugger, Leonard, Refsnider and Dolan2021) 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, Reference Munroe and Laabs2013), 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., Reference Winnick, Carroll, Williams, Maxwell, Dong and Maher2017; Wan et al., Reference Wan, Tokunaga, Williams, Dong, Brown, Henderson, Newman and Hubbard2019, Reference Wan, Tokunaga, Brown, Newman, Dong, Bill and Beutler2021). 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., Reference Wan, Tokunaga, Williams, Dong, Brown, Henderson, Newman and Hubbard2019, Reference Wan, Tokunaga, Brown, Newman, Dong, Bill and Beutler2021). 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−1. 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., Reference Uhlemann, Dafflon, Wainwright, Williams, Minsley, Zamudio, Carr, Falco, Ulrich and Hubbard2022), 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
10Be 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 10Be 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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/qua.2024.5
Acknowledgments
This work was supported by U.S. Department of Energy grant DE-SC0020354 to Larsen. We thank J. Reithel at the Rocky Mountain Biological Laboratory for securing access to study sites. We also thank two reviewers and the associate editor for providing helpful comments. Prepared in part by LLNL under Contract DE-AC52-07NA27344; this is LLNL-JRNL-842938. Data from this study are archived at the Environmental System Science Data Infrastructure for a Virtual Ecosystem (ESS-DIVE) data repository: doi:10.15485/2316162.
Open access
