Introduction

Big Timber Canyon in the Crazy Mountains of Montana at the western margin of the Great Plains ( Fig.1) is a prime example of a glacially sculpted valley with its U-shaped valley and well-marked trimlines. Because of these well-defined features and its straight channel, the canyon provides an excellent location in which to study glacier dynamics using flow-law theory. By reconstructing the former glacier profile and areal extent, calculations of velocity and mass flux at the equilibrium-line altitude can be made; from these, the glacier’s annual mass balance can be calculated.

The mass balance of a glacier is controlled by the climate at the glacier’s location. Therefore, if the mass balance of a paleoglacier can be calculated, a reconstruction of the paleoclimate at the time when the glacier existed may be possible. Mass balance is usually thought of as a balance between winter accumulation and summer ablation, which are usually investigated as functions of winter precipitation and summer temperature, respectively. However, there are many variables in the microclimate of a glacier that influence these factors including accumulation-area aspect (shading), avalanches, and amount of wind-blown accumulation. If these variables can be held constant by selecting valleys in a region with similar aspect and topography, a model that calculates mass balance of paleoglaciers can be used to produce a regional estimation of paleoclimate.

Fig.1. Location, map view, and longitudinal profile of the late Pleistocene Big Timber glacier.

This study is the test case for such a model using the late Pleistocene Big Timber Canyon glacier. The following steps were undertaken for this process: (1) reconstruction of the glacial profile and areal extent from morphologic evidence and published data, (2) calculation of basal shear stresses along the longitudinal profile, (3) calculation of flow rates from ice-flow equations, and (4) estimation of mass balance from the calculated mass flux through the equilibrium-line altitude (ELA). During steps 2, 3, and 4, sensitivity and error analysis was performed by varying parameters such as glacier thickness and slope.

Areal Extent and Longitudinal Profile

The areal extent of the paleoglacier ( Fig.1) was determined from USGS topographic maps (scale 1 : 24 000), aerial photographs, and published documentation (Reference AldenAlden, 1932; Aten, unpublished). Geomorphic evidence used to estimate ice extent included terminal moraine locations, lateral moraine elevations, truncated spurs, elevations of hanging tributary valleys and glacial trimlines ( Fig.1). Only deposits and scouring features from the Pinedale I advance identified by Aten (unpublished) were used in this study, since this was the last major glaciation of the valley and these features are best preserved. Since only one maximum Pinedale moraine was identified (Aten, unpublished), it is thought to correlate with the last glacial maximum approximately 20 000 to 15 000 years ago (Reference Porter, Pierce, Hamilton and PorterPorter and others, 1983). Prior glaciation limits have been obscured by this advance and subsequent glaciations were minor.

From this evidence, the longitudinal profile along the axis of the glacier was drawn ( Fig.1). In the upper reaches of the valley, depositional and scouring features were not as prevalent. In this area, the profile was determined from a model (Reference Schilling, Hollin, Denton and HughesSchilling and Hollin, 1981) which calculates ice elevations using a theoretical average basal shear stress (see next section). In the middle reaches of the valley, there was good correlation between the geomorphic limits and the theoretical limits using an average basal shear stress of 1.0 bar (Table I). It is thus considered reasonable to extend this correlation into the upper reaches.

In the lower reaches of the valley, ice-surface elevations were assumed to be indistinguishable from the elevations of the ice-marginal features. Comparison of the surface contours of several modern-day glaciers in Alaska (Anonymous, 1960) which are similar in plan form to the Big Timber glacier showed a mean difference of only +3.9 m between the ice center line and ice-marginal elevations below the equilibrium line ( Fig.2). Because of this minimal variance, no correction factor was added to the ice-surface elevations below the equilibrium line. Ice center-line elevations are lower than the ice margin above the equilibrium line of modern glaciers. Since the theoretical ice elevations were lower than the ice-marginal features on the upper parts of Big Timber glacier (Table I), the theoretical ice elevations were used where the trimlines were not well defined.

