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 23) 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.

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:

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

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).

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

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)

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.