Introduction

“Fox Glacier”Footnote * is a small valley glacier in the Yukon Territory. Canada; it is located between Steele and Hodgson Glaciers and shares a common drainage basin with “Jackal” and “Hyena Glaciers” (Fig. 1). All five of these glaciers have surged and only “Fox” and “Hyena Glaciers” are at present inactive. The fragmented condition of an actively surging glacier makes surface measurements impractical. In anticipation of a future surge of “Fox Glacier”, an extensive study of the glacier, co-ordinated by the Icefield Ranges Research Project, Footnote † was begun in 1967. In 1968 the reported gravity survey was carried out to find the ice thickness; near-surface temperature measurements were also made.

The earliest work in the “Fox Glacier” area was by the Wood Yukon Expedition in 1935, which conducted a program of reconnaissance aerial survey during the ascent of Mount Steele (Reference WoodWood, 1936). Reference SharpSharp (1943) reported a geological study of the Steele Glacier valley which included some of the terrain to the north-east of “Fox Glacier”.

Recent scientific investigations include surface-flow surveys by S. Collins of the American Geographical Society, mass-balance studies by T. Brewer of Boston University, and unsuccessful attempts at seismic and radio depth soundings by the authors. In view of the gravity results, the failure to obtain seismic and radio echoes is probably due to the thinness of “Fox Glacier”. The hydrology associated with the melt-water drainage of “Fox”, “Jackal” and “Hyena Glaciers” was investigated by T. Faber of the Canadian Department of Energy. Mines and Resources (Inland Waters Branch). K. West of the University of Alberta has analysed oxygen isotope ratios of ice cores. Morainal geology has been examined by G. Denton of the American Geographical Society.

In 1968, temperature measurements were taken to depths of 3 m at stakes 2, 16 and 34, and to 8 m at stake 7 (for location of these stakes see Figure 2). Mechanical drilling and the absence of anti-freeze solutions minimized thermal contamination. At stake 7 the minimum temperature obtained was —5-5°C at 3 m on 10 July; by 16 August the temperature at this point had risen to — 4.1°C. Similar temperatures were obtained at the other stakes. An analysis of the effect of seasonal variations in the temperature at the glacier surface indicated the low temperature is not due to the annual cold wave; for this reason “Fox Glacier” was classified as a sub-polar glacier.

The gravity survey

The instrument used in the survey was a Sharpe gravity meter (No. C132) with a factory calibration constant of 0.101 35 mgal /division.

Figure 2 is a sketch of the glacier showing the positions of 66 stakes which had been drilled into the ice during the summer of 1967. These stakes served as the basis for the gravity network. This network consisted of 22 closed loops which included, in addition to the stakes, nine survey stations situated on elevated rock outcrops around the edge of the glacier, and 14 rock cairns set up on bedrock at the ends of glacier traverses.

Fig. 1. Locution of “Fox Glacier”, Yukon Territory.

As the gravimeter was a short-range instrument with a scale width of 1000 divisions (about 100 mgal), several re-scaling operations were necessary to cover the elevation difference between the lowest and highest stations. Insufficient time prevented the tying of the network to an absolute pendulum measurement in the vicinity.

The UTM (Universal Transverse Mercator) coordinates and elevation above sea-level of the glacier stations were determined from a theodolite survey in 1968 to an absolute accuracy of ±1 m for the UTM coordinates and about ±0.02 m for the elevations. The coordinates and elevations of the edge stations were obtained from a less precise independent survey, within ±1 rn and an estimated ±0.15 m, respectively.

Since the main survey determined the elevation to the top of the stakes, a correction was applied to allow for the height of the stake above the snow or ice surface. Uncertainties in the stake heights were due to measurement errors and fluctuations in the ablation rate; the lowering of the glacier surface between gravity readings on different loops was ignored. Combined errors in the station elevations amounted to ±0.05 m for the glacier stations and ±0.18 m for the edge stations. The lack of any obvious correlation between loop misclosures and elapsed time indicated that drift was small compared to random measurement errors. Misclosures were adjusted by the weighted least-squares method of Reference GibsonGibson (1941). By taking the r.m.s. deviation of 48 calculated minus measured differences between junctions, it was estimated that the least-squares values were not in error by more than ±0.68 divisions (±0.07 mgal).

Fig. 2. Sketch map of “Fox Glacier”, showing Bouguer anomalies. Terrain correction removed.

As noted by Reference ClarkeClarke. (1967) and others, rock density in the vicinity of glacierized areas is difficult to obtain accurately, since accessible rock outcrops tend to be more resistant to erosion and hence are not typical of the subglacier bedrock. For this reason, no sampling of rock was attempted for density and the commonly quoted value of 2.67 Mg/m³ was used in the calculations. For the type of rock described by Sharp for the “Fox Glacier” area, this density does not seem at all unreasonable. Outcrops on “Fox Glacier” were found to be predominantly basalt and andesite; for these rock types a typical density range is 2.6-2.9 Mg/m³. Unconsolidated morainal material with lower density was also present. Measurements of ice density from the surface of “Fox Glacier” yielded an averaged value of 0.89 Mg/m³, but allowing for ice compaction at depth a value of 0.90 Mg/m³ seems better. The resulting density contrast of 1.77 Mg/m³ is probably accurate to ±0.20 Mg/m³.