Ice thicknesses were calculated from the difference of ice-surface elevation and bedrock elevation. One of the underlying assumptions is that the bedrock elevation has changed little since late Pleistocene time. Because there have been subsequent glaciations (Aten, unpublished), there has been minor scouring of the bedrock in the upper reaches. Thus, the ice thicknesses in this part of the valley are maximal. In the lower reaches of the valley, re-advance and recessional till and post-glacial fill make the estimated ice-thickness values minimal. From the mid-point of the glacier on up the valley, the present stream runs along ice-scoured bedrock, so fill does not present a problem in estimating ice thicknesses.

Table.1. Paleoglacier morphology and rheology interpreted from topographic maps and comparison with theoretical values

Fig.2. Difference in elevation between ice center line and ice-marginal features of modern glaciers plotted against elevation relative to the ELA (after Anonymous. 1960).

Reconstruction based on this evidence shows that the glacier was approximately 18 km long. Approximately 14 km up the valley from the terminus, the main glacier split into two branches which headed in mainly north-east-facing cirques. Three small cirque glaciers fed the main glacier between 10 and 13 km from the terminus. Maximum thicknesses were reached between 9 and IOkm up the valley from the terminus (Table I), where the glacier was 325–350 m thick.

Basal Shear Stresses

Where ice thicknesses were known from geomorphic evidence, average basal shear stresses were calculated along the glacier center line. For a valley glacier, the average basal shear stress (τ_b) can be calculated from the relation:

(1)

(Reference PatersonPaterson, 1981) where ρ is specific gravity of ice (910 kg m⁻³); g is acceleration due to gravity (9.81 m s⁻³); H is center-line ice thickness (m); F is shape factor to account for drag on the valley sides (Reference NyeNye, 1965); and α is ice-surface slope (averaged over 2.4 km). From laboratory and field data, the theoretical values for τ_b should range from 0.5 to 1.5 bar (Reference PiercePierce, 1979; Reference PatersonPaterson, 1981). Values for τ_b under Big Timber glacier where ice-marginal features were well defined (4.0–12.0 km from the terminus) ranged from 0.88 to 1.15 bar (Fig. 3; Table I). In the lowest reaches of the glacier (0–4.0 km from the terminus), unusually low values of basal shear stress (0.31–0.51 bar) were calculated. The addition of melt water to the glacier bed in this region from Amelong Creek could have significantly increased basal sliding. In this situation, the effective basal shear stress would be very low and the glacier-surface gradient and thickness would have been lowered. Ice-contact features have also been modified by stream erosion and mass wasting since they were deposited. Therefore, reconstruction of the glacier in this region is tentative.

A theoretical profile of the glacier was calculated using the relation:

(2)

which was derived from Equation (1) (Reference Schilling, Hollin, Denton and HughesSchilling and Hollin, 1981). This iterative approach was used to calculate glacier thickness (H) at 1000 ft [305 m] intervals (x) along the glacier profile. Here e_i and e_i+1 are the ice-surface elevations at steps x_i and x_i+1 respectively, and the other variables are the same as in Equation (I). The profile calculated using an average basal shear stress of 1.0 bar provided good agreement with the ice-marginal features (Table I).

Fig.3. Calculated effective basal shear stress and mass flux vs distance from the termius. Slip flux is the difference between deformation and continuity fluxes; see text for discussion.

Location of the Ela

Once the glacial extent and longitudinal profile were determined, the location of the equilibrium-line altitude (ELA) could also be determined. The ELA is an important descriptor of any glacial system, since it is the point where mass balance changes from net accumulation to net ablation. At this point, flow lines are also parallel to the glacier surface. For a steady-state glacier, the amount of mass that accumulates should equal the amount of mass that ablates, and this should also equal the amount of mass that flows through the ELA during the same time period (Reference AndrewsAndrews, 1975). If the ELA of a glacier can be determined, glacial mass balance can be determined by calculating the glacier velocity and multiplying it by the cross-sectional area at the ELA.

