Reconstructing the Glacial History of Green Lakes Valley, North Boulder Creek, Colorado Front Range

Abstract

We document the post–Last Glacial Maximum (LGM) deglaciation history of a small catchment in Colorado’s Front Range and model the glacial history responsible. We combine cosmogenic exposure dating of moraine boulders and glacial polish from 10 sites down the valley axis with a set of 1-D and 2-D numerical glacier models in which climate is represented by prescribed winter and summer mass balance profiles. Moraine ages of 24–18 ka constrain the LGM maximum glacial extent at 15 km downvalley from the range crest. Glacial polish sites decline in age along the center portion of the valley, dropping to roughly 12 ka by 5 km from the crest. This supports a monotonic but non-steady decline in glacier length. Modeling of the deglaciation history reveals that the equilibrium line altitude (ELA) rise from ∼18 ka to ∼14 ka was between 250 and 350 m, and that from ∼14 ka to ∼12 ka was between 100 and 150 m. Complete deglaciation requires at least another 250–300 m of ELA rise. A transient model run employing this ELA history can both reproduce the spatio-temporal pattern of glacial polish ages, and honor the constraints on maximum glacial extent. The full deglaciation between 18 ka and 10 ka can be modeled successfully with a 4.5–6 °C temperature rise and little to no change in precipitation.

Introduction

The paleoclimate community has acknowledged that alpine glaciers are sensitive indicators of climate variations and has employed glacial geologic features to reconstruct the timing of glacier advances and retreats, which may then be used to reconstruct the climate history. The absolute dating of moraine boulders, using cosmogenic ¹⁰Be exposure dating, has evolved as a preferred method for determining the absolute timing of the formation of terminal and lateral moraines. While dated moraines provide strong control on the timing of the maximum glacier extent and the onset of glacial retreat, reconstruction of the full glacier retreat history has been more challenging, in particular if no recessional moraines with boulders are available for dating. Dating of glacially polished bedrock surfaces (Nishiizumi et al., 1989; Briner and Swanson, 1998; Colgan et al., 2002; Miller et al., 2006) along a longitudinal valley transect has recently emerged as an alternative strategy to develop retreat histories (Guido et al., 2007; Ward et al., 2009; Dühnforth et al., 2010). The underlying assumption when using this approach is that glacial erosion during the last glacial occupation of the valley entirely removes the cosmogenic production zone, leaving behind a bedrock surface with a cosmogenic inventory that is fully reset to zero, an assumption that is often reasonable in an alpine glacial environment. The inventory is not contaminated with nuclides from prior interglacial exposure. Hence, as the ice retreats and bedrock surfaces are again exposed to cosmic rays, the surfaces accumulate cosmogenic atoms, and the number of atoms measured in a sample directly corresponds to the time of deglaciation at the sample location. However, bedrock surfaces may contain inheritance if glacial erosion does not completely remove the cosmogenic production zone, which is about 2–3 m deep (e.g. Briner and Swanson, 1998). Potential causes for incomplete removal of the cosmogenic production zone include wide fracture spacing (Dühnforth et al., 2010), or extreme hardness of the underlying bedrock (e.g. Guido et al., 2007), or the presence of cold-based ice, which disallows glacial erosion (Sugden, 1978).

Alpine glaciers in the western United States reached their maximum glacier extent asynchronously but deglaciation appears to have occurred nearly synchronously (Young et al., 2011). While terminal moraine ages from the last glacial maximum (LGM) in the Wind River Range, Wyoming, and in parts of Colorado, span ∼24 to ∼16 ka, the ages from terminal moraines in Montana, Idaho, and northwestern Wyoming are generally younger (e.g. Gosse et al., 1995; Schildgen et al., 2002; Licciardi et al., 2001, 2004; Thackray et al., 2004; Benson et al., 2005; Munroe et al., 2006; Brugger, 2007; Licciardi and Pierce, 2008; Young et al., 2011). Various hypotheses have been put forth to explain the discrepancy in the terminal moraines ages. These include appeals to the influence of the Laurentide ice sheet on westerly atmospheric flow and the pattern of moisture delivery (e.g., Kutzbach and Wright, 1985; Hostetler and Clark, 1997; Licciardi et al., 2004; Thackray et al., 2004; Thackray, 2008), or to the Great Basin paleolakes as important moisture sources (Laabs et al., 2006, 2009; Munroe et al., 2006), or to non-climatic, local control mechanisms such as basin hypsometry (Licciardi and Pierce, 2008; Ward et al., 2009). Additional chronologies, in particular from areas with sparse or even no absolute age control, are required to discriminate among these possible explanations. In the Colorado Front Range, only one absolute deglaciation history has been documented, that from Middle Boulder Creek near Nederland (Ward et al., 2009). This record shows that deglaciation initiated around 20 ka and was completed by 12–13 ka, with slight interruption by a stillstand or slight re-advance between 16 and 14 ka. Additional data from other valleys in the Front Range would be useful to compare with Ward and others' chronology and inferred climate history, and to assess the degree to which the chronologies vary from one valley to the next. Such information is typically sparse, such that the glacial history from a range commonly leans heavily on that of a single valley. Thus, it is all the more important to document the fidelity of such records by assessing the spatial variability within a range.

In this study, we focus on the deglaciation history in Green Lakes Valley (GLV), which is part of the North Boulder Creek (NBC) drainage, directly adjacent to Middle Boulder Creek Valley (MBCV) (Figs. 1 and 2). We employ cosmogenic ¹⁰Be exposure dating of glacial polish to document the terminus history as the last glacial (Pinedale) waned, followed by 1D and 2D glacier models to explore the climate responsible for this history. The goals are to (a) place absolute age control on the pattern of glacial retreat, (b) reconstruct the climate history in GLV from the LGM to complete deglaciation using the equilibrium line altitude (ELA) as a climate reference, and (c) compare results from North Boulder Creek with the glacial history deduced from nearby Middle Boulder Creek using a similar strategy.

FIGURE 1 Shaded relief map of the Colorado Front Range showing North and Middle Boulder Creeks. Black box outlines the location of our study site in Green Lakes and Arapaho Valleys (see Fig. 2).

FIGURE 2 Sample sites and ages (Tables 1 and 2) shown on 1 m Digital Elevation Model (DEM) based upon LiDAR data for upper Green Lakes Valley (GLV). The mapped Last Glacial Maximum (LGM) glacier extent (light blue) in the North Boulder Creek drainage is based on Madole (1986) and LiDAR. DEM shows coalescence of ice from Arapaho and GLV tributaries, and the headwaters of Middle Boulder Creek valley to the south (also called 4th of July Valley), in which the deglaciation history was documented by Ward et al. (2009). The spatial pattern of ¹⁰Be ages is replotted in Figure 5 as a function of distance downvalley. Inset DEM shows the high resolution topography of the LiDAR data emphasizing the distinct boundaries of the lateral moraines.