Conventional corrections for latitude and elevation were applied. Details of these corrections and tables of results can be found in Crossley (unpublished).

Topographic corrections

In any survey where the topography near a gravity station differs from a horizontal plane through the station, a correction has to be made for the attraction of the material above and below that plane. In mountainous terrain, this attraction can be of the order of several milligals and can become the least certain addition to the gravity readings. An initial inspection of the simple Bouguer anomalies at stations on the edge and middle of “Fox Glacier” suggested that the difference of 2.5 mgal between these stations was the same magnitude as the expected terrain corrections.

The topography surrounding “Fox Glacier” consists of high lateral moraines, ice-cored to uncertain depth, and glaciers to the south and east occupying deep valleys. The lack of symmetry in the geometric shapes of the terrain precluded two-dimensional integration for the attraction. Due to the limited extent of the survey, hand computations using the zone system of Reference HammerHammer (1939) were considered suitable. The attractions of the zones were computed for both rock and ice densities to allow for compartments which included both materials. There are several ways of making these corrections and the method adopted is given in more detail.

Figure 3 shows a gravity station P situated on uneven terrain, (i) is at the edge of the glacier and (ii) is on glacier ice. In both cases the Bouguer correction compensates for rock material up to the datum through P and the terrain correction must not remove the anomaly due to the ice. In case (i) the effect of the rock above P and the deficiency of the material between the datum through P and the ice surface were both removed as rock. The same correction for the deficiency below P still applies for case (ii). The material above p now consists of ice and rock and was removed according to the different densities and fractions lying within a single compartment. It is clear that corrections made in this way cannot be made without an initial assumption about the depth of ice to be removed in case (ii), and for this reason 50 m was assumed to be the depth of all glacial ice. Naturally, the depth of “Fox Glacier” will differ from the 50 m allowed in the terrain corrections; however, glacial ice was present on much of the surrounding terrain, Since the depth of this ice was indeterminate, it was decided to take 50 m, obtained by considering the anomaly clue to an infinite slab, as a first approximation for all visible ice. Corrections for distant zones were not treated using the zone scheme due to the lack of a suitably scaled topographic map, but instead were incorporated as an additional slowly varying component into the regional gradient. Within 90% confidence limits the error in the terrain correction was ±0.37 mgal for a station.

Fig. 3. Topographic correction detail.

Fig. 4. Residual anomalies (mgal). Regional gradient removed.