The paleo-ELA of the late Pleistocene Big Timber glacier was estimated from (1) highest lateral moraines, (2) lowest cirque-floor elevation, (3) toe-headwell area ratio (THAR), and (4) accumulation-area ratio (AAR) after Reference MeierdingMeierding (1982). The elevations estimated using these methods ranged from 2163 to 2432 m with a mean of 2252 m.

Lateral moraines

Lateral moraines develop in the ablation zone where ice flows outward towards the margins of the glacier and deposits debris. Since ice-flow lines are descending in the accumulation area, the highest point on the lateral moraine would indicate a change from descending to ascending flow and would therefore be an approximation of the location of the ELA (Reference AndrewsAndrews, 1975; Reference MeierdingMeierding, 1982). Lateral moraines are well developed in the lower reaches of Big Timber glacier and extend about 8.8 km up the valley. The elevation of the highest of these moraines is 2170 m, and the ice thickness at this point is 323 m. Mass wasting of the moraines since they were deposited makes this estimate minimal.

Lowest cirque elevation

Cirque-floor elevations have been widely used as a measure of the ELA of former cirque glaciers (Reference AndrewsAndrews, 1975; Reference MeierdingMeierding, 1982). Reference MeierdingMeierding (1982) showed that, although this method is rapid, it is also highly subjective because cirques floors are not always easily identifiable. Reference AndrewsAndrews (1975) pointed out that this method is better for a regional approximation of ELA if a trend surface connecting the lowest north-facing cirques is constructed. Also, since valley glaciers extended outside the cirques, ELAs using this method on valley glaciers should provide a maximum estimate. In Big Timber valley, there are many well-defined cirques. The lowest occurs in the north branch of the valley at an elevation of 2432 m.

Toe-headwall altitude ratio (THAR)

An empirical relation between the highest and lowest ice limits has been used extensively for the rapid determination of ELA. Meierding (1982) found ratios of 0.35 and 0.40 produced the best results for the Colorado Front Range. Highest ice limits were found on the north face of Crazy Peak at an altitude of 3170 m. Since the bedrock surface at the terminus was covered with till, it was estimated to be at 1622 m. This was inferred from the elevation of the sub-alluvial terrace upon which the till rests (Aten, unpublished). ELAs using ratios of 0.35, 0.40, and 0.45 ranged from 2163 to 2318 m (Table II).

Table.2. Ela locations using each of the methods listed in the text

Accumulation-area ratios (AAR)

Studies of modern glaciers have shown that the accumulation area of a glacier is about 0.65 of the total area. However, this percentage may vary between 0.60 and 0.70 (Reference AndrewsAndrews, 1975). Several studies of paleoglaciers have used the 0.65 value (Reference PorterPorter, 1975; Reference MeierdingMeierding, 1982; Reference LeonardLeonard, 1984) to locate the ELA. Meierding (1982) calculated ELAs using values ranging from 0.50 to 0.75 and found 0.65 to have the least error. In this study, the area between successive 400 ft [122 m] ( Fig.1) contours was planimetered and a cumulative total was plotted ( Fig.4). From this graph, ELAs using AAR values of 0.60, 0.65, and 0.70 were determined. ELAs using this method ranged from 2170 to 2286 m.

ELAs using a THAR of 0.35, the highest lateral moraine, and an AAR of 0.70 were within 7 m of each other (Table II). There was also good agreement between the THAR of 0.40 and the AAR of 0.65, both of which Reference MeierdingMeierding (1982) found to give the best result. Since mass wasting of the lateral moraines would produce a low estimate for the ELA, and since Meierding found the latter methods produced the best estimate, an ELA of 2240 m was used in this study. By comparison, ELAs on the reconstructed Amelong Creek and Middle Fork Big Timber glacier ( Fig.1) were approximately 2430 m, using THAR = 0.40 and AAR = 0.65.

Fig.4. Hypsometric cune of Big Timber glacier constructed by planimetry of the 400 ft [122 m] contours of the reconstructed glacier ( Fig.1). Elevations of ELAs for various AARs (0.65 shown) were calculated from this curve.

