Climate During the Last Glacial Maximum in the Southern Sangre de Cristo Mountains, Colorado, U.S.A.

1. Introduction

Knowledge of climate during the Last Glacial Maximum (LGM), an extreme climate state in Earth’s past, is fundamental for understanding climate dynamics that ultimately informs model projections of future climate change. This knowledge has largely relied on pollen spectra and other microfossils, speleothem records, lake level reconstructions, and so forth. These proxies in-and-of-themselves are important in understanding paleoclimate dynamics, but also because they constrain regional scale climate modeling, particularly in the western United States [1,2,3,4]. Yet while forming the core of available climate proxies, these often do not yield quantitative measures of climatic parameters. More significantly, some proxies begin some time after, and thus do not strictly record LGM climate.

Alternatively, modeling of mass balances of paleoglaciers has provided quantitative estimates of LGM climate in the Rocky Mountains. In general, temperatures and/or precipitation are sought that maintain steady-state mass balances of paleoglaciers at their LGM extents using difference approaches. These include more physically-based coupled energy/mass balance and ice flow models [5,6,7,8,9,10,11] and more empirical temperature-index modeling [12,13,14,15,16]. While these approaches have in some instances differed on the precise magnitude, they suggest LGM temperature depression of ~5 to 9 °C cooler than present. However, most studies have focused on glaciated ranges in Utah, Wyoming, and northern and central Colorado; few studies [8,10,16] have examined the most southerly regions of the Rocky Mountains region. Moreover, temperature depression in the southernmost Sangre de Cristo Mountains in New Mexico appears to be cooler than expected [10] from that suggested by climate modeling (e.g. Oster et al., 2015 [1]; Lora and Ibarra [4]; Shafer et al., 2021[17]) that show a consistent north-to-south decrease in cooling along the Rocky Mountains.

Based on the foregoing, it is clear that our understanding of paleoclimate over the southern extent of the Rocky Mountains would benefit from additional estimates of LGM temperature depression. These would add to the existing sparse record of late glacial climate change and thus better constrain regional climate simulations. Accordingly, in this paper we provide another estimate of LGM temperature depression in the southern Sangre de Cristo Mountains in Colorado based on glacier reconstruction and temperature-index modeling.

2. Regional Setting

2.1. Geologic and Geomorphic Setting

Forming a prominent physiographic barrier, the fault-bounded Sangre de Cristo Mountains extend from south central Colorado south-southeast into northern New Mexico (Figure 1). Rocks of the region include Proterozoic igneous and metamorphic basement and Paleozoic through Cretaceous sedimentary sequences. Overall the structural history of the range is complex, involving deformation through the Late Paleozoic to the present. The range is bounded on the west by the Rio Grande Rift that became active as early as the Miocene [18]. To the east lie the Wet Mountain Valley, a down-faulted block in the north, and southward the downwarped trough of the Raton Basin.

Rising some 1000 to over 1500 meters above the adjacent valleys, the crest of the Sangre de Cristo consists of numerous peaks exceeding 4000 meters. Ongoing tectonism has resulted in deeply incised and steep valleys. Many valleys exhibit features of alpine glacial erosion and deposition of the last, or Pinedale glaciation (corresponding to the LGM and Marine Isotope Stage (MIS) 2) and to some extent the penultimate Bull Lake glaciation (MIS 6). Scant evidence exists for a pre-Bull Lake glaciation [19]. Rock glaciers, possibly still active, and other periglacial landforms are common in cirque basins. Valley mouths are deeply incised owing to continued tectonism of the Rio Grande Rift and are fronted by well-developed alluvial fans with surfaces indicating multiple stages of development [20,21].

2.2. Modern Climatology

Existing meteorological stations within the study area are limited to elevations (<~2700 m) that are below those of glacier extents. Modern climate is therefore characterized in the model by the PRISM gridded climatology (Parameter-elevation Regressions on Independent Slopes Model; http:/www.prism.oregonstate.edu; [22]) augmented with data from three SNOTEL (Snow Telemetry) stations all located at similar higher elevations (Figure 1 and Figure 2). Median values of 1990-2020 normals are used here. Mean annual precipitation at the highest elevations (~4000 m) is ~100 cm decreasing to about half that at valley mouths (~2900 m).