The final correction to the observed anomalies was an estimate of regional variations in bedrock density, geologic structure and distant terrain effects in the form of a regional gradient. The mathematical form of (his gradient can be approximated by a three-dimensional geometric surface of any degree, but to describe non-local effects low-degree functions are normally used. In theory a good estimate of the regional gradient can be obtained by examining the anomalies at stations distant from the ice body. In the present survey, none of the stations was free of the effect of the glacier and hence an unbiased estimate of the regional gradient was not immediately possible.

The regional gradient was obtained by least-squares fitting of a polynomial to stations situated near the edge of the glacier. By subsequently allowing for the attraction of the ice at these stations, it would be possible to approximate by iteration an unbiased estimate of the gradient. A third-degree polynomial gave a good fit to the edge stations; for a higher-degree polynomial, with all possible terms, there would have been more least-squares parameters than edge stations. Figure 2 shows the Bouguer anomalies with the terrain correction applied and Figure 4 shows the residual anomalies

Here g_t is the Bouguer anomaly at station i and F(x_i,y_i) is the regional gradient.

Fig. 5. Part of integration scheme. Stations created for network are shown with a cross.

Intrrpretation

The residual anomalies represent the local density variations due to glacier ice thickness and fluctuations in bedrock density. For practical reasons these two effects become indistinguishable in any interpretation and the accuracy of the agreement between gravity and other methods of depth estimation will reflect variations in the bedrock density.

The estimated uncertainties in the Bouguer anomalies include errors accumulated from least-squares fitting of the initial data, standard reduction errors, terrain correction errors and regional gradient errors. Totals for the first three of these are 0.39 and 0.40 mgal for the glacier and edge stations, respectively (90% confidence limits). Errors arising from the regional gradient fitting are difficult to estimate; an r.m.s. deviation of 0.507 mgal reflects random errors and the simplicity of the least-squares function. Assuming that these two causes equally contribute to the r.m.s. deviation, 0.25 mgal might reasonably be combined with the random errors. Assuming a Gaussian error distribution, 95% confidence limits on the mean error are found to be ±0.67 mgal.

Fig. 6. Glacier depth (m).

It is clear from the glacier sketch (Fig. 2) that a two-dimensional model of “Fox Glacier” would be a poor approximation and therefore a three-dimensional model was sought. Due to errors in the data, the problem of finding a satisfactory solution becomes one of adjusting some model to the residuals to within the accuracy desired. For a three-dimensional body of non-geometric shape, the integrated mass effect requires dividing the body into simpler shapes which can be handled analytically. The computerized sum of a number of horizontal laminae formulated by Reference Talwani and EwingTalwani and Ewing (1960) has been widely used, and once the coordinates of the lamina have been obtained the calculation is fairly rapid. This method has the disadvantage that each time an approximation is made to the shape of the body these coordinates have to be re-determined to allow integration to proceed.

Fig. 7. a. Transverse depth profiles (rows 12 and 16). Rock cairns are denoted by Δ.Elevation h is above mean sea-level.b. Transverse depth profiles (rows 20 and 24). Doubtful lines are dotted.

Fig. 7.1. c. Transverse depth profiles (rows 23 and 32). d. Transverse depth profiles [raws 34 and 36). Glacier edge is shown by ↓. e. Longitudinal depth profile. Vertical exaggeration 2.54.

For “Fox Glacier” an alternative method of integrating a depth distribution was attempted, based on the attraction of a number of triangular vertical prisms. The vertices of each prism at the surface of the glacier are the known coordinates of the gravity stations and remain fixed. The z coordinate of a vertex at the glacier bed is obtained from an initial model, with the depth of a prism taken as the average of three vertices. For steeply dipping bedrock this averaging contributes the main source of error.

Dividing the glacier in this manner resulted in 173 prisms; the main pre-computational task was relating the prism vertices to the gravity stations. Figure 5 illustrates a favourable aspect of the system; any number of extra stations can be placed around the glacier edge without adding unknown depths. Thirty stations were added in the manner shown in Figure 5. The calculation of the attraction of a single prism is outlined in the Appendix and this is based on the Talwani and Ewing formula for a triangular lamina. For stations far enough away from a prism, a center-of-mass approximation was used; the minimum radius for which this approximation is valid was determined by trial and error to be 0.7 km.

Integration of the glacier in this way required an initial depth distribution and this was obtained by an infinite slab assumption for each station. Figure 6 shows the resulting depth distribution; cross-sections at the glacier traverses (Fig. 7a-d) and a longitudinal profile (Fig. 7e) are also shown. The mean deviation between observed and calculated anomalies for the glacier stations is 0.44 mgal which is within the calculated confidence limits (0.67 mgal),

An adjustment to the depth distribution would normally have been made to reduce the r.m.s. deviation but, since the infinite slab depth distribution gives calculated anomalies which fit the data, no further adjustment was considered necessary. Adjustment of the depths by two methods was attempted to test the efficiency of the prism integration procedure. The first adjustment was obtained by replacing g_i in the original infinite slab anomaly calculation by the error Δ g _i obtained from observed minus calculated anomalies. Although an improvement was readily obtained, continued application showed thai the convergence was slow.

Reference CorbatóCorbató (1965[a]), noting this slowness, has given an application of the least-squares approach for improving the fit. Due to the lack of a manageable integration method and the large number of simultaneous equations involved, he remarked that a three-dimensional least-squares solution is cumbersome (Reference CorbatóCorbató, 1965[b]). An advantage of a least-squares technique is that the desired solution, if it exists, is obtained in a single step. The coefficients in the resulting set of simultaneous equations are partial derivatives of the residuals with respect to the depth at each station (Reference CorbatóCorbató, 1965[a], p. 229, equation (5)). In the prism integration these partial derivatives can be obtained rather easily since they are related to the quantity V₃ , the gravity anomaly per unit thickness at the base of a prism.

Fig. 8. Sample of least-squares adjustment (pecked).

The application of this procedure required the solution of 70 simultaneous equations but the results were unstable and physically unacceptable (Fig. 8). This was taken to indicate that either the data were inconsistent or that the least-squares parameters had too many degrees of freedom.

The depths obtained for the tongue region of the glacier are unreliable, but this is due to the small physical dimensions of the tongue and the relatively large terrain corrections, There is some doubt about the depths below the ice falls on the west side of the glacier, and little confidence should be placed in the results for this area.

Preliminary drilling results became available (personal communication from Mr Classen) from work completed during the summer of 1969. The predicted and actual depths are compared below. Errors in the gravity depths were estimated from the 95% confidence limits on the mean error in the anomalies combined with the uncertainty in the density contrast.

Acknowledgements

We wish to thank the following people and organizations for their assistance: the seismic apparatus was provided by Dr W. A. Weeks and the U.S. Army Cold Regions Research and Engineering Laboratory, Hanover, New Hampshire; Dr J. Tanner and the Gravity Division, Dominion Observatory, Ottawa, generously loaned a gravimeter at short notice. The Icefield Ranges Research Project gave invaluable logistic support; special thanks go to Dr M. Marcus the Director, and Phil Upton the pilot. Sam Collins was an excellent leader at Fox Camp; David Harrison and Patrick Powell were commendable field assistants. The help given both by Dr A. Stanley, Inland Waters Branch, Department of Energy, Mines and Resources, Ottawa, and Dick Ragle was greatly appreciated. Financial support was obtained from the Arctic Institute of North America and the National Research Council of Canada.