Glacier-flow Velocities

By calculating the glacier-flow velocity at the ELA and multiplying by the cross-sectional area of the glacier at that point, an estimate of mass flux (therefore ablation and accumulation on this steady-state glacier) can be made. Glacier movement can be broken down into components of ice deformation and basal slip. By integrating Glen’s flow law of ice (Reference GlenGlen, 1952) over small increments of glacial thickness, Reference NyeNye (1952) was able to calculate a profile of the velocity due to deformation in a glacier. Assuming that the basal velocity is zero, the center-line surface velocity due to ice deformation (V_c) can be calculated from:

(3)

where A is the temperature-dependent constant of the flow law (0.167 bar⁻³ a⁻¹ at 0°C; Reference PatersonPaterson, 1981); H is ice thickness at the center line (m); τ_b is average basal shear stress (bar; calculated from Equation (1)); and « is exponential constant of the flow law (Equation (3); Reference PatersonPaterson, 1981). Longitudinal variations in velocity are minimized by averaging ice-surface slope over a distance of 8–20Я (Reference Raymond and ColbeckRaymond, 1980). Local values of H are used in this study because they were the same as the average values in most cases. Substituting Equation (1) for τ_b shows that velocity due to deformation varies as H⁴ and sin α³. Therefore, estimation of these values must be as precise as possible.

The constant A is a function of the ice temperature and thus should vary through the glacier thickness and also with altitude. In a temperate or sub-polar glacier, the basal temperature should be at the pressure-melting point, approximately 0°C. Since most of the motion due to ice deformation occurs near the base of the glacier, and since the temperature in this region should be near the pressure-melting point (0 to −1°C), A is assumed constant in this study.

The other component of glacier flow is basal slip. This has been measured to account for 10–90% of total velocity (Reference AndrewsAndrews, 1975; Reference PatersonPaterson, 1981). Sliding velocity in present glaciers can be directly observed in bore holes or calculated by subtracting the velocity due to deformation from the observed surface velocity. Several models have been developed (summarized in Reference WeertmanWeertman, 1979; Reference Raymond and ColbeckRaymond, 1980; Reference PatersonPaterson, 1981) to explain basal slip but comparisons with observations do not usually agree with theory. Since actual surface velocity cannot be compared with calculated creep velocity on paleoglaciers, estimation of basal slip for these glaciers presents a problem.

Basal slip may occur to varying degrees along a glacier depending on the thermal and moisture regime. Even sub-refreezing conditions do not prevent basal slip (Reference Echelmeyer and ZhongxiangEchelmeyer and Zhongxiang, 1987). However, any glacier must have a point at which basal slip is at a minimum. This point can be identified by calculation of effective basal shear stress, which will generally be the greatest where slip is least, and deformation mass flux, which will be greatest where slip flux is least.

Because shear stresses may vary due to extending and compressing flow arising from irregularities in the bedrock slope, this theory will not work in a multi-stepped glacial valley. In Big Timber Canyon, basal slope does not vary greatly up the valley ( Fig.1). Therefore, the place where the highest effective shear stress occurs should be the place where slip is at a minimum. For Big Timber glacier, the highest basal shear stress coincidcntally occurs at the ELA ( Fig.3).

Assuming no slip at the ELA, a calculated effective basal shear stress of 1.15 bar and an ice thickness of 351 m, the center-line velocity at the ELA was 44.5 m a⁻¹. This is well within the range of velocities observed on modern glaciers (Reference PatersonPaterson, 1981). Since the velocity is sensitive to changes in ice-surface slope and thickness (which were estimated from topographic maps), the actual creep velocity may vary by as much as 40% due to errors in the calculation of these variables. However, the probability of such a large error is very low (p = 0.00006) and the value for velocity is considered accurate (see Appendix).

Mass-balance Estimates

