Introduction
Athabasca Glacier is located at 52°12ʹN, 117°14ʹW, inside Jasper National Park on the Alberta—British Columbia border (see Fig. 1). It flows in a north-northeast direction off the Columbia Ice Field over three icefalls into an alpine valley. Debris covers approximately one-third of the western margin of the ablation area. Athabasca Glacier was mapped many times using different methods between 1919 and 1979. These maps were converted into raster (grid square) digital elevation models (DEMs) and used to calculate the surface and volumetric change of the glacier over that period of time.
The main part of the ablation area lying between the terminus (at ≈1950 m) and 2400 m a.s.l. was studied. It is in this zone that the greatest changes have occurred and that the precision of measurement by the methods of mapping used is highest.
Glacier dynamics, the influence of bedrock configuration on glacier motion, and the relative importance of surface melt and emergence velocities are not included within the scope of this paper. It was initially hoped that the results would be of relevance in the study of kinematic waves. However, it was found that the original data are not sufficiently precise to reliably detect surficial changes associated with kinematic waves.
Data Sources
Using DEMs to calculate and show how the glacier has changed over time requires accurate, quantitative information portraying the glacier surface at different times. For Athabasca Glacier, this is stored in a series of contour maps. Table 1 summarizes the information concerning the production of the maps used for this study.
Glacier margins
The studied part of the glacier consists of flowing debris-covered and clear ice, with stagnant debris-covered ice in lateral moraines surrounding the glacier. Changes in stagnant ice were not measured. Therefore, the lateral extents of the flowing ice were estimated for each of the maps used in this study using both the maps themselves and aerial photography from the nearest year. figure 2a, d and e demonstrate that the method of edge estimation used was not very successful since the borders recorded do not display a logical progression of change between years.
Georeferencing
On each map there are printed points which state explicitly their location in Universal Transverse Mercator coordinates or latitude and longitude. If georeferencing is the same on all maps then they can be compared without difficulty. Unfortunately, the georeferencing of the small-scale maps did not agree with that of the large-scale maps. Due to its superior labelling of georeferencing points and since its referencing is determined in relation to the national georeferencing network, the 1955 map is assumed to be the most positionally accurate.
The small-scale 1919 map required a slight stretch in the north–south direction and a small displacement to bring it
Fig. 1. Athabasca Glacier, (a) Location (from Reid and Charbonneau, 1981). (b) Ablation area (Parks Canada, 1980).
into agreement with the 1955 map. The large-scale maps shared an internal georeferencing but were displaced slightly with respect to the 1955 map. To be brought into agreement with the 1955 map, they were displaced 161.4 m in the direction 66.7° east of south. Neither the 1919 map nor the large-scale maps displayed any rotational errors with respect to true north.
Generating The Surfaces
The contour lines of the source maps were digitized, producing a series of (x, y, z) elevation points from each map. The number of points digitized varied with the scale of the map, having 4000–7000 points for small-scale maps and greater than 11 000 points for each of the large-scale maps. For each of the maps, the surface of the glacier and the valley walls surrounding it had the greatest density of digitized points, while the surrounding area had few data points, widely dispersed. The effect of this was that the extra-glacier surface interpolation was less accurate. In a high-relief area, a slight horizontal displacement results in a large vertical non-coincidence.
The package Surfer 3.01 for Windows was used to generate raster surfaces from the (x, y, z) data files. Surfaces were created by overlaying a regular rectangular grid having spacing of 5 m over the mapped area of the glacier. The elevation at each of the grid nodes was interpolated from surrounding data points using universal Kriging with a linear semi-variogram. This method was relatively fast and produced surfaces that when compared with the original data points had residuals ranging in magnitude from 0.09 to 0.34 in over 50 m elevation zones of the glacier.
Change And Comparison: Results
The results reported in this paper do not measure mass balance, since information was neither collected on a yearly basis nor acquired at the end of the ablation season. Thus the values of surface and volumetric change derived do not reflect changes between one hydrological year and the next: change is measured between one map and the next.
The results of surface change in terminal position will be considered first, followed by a discussion of volumetric change.
Terminal position
figure 2a illustrates a retreat of the central part of the glacier terminus of ≈ 1km between 1919 and 1979. The former extent of the glacier and its rates of recession were determined by Luckman (1988) using old photographs, dendrochronology and the dating of recessional moraines.
The positions of the clear-ice terminus were determined solely from the terminus positions marked on the maps used to produce DEMs. The terminus of the debris-covered part was marked on no map. Air photo interpretation was used to estimate the terminus of the flowing debris-covered glacier in each of the years of mapping.
Additionally, terminus positions were not located with sufficient precision on the large-scale maps. Field observations at the time of record (Luckman, 1988) indicated that no readvances occurred, but some of the maps show slight readvances over limited parts of the terminus.
Area
As the glacier receded from 1919 to 1979, its area below 2400 m a.s.l. also decreased. One of the most notable trends is its great loss of area below 1950 m. In 1919, more than 100 000 m2 of its planar area was below 1950 m a.s.l. In 1979, less than 1500 m2 of its area was below that elevation. The elevation zone 1950–2000 m a.s.l. records a similar but not as precipitous decrease in area.
When 50 m elevation zones between 2150–2400 m a.s.l. are considered, measured areas in given zones arc fairly constant. The greatest yearly variability of area occurred in the elevation zone 2200–2300 m, the flattest section of the glacier.
Table 1. Summary of data sources
Fig. 2. Surface change of Athabasca Glacier. (a) Borders. Note non-coincident sections, especially on western margin. (b) Elevation change, 1919–55. Note apparent upgrading on terminus of glacier, (c) Elevation change, 1955–70. (d and e) Comparison of 1919, 1955 and 1979 cross-sections. Thick lines represent surface as digitized. Thin lines surrounding represent photogrammetric margin of error. Italicized years attached to arrows represent the edge of the glacier as interpreted for that map. Note non-coincidence of valley walls on the cross-sections. Since valley walls are steep, a small horizontal displacement results in a large vertical non-coincidence.
Elevation change
To calculate elevation change, nodes of identical (x, y) location are compared between DEMs. Figures were prepared showing the distribution of elevation change over all of the nodes of a glacier surface. figure 2b and c show total elevation change of the glacier over the intervals 1919–55 and 1955–79 and indicate that elevation change is heterogenous over the glacier surface. The elevation change that the glacier displayed between 1919–55 (fig. 2b) records an anomalous rise at the terminus, despite substantial volumetric loss in that interval. The rise as recorded is spurious, an artifact of the horizontal and vertical inaccuracies associated with the mapping, detailed below.
The amount of change registered at each of the nodes was recorded and averaged to find the mean elevation change of the surface over a given interval (table 2).
Volumetric change
Total volumetric change between successive mappings is summarized in Table 2. The change between the successive maps is calculated within the borders of the earlier of the paired maps. Change has not been consistently negative. There have been periods of mass increase within an overall pattern of glacier shrinkage. The annual rate of volumetric loss between 1919–55 is −6.25 × 106 m3 year−1, while the rate of volumetric change between 1955–79 is less, being −3.92 × 105 m3 year−1. This is in agreement with the record of terminal recession discussed in Luckman (1988).
The results of volume and elevation change recorded in the intervals 1969-71 and 1971–73 show a great decline in volume and elevation between 1969–71 followed by a sharp gain between 1971–73. This is believed to be due to a systematic error in the production of the 1971 map, since it is
Table 2. Total volumetric change of the complete ice surface, terminus–2400 m, 1919–79
difficult to conceive of a process that could cause this sharp decline and rebound effect.
Magnitudes Of Uncertainty
To quantify measures of uncertainty for both elevation and volumetric change, the precision of measurement of the source maps must be considered for the corresponding DEMs. Several different methods (listed in Table 1) were used to prepare the contour maps. Calculations of vertical uncertainty assume that the effects of horizontal uncertainty are small and that the uncertainty of vertical estimation at a point is primarily due to the distance from the point to the photographic station, given consistency in the optics of the cameras used. The effects of horizontal uncertainty are estimated with respect to georefereneing corrections.
Vertical uncertainty
Since distance to the photographic station is the primary factor in determining the precision of vertical estimation, a series of equations that describe the magnitude of uncertainty at given points on the glacier can be produced. Maps produced using aerial and terrestrial photogrammetry have different equations of uncertainty (table 3).
The two varieties of terrestrial photogrammetry used have different precisions of angular measurement Readings from the 1965–79 terrestrial photogrammetry (cf. Reid and Charbonneau, 1981) can be made to 1ʹ of arc (Young and others, 1978). The precision of angular measurement made from the 1919 photography was not stated in Cautley and others (1924), but was assumed to be 0.0846° of arc from the focal length of the camera and the value estimated for the precision of measurement of the equipment used in the early 1920s.
Of the measurements made using aerial photogrammetry, the 1959 map was studied by Paterson (1966), whose findings suggested that its precision was 2.0 × 10−4 of flying height. No special verification was carried out for the 1955 map; thus, its precision of vertical estimation was taken to be the mean of the range of other empirical values produced for similar measurements as reported in Young and Arnold (1978): 4.5 × 10−4 of flying height. Camera altitudes were markedly dissimilar: the 1955 photography was taken at a height of 10 700 m a.s.l. and the 1959 photography was taken at 4110 m a.s.l.
These equations of uncertainty are applied to the entire area of the DEM representation of the glacier. Every raster grid node on the glacier surface was treated by the appropriate equations. This produced surfaces of uncertainty for each map. The surfaces of uncertainty produced were added to and subtracted from surfaces representing the elevation of the glacier for every year of measurement. Thus, the glacier for each year of mapping is represented by three surfaces. When different mappings are compared, values of elevation and volumetric change are calculated for each of the nine combinations of the six surfaces. The uncertainty of measurement, a measure of dispersal, is the standard deviation of the differences. figure 2d and e show the surface of the glacier for the cross-sections AB and CD (shown in figure 2a–c) for the years 1919, 1955 and 1979. The surfaces of uncertainty associated with each year of mapping are displayed as thin lines surrounding the surface as digitized.
With quantities representing magnitudes of uncertainty, the imprecision of measurement can be compared to the values of elevation and volumetric change. This is shown in Table 2, which sums the values of uncertainty of measurement from the terminus to 2400 m for each of the intervals between successive mappings. In most intervals the uncertainty
Table 3. Accuracy of vertical estimation for different maps or map series
Table 4. elative magnitudes of error sources
Correction of geo-referencing Limited number of points to interpolate from Not dealt with “Minimizing differences” rather than rubber sheeting Estimated the effects of a 10 m displacement Area and volumetric calculations Depends on fineness of grid Finer grid spacing produces more precise measures of area and volumetric change of volumetric change was greater than the calculated volumetric change. This is due to heterogeneity of change over the surface: part of the surface may be downwasting while another part upgrades. When change is measured over intervals finer than the entire studied surface, such as 10 or 50 m elevation zones, there is a much higher incidence of measured change being greater than the uncertainty of measurement.
Horizontal uncertainty
If maps do not overlay precisely, measures of change between them will be affected. These problems have been rectified, thus eliminating most of this imprecision.
Tests were performed to determine the magnitude of imprecision caused by a slight displacement in the horizontal plane such as that which could exist between different map series. Each of the digitized data points representing the 1955 surface had their georeferencing changed by 10 m horizontally in each of the eight cardinal and diagonal directions. Surfaces were generated again for each displacement, and volumetric change calculations for the 1955–59 interval were repeated.
The standard deviation of the eight recalculations and the undisplaced value was 2.8 × 106 m3. This is less than a quarter of the value of volumetric uncertainty associated with the limits of photogrammetric vertical estimation for that interval.
Other sources of uncertainty
The estimates of vertical and horizontal uncertainly discussed above account for only two of the potential sources of error. Table 4 lists other error sources involved in this method and, where they were quantified, reports the values produced.
Conclusion
The large-scale terrestrial photogrammetric maps were made at 2 year intervals. When the uncertainty of vertical estimation was calculated, it was apparent that large parts of the surface of the glacier often had a magnitude of elevation change that was less than the uncertainty of vertical estimation. When similar calculations were made for longer intervals, the amount of such areas declined as the time interval between maps increased. In a time of recession, this is a result to be expected: the uncertainty remains constant, but the magnitude of change increases. A greater interval between maps would allow a more effective measurement of change.
When a series of maps is created to measure change, the interval between them should be determined by a comparison of the expected magnitude of elevation change and the precision of measurement of mapping. Despite the problem of high imprecision in the 1919 and 1955 maps, their comparison (table 2) reveals that very little of the change recorded between those two DEMs was less than the uncertainty of volumetric estimation. This is in contrast with many of the higher-precision DEMs which show substantial parts of the glacier having volumetric change that is less than the uncertainly of measurement.
This study demonstrates that since there are many sources of error in map-making and even more when maps are compared, extreme caution is needed when quantitative results are produced from the maps. In a study of this nature, it would be best to return to the original data sources and, using current software, produce DEMs directly from the photography which was used to make the maps.
Acknowledgements
The authors wish to acknowledge the financial support of Alberta Environment; the Atmospheric Environmental Service of Environment Canada; and U.S. Natural Sciences and Engineering Research Council Operating Grant No. OGP0038355.