LOCAL SETTING: BOULDER CREEK WATERSHED

The headwaters of the Boulder Creek drainage basin are presently occupied by the very small Arikaree and Arapaho Glaciers, which have surface areas of 0.06 km² and 0.24 km², respectively (Fig. 2). The presence of these diminutive glaciers indicates that the modern interglacial climate conditions are not sufficient to allow larger glaciers to occupy a significant portion of the upper Boulder Creek headwaters. The combination of topographic shading, wind-blown snow accumulation, and cold high-elevation temperatures allow presence of small fragments of glaciers in this and nearby catchment headwaters in the modern climate. Their elevations imply that the modern ELA lies not too far above the present topography—perhaps around 4100 m. The Arapaho and Arikaree Glaciers are most likely preserved due to their shaded location in north- and east-facing cirques near the range crest. Mass balance measurements on Arikaree glacier from 1968 to 2005 showed a balanced net budget for the period from 1968 to 1998, followed by a decade in which the net balance has been consistently negative. As a consequence, the Arikaree Glacier surface has lowered by up to 6.8 m during the 1968–2005 period (Caine, 2006). In addition, Haugen et al. (2010) show that rapid loss of area of the Arapaho Glacier slowed in the latter half of the century, presumably as it found its protected niche.

During the Pleistocene, large glaciers occupied the upper Boulder Creek watershed, originating at the modern continental divide and flowing ∼15–20 km down their respective valleys. Madole (1986) mapped terminal moraines and glacial deposits above the town of Nederland and correlated these moraines with glacial advances in the Middle and North Boulder Creek valleys during the Bull Lake and Pinedale glaciations.

Absolute ages from boulders on the Middle Boulder Creek (MBC) Pinedale terminal moraine using cosmogenic radionuclide exposure dating have provided the first absolute constraints on the timing of the local LGM (Schildgen, 2000; Benson et al., 2005; Ward et al., 2009). Incorporating previous results from Schildgen (2000) and Benson et al. (2005), Ward et al. (2009) calculated a range of terminal moraine ages from 18.2 ± 2.9 ka to 22.7 ± 2.0 ka with a weighted mean age of 20.0 ± 1.5 ka. Additional age control is available along a longitudinal valley transect in the MBC valley, where Ward et al. (2009) dated glacially polished bedrock surfaces to reconstruct the deglaciation history. They measured ages between 17.7 ± 2.7 ka close to the terminal moraine and 12.0 ± 1.1 to 13.5 ± 1.3 ka in the upper headwaters of MBC, suggesting that the deglaciation in MBC was completed by 12–13 ka (Ward et al., 2009). Independent ¹⁴C-age controls from Lake Dorothy, a tarn in one of the cirques in upper MBCV, supports the cosmogenic data set (Davis et al., 1992; calibrated using CALIB v. 5.0.1 (Stuiver and Reimer, 1993) with the calibration data set of Reimer et al., 2004). Ward et al. (2009) further modeled the deglaciation using a 2D glacier model. They showed that the climate responsible for the pattern of deglaciation most likely consisted of two main pulses of ELA rise, with an intervening stillstand or slight re-advance around 16 ka. Once the second pulse of deglaciation was initiated, the final demise of the MBC glacier appears to have been swift, fully deglaciating the valley by roughly 12 ka. This final pulse corresponds closely with the Bølling-Allerød interstadial climate event.

In this study, we explore the degree to which the deglaciation history documented in the MBC is matched in the adjacent North Boulder Creek valley headwaters. The valleys are similar in elevation, and are within a few kilometers of each other along the crest of the Front Range. The footprint of the MBC glacier was roughly twice as large as the NBC glacier, reflecting a larger number of contributing headwater basins.

STUDY SITE: GREEN LAKES VALLEY

GLV is the northern tributary of North Boulder Creek (Figs. 1–^{FIGURE 2} 3) whereas the Arapaho Valley is its southern tributary. High peaks reaching 4010 m (Arikaree Peak) and 4088 m (Navajo Peak) and ridgelines at 3960 m in the headwall of GLV are a part of the Continental Divide. GLV displays several features characteristic of an alpine glacial landscape: a U-shaped valley cross section, a valley floor dotted with roches moutonnées and punctuated by several rock basins now occupied by lakes, polished and striated bedrock surfaces, and erratics left behind after the last glaciation. All of these characteristic features suggest that GLV was significantly modified by glaciers during previous glacial episodes.

FIGURE 3 (A) Aerial photo of the Green Lakes Valley (GLV) and Front Range crest. GLV is nestled between Kiowa Peak and Niwot Ridge. The moraine complex extends significantly downvalley into a terminal moraine complex that reflects the Last Glacial Maximum (LGM) extent of the coalesced GLV and Arapaho Valley glaciers. Arapaho Glacier and Arikaree Glacier (not visible) are small glaciers in the headwaters of North Boulder Creek at present. (B) White arrow and point indicate location and direction of view of the upper Green Lakes Valley with its knobby bedrock terrain and interspersed lakes. Photos by R. S. Anderson (A), M. Dühnforth (B).

The coalesced NBC glaciers (from Arapaho and Green Lakes Valleys) flowed approximately 15 km downvalley from their respective headwaters, and generated moraines that occupy a large swath of relatively open landscape that is not well confined by steep valley walls (Fig. 3, a). The moraines are geometrically complex, consisting of several tree-covered ridges. Those from the LGM (Pinedale) glacial advances are distinctly more complex and sharply defined than those of older glacials, as revealed by a recently collected LiDAR data set that forms the backdrop to Figure 2. The right lateral moraine from this coalesced glacier dammed a small tributary to generate what has been dubbed Lake Devlin (Madole, 1986). Stratigraphic (Madole, 1986, 2010) and geophysical (Leopold and Dethier, 2007) studies of the lake deposits reveal that the lake was in place from roughly 31 ka (Madole, 2010) to 14 ka, and was drained by a rapid release when the lateral moraine dam was breached.

Methods

COSMOGENIC ¹⁰Be EXPOSURE DATING

In order to document the timing of glaciation and glacial retreat in GLV, we collected glacial polish samples along a longitudinal valley transect, and samples from the tops of boulders from the terminal moraine (Fig. 2). The glacial polish samples were taken from bedrock knobs in the valley center in order to prevent or at least to minimize the influence of snow and sediment shielding. On the moraine, we sampled boulders along the moraine crest that stood at least 1 m above the ground and appeared to have remained in a stable position since moraine formation. Based on the moraine height of 35 m above the valley floor and the expected moraine age of ∼20 ka, we followed sampling guidelines by Putkonen and Swanson (2003) and collected samples from 5 moraine boulders.

