Introduction
In the summer of 1959 the University of Alberta and the University of British Columbia undertook an expedition to the Athabaska Glacier to study particular problems in the geophysics of glaciers. The program was sponsored by the National Research Council of Canada and was under the supervision of Professors J. A. Jacobs and G. D. Garland. The major factor determining the selection of the glacier to be studied was ease of access. The Athabaska Glacier was an obvious choice as it is a smooth valley glacier with a well-defined terminus which is only 0.5 miles (0.8 km.) from a major paved route. It is also 65 miles (105 km.) from the town of Jasper on the main railroad line between Edmonton and Vancouver.
The major part of the program was designed to determine the mechanics of glacier flow. In this connection 80 stakes were set out to study surface movement and a series of bore holes was drilled. During the 1960 season more bore holes were drilled and several were cased for inclinometer measurements. The thickness of the ice and the depth of the bottom were determined at 32 points by reflection seismic methods. The thickness of the ice and the shape of the bedrock along transverse profiles were to be obtained by a gravity investigation using a Worden gravity meter. The gravity and seismic studies were planned to provide a general description of the thickness of the glacier so that movement studies could be properly evaluated and the average shear stress on the bedrock calculated. In addition, the gravity data checked areas with poor seismic reflections and depth points determined from events suspected of being multiple in origin.
In addition to the above studies, Professor G. D. Garland and D. Lennox made temperature measurements in uncased bore holes. Independent investigations were carried out by other parties on the same glacier during the 1959 season. These included a detailed aerial photographic and topographic survey by the Water Resources Branch of the Dominion Government. Various experiments to determine electrical properties of an ice sheet by electrical and electromagnetic methods were conducted by A. Don Watt and Gene Maxwell of the United States Bureau of Standards. Dr. G. Keller from the United States Geological Survey made depth observations using electromagnetic methods.
Physical Description of the Area
The Athabaska Glacier is the second largest ice stream issuing from the Columbia Icefield, an ice cap covering an area of 110 square miles (285 km.2) and straddling the continental divide which marks the boundary between the provinces of Alberta and British Columbia. The glacier issues from the Columbia Icefield at an elevation of about 9,000 ft. (2,740 m.) and flows 4 miles (6.4 km.) north-eastwards to a well-defined terminus at 6,320 ft. (1,926 m.) above sea-level. The glacier is the source of the Sunwapta River which joins the Athabaska, Slave and Mackenzie Rivers to flow 2,200 miles (3,540 km.) into the Arctic Ocean. The terminus of the glacier has shown a slow general recession and an estimate of the rate of recession can be made from topographic maps and photographs made by Cautley and Wheeler in 1919. These indicate a recession of 2,700 ± 200 ft. (823 ± 61 m.) in 40 yr. or an annual mean recession of 68 ± 5 ft. (20.7 ± 1.5 m.).
The Athabaska Glacier flows across the axis of a long anticline. The anticlinal axis crosses the glacier at line B just below the first step and the limbs dip away from the axis at 10–15°. The predominant rock type which outcrops in this area is limestone of Middle and Upper Cambrian age. From measurements on 28 rock samples it was concluded that a density of 2.67 g./cm.3 could be used in making gravity corrections. An examination of the individual measurements indicates that this figure is correct within about 2 per cent.
The Gravity Measurements and their Reduction
The observations on lines A to F were made with the University of British Columbia’s Worden gravity meter No. 35. Later in the season the University of Alberta’s Worden gravity meter XPO arrived from factory servicing and observations were made with it on lines B to H inclusive. There are, therefore, two independent readings on 85 of the 127 stations occupied. A comparison of the results from the two gravity meters indicates that errors between stations on one line do not exceed 0.1 mgal.
Standard procedures were used in surveying and in the reduction of the gravity measurements, and these will not be discussed here. To correct for the topography Hayford’s values for the radii were used together with tables given by Reference BullardBullard (1935). These tables were not accurate enough for modern surveys and new tables were calculated following Bullard’s method. The topographic effect of mountains and valleys was corrected out to a distance of 12.4 km. (zones A to J). Estimating the heights in each of the 96 sectors surrounding a station involved the greatest amount of labour. Care had to be taken to distinguish between sectors containing ice which has a density of 0.90 g./cm.3 and sectors containing rock with a density of 2.67 g./cm.3. The reductions of terrain corrections for 125 stations required about one month to complete. The corrections were very large because of the proximity of Mt. Athabaska (3,491 m.), Snow Dome (3,456 m.), Mt. Kitchener (3,505 m.) and the Sunwapta valley (1,920 m.). The maximum correction was 14.9 mgal and the minimum was 6.31 mgal.
Accuracy of the data
The relative field measurements on each line are accurate to within 0.1 mgal. The elevations are known to a precision of ± 1 ft. (0.3 m.) so that the combined free air and Bouguer corrections are accurate to ± 0.06 mgal. The principal error and the one that is most difficult to evaluate is that due to the terrain corrections. The uncertainty in the density contributes an error of about 2 per cent. This produces an error of 0.1 mgal between a station in the middle and a station at the sides of the glacier. For each station there are 96 sectors each of which is subject to an error of about 0.1 mgal. These reading errors are in the nature of a one-dimensional “random walk”. In the measured quantity there are 96 causes of error, each of which will add the amount + 0.01 or − 0.01 mgal with equal probability. The probability of an error of ± 0.04 mgal is only 0.2, while the probability of an error as large as ± 0.1 mgal is 0.004.
In the preceding discussion it has been tacitly assumed that the topographic maps are perfect. In some of the zones a determinate error may be introduced due to a systematic distortion of the contours. A large-scale map on the scale of 1 : 4,800 was constructed and contoured in the field. Zones B and C, which cover an area with a radius of 230 m., were entirely corrected using this map and it is unlikely that any significant error was introduced here. Zones F to J, which cover an annular area 2.29 to 12.40 km. from a station, were also corrected on maps of adequate scale and accuracy. The topographic maps on a scale of 1 : 62,500 were not adequate for correcting zones D, E and F. The magnitude of the error can only be roughly determined from plausible changes in the shape of the contours. It is estimated that the error is no more than ± 0.2 mgal. The terrain corrections are, therefore, thought to have a precision of at least ± 0.4 mgal and the precision of the Bouguer anomaly is ± 0.5 mgal.
The Interpretation of Gravity Anomalies over a Glacier
The interpretation of gravity anomalies involves an evaluation of the mass distribution causing the anomaly. A unique interpretation cannot be found as this is a problem in potential theory, and one can rearrange the sub-surface mass distribution in an infinite number of ways to reproduce the same gravitational field. In the present problem there is additional information which places restrictions on the interpretations possible. The density difference producing the anomaly is known within 2 per cent. If the shape of the bedrock profile were known, the central depth to bedrock would be the only unknown variable and it could be determined exactly from a knowledge of the gravitational field. Sometimes the depth of a glacier will be known from a bore hole or from seismic work. In this case the shape of the bedrock profile can be determined approximately. The detection of small variations on the bedrock would be limited by the accuracy of the gravity data.
Previous determinations of the thickness of a glacier with gravity data have been made by Reference Bull and HardyBull and Hardy (1956), Reference ThielThiel and others (1957), and Reference RussellRussell and others (1960). Topographic maps were not available so the terrain corrections had to be estimated with two-dimensional models. The depth determinations were limited to the smooth, central parts of the glaciers, since a two-dimensional technique was used in the computations.
A two-dimensional model of the glacier
Geomorphologists have described glacier valleys as being U-shaped or parabolic. Reference SvenssonSvensson (1959) has demonstrated that the glacial valley Lapporten in northern Sweden is very nearly parabolic in shape. It was, therefore, decided that a good first approximation to the shape of the bedrock beneath the Athabaska Glacier would be parabolic in form. An elliptical shape was also modelled and found useful when the side walls were very steep.
The process of finding the depth and shape of a glacier is essentially a trial and error process. The gravity anomaly of a two-dimensional model with arbitrary cross-section can be found by a method of numerical integration using a graticule chart. This method is not particularly accurate and is quite tedious if a large number of shapes, depths and observation points are used. The theoretical gravity anomaly of a two-dimensional body with a parabolic cross-section can be approximated as closely as required by the repeated use of the formula for the gravity of a rectangular prism (Reference HeilandHeiland, 1940, p. 151). The vertical component of gravitational acceleration at the point P due to the shaded rectangle in Figure 1 is given by
Fig. 1. Symbols used in the calculation of the gravitational attraction of a two-dimensional body at point P. The parabolic prism is approximated by a series of rectangular prisms. D and d are distances to the top and bottom of a particular prism and b is its width
The gravitational constant is denoted by the symbol k and the density difference between bedrock and ice by the symbol δ. The expression given above was programmed for the Royal McBee LGP–30 digital computer. The program was designed to choose a given depth and observation point, calculate the equation of the parabola and fit a series of rectangular prisms, each 50 ft. (15.2 m.) thick, to the parabola. The computer then calculated equation (1) for each of the prisms, summed these and printed out the answer. The process was repeated for all observation points and for each maximum depth.
Figure 2 gives a comparison of the shape and theoretical gravity anomaly of a rectangular, elliptical and parabolic cross-section for the glacier at the approximate position of the bore hole near line C. One of the points to notice is that the gravity anomaly is still 4 to 6 mgal where the glacier has zero thickness. The zero-thickness line can seldom be determined very accurately and it is often so close to the steep mountain sides that a gravity station cannot be located on it because the terrain correction will be uncertain by several milligals. Figure 2 also indicates that the gravity anomaly has the greatest rate of change near the edge of the glacier. These considerations will usually make it unwise to tie the Bouguer anomaly from the field data to the end points. The Bouguer anomaly does not give information about the value of the zero ordinate for the glacier anomaly. There is merely a relative anomaly from the centre to the end station. The centre points have been tied together in any comparison of gravity anomalies from field or model studies. Figure 2 indicates that, if the depth to the centre is known, the gravity data can predict the shape with considerable success.
Fig. 2. (a) Models used for a comparison of the theoretical gravity anomalies of three two-dimensional bodies. The bodies are 1,305 m. wide along the surface, 312 m. deep and extend to infinity along the third direction. The density contrast between the ice and the underlying bedrock is 1.77 g./cm.3. (b) Theoretical gravity anomalies for the three bodies illustrated in Figure 2a
Table I gives the results of a two-dimensional computer program for determining the depth of the glacier along line C with the assumption that the cross-section is parabolic in shape. For a determination of the best fit, the Bouguer anomaly of the centre station was tied to the theoretical gravity value. The root mean square of the deviations from the theoretical values was computed for each of the maximum depths using the 15 observation points on the line.
Table I Results of Two-dimensional Computer Program Line C—15 observation points
From a seismic program and from two nearby bore holes which reached bedrock, the depth is known to be 1,050 ± 30 ft. (320 ± 9 m.). On the basis of the above analysis a depth of 1,100 ft. (335 m.) would have been chosen.
The procedure illustrated above can be used in a systematic way on the central part of any glacier for which the observed gravity anomaly is known. The expression for the gravity of a rectangular prism can be adapted to an elliptical cross-section with equal ease. If the Bouguer anomalies are asymmetrical, a similar degree of asymmetry can be programmed into the basic parabolic or elliptical cross-section and the maximum depth determined from the least root mean square of the deviations.
Computer program for a body with finite dimensions
Once a two-dimensional approximation has been found for the glacier, the body of ice can be treated as an irregular, three-dimensional “ore” body whose gravity anomaly can be computed by a method developed by Reference Talwani and EwingTalwani and Ewing (1960). A cylindrical coordinate system (r, ψ, z) is employed with the z-axis vertically downward and the origin at the gravity station. The gravity anomaly due to a horizontal lamina of infinitesmal thickness dz is
where
In Cartesian coordinates, r 2 = x 2 + y 2.
Equation (3) was programmed for the LGP–30 digital computer to obtain the gravitational attraction of the glacier at any station following the method of Reference Talwani and EwingTalwani and Ewing (1960). These equations also form the analytic basis of a mechanical integrator for the computation of gravity anomalies as developed by Reference SiegertSiegert (1942).
A contour map of the model of the glacier is drawn up and each contour is replaced by a horizontal irregular polygonal lamina as shown in Figure 3b. The contour interval used was 100 ft. (30.5 m.). The polygon should fit the contour map closely if an observation point is nearby. A glacier is a special case of an “ore” body, because the observation points straddle the “ore” body itself. Therefore, only the laminae below the observation point need be computed. The material above the level of the observation point has been accounted for by the terrain correction.
Fig. 3. (a) Three-dimensional geometry for the evaluation of the gravity anomaly caused by a density contrast between ice and bedrock. (b) Two-dimensional polygonal lamina for particular elevations taken from the glacier model in Figure 3a. (c) Numerical integration necessary to obtain the gravity anomaly at a point A-5 due to a series of horizontal laminae at depth z. Data for V(z) were obtained from the computer program
The input data to the computer consist of the coordinates and elevations of all the observation points and the coordinates and elevations of each of the corners of the polygonal lamina. The computer program solved for the function V(z), which is the gravity anomaly at a point due to a horizontal lamina at a depth z. The gravity anomaly caused by the whole glacier at any observation point was determined by integrating equation (2) using the graphical procedure illustrated in Figure 3c.
Two models were computed for each of the eight lines. For the majority of the lines, the second model was 20 per cent shallower up-dip from the line in question. In addition, two extra models were determined for lines B and D to provide a better fit. The total computer time for 45 observation points for the two models was
Results
The Bouguer gravity map of the Athabaska Glacier is shown in Figure 4. The results of the computer program are presented in Figures 5, 6 and 7. In each of these figures the top diagram is a transverse cross-section of the glacier, whereas the various other diagrams are models computed for the purposes of comparison. The lower two graphs show a comparison of the observed values with the theoretical gravity anomaly of the model. The observed Bouguer gravity values did not have to be corrected for a regional gradient because the cross-sections are along the strike of the geologic structures. The final model was generally interpolated from the computed model to give a better fit than is illustrated by the model anomaly curve. Figure 8 is a longitudinal profile of the glacier along the central line.
Fig. 4. Bouguer anomaly map of the Athabaska Glacier, Alberta, Canada
Fig. 5. (a) Cross-section of the Athabaska Glacier along line C. Model I is the best fit according to the gravity data. Horizontal and vertical scales are equal. A comparison of bore hole and seismic results is also shown. (b) A comparison of theoretical and observed gravity anomalies along line C for Model 1. (c) A comparison of theoretical and observed gravity anomalies for Model 11. This model was rejected because of the poor, fit
Fig. 6. (a) Cross-section of the Athabaska Glacier along line E. Model III was determined as the best fit from a study of gravity results on Models I and II. Horizontal and vertical scales are equal. (b) Theoretical and observed gravity anomalies for line E. (c) Cross-section of the Athabaska Glacier along line H. Model III was determined as the best fit. (d) Theoretical and observed gravity anomalies for line H
Fig. 7. (a) Cross-section of the Athabaska Glacier along line F. Horizontal and vertical scales are equal. (b) Theoretical and observed anomaly along line F. Model II provides the best fit. (c) Cross-section of the Athabaska Glacier along line G. Horizontal and vertical scales are equal. (d) Theoretical and observed gravity anomaly along line G. Model IV is a slightly modified version of Model III to provide the best fit between theoretical and observed gravity results
Fig. 8. (a) Gravity anomalies along line L of the Athabaska Glacier, The Bouguer anomalies along the glacier are tied into Dominion Observatory gravity station 459 (lat. 52° 12.9′ N., long. 117° 12.3′ W.). If the glacier were entirely removed, the calculated gravity anomaly remaining would be along the curve marked “Residual Bouguer Effect”. (b) Longitudinal profile (along line L) of the Athabaska Glacier. A comparison of gravity, seismic and bore-hole results is shown. Vertical and horizontal scales are equal
For a comparison of the accuracy, the figures include results from the seismic and borehole program carried out by J. C. Savage. The gravity results appear to be accurate within −10 to + 15 per cent. Lines E and H are 10 per cent shallower than the seismic data indicate. The longitudinal seismic section indicates that the bedrock is quite rough and, as expected, gravity data tend to smooth out the irregularities. Only bore holes which are thought to have reached bedrock are shown in the figures.
The value of three-dimensional computations is best illustrated on lines near the terminus of the glacier. It was possible to assign an absolute value to the anomaly produced by the field data by means of a station near the snout of the glacier. The terrain corrections here are certainly accurate within 0.2 mgal and the absolute depth should be correct within 10 per cent. The seismic record near G (Fig. 8b) has a poor reflection at a depth of 490 ft. (149 m.). As illustrated in Figure 8b, this is incompatible with the gravity results. Models were programmed for a depth of one-half and one-third of the seismic depth to check on multiple reflections. On the basis of the gravity study the reflection appearing on the record is a second multiple. The primary reflection should appear at 0.027 sec. (plus phase lag due to filtering). The primary reflection and its first multiple are probably obscured by surface waves.
Figure 8b represents a longitudinal cross-section of the glacier without any vertical exaggeration. Figure 8a shows the Bouguer anomaly along line L. It gives a qualitative representation of the thickness of the glacier. It is now possible to subtract the gravity anomaly of the glacier from the Bouguer anomaly and obtain a residual Bouguer anomaly which should reflect the underlying geology. If the glacier computations on each transverse line have been properly made, the remaining curve on the longitudinal profile should be very smooth and become more negative to the south-west. This negative trend into the Rocky Mountain Trench is shown on a Bouguer gravity map of British Columbia and Alberta (Reference Garland and TannerGarland and Tanner, 1957). The residual anomaly is indeed as described above and reflects the structural geology of the underlying section.
Application of the Results to Glaciology
In a glacier of arbitrary but constant cross-section, Reference NyeNye (1952[b]) has shown that the average value of shear stress, τ a , at the bed is given by the relations
and
where A is the area of the cross-section perpendicular to the bed, and P is the perimeter of this cross-section. The equation is approximately true if the cross-section varies slowly and the current value of R and the surface slope, a, are used. The gravitational constant, g, was taken to be 981 cm./sec.2 and the density, ρ, was 0.90 g./cm.3. The average shear stress on the bedrock is given in Table II.
Table II. Average Shear Stress on the Bedrock
The bedrock on line F is rough and the line is close to the terminus, so the conditions of the theory may not be satisfied. If this low value is omitted, the average shear stress exerted by the bed on the ice for the five upper lines becomes 1.0 bar. Reference NyeNye (1952[a]) has made calculations for sixteen alpine glaciers and finds values ranging from 0.49 to 1.51 bar. The values are in the middle of Nye’s range and show remarkably small variation.
Conclusions
The thickness of the Athabaska Glacier has been obtained along eight transverse lines by an investigation of the gravity anomalies. Previous workers obtained the depth with a precision that varied from − 10 to + 25 per cent or in some cases + 40 per cent. A higher degree of accuracy was attained in this study by making detailed terrain corrections and by treating the glacier as a three-dimensional body. The Bouguer anomalies were determined within 0.5 mgal or better and the depth values are believed to be accurate within − 10 to + 15 per cent. Previous depth determinations were limited to the smooth, central part of a glacier since two-dimensional calculations were used. Three-dimensional computations with a low-speed computer were made in this study and the depths near the terminus of the glacier and in the upper elevated section were obtained.
The major disadvantage of the method is that it requires a large-scale contour map of the area to reduce the terrain corrections. This problem is rapidly being overcome for most mountainous regions by the use of aerial photogrammetry. Making the terrain corrections is a long and tedious operation but this can be overcome by the use of Siegert’s gravity integrator or by the use of a digital computer with a large storage capacity.
In making depth computations, a useful first approximation to the cross-section of the valley of a glacier is a parabola. The cross-section of the valley through which the Athabaska Glacier flows has been found to approximate a parabola along several of the lines.
A basin has been formed in the centre of the glacier, probably by some form of erosion, so that the glacier now flows along a surface that is flat or slopes slightly upwards toward the terminus. Combined gravity and seismic data suggest a similar up-valley slope on the first step of the glacial stairway. If the glacier should disappear entirely, the chain of basins would form “paternoster” lakes. The glacier has thickened to its maximum extent in the central region, so that it now has a surface slope which allows it to flow along the nearly level valley floor.
From a knowledge of the depth, shape and surface slope of the glacier, the average shear stress exerted by the bed on the ice was found to be 1.0 bar.
Acknowledgements
The author wishes to thank Professor G. D. Garland who directed this research. Dr. J. C. Savage of the University of British Columbia arranged for the use of that institution’s gravity meter and kindly allowed me to use his unpublished seismic data. Messrs J. Stacey, S. Paterson, H. Vetter and J. Fairley assisted in many phases of the field work. W. S. Adams of the University Computing Centre assisted in programming the material for the digital computer. The work on which this paper is based was done at the University of Alberta, Edmonton.
Appendix
Note on Seismic Results over the Athaisaska Glacier by J. C. Savage
The seismic work on the Athabaska Glacier was carried out using standard reflection techniques. A twelve-trace high-resolution seismograph manufactured by Houston Technical Laboratories was used. A 200 m. geophone spread was employed and three geophones were used for each trace. The standard shot was a half stick of 60 per cent dynamite set in a water-filled drill hole at a depth of 3.5 m.
Thirty-six reflection sites were occupied. At each of these sites at least one spread transverse to the glacier flow and one spread along the direction of glacier flow were shot. Usable records were obtained at all but four of these sites. As a general rule high-quality records could be obtained for ice depths in excess of 200 m. The records became marginal at ice depths less than 130 m.
A 2 km. refraction line was also shot. The velocity of the P-wave was 3,610 m./sec. (± 10 m./sec. standard error). The velocity in the glacier bed was determined to be about 4,500 m./sec.