SNOTEL data are representative of the monthly variation in precipitation (Figure 2). Stations on the eastern slopes (North Castilla and Whiskey Creek) show fairly distinct peaks in precipitation during mid-spring and summer, the latter being more pronounced. This peak is not present in the Culebra SNOTEL record on the western slope. This suggests that the mid-spring peak is possibly related to either late winter/early spring upslope precipitation events over the eastern slopes [23], or a result of snow advection by the prevailing westerlies augmenting snow accumulation on eastern slopes [24]. The summer peak in precipitation, evident in all three records, presumably reflects the influence of the summer monsoon [25] that brings subtropical moisture from the Pacific Ocean and Gulfs of California and Mexico. Regionally within the study area, the PRISM data (1991-2020 normals) indicate that eastern slopes annually receive ~15-20 cm more precipitation over elevations between 3000 and 4000 m, this difference increasing with elevation. Precipitation dependence on elevation, a critical aspect of the modeling pursued here, is statistically significant; all correlation coefficients (r²) exceed 0.89 with the number of PRISM gridpoints n used to develop the regressions varying between 45 and 187 depending on the valley or location of interest.

Maximum summer and minimum winter temperatures are typically ~15 and –3 °C, respectively, in valleys, and about ~8 and –11 °C at the highest elevations in the study area. Temperatures strongly vary linearly with elevation, with r² values all > 0.95 (again, n = 45 to 187). Seasonal temperature variations (January to July) are on the order 18 to 20 °C and relatively constant, but in some valleys show some dependence on elevation. No appreciable differences in seasonal temperature variations are evident over eastern and western slopes as exemplified by the SNOTEL records (Figure 2).

2.3. Timing of the Local LGM

A glacial chronology for the study area has not been established but ¹⁰Be exposure ages of boulders on terminal moraines in other areas of the Sangre de Cristo Mountains provide a temporal setting for the local LGM. Preliminary exposure ages (Feldman et al., in prep.) from moraines in Long Valley, New Mexico, about 25 km south-southwest of the study area, suggest that glaciers in that valley reached their maximum extent at ~24.5 ± 1.6 ka. In the Winsor Creek valley, also in New Mexico, about 150 km south-southwest, the LGM is dated to 21.2 ± 2.0 ka [10]. In the South Crestone and South Colony valleys of Colorado, approximately 100 km north-northwest of the study area, glaciers advanced to their maximum ~21 ka [8]. Considering these ages in the context of those in other ranges in Colorado [14,15,26,27,28], the timing of the local LGM in the study area was likely ~20-25 ka.

3. Materials and Methods

3.1. Glacier Reconstruction

The LGM ice extents of seven paleoglaciers (Figure 3) were reconstructed using varying combinations of field mapping of glacial features, high-resolution lidar imagery and derived digital elevation models, geological [29,30,31] and topographic maps, and Google Earth® imagery. Private land ownership limited field mapping in some areas, therefore the reconstructions presented here are more reliant on “remote approaches,” especially lidar. Additionally, detailed reconstructions were not undertaken in valleys where ice limits were ambiguous owing to mass movements, possible remobilization of till in terminal and lateral moraines, among other complexities that would increase uncertainties in model results. Ice in these valleys are shown schematically in Figure 3.

Although software exists for paleoglacier reconstructions, ice surface contours were drawn manually with due consideration of ice extent, local flow directions indicated by erosional forms – often subtle on imagery, and general convergent flow in accumulation zones and divergent flow in ablation zones. Contours were adjusted iteratively so that reconstructed ice surface slopes were sub-parallel to those of the valley and to ensure driving stresses τ were between 50 and 200 kPa, values typically associated with modern glaciers [32]. Stresses were calculated using:

$$\tau = S_{ƒ}\rho gh sin\alpha$$

(1)

where ρ is the density of ice, g is gravitational acceleration, h is ice thickness, α is the slope of the ice surface, and S_ƒ is a shape factor to account for drag of the valley sides [33]. Surface slopes were averaged over distances of 10h to account for longitudinal stress gradients [32,34].

3.2. The Temperature-Index Model and Its Implementation

The temperature-index model (TM) and the rationale for using this approach (as opposed to energy-balance modeling) has been discussed in greater detail elsewhere [14,15,16], therefore a brief overview is presented here. The net mass balance B_n of a glacier in steady-state with climate will be equal to zero; that is averaged over some years and over its surface area the glacier is neither gaining nor losing mass. The annual variation of the specific net mass balance, here at a specific elevation z, is