Since yearly accumulation, yearly ablation, and yearly mass flux through the ELA are all equal for a steady-state glacier, mass balance can be estimated if the velocity at the ELA and the cross-sectional area of the ELA are known. Minimum center-line velocity has already been calculated; therefore, multiplying this value by a factor to average velocity over the cross-section and then multiplying by the cross-sectional area provides a calculation of minimum mass balance.

Due to frictional drag on the valley walls, the mean velocity of a glacier through a cross-section is less than the center-line surface velocity. In a no-slip environment, the ratio of average velocity through the cross-section to the surface velocity is 0.63 (Reference NyeNye, 1963; Reference Raymond and ColbeckRaymond, 1980). Since it is assumed that there is no slip at the point of highest effective basal shear stress (the ELA on this glacier), multiplying 0.63V_C by the cross-sectional area at the ELA (determined from topographic maps) yields an approximation of yearly mass balance.

Calculated mass flux (mass balance) for Big Timber glacier is 8.8 × 10⁶ m³ a⁻¹. Mass flux using the theoretical profile and a shear stress of 1.0 bar yields a mass flux of 5.9 x 10⁶ m³ a⁻¹. Since H only differed by 1.5 m, the major difference is in slope. A constant effective basal shear stress of 1.0 bar implies a lower slope (2.7°) than the observed values (3.1°).

Ablation/Accumulation Gradients

Mass balance provides a key to the activity of a glacier. A measure of a glacier’s activity is the ablation gradient, which can be defined as the relationship of summer balance to elevation (Reference AndrewsAndrews, 1975). High ablation gradients (>10mm/m) are typical of glaciers in maritime environments, while these gradients decrease inland towards more continental environments (Reference Meier, Tangborn, Mayo and PostMeier and others, 1971; Reference AndrewsAndrews, 1975). An exception to this trend is the occurrence of high gradients on small cirque glaciers in the Rocky Mountains, which can be accounted for by the local microclimates that produce large amounts of snow accumulation from wind drifting and orographic precipitation (Reference Meier, Tangborn, Mayo and PostMeier and others, 1971). However, on moderate-sized glaciers ablation gradients should be indicative of climate.

Since the mass that moves through the ELA (8.8 x 10⁶m³ a⁻¹) must be ablated on this steady-state glacier and since the area between successive contours is known, an ablation gradient can be determined (Fig. 5 and Table III; after Reference PiercePierce, 1979). For Big Timber glacier, an ablation gradient of 3.0 ± 0.6 mm/m produced this balance. An accumulation gradient of 1.0 ± 0.2 mm/m above the ELA was calculated using the same method. These low gradients are typical of cold, low-moisture environments (e.g. McCall Glacier, 2 mm/m; Meier and others, 1971) showing that the glacier at its peak was sustained more by low temperatures than high precipitation. This is in contrast to the higher ablation gradient (9.0 mm/m) that Reference PiercePierce (1979) determined for Yellowstone ice cap which lay just to the south. The addition of moisture from Snake River Plain to Yellowstone Plateau could account for this difference (Reference Porter, Pierce, Hamilton and PorterPorter and others, 1983; Reference Locke and KemphLocke and Kemph, 1987). Using a model similar to the one in this study, Reference Leonard, Huston and ManleyLeonard and others (1986) calculated very low gradients for late Pleistocene glaciers in the Colorado Front Range, showing that those glaciers were also sustained by a cold and dry climate during that time.

Fig.5. Mass balance of Ihe late Pleistocene Big Timber glacier. The reconstructed surface area (left) and the estimated mass flux at the ELA allow the calculation of average mass-balance gradients (center), thus specific (center) and net balance (right).

These gradients also allow an estimation of basal slip around the point of assumed no-slip, which in this case occurs at the ELA. Using the theory of continuity, the difference between the mass flux at one cross-section and a down-ice cross-section should be equal to the amount of mass that is ablated or accumulated over the glacier surface area between the two cross-sections. If the mass flux resulting from deformation through a cross-section down-ice from the ELA • is calculated assuming no slip, and the amount of ablation calculated from the gradient and surface area is less than the differences between the mass flux atthis point and the ELA, then the remaining difference is assumed to be from mass flux due to slip.