We sampled glacial polish and 2 glacial erratics at 11 locations in the valley and 5 boulders on a moraine 15 km downvalley from the drainage divide (Fig. 2). Bedrock exposures in the valley are common within the first 5–6 km from the drainage divide but are absent further downvalley where bedrock is blanketed by sediment cover. We therefore lack constraint on deglaciation history in the reach between the moraine and our lowest glacial polish sample ∼5.5 km downvalley from the drainage divide.

We chiseled 1-cm-thick bedrock flakes of glacial polish from the bedrock surfaces. The boulder samples were cut into 1-cm-thick samples. The samples were processed according to standard procedures at the University of Colorado cosmogenic isotope lab, and the samples were measured at the PRIME Lab AMS facility (Table 1). The exposure ages were determined from ¹⁰Be concentrations using the CRONUS-Earth online calculator (http://hess.ess.washington.edu/math/; accessed in November 2010) (Balco et al., 2008). Ages were calculated based on a sea level, high latitude production (SLHL) rate of 4.49 ± 0.39 atoms ¹⁰Be gram quartz⁻¹ year⁻¹ by neutron spallation and a very small muogenic component using the constant production rate scaling routine of Lal (1991) and Stone (2000) to obtain local production rates (Table 1). The ages we report (Table 2) are calculated based on the Lal (1991) and Stone (2000) constant production rate model. We used the constant production rate model in order to be consistent with other publications on western United States deglaciation chronologies. Using the time-dependent production rate model results in ages that are 1.4–3.0% younger. We note that regionally calibrated ¹⁰Be production rates are increasingly being used for selected sites in northeastern North America (NENA) (Balco et al., 2009) and the Southern Alps in New Zealand (Putnam et al., 2010). The calibrated NENA production rates are about 6–12% lower compared to the global production rates. So far there is no calibration site available in the Rocky Mountains that could inform us about the possible deviation from the global production rate. Use of the NENA calibration of Balco et al. (2009) and the scaling scheme of Lifton et al. (2005) that explicitly accounts for the effect of elevation may yield the best current approximation for regional production rates. Using this scheme, the ages of our sampling sites are up to +7% younger when compared with those using the global production rate and the time-independent scaling by Lal (1991) and Stone (2000).

TABLE 1 Sample names, locations, sample background information, and calculated production rates using CRONUS-Earth calculator (Version 2.2).

TABLE 2 Sample names, AMS-measured ¹⁰Be concentrations, and calculated exposure ages and uncertainties using CRONUS-Earth calculator (Version 2.2).

While the boulder ages we report are corrected for erosion, the glacial polish samples do not include an erosion correction, as the erosion is zero or at least so close to zero that the influence on the age is negligible. We do not apply a snow shielding correction for reasons outlined below.

CONSTRAINTS ON PRESENT-DAY SNOW SHIELDING IN GLV

In order to test whether the neglect of snow shielding is justified, we analyzed a 13-yr (1997–2009) record of snow depth data from GLV for the pattern of snow thicknesses in the upper part of GLV extending from the drainage divide to Green Lake 3 (Fig. 4). The data were collected during annual spring snow surveys by the Niwot Ridge Long-Term Ecological Research (LTER) project and the University of Colorado Mountain Research Station and can be downloaded from the LTER website (http://culter.colorado.edu).

FIGURE 4 Shaded relief map showing a 13-year record of snow depth measurements in Green Lakes Valley (GLV) (snow data from the Niwot Ridge LTER website, http://culter.colorado.edu). Snow depths vary non-uniformly and appear to highly depend on the degree of wind exposure. Greater snow depths are generally found in the lee of bedrock bumps or in sheltered areas in the cirque where we find the little fragments of the modern Arikaree and Arapaho glaciers. Circles A–D mark the locations of the Tidbit temperature data loggers (see also Fig. 7, temperature plots).

As the snow survey data only provide information about the depth and distribution of snow at the time of the survey, we used HOBO Tidbit v2 temperature data loggers (accuracy ±0.2 °C) to record temperature over the course of an entire calendar year. In this exercise we employ the temperature data as a proxy for snow coverage, which allows us to reconstruct the duration of snow cover at each individual logger site. We anticipated that when the Tidbit was beneath any significant snowpack the amplitude of the temperature variations would be significantly muted. We strategically placed Tidbits on two bedrock knobs that stand at least 10 m above the valley floor (Fig. 4; locations B and C). These bedrock knobs are potentially good sites for cosmogenic sampling, one of which we indeed sampled for glacial polish, as we expect that high wind speeds and high negative curvature of the knobs prevent significant snow accumulation during winter. Two additional Tidbits were placed on the north- and south-facing valley hillslopes, approximately 50 m above the valley floor, to record the importance of aspect and shading on diurnal and seasonal temperature variations. The Tidbits were covered by approximately 3- to 5-cm-thin rock slabs in order to prevent direct insolation. The recorded temperatures should therefore represent near subsurface rock temperatures.

1D AND 2D NUMERICAL GLACIER MODELS

In order to extract a climate history from the history of glacier terminus positions we have deduced from ¹⁰Be concentrations in glacial polish samples, we turned to models of the GLV glacier. We first simulated glaciers with a 1D flow line model and then employed a 2D map plane model. While the 1D models are very efficient, allowing us to explore the sensitivity of the model results to the various rules employed in the model, the 2D models honor more faithfully the geometric complexity of the valleys that the glaciers occupy. The models have been described extensively in earlier papers. The 1D models were built upon those described in MacGregor et al. (2000, 2009), while the 2D model were adapted from that reported in Kessler et al. (2006). We focus here only on the modifications we have made.

All models start with conservation of ice mass, which includes the net mass balance at the ice surface, b, and the divergence of ice discharge downvalley.

In the 1D flowline simulation, variations in width of the valley, W, were taken into account, as discussed in Anderson et al. (2004).

The models therefore require the meteorological forcing, here captured by a net mass balance profile, b(z). In past models (e.g. MacGregor et al., 2000) we prescribed an equilibrium line elevation (ELA), a maximum positive balance, b_max , and a net mass balance gradient, db/dz, to produce a pattern of net mass balance. In the present 1D simulations, we attempted to disentangle the winter and summer balances, b_w and b_s , respectively, so that they can be independently varied. Winter balance was taken to be controlled by a maximum snow accumulation, b_{w_max} , an elevation scale over which snow depth declines significantly, σ_bw , and an elevation around which this snow pattern pivots, z_bw :

The winter balance is always positive. The summer balance was captured by employing a simple positive degree day (PDD) melt algorithm (e.g. Kessler et al., 2006), in which melt rate is a function of the temperature above 0 °C. We prescribed a temperature field in which the mean annual temperature declines with elevation above a reference elevation, z_bs , at a rate governed by the local lapse rate, Γ, around which oscillation occurs on both annual (P_a  =  1 year) and daily (P_d  =  1 day) cycles.

The summer balance was therefore calculated from the sum of melt from all time increments, i, in which T_i > 0: where the positive degree factor, pdd, was taken to be 6 mm s.w.e./positive degree day (see for example compilations in Hock, 2003; and Brugger, 2010). The summer balance is always negative. An ELA was not prescribed, but instead fell out naturally by summing the winter and summer balances to obtain the net balance profile

We determined the ELA after the fact by noting the elevation at which the net balance profile b(z) crosses zero. As discussed below, this method allowed us to both honor constraints from the modern meteorology, and control separately the snow accumulation and temperature changes between LGM and present.

Ice discharge was accomplished by both internal deformation and by sliding, such that where U_s is the local sliding speed, and the vertically averaged horizontal speed associated with internal deformation is captured by where A is the flow law constant (A  =  6.77 × 10⁻¹⁷ Pa⁻³ yr⁻¹ or A  =  2.15 × 10⁻²⁴ Pa⁻³ s⁻¹), τ_b is the local basal shear stress, and H_i is the local ice thickness. This equation honors the nonlinear rheology of ice through Glen's flow law. Sliding is much less well understood, but is known to be strongly influenced by the hydrology of the glacier (e.g., Marshall, 2005; Bartholomaus et al., 2007; Cuffey and Paterson, 2010). Following Marshall (2005), we employed a sliding rule in which the local glacier sliding speed is governed by the product of the local driving stress, τ_b , with the flotation fraction, τ_b , taken to a power, P.

We modified this approach by adding two additional ratios, so that the dimensional constant U_s0 can be set by choosing the sliding speed for a given basal shear stress, τ_b  =  τ_b0 of 1 bar ( =  10⁵ Pa), a given ice thickness of H_i  = H_* (150 m in the models shown), and a given flotation fraction (we used 0.8). We employed P  =  2, and chose U_s0 to result in annual sliding speeds of 20–30 m yr⁻¹. We did not explicitly calculate the evolving water pressure field, but instead attempted to mimic the annual average flotation fraction along the glacier in the following way. While the melt rate is highest near the glacier terminus, the conduit system that serves to bleed off high water pressures evolves most rapidly at the terminus (Kessler and Anderson, 2004). Conversely, while the conduit system is less likely to reach the upper glacier, the melt rates there are lowest. We therefore suspected that the annual average flotation fraction is highest near the ELA and declines both up and down-glacier. We mimicked this expected pattern by prescribing a Gaussian pattern of flotation fraction, with maximum of 0.8 at the ELA and a length scale over which this decays to either side, L^* . We took L^* to be a given fraction of the downvalley location of the intersection of the ELA with the glacier surface, or L^*  =  σ_f *x_ela . In the models shown, we employed σ_f  =  0.8 to 1.2. It therefore scales with the length of the glacier.

While the surface mass balance and the characterization of ice internal deformation are the same in both 1D and 2D glacier models, the 2D model does not include the proxy for hydrology in the basal sliding parameterization. Instead, sliding was formulated as an attractor sliding rule in which basal sliding is controlled by the tendency of a glacier to achieve and maintain a basal shear stress of 1 bar (Kessler et al., 2006). The geometric representation of the glacier in two dimensions is one apparent advantage of using the 2D over the 1D model. Hence, glacial deposits and landforms can be used to calibrate the extent of the glacier, which in turn provides some leverage on the mass balance of the glacier and the position of the ELA. A field-based calibration of the glacier extent in the 2D glacier model is particularly important if the model is applied to areas that are strongly affected by wind-blown snow and the redistribution of snow in the landscape. Test simulations in GLV show that the 2D glacier model accumulates ice on broad flat areas in high wind-exposed terrain, which are known to have been ice-free during the last two glaciations. In order to reconstruct the ELA history during glacier retreat as precisely as possible, ice should only be allowed to accumulate where there is clear evidence for glacial occupation. In the absence of a windblown snow algorithm, we created a mask in the 2D simulation that effectively prevents ice from accumulating outside of a pre-defined maximum glacier outline. We used the glacier extent map by Madole (1986) in addition to topographic information from Light Detection and Ranging (LiDAR) data in order to delineate the boundary of the LGM glacier (Fig. 2). The LiDAR data were collected and processed by the National Center for Airborne Laser Mapping (NCALM) in the spring of 2010 for the Boulder Creek Critical Zone Observatory (CZO). As the LGM glacier deposited distinct glacial moraines (Fig. 2), the footprint of the LGM glacier can be easily mapped, in particular below the confluence of the GLV glacier and the Arapaho Glacier. A sharp boundary of the LGM deposits, identified by sharp-crested sinuous longitudinal ridges and kettle ponds, delineates the latest glacial extent in GLV (Fig. 2).

Results

The five terminal moraine boulders show ¹⁰Be ages between 18.3 ± 0.8 ka and 23.7 ± 1.2 ka, resulting in a mean moraine age of 20.9 ± 2.5 ka (Tables 1 and 2; Figs. 2 and 5). Our reported ¹⁰Be ages upstream of the terminal moraine are significantly younger than the moraine boulder ages; they range from 16.1 ± 0.3 ka to 10.3 ± 0.4 ka (Table 2). In the lower part of GLV, the cosmogenic ages decrease with distance upvalley from the terminal moraine, declining from 14.4 ± 0.4 ka, 14.3 ± 0.5 ka, 13.6 ± 0.9 ka, 12.8 ± 0.5 ka, and 12.4 ± 0.7 ka to 12.0 ± 0.4 ka near Green Lake 4. From there, the sample ages increase upvalley from 14.3 ± 0.5 ka to 14.8 ± 0.6 ka to 16.1 ± 0.3 ka (Fig. 2). Our youngest age (10.3 ± 0.4 ka) is obtained from a glacial erratic on a bench that sits about 200 m above the modern valley floor at Green Lake 3 (Fig. 6).

FIGURE 5 Terminus history of Green Lakes Valley (GLV) deglaciation constrained by ¹⁰Be ages and corresponding modeled equilibrium line altitude (ELA) history. Error bars on ¹⁰Be data depict internal (thick black) and external (thin black) uncertainty. Dotted circles represent ages from boulder samples; open circles are ages from glacial polish samples. Terminus history (gray) resulting from a transient 1D model in which the ELA history (black) implied from the set of 1D models is prescribed. Mean annual temperatures that would be measured at 3700 m are noted beside the resulting ELA (based upon modern values of −3.7 °C at this elevation). The model run focuses on period from 20 to 10 ka during which deglaciation occurs.

FIGURE 6 Large boulder sampled for ¹⁰Be. Note 1-m-diameter rounded boulder sitting atop larger boulder being sampled (sample # GLV5) at roughly 200 m above the floor of the adjacent valley. View up Green Lakes Valley, showing lateral bench on northern wall of the valley at this location. Photo by R. S. Anderson.

Snow depths measured during the annual GLV snow survey vary locally from 0 to up to 10 m thickness. The spatial distribution of snow depths is well correlated with local topography and reflects significant redistribution of snow by wind (Fig. 4). Winds here are dominantly from the west. Figure 4 demonstrates that deep snow preferentially occurs in wind-sheltered areas. For example, the greatest snow depths occur in the cirque of GLV, coinciding with the small Arikaree Glacier. Another zone of high snow accumulation occurs immediately downvalley of a prominent bedrock ridge protruding from the north wall of the valley, where divergence of wind should promote deposition. Typically, the tops of bedrock knobs remain free of snow, as the wind sweeps these areas clean. On the lee of these obstacles, wind speeds decrease such that snow can accumulate directly adjacent to the knobs. Snow depth measurements from adjacent sites, separated by less than 2 m horizontal distance, can therefore vary by a few meters. These local differences in the snow depth pattern highlight the influence of the highly non-uniform valley topography.

The Tidbit data support our contention that the tops of bedrock knobs remain mostly snow-free during the winter season. The Tidbits placed on knobs recorded large temperature fluctuations over the entire year, with shielded seasonal temperature maxima of nearly 40 °C and minima of less than −20 °C (Fig. 7). Diurnal temperature variations show a distinct trend of decreasing variability from south- to north-facing slopes (Fig. 7), with up to 40 °C temperature variability on the south-facing hillslope (location A in Fig. 4), up to 25 °C on the bedrock knobs in the valley center (locations B and C in Fig. 4), and up to 10 °C on the north-facing hillslope (location D in Fig. 4). During the snowfall season between October and the beginning of July, the recorded temperatures at sites B to D fluctuate above and below 0 °C, implying that these locations remained snow-free. Only site A on the south-facing slope shows steady temperatures of 0 °C for an extended period, suggesting that this site was covered by snow for a significant portion of the year. As we sampled bedrock knobs in the valley centerline, the reduction of cosmogenic production at these sites by snow shielding should be minimal.

FIGURE 7 Shielded rock surface temperature data from four locations in Green Lakes Valley (see Fig. 4 for locations) showing the temperature variability on (A) the south-facing hillslope, (B, C) the valley center, and (D) the north-facing hillslope. Red dashed line equals 0 °C temperature. The temperature data indicate that sites (B)–(D) remain snow-free during the winter months. Only site (A) experiences an extended period of snow coverage, as indicated by the time spent at low amplitude oscillations and in the “zero-curtain” associated with the phase change.

In our glacier simulations, we must impose a snowfall pattern, or winter balance profile, and a melt pattern that reflects a drop in temperature from present day conditions. Modern snowfall results in roughly 1 m of runoff from the GLV (Caine, 1995). The mean annual temperature measured at 3739 m at the D1 met station on the ridge north of GLV is −3.7 °C. The modern fragments of glaciers in the headwaters of GLV and the adjacent Arapaho Valley indicate that the present day ELA is slightly above the ridgetops (say 4000–4100 m). Only topographic depressions that promote low radiative forcing and high avalanche-derived accumulation support these small present-day glaciers. Recent work of Brugger (2010) and of others suggests that the extents of LGM glaciers in central Colorado required a decrease in mean annual temperature of roughly 6–8 °C from present-day values, and little change in precipitation.

In Figure 8 we show the modeled steady-state GLV glacier with a 15 km length that matches the location of the terminal moraine complex. The mass balance profiles in the 1D and 2D simulations that result in this match employ a maximum positive balance of 1 m (modern values), and a 6 °C drop in mean annual temperature from the −3.7 °C measured in the present GLV landscape. The mass balance profile achieves a gradient of 0.007 m yr⁻¹m⁻¹ in the lower ablation zone. The best fitting ELA resulting from these prescribed meteorological parameters in the 1D simulation ranges between 3280 and 3320 m, depending on the parameters chosen in the sliding velocity algorithm. Using the same climate forcing parameters in the 2D run, the best fitting ELA is 3300 m. The ELA drop between present and LGM must therefore have exceeded 700 m (4000–3300 m).

FIGURE 8 (A) 1D model of the Green Lakes Valley (GLV) glacier for an imposed steady-state climate run-out to 1000 years. Ice thickness profiles at 50-year intervals are shown in the bottom panel. Model terminus position coincides with the mapped terminal moraines. Mass balance profile for 1D model is shown on the right, where b is net mass balance, b_s summer balance, and b_w winter balance. The resulting equilibrium line altitude (ELA)  =  3320 m, where b  =  0. (B) 2D model of the GLV glacier reaches the terminus position during an ELA lowering experiment when ELA  =  3300 m (red solid line), which coincides well the 1D glacier results. The 2D simulation is run with a mask based on mapped glacial deposits (Madole, 1986), for example on the southern rim of the Arapaho Valley. Glacier thicknesses range from a minimum of 10 m to 320 m (see color bar). The yellow star indicates the position of the sample GLV-2 with significant ¹⁰Be inheritance. The relatively thin ice cover at this site, and the correspondingly slow sliding speeds there, suggest that the glacier at this location should not have eroded efficiently.

The calculated ice thickness at the location of the upper ¹⁰Be sampling site is roughly 100 m, and the calculated sliding speed there is only a few m/yr (Fig. 8). This should result in only minor glacial erosion at upper valley sites. In addition, the ice thickness at the location of the isolated sampled boulder is roughly 180 m. Its height above the present valley floor is roughly 200 m. Minor increases in ice thickness, caused by an increase in mass balance or a decrease in sliding speeds, could lead to a temporary thickening of the ice with subsequent deposition of the erratic boulder on the bench above Green Lake 4.

The response time of these small valley glaciers to fluctuations in climate is relatively short. We have verified this in model runs in which we prescribe an oscillating climate: the lag between ELA and terminus positions is on the order of 100–150 years for our model GLV glacier. Scaling of expected response times is commonly taken to be the maximum ice thickness divided by the ablation rate near the terminus (e.g., Cuffey and Paterson, 2010). In our simulations, the LGM glacier has maximum ice thickness of about 300 m and ablation rates near the terminus of the order of 3 m/yr. The division therefore yields response times of 300/3  =  100 years. We may therefore deduce the ELA responsible for each glacier terminus position by running the glacier model out to steady state for a series of chosen ELAs (Fig. 9). The relationship between glacier length and ELA is nonlinear, reflecting the geometries of the valley profile and the valley width. Knowing the age of each of these glacier terminus positions, we may now transform the history of glacier terminus into a history of ELA (Figs. 5 and 10). Following the results from the 1D model, the estimated ELA history displays a 250–300 m rise of ELA from 3300 to 3550–3600 m between ∼18 ka and ∼14 ka, with a subsequent rise of another 100–150 m to 3700 m by ∼12 ka. The history thereafter is unconstrained by ¹⁰Be; but a rise of another 300 m is required to achieve the present minimal glacier cover. The ELA history derived from the 2D model corresponds well with the 1D results: the ELA rose about 350 m between ∼18 ka and ∼14 ka, followed by a 100 m rise from ∼14 ka to ∼12 ka, and another 300 m required to restrict ice to the modern glacial footprint.

FIGURE 9 Modeled glacial terminus position as a function of imposed equilibrium line altitude (ELA). Results from 1D and 2D glacier models for a set of fixed ELA positions. Each run is modeled out to steady state, many hundred years. The nonlinear relationship reflects the valley hypsometry and the various steps in the valley profile. Two sets of model runs were performed in 1D. The difference in the results reflects choice of the parameters governing the sliding rates (labeled). The mean annual temperatures (MAT) used in the 1D mass balance routine that lead to the reported ELAs are also shown.

FIGURE 10 Equilibrium line altitude (ELA) history in Green Lakes Valley (GLV) resulting from combination of dated terminus positions and modeled relationship between ELA and terminus position (Fig. 9). Initiation of glacier retreat (left shaded bar) coincides with amelioration of climate deduced from δ¹⁸O record in GRIP ice core (Johnson et al., 1997). Reduced winter temperatures resulting from extensive sea ice cover in the North Atlantic appear to have caused an extended cold period (unshaded) between 17.5 and 14.7 ka (Denton et al., 2005; Schaefer et al., 2006). The near-synchronous retreat of glaciers in the western United States at ∼16 ka (Young et al., 2011) overlaps with the warming trend that starts at ∼16 ka and culminates in the Bølling-Allerød interstadial climate event (right shaded bar).

We show in Figure 5 a 1D model run in which the ELA history deduced from the set of steady-state glacial runs is imposed over the period from 20 to 10 ka. The history of glacial terminus positions constrained using ¹⁰Be is indeed well mimicked by the model. This supports the contention that the time scale of glacial response to climate change is short enough to warrant use of steady-state models. The transience in the calculation is shown only in the rounding of the edges of the terminus history, showing that the time scale for response is of the order of 100 years.

Discussion

In order to reconstruct the timing and pattern of deglaciation in glacial valleys, we expect a trend of ages that decreases monotonically with distance upvalley from the terminal moraine. As seen in Figures 2 and 5, this is indeed broadly the case. There are, however, significant departures from this simple trend. We discuss these first, before addressing the broader pattern and its implications.

DISCUSSION OF THE COSMOGENIC ¹⁰Be AGES

Two types of outliers clearly do not follow the trend of decreasing age with distance upvalley: the three glacial polish ages downvalley from the uppermost sample in the transect near the cirque, and the erratic boulder age from the bench above Green Lake 4 (Fig. 2; Tables 1 and 2). The glacial polish age of 16.1 ± 0.3 ka in the uppermost part of the valley is the most obvious exception, although the ∼14 ka ages of the next two samples downvalley also disobey the expected upvalley decline in age. We argue that the sampled bedrock surfaces were not completely reset during the last glaciation (Briner and Swanson, 1998; Fabel et al., 2004; Guido et al., 2007; Ward et al., 2009; Dühnforth et al., 2010), and that the incomplete removal of the cosmogenic production zone left a component of inherited cosmogenic ¹⁰Be atoms behind. At least three scenarios could explain the inheritance at the sample sites: (1) the rock here is less erodible, (2) the glacier was cold-based and less erosive, or (3) the glacier was less erosive due to low sliding speeds.

When comparing the nature of the bedrock at the uppermost sample site with that at other localities along GLV, it is clear that a compact, massive area of bedrock, standing at least 10 m above the valley floor, is exposed across the entire valley width in this reach (Fig. 11). The presence of this massive bedrock bench may therefore reflect slower rates of glacial erosion due to higher resistance to erosion. It is also possible that the glacier in this upper valley was insufficiently erosive to remove the full production profile at these upper sites. The presence of cold-based ice, leading to reduced or absent glacial erosion (Sugden, 1978; Briner et al., 2006) could be another possible explanation for the inherited ¹⁰Be. In the absence of a thermally coupled glacial flowline model, we are unable to evaluate the likelihood of cold-based glacier conditions in the upper part of GLV. However, a lowering of basal glacier sliding speeds and hence glacial erosion rates can also occur when the glacier is relatively thin. It is indeed striking that the uppermost sample, with ¹⁰Be highest inheritance, occurs exactly where the 1D model suggests ice thickness is significantly thinner (<100 m) compared to all other sampling sites (Fig. 8).

FIGURE 11 View from edge of Green Lake 5, showing massive bedrock band from which sample GLV2 was chiseled. Arikaree Peak in background. Photo by M. Dühnforth.

Here we first estimate how much more erosion would have been needed to reset these upper three sites. Consider the uppermost site. Given the discrepancy between the measured 16.1 ka age and the age one would project from the sample sites immediately downvalley (∼12 ka), we can estimate how much erosion would have been required to lower this 4 ka worth of inheritance to levels that are within an error of zero. If the site were exposed to ¹⁰Be production for roughly 50 ka prior to reoccupation by glacier ice, it is this inventory that must be removed. Of 50 ka only 4 ka now remains. The production profile is well described by where P₀ is the local surface production rate, and z* is the depth scale for decay of production ( =  Λ/ρ, where Λ is the attenuation length (g cm⁻²), and ρ is the density of rock; z*  =  160 g cm⁻²/2.7 g/m³  =  60 cm  =  0.6 m). The inheritance is simply the product of the duration of exposure during the interglacial, t, with this production profile. Given the length scale for decay of production with depth is about 0.6 m, we can solve Equation 9 for the depth of erosion required to retain only 4/50 of the production profile.

This calculation suggests that 1.5 m of erosion, z_e , must have occurred. To reduce the 4 ka of apparent inheritance to the 0.3 ka that we find as uncertainty in the ¹⁰Be measurements of this and nearby sites, another 1.5 m of erosion would have been required. The site was roughly halfway toward being reset. This calculation obviously hinges upon the assumed duration of interglacial exposure, but illustrates that (i) the depth of erosion required to remove enough of the ¹⁰Be inventory for the remaining concentration to lie within the errors of measurement in this setting is roughly 3 m (1.5 ± 1.5 m), and that (ii) this uppermost site retains a significant portion of this inheritance. Only half of the required erosion occurred here. The same calculation suggests that 1.9 m of erosion had occurred at the other upper sites with ages of ∼14 ka; another 1.1 m would be required to reset them fully (Fig. 12).

FIGURE 12 ¹⁰Be production profile showing truncation needed to retain 4 ka, 2 ka, and 0.3 ka out of 50 ka. Reduction of the production rate to only a few percent of the surface production rate requires erosion of about 3 m, given the length scale for the decay of production with depth of about 0.65 m. Smaller amounts of erosion result in more retention of the signal from nuclides inherited in prior interglacials (e.g., 8% (4/50) or 16% (4/50) of the inherited signal remains for erosion of 1.9 or 1.5 m, respectively).

The second type of outlier, the erratic boulder with an age of 10.3 ± 0.4 ka, is clearly younger than the timing of deglaciation. We are certain that this boulder is a glacial erratic deposited during the Pinedale deglaciation, as another rounded boulder (∼1 m in diameter) of different lithology is deposited on top of our sample boulder (Fig. 6). We therefore expect an age between 12 and 14 ka, based upon its location in the GLV profile and the other glacial polish ages nearby. Possible explanations for the several ka discrepancy between these ages and our 10.3 ka ¹⁰Be measured age include erosion of the boulder surface and shielding of the sample site by snow or sediment. It seems rather unlikely that the postglacial erosion was great enough to accomplish a minimum of 20–30 cm rock removal that would be required in order to obtain a deglaciation age of 12–13 ka. An alternative explanation of the unusually young age could be shielding by stagnant ice after the boulder was deposited along the LGM glacier margin. As the boulder sits on a bench in a position comparable to the Martinelli snowfield under present-day conditions (Williams et al., 2000; Caine, 1992), snow and ice may have lingered longer there than elsewhere in the valley, effectively shielding the boulder from the full cosmic ray flux.

SNOW SHIELDING

Snow shielding can potentially influence the measured cosmogenic exposure ages; it results in ages that are younger than the true age. Here we discuss the possible influence of snow shielding on the GLV sample results. In contrast to the discussion above, where we treated the 16.1 ± 0.3 ka and the two ∼14 ka ages as exceptions, we treat now the three ∼12 ka ages in the middle section of our longitudinal transect as potential exceptions, implying that these measured samples experienced snow shielding in the past that lowered their ages. In the following exercise we calculate how much snow shielding is needed in order to increase the ∼12 ka ages to ∼14 ka ages such that they are in line with the ages in the lower and upper part of the valley, and we discuss how realistic this scenario is. For this exercise we use the measured range of snow thicknesses from the Niwot Ridge LTER snow survey data set and assume that measured snow thickness shields the sample surface for two end-member durations of 6 and 3 months per year. We estimated the snow water equivalent (SWE) by multiplying snow depth by snow density. The SWE is needed to calculate the snow-shielding factor (S _snow) which in turn is used to adjust the cosmogenic ¹⁰Be production rate (Gosse and Phillips, 2001; Schildgen et al., 2005).

Figure 13 illustrates that an annual mean snow cover of 0.75 to 4.5 m, depending on the snow density and duration of snow cover, are required to raise the ∼12 ka exposure ages to 14 ka. These snow thicknesses are within the range of snow depths that were measured by the LTER snow survey and therefore the young ages could potentially be explained by snow shielding. However, the fact that the bedrock knobs in the upper part of the valley appear to be unaffected by snow shielding under present-day conditions suggests that snow shielding most likely does not occur at the sampling sites in the lower part of GLV, where we used the same sampling strategy targeting high bedrock knobs to avoid shielding effects. Neither snow survey nor temperature data are available to show the true snow thickness in this lower part of GLV. One possible explanation for this pattern of cosmogenic exposure ages is that stagnant ice or sediment shielding could have led to the younger ages.

FIGURE 13 Mean annual snow cover plotted against the snow shielding factor used for correction of the local cosmogenic ¹⁰Be production rate. Two scenarios, one with an annual snow cover of 3 months (dashed lines) and the other with 6 months (solid lines) are plotted using different snow densities. The black dotted lines show the range of different snow depths that would be needed in order to increase the ages of the four outliers GLV-1, GLV-5, GLV-6, and GLV-7 to 14 ka such that they are in line with the samples up- and downvalley.

Between these two possible explanations of the outliers presented here, we adopted the first scenario in which the upper valley samples contain inherited nuclides due to incomplete glacial erosion of the bedrock surfaces. The incomplete removal of rock led to ages that are older than expected based on the succession of ages upvalley from the terminal moraine. The exact reason for the incomplete removal of the production zone is difficult to determine, and it may even be that more than one of the mechanisms we outlined above is responsible for the inheritance in the upper part of the valley. As the temperature data show that bedrock knobs in the upper GLV remain snow free during the winter, we believe that snow shielding exerts no significant influence on the measured ages in the valley. We argue that it is more plausible that the erosional efficiency of the LGM glacier was insufficient to completely reset the three cosmogenic sampling sites in the upper valley. This is reflected in ages that are too old to represent true deglaciation ages.

COMPARISON WITH GLACIAL HISTORIES DEDUCED FROM OTHER LOCAL CONSTRAINTS

Our results indicate that the deglaciation of GLV occurred between 18 ka and 12 ka and that GLV was largely ice-free by 12–13 ka. The timing of glacier retreat may additionally be supported by radiocarbon and Optical Stimulated Luminescence (OSL) data from sediments from Lake Devlin, an ice-marginal lake impounded by the North Boulder Creek Valley Pinedale glacier 3 km upvalley of the Pinedale (LGM) moraine. Lake Devlin existed between about 31 and 14 ka (Madole, 1986, 2010; Leopold and Dethier, 2007). The younger age constraint can be interpreted as when the lake breached catastrophically after the glacier had receded from the position of the Lake Devlin outflow (Fig. 2). The ponding of Lake Devlin was accomplished by the right lateral moraine of the coalesced glacier from both GLV and Arapaho Valley tributaries. As such, its history is likely more governed by the larger of these valleys (Arapaho) than by GLV, and therefore represents only a soft constraint on the GLV glacier. The ¹⁰Be-based deglaciation ages suggest that the GLV glacier terminus was already several kilometers upvalley from the Lake Devlin breach site when the moraine-dammed lake drained. This implies that either the water balance in the lake prevented it from overtopping the moraine until well after the glacier terminus receded, that a portion of the Arapaho Glacier remained to dam the lake, or that it simply took that long for the moraine dam to fail. Whereas the initial ponding of the lake requires that the glacier terminus be downvalley of the lake site, the recession of the glacier back past that location is not a sufficient condition for failure of the lake.

The history of deglaciation in GLV is in good agreement with the deglaciation history in MBCV (Ward et al., 2009; Davis et al., 1992) in the next valley network to the south. Deglaciation appears to have initiated around 18 ka in both valleys. In MBCV the best-fitting numerical models suggest that the deglaciation stalled at around 16 ka, with some models favoring slight readvance. Final deglaciation of MBCV appears to correspond well with the Bølling-Allerød interstadial climate event, captured in the GRIP δ¹⁸O record (Fig. 10), and lasting formally from 14.7 to 12.7 ka. We note, however, that the use of a lower production rate such as the NENA scaled with the Lifton et al. (2005) scheme, results in older ages that overlap with both the Bølling-Allerød and the Younger Dryas event. In neither the MBCV nor GLV is there a field evidence for a readvance in the subsequent Younger Dryas (∼12.8–11.5 ka) cold event.

DISCUSSION OF MODEL RESULTS

The calculated ELA required to reach the LGM terminal moraines in both the North Boulder Creek and Middle Boulder Creek valleys is the same: 3300 m. It therefore appears that the same climate history can indeed be invoked in these adjacent valleys. This is encouraging. The rise of ELA needed to deglaciate the GLV valley is at least 700 m (assuming present ELA  =  4000 m), but may be as high as 800–900m, depending upon the present ELA. As 4100 m is the top of Navajo Peak at the top of the present drainage divide of GLV, 4200 m would be 100 m above it. If we employ a lapse rate of 6.5 °C km⁻¹, this implies a minimum temperature drop of 4.5 °C, which is less than most temperature depressions suggested for other locations in the Rocky Mountains. If the present ELA is indeed at 4200 m, the temperature drop would be roughly 5.8 °C. In the Sawatch Range and Elk Mountains in southwestern Colorado, the reconstructed mean annual temperature depression during the LGM was around 7 °C (Brugger, 2010). Even changes in the precipitation rates of ±25% appear to require only minor (0.5 °C) reduction in the estimated temperature drops. Our model suggests that in order to match the LGM moraines, a 50% increase in snow accumulation from modern rates would require a smaller drop in temperature, down from ∼6 °C to 4.5 °C. On the other hand, LGM snow accumulation of 50% of the modern rate would require a greater temperature decline, up from ∼6 °C to 8 °C. Leonard (2007) suggested that mean annual temperatures were 4 to 9 °C lower during the LGM compared to today. In his calculations, LGM glacial lengths could be explained by anything between a 100% increase in precipitation combined with a 4–5 °C temperature depression, and modern precipitation combined with a temperature depression of 6–8 °C. Significantly larger temperature depressions have been suggested by Leonard (1989), Elias (1996), Hostetler and Clark (1997), and Kaufman (2003), who proposed that LGM summer temperatures in the central and southern Rockies were 7 to 13 °C lower than today. Our estimate of a 4.5–6 °C temperature depression at the LGM in the Front Range may be compared with estimates from other ranges in Colorado. Brugger (2010) showed that LGM ELAs in central Colorado display a regional trend of increasing elevation toward the east, at ∼4.5 m km⁻¹. Projection of his LGM ELA trend to the longitude of the GLV (105.65°W, 31 km east of the easternmost point in Figure 7 in Brugger (2010) would yield an ELA estimate of 3600 ± 4.5m/km*31 km)  =  3740 m. Similarly, projection of the calculated modern snowline trend from the Sawatch Range to the longitude of the Front Range would yield a snowline estimate of 4850 m. Neither of these projections is supported by our findings. The LGM ELA in GLV, as in the adjacent Middle Boulder Creek Valley, is 3300 m, which is similar to many of the ELA estimates in the Sawatch Range (Brugger, 2010); no significant gradient in ELA exists at this range-to-range scale. In addition, the fact that the Front Range preserves numerous small glaciers even under modern climate conditions, in contrast to the Sawatch Range, with its many tall peaks exceeding 4250 m but absence of glaciers, suggests that a regional snowline gradient of the opposite sign exists today between the Sawatch Range and the Colorado Front Range. We therefore argue that the range-wide climate governing glaciation in the Front Range cannot be simply extrapolated from trends in ELAs estimated from other ranges in Colorado. While trends within a single range or nearby ranges may reflect local patterns of precipitation and melt, as shown by Brugger (2010), the meteorological settings of ranges vary significantly within the complex topography of the Rockies. The temperature drop of roughly 4.5–6 °C that we infer from the estimated drop on ELA and the simple lapse rate scaling is supported by the 1D and 2D modeling in the absence of significant changes in precipitation.

Conclusions

Absolute ages on the deglaciation history in GLV, the northern tributary of the NBCV, show that the onset of post-LGM glacier retreat occurred at 18 ka, and was followed by a rapid deglaciation that was completed by ∼12 ka. The timing of deglaciation correlates with the reconstructed retreat history in the adjacent Middle Boulder Creek Valley to the south. While our results appear to suggest a monotonic retreat, the lack of data for nearly 10 km between the terminal moraine and our lowest samples in the upper valley prevents us from detecting the history of NBCV between between 16 and 14 ka, although a readvance or stillstand has been suggested for the glacial history in the adjacent MBCV. The ELA history deduced from modeling the glacial retreat history displays a rise of at least 700 m over the 6 ka in which the LGM glacier met its demise. Our modeling suggests a mean annual temperature drop of 4.5–5.8 °C during the LGM. The absolute ages from GLV confirm the overall pattern of LGM ages across the western United States with older terminal moraine ages in Colorado and parts of Wyoming and generally younger local LGM ages in Montana, Idaho, and northwestern Wyoming. In addition, the initiation of deglaciation in GLV corresponds with the nearly synchronous onset of deglaciation across the western United States at 18–16 ka (Fig. 3 in Young et al., 2011). As in the adjacent MBC drainage, it appears that the final pulse of deglaciation correlates with the Bølling-Allerød event.

Acknowledgments

This research was funded by a grant from the National Science Foundation (Boulder Creek CZO; NSF-0724960). Logistical support and snow survey data were provided by the NSF-supported Niwot Ridge Long-Term Ecological Research (LTER) project and the University of Colorado Mountain Research Station. We thank Craig Lindel for access to the Boulder Creek Watershed in which the sampling sites sit. We benefited from discussions with Nel Caine, and with Dylan Ward, whose work in the adjacent Middle Boulder Creek watershed was inspirational. Critical reviews from two anonymous reviewers helped hone the manuscript significantly.