$$b_n\left(z\right)={\int }_{t_1}^{t_2}\left(P_s\left(t,z\right)+M\left(t,z\right)\right)dt$$

(2)

where P_s(t,z) is the rate of snow accumulation, M(t,z) the rate of snow or ice melt (ablation) over the glacier’s surface during the interval t₁ to t₂, taken as the hydrologic year. In practice, Equation (2) is numerically integrated over a monthly time-scale to yield monthly melt that is combined with monthly precipitation. This sum is then integrated over the hydrologic year to obtain the net balance at z.

Monthly melt is calculated by summing daily melt using a melt factor m_ƒ that empirically relates ablation to mean daily air temperature T_d(t,z):

$$M\left(z,t\right)=\left\{\begin{array}{c}{m}_{ƒ}{T}_d\left(t,z\right) T_d\left(t,z\right)>T_m \\ 0 T_d\left(t,z\right)\le T_m\end{array}\right.$$

(3)

where T_m is a threshold temperature above which melting occurs. A compilation of measured melt factors for ice (supplementary material in [16]) yields mean (and standard deviation) and median values of 0.73 ± 0.27 and 0.67 cm w.e. (water equivalent) °C^–1 d^–1 respectively (n = 92). Similarly, for snow the values are 0.43 ± 0.2 and 0.4 cm w.e. °C^–1 d^–1 (n = 61). For the sake of modeling therefore, the values of m_ƒ for ice and snow are taken as 0.70 and 0.45 cm w.e. °C^–1 d^–1 but in subsequent sensitivity analyses these values were allowed to vary by 0.2 cm w.e. °C^–1 d^–1. In the model, m_ƒ is initially set for snow until all snow is melted after which that for ice is used. The threshold for melt was set to +1 °C; additional modeling runs used 0° C.

Daily air temperature as a function of time of year and elevation is

$$T_d\left(z,t\right)=\left[H\left(z\right){\left[{\displaystyle \frac{1-\cos \left(\left(\frac{2\pi d}{365}\right)-\varphi \right)}{2}}\right]}^{k\left(z\right)}-T_{jan}\left(z\right)\right]-∆T$$

(4)

where H(z) is the difference between mean July and January temperatures (amplitude of the sinusoid), d is the day of the year, ϕ is the phase lag (= 0.359 rads), and T_jan(z) is the mean January temperature at elevation z based on the January linear lapse rate obtained by PRISM climate data (r² ~0.99). In implementing the model, ∆T is the prescribed LGM temperature perturbation. A tuning parameter k is used to adjust the sharpness of the temperature sinusoid so as to minimize the difference between modeled monthly temperatures during the ablation season (May through September) and those from the PRISM data. Because of valley-to-valley differences in climate, especially precipitation, subsets of the PRISM data vary accordingly. H(z), T_jan(z), and k(z) are also specific for each glaciated valley and in the case of SNOTEL stations location (see Section 4.2). Note that H and k are given here as functions of elevation, varying linearly, but in some valleys or regions the correlations are weak and therefore modeling was done with constant (mean) values of these parameters.

Monthly snow precipitation P_s(t, z) is partitioned from modern monthly precipitation P_mod(t,z) and the prescribed change in LGM precipitation ∆P via a continuous temperature-dependent function determined from the SNOTEL records:

P_s(t,z) = ƒ(T_month) (P_mod(t,z) + ∆P).

(5)

In detail, ƒ(T_month) has the form of a logistic curve relating the snow-to-rain ratio to mean monthly temperature and ensures that in some months a fraction of precipitation falls as snow, some as rain rather than all snow or all rain according to a threshold temperature. Modeling runs used the partitioning function based on proximity to a SNOTEL site, but the effect of using different functions was investigated and discussed subsequently.

4. Results

4.1. Glacier Reconstruction and Equilibrium-Line Altitudes

Attributes of the reconstructed paleoglaciers shown in Figure 3 are summarized in Table 1. Driving stresses are within limits of what is observed in modern glaciers. In general, larger stresses are associated with greater ice thicknesses and/or steeper ice surface slopes, the latter following the steepness of subglacial topography in the valleys. Conversely lower driving stresses are due to thin ice and/or more moderate slopes, and are often found near glacier termini.

Of particular note here are the equilibrium-line altitudes (ELAs) derived from area-altitude ratios (AARs). Regardless of the precise value of the AAR, the reconstructions suggest LGM ELAs on the eastern slopes in the study area are slightly more than 100 m lower than those on the west-facing slopes. Similar east-west asymmetries were recognized elsewhere in the Sangre de Cristo Mountains in Colorado [23] and in the Mosquito Range in central Colorado [14] that, like the Sangre de Cristo, is the most prominent topographic barrier for southeasterly-derived moisture transport. This apparent difference in ELAs on opposing slopes has been attributed to wind transport of snow from western to eastern slopes, either by deflation on the former [35] or advection prior to deposition [16,24], and/or by virtue of difference in net radiation.

4.2. Temperature-Index Modeling: Model Skill

Evaluation of model skill follows the criteria previously used [12,13,14,15,16]: how well does the model (1) capture modern climate, and (2) reproduce snowpack evolution at SNOTEL stations using parameterization that is unique to the SNOTEL sites. Modeled modern temperatures and precipitation at SNOTEL sites are graphically compared with those observed in Figure 4a and b. For the purpose of modeling, accurate temperatures are most critical during the ablation season (taken here as May through September) and precipitation during the accumulation season (November through April). Cumulative temperature differences during the former are ~0 °C, which is not surprising in view of the fact that the model is tuned (via H(z) and k(z) in Equation 4) specifically for the SNOTEL locations. Cumulative precipitation differences (no tuning) are very small and negative for all three SNOTEL sites, meaning accumulation season precipitation is very slightly overestimated in the model.

Model-data agreement is also quantified using a Nash-Sutcliffe Efficiency criteria commonly used in hydrologic modeling

$$NSE=1-\left[{\displaystyle \frac{\sum _{i=1}^n{(O_i-p_i)}^2}{\sum _{i=1}^n{(O_i-\overline{O})}^2}}\right]$$

(6)

where O_i are the observed values, p_i are those modeled, and the overbar denotes the mean of n number of values. A value of 1.0 indicates perfect agreement between modeled and observed values, whereas values approaching or less than 0.0 suggest a poor fit to the observed values. NSE values (Figure 4a,b) are all close to one and thus suggest model results faithfully match the observed SNOTEL data.

Similar comparisons between modeled modern temperatures and precipitation and PRISM values for the valleys of interest are summarized in Table 2. These comparisons are perhaps more significant in that they are more regional in nature and also account for elevational variations in differences. With one exception, modeled temperatures differ from PRISM values by less than 1°C and on average less than 0.5 °C. NSE statistics for temperature are ≥ 0.97 indicating nearly perfect agreement. Modeled precipitation values differ by ~1-2 cm or less than 5% of total precipitation. NSE scores for precipitation, with the exception of two, exceed 0.84 indicating good correspondence between modeled and PRISM values. Table 2 also shows that there is no systematic differences, that is the TM approach neither consistently over- or underestimates temperature or precipitation.

Finally, modeled snowpack evolution is compared with median 1991-2020 snowpack depth (in snow water equivalent) recorded at the three SNOTEL sites over a hydrologic year (Figure 4c). Here again, PRISM data, the H(z) and k(z) parameters, and snow-to-rain ratio used in the model are specific to each site. Despite difference in temporal resolution (monthly in the model, daily in the SNOTEL records), the agreement is quite good, particularly for the North Castilla and Whiskey Creek sites.

4.3. Steady-State Mass Balances of Paleoglaciers and Implications for LGM Climate

Estimates of LGM climate in the study area were determined by finding temperature and/or precipitation combinations that maintain steady-state mass balances B_n of the reconstructed glaciers, i.e. those that satisfy

$$B_n={\int }_A{b}_{n}dA\approx {\sum }_{i=1}^j{b}_{n_i}{A}_i=0$$

(7)

where b_n is the annual specific net mass balance given by Equation (2), A is glacier area composed of j number of elevation increments, and i denotes the value at an elevation increment, thus glacier hypsometry is explicitly included.

As there are an infinite number of solutions for Equation (7), limits must be placed on possible changes in precipitation in order to obtain estimates of LGM temperature depression. Proxy records of precipitation would serve as such limits, however such records are sparse in the region and typically do not provide the magnitude of that change. Nevertheless, pluvial lake levels and speleothem records, as well as climate modeling [1,4,36,37,38,39,40,41], point to wetter conditions in the American Southwest during the LGM. Therefore, modeling examined the required temperature depression under imposed increases in precipitation of 10 and 20 cm yr^–1 corresponding to ~10-15% and 20-30% increase over modern precipitation respectively, depending on elevation. It should be noted that paleohydrology of Lake Estancia, about 280 km south-southwest from the study area, suggest possible increases in precipitation at times during the LGM could have been as large as 50 to 100% over modern [36,39]. While a number of contrasting mechanisms have been invoked to explain apparent increases in LGM precipitation in the Southwest [1,4,42,43,44,45,46], the question remains whether these would have affected the study area. In particular, climate models [1,3,4] place the study area at or near the boundary between wetter conditions to the south and drier conditions to the north (the so-called “precipitation dipole”). Thus in the interest of completeness, addition modeling runs investigate the effect of reducing LGM precipitation by similar amounts.

Table 3 and Figure 5a summarize temperature depression required to maintain steady-state mass balance of the seven paleoglaciers under different scenarios. These include assuming LGM precipitation was similar to that today, increased or decreased precipitation, and differences in melt factors. Under the assumption of no change in precipitation and average melt factors, temperature depressions range from 8.7 to 9.8 °C. Figure 5a also shows the associated ELAs from modeling. The individual differences notwithstanding, mean temperature depression appears to be slightly greater, by ~0.5 °C, on the eastern slopes of the study area. Paralleling this is that ELAs tend to be ~110 m lower on the eastern slopes, consistent with what was noted above for the AAR-derived ELAs. Moreover, Figure 5b shows ELAs determined by modeling closely agree with using the AAR method. However, given these difference generally fall within estimates of uncertainties (discussed subsequently) and in the interest of determining regional LGM climate, the emphasis here is placed on the (grand) mean temperature depressions of ~9.6 °C for the standard melt scenario, and ~8.9 and 10.0 °C for the reduced and enhanced melt scenarios, respectively.

As discussed above, considering the possibility of increased LGM precipitation of 10 cm yr^–1 the mean temperature depression to maintain steady-state mass balances is ~9.4 °C, and for 20 cm yr^–1 ~9.2 °C. In combination with lower melt factors, these temperature depressions could be on the order ~1 °C less; conversely with higher melt factors ~0.5 °C greater (Table 3) and minimally dependent on these prescribed changes in precipitation (Figure 5c). Decreases in LGM precipitation of 10 cm yr^–1 would, on average, require temperatures to be ~9.8 °C cooler than present. A decrease of 20 cm yr^–1 would require temperature depressions of ~10.0 °C. Under the reduced and enhance melt scenarios temperature depressions would be ~0.7 less and ~0.4 °C greater, respectively, and again not strongly dependent on the prescribed changes in precipitation.

4.4. Uncertainties in the Estimates of LGM Temperature Depression

Uncertainties in the TM results were evaluated by performing sensitivity analyses to different model parameters. Specifically, these were parameters in Equation (3): the melt factors m_ƒ and melt threshold T_m; in Equation (4): k(z) – constant (average value) or non-constant, H(z) – constant or non-constant, the lapse rate T_jan(z); and in Equation (5): ƒ(T_month). Additional analyses examine the effect of the melt threshold and glacier hypsometry. Individual uncertainties were added in quadrature. Maximum uncertainties associated with these parameters from each of the seven paleoglaciers, where appropriate, were used to calculate total uncertainty.

Variation in the melt factors for snow and ice yield uncertainties in the estimates of temperature depression of about +0.5/–0.8 °C, corresponding to more or less melt. If the melt threshold is assigned a value of 0 °C instead of +1 °C, required temperature depressions increase by ~0.2 °C. Assuming the tuning parameter k(z) is constant (i.e. a mean value) or not (i.e. a function of elevation) introduces a very small uncertainty less than ± 0.05 °C. Similarly, whether the yearly amplitude in temperature H(z) is assumed constant with respect to elevation or not results in a uncertainties that do not exceed ~ ±0.3 °C. Reasonable variations in the lapse rate T_jan(z) based on the regressions introduces an uncertainty of ~ ±0.3 °C. Temperature depressions using different partitioning functions defined at each of the three SNOTEL sites differ by a maximum of ~0.3 °C.

The uncertainties associated with glacier hypsometry are more challenging to assess. Following Brugger et al. [15] a Monte Carlo scheme was used wherein glacier area within elevation intervals was allowed to increase or decrease randomly by 15%; this essentially addresses uncertainties in ice contours as reconstructed and thereby mass balance via Equation (7). This approach suggests uncertainties associated with the hypsometries of the paleoglaciers does not exceed ±0.2 °C. Based on the foregoing the total uncertainty in the estimates of temperature depression are on the order of ±1.0 °C.

5. Discussion

The results of temperature-index modeling suggest that regional LGM temperatures in the southernmost Sangre de Cristo Mountains in Colorado were on the order of 9 to 10 °C cooler even allowing for a modest increase in precipitation. With due consideration of uncertainties, the apparent difference in temperature depression of ~0.5 °C on east- and west-facing slopes (Figure 5a, Table 3) is likely not significant. Nevertheless, in view of the fact that ELAs, determined by both the modeling and using the AAR method, are significantly lower on the eastern slope warrants some attention. As noted above this has been observed elsewhere and possibly reflects less incoming solar radiation hence less ablation, and/or any mechanism whereby accumulation is enhanced [16,23,35] on the eastern slope. Under the assumption that precipitation differences existed during the LGM as they do today, say +10 cm yr^–1 on the eastern slope, this difference is reduced, and when combined with the assumption of less ablation (reduced melt scenario), LGM temperatures on the eastern slope would be ~9 °C cooler, again within uncertainties. Thus the estimates of LGM temperature depression presented here appear to be robust. It bears mentioning that under these assumptions, TM-derived ELAs on the eastern slope are at most ~20 m higher, thus maintaining accord with those using the AAR approach.

Perhaps the most intriguing result is that of the magnitude of LGM temperature depression. Estimates using a similar TM approach [12,13,14,15,16] and coupled energy/mass balance and flow models [7,8,11] in Colorado show neither a consistency in magnitude nor a coherent spatial pattern (Figure 6). (It should be noted that previous TM modeling referred to here used slightly different melt factors, but this accounts for a difference of ~0.1 °C.) For example, in the Front Range of northern Colorado, temperature depression has been estimated to be 4.5 °C assuming precipitation was similar to that today [7]. In contrast, temperature depression over the Blanca Massif in the Sangre de Cristo Mountains may have been nearly twice as large [16] with the same assumption regarding precipitation. This contrast runs counter to the southward decrease in LGM cooling evident in climate models [1,4,17], presumably due to less climatic influence of the Laurentide Ice Sheet. While our estimates of ~9-10 °C appear to support a more severe cooling during the LGM, at least in the Sangre de Cristo Mountains, a temperature depression of ~5 °C has been suggested for the Crestone Peaks region of the range, about 100 km north of the study area [8]. Moreover, recent work [48] using lipid biomarkers in bog sediment in the nearby San Juan Mountains, about 110 km west (Figure 6), yielded an LGM temperature 6.5 °C colder than today. However, (1) basal sediment age was estimated to be ~19.2 ka, potentially predating the LGM in the study area by several thousand years; and (2) this value is the departure from pre-industrial temperatures and therefore probably ~1 °C greater if comparing to modern conditions. Additional support for a colder LGM in the region comes from the southernmost Sangre de Cristo Mountains in New Mexico, where temperature depression was estimated to be 8.6-9.0 °C providing precipitation was similar to the present, or ~8.0 °C cooler if precipitation was greater than modern [10].

The disparate estimates of temperature depression in Colorado are likely due in part to local modulation of regional climate (i.e. the role(s) of microclimate), differences in the precise timing of local LGM (comparing climate at different times), and/or the particular methodology used to obtain the estimate. But these in-and-of-themselves are not likely to account for such large differences in the measure of LGM climate change. What is not known and therefore somewhat unconstrained in modeling, are changes in precipitation that could significantly affect glacier mass balance and consequently temperature depression. This further complicates the picture of regional climate change during the LGM.

Although the biomarker proxy in the San Juan Mountains [48] is indicative of drier conditions during the LGM, several climate models suggest possible increases in precipitation. An ensemble of nine PMIP3 models (Paleoclimate Modeling Intercomparison Project Phase 3 [1]) suggests a mean increase in precipitation within the study area of ~20%, or about 10-15 cm yr^–1. However, five of the models indicate a reduction of precipitation of similar magnitude, therefore rendering any conclusion with regard to precipitation changes based on PMIP3 models equivocal. More recent TraCE21k modeling (Transient Climate Evolution of the last 21k years) using the CHELSA downscaling (Climatologies at High resolution for the Earth’s Land Surface Area, [49]) suggests substantial increases in LGM precipitation in many mountain ranges in Colorado (Figure 7a). In the study area (Figure 7b) more modest increases in precipitation are indicated, on the order of 3 to 7 cm annually. The CHELSA-TraCE21k output also implies a disproportionate increase in precipitation on the western slopes, furthermore suggesting LGM precipitation was more uniform over the study area. This appears to be at odds with the ELA differences on eastern- and western-facing slopes discussed previously. If, however, the greater increase in precipitation on the western slopes resulted from intensification of the prevailing westerlies during the LGM [1,2,3,46,47], this could have also resulted in greater wind-drift and/or advection of snow that enhanced accumulation on the eastern slopes. (Obviously, differences in insolation and its affect on local energy balances would factor in, at least to some extent, to any explanation of these differences.)

Increased precipitation would have the effect of decreasing the magnitudes of LGM temperature depression, particularly in the northern mountain ranges in Colorado. But in the study area the small increase would not significantly affect our results, making it seem more of an outlier. Yet in examining the CHELSE-TraCE21k temperature output, regionally LGM temperatures (Figure 7c) might have been ~8 °C cooler in southwestern Colorado to ~14 °C cooler in north-central area of the state. In the Sangre de Cristo Mountains (Figure 7d), temperatures appear to have been ~8.5 to ~10 °C colder than present, the change perhaps being more pronounced at lower elevation possibly due to changes in lapse rates. Such departures from modern temperatures are also corroborated by PMIP3 models that show about an 8.7 ± 2.5 °C cooling in the Sangre de Cristo [1]. Furthermore, using the same approach as that used here, temperature depression on the Blanca Massif just 50 km to the north of the study area averages ~8.6 °C assuming no change, or ~8.3 °C with a slight increase in precipitation [16].

It is difficult to reconcile this magnitude of temperature depression with other estimates in the region based on glaciological approaches, particularly ~5 °C in the Crestone Peaks area of the Sangre de Cristo Mountains [8] (Figure 6) about 100 km away. Under reasonable optimal conditions for glaciers in the study area, i.e. modest increases in LGM precipitation and reduced melt (Table 3), and allowing for uncertainties, temperature depression might have been a small as ~7.6 °C. Conversely, in the absence of these conditions, our estimate of temperature depression is consistent with the CHELA-TraCE21k downscaling and PMIP3 models, and that for the southernmost Sangre de Cristo [10]. Lacking more definite temperature and precipitation proxies that would constrain modeling, the LGM in the study area was ~9.6 to 8.6 °C cooler than present.

6. Conclusions

The salient points based on the reconstruction of glaciers in the southern Sangre de Cristo Mountains of Colorado are:

(1) At the glaciers maximum LGM extents, ELAs on east-facing slopes were on average ~100 m lower than those on west-facing slopes. This east-west asymmetry in ELAs has been documented elsewhere in the Sangre de Cristo [23] and in the Mosquito Range farther north [14]. This is possibly a result of less insolation, hence ablation, and/or mechanisms through which accumulation might have been enhanced via wind drift or snow advection [16,24,35] or late winter/early spring upslope precipitation events [23].

(2) Temperature-index modeling was used to determine the temperatures and precipitation necessary to maintain steady-state mass balances of the glacier. Modeling suggests mean temperature depressions during the LGM of between 8.9 and 10.0 °C, this range being associated with different melt factors used in the model, and if LGM precipitation did not differ from that today. Considering modest increases in precipitation, the magnitude(s) of temperature depressions is slightly reduced by ~0.3 °C, well within estimated uncertainties of ±1 °C. Under the most reasonably favorable glacial conditions, temperature depression could have been as small as 8.6 or even 7.6 °C incorporating uncertainty.

(3) While the temperature depression obtained for the study area is significantly greater than some other estimates elsewhere in Colorado, it is consistent with those suggested CHELSA-TraCE21k downscaling [49] and PMIP3 ensembles [1]. With due consideration of the associated uncertainties, it is also in agreement with the mean temperature depression determined for the nearby Blanca Massif [16] and the southernmost Sangre de Cristo Mountains in New Mexico [10].

(4) Because of their dependence on knowledge of paleoprecipitation, models seeking to understand the nature of LGM temperature change in the Colorado Rocky Mountains would be better constrained by, and benefit from new and/or improved precipitation proxies.