Table.3. Average slip velocities below the ela calculated from continuity

Table III and Figure 3 show the estimated slip up- and down-ice from the ELA on Big Timber glacier. Slip becomes the dominant part (>90%) of the flow in the lower 4.0 km of the glacier. Since Amelong Creek was supplying melt water to the glacier bed at this point, a substantial increase in basal slip is logical. Even with the accepted error in the deformation mass flux on this glacier ( Fig.3), continuity theory shows that basal slip had to be the dominant part of flow in this part of the glacier. Slip is also required above the ELA, and was responsible for at least 50% of the mass flux.

Conclusions

A reliable reconstruction of the areal extent of late Pleistocene glaciers can be made using topographic maps and aerial photographs if evidence such as moraines and trimlines is preserved. From this reconstruction, calculations of basal shear stress, ice velocity, and mass balance can be made. Allowable error in the topographic maps affects the precision of these calculations, but the probability of a 20% error in calculated velocity (thus mass flux) is only 1% (see Appendix). Field survey only needs to be performed to check that the valley floor is bedrock and not fill.

Calculations of effective basal shear stress using the data derived from the maps can be used as a check on the validity of the reconstruction. Normal effective basal shear stresses should range from 0.5 to 1.5 bar if the reconstruction is valid. For Big Timber glacier, a 1.0 bar effective basal shear stress matched observed ice-margin indicators well across most of the glacier. The lowermost 4.0 km of the glacier were apparently slip-dominated.

Sensitivity of the flow-law equations to changes in ice slope can be minimized using a slope averaged over a distance of 8–20 times the thickness. Local variations in surface slope do not affect glacier velocity over short distances. This averaging also reduces the effect of errors in the topographic maps by making the error a smaller percentage of the distance over which slope is averaged. For Big Timber glacier, slope was averaged over a distance of 8000 ft [2.4 km].

For a valley with a floor of constant or slowly varying, gentle slope, mass balances can be calculated assuming no slip at the point where the highest effective, basal shear stress occurs. This point should occur above or near the ELA where lower temperatures should reduce the amount of basal melt water. If this point occurs above the ELA, the amount of slip at the ELA can be calculated using the continuity equation and an appropriate accumulation gradient. Mass flux due to slip and ice deformation at the ELA can be calculated and then an ablation gradient can be determined for the paleoglacier. For Big Timber glacier, the estimated accumulation gradient was 1.0 ± 0.2 mm/m and the ablation gradient was 3.0 ± 0.6 mm/m. These estimates are minima, as slip may occur even at the point of maximum deformation mass flux.

Comparison of calculated ablation gradients to modern analogs provides an estimation of the climate that occurred when the glacier was at its peak. For the late Pleistocene Big Timber glacier, a modern analog is McCall Glacier in Alaska with an ablation gradient of 2 mm/m (Meier and others, 1971). Higher gradients would show areas where more moisture is available (Reference Haeberli and PenzHaeberli and Penz, 1985). By selecting valleys of comparable aspect, size, and relief, microclimatic variations can be minimized. Using valleys with relatively constant, gentle slope and few tributary cirques, variations in basal shear stress due to extending and compressing flow are also minimized. If this method is applied to a region using valleys of similar aspect and topography, regional variations in paleoclimate may become apparent.

The model developed on Big Timber glacier is being used as an estimator of paleoclimate for the northern Rocky Mountains of Montana and Idaho. Careful selection of study valleys will reduce the errors associated with the model and therefore increase the accuracy of the climatic reconstruction, which can be compared to other independent estimates of late Pleistocene paleoclimate (Reference Barry and PorterBarry, 1983; Reference Locke and KemphLocke and Kemph, 1987).

Acknowledgements

The concept on which this paper is based was developed in unpublished manuscripts by E. Leonard (Colorado College) and I. McCalpin (Utah State University). This paper was inspired by those men and improved by their unselfish assistance.