Introduction
Small valley glaciers are believed to be sensitive indicators of climatic change, and many attempts have been made to establish the relation between glacier extent, volume and climate (e.g. Reference Johannesson, Raymond and WaddingtonJohannesson and others, 1989; Oerlcmans. 1994). Mountain glaciers and small ice caps are also believed to account for between one-third and one-half of observed rise in sea level, with the largest contribution From glaciers in the coastal mountains bordering the Gulf of Alaska (Reference MeierMeier, 1984, 1990; Reference Dyurgerov and MeierDyurgerov and Meier, 1997). However, these analyses were based on very limited data-sets, which, in many cases, entailed the assumption that the mass-balance record of one glacier is representative for a large region. Tests of this assumption arc limited (Rabus and Echelmeyer, in press), and it remains of interest to determine how glaciers of different geometries and types (tidewater vs land-terminating) in one region respond to a similar climate change.
A new method has recently been developed for measuring profiles of mountain glaciers relatively quickly and accurately, that of airborne surface-elevation profiling (Reference EchelmeyerEchelmeyer and others, 1996; Reference Sapiano, Harrison and EchelmeyerSapiano and others, 1998). As part of these studies, we profiled the upper region of the Harding icefield, Alaska, and 13 outlet glaciers emanating from the icefield. Six of the outlet glaciers were profiled in 1994, seven in 1996 and three in both 1994 and 1996. The profiles were compared to the U.S. Geological Survey (USGS) 1:63 360 topographic maps constructed from aerial photographs taken in 1950-52, Here we present their elevation and volume changes over this time interval. The changes of different types of glaciers, tidewater on the east side of the icefield and land-terminating glaciers on the west, all of which drain the same ice cap, arc compared to reveal differences in their response to the same large-scale climatic change. From the repeated profiles we are able to calculate the short-term elevation changes for three glaciers and the upper area of the icefield from 1994 to 1996. We also present terminus and volume changes for each glacier, as well as an estimated total volume change for the entire icefield. The repeatability of the profiling system and the quality of the maps are also discussed.
The Harding Icefield
The Harding Icefield is located on the Kenai Peninsula in south-central Alaska (Fig. 1). It is the largest icefield completely contained within the boundaries of the United States. The icefield is about 80 km long (northeast-southwest) and 50 km across. Including the outlet glaciers, it covers an area of about 1800km2. Slightly more than half of the icefield lies within the present boundary of Kenai Fjords National Park; the remainder lies within the Kenai National Wildlife Refuge. At least 38 glaciers of different sizes and types flow from the Harding Icefield. Seven of them are presently tidewater glaciers.
Present climate may be described as sub-Arctic maritime along the south and east sides of the Kenai Mountains, and sub-Arctic continental on the northwest. The east side is open to the ocean and receives copious precipitation in the form of maritime snow (Reference BensonBenson, 1980). This precipitation decreases towards the west, as can be seen in the long-term average annual precipitation records: Seward, on the east, receives 1.7 m, while Homer and Kenai, to the west, receive only 0.6 and 0.5 m, respectively. These latter two stations are in the precipitation shadow of the Kenai Mountains (S. Bowling, ftp: / /climate.gi. alaska.edu/publ ici Mo m hly.P). The icefield receives substantially more precipitation than Seward. For example, during one storm almost three times as much precipitation fell at 1275ma.s.l. on the Harding Icefield as in Seward (Reference RiceRice, 1987).
Fig. 1. (a) Harding Icefield, shouting the profiles flown in 1994 and 1996. Coordinates are in UΓ?, zone 5. The circles are the icefield crossing points between 1996 and 1994 profiles, and the diamonds are localions of snow-depth measurements and the snowline as measured in May 1996. (b) Boundaries of the different glaciers on the Hording icefield. The 1000 m contour on the main icefield and the 1200 m contour on the south icefield are shown as solid black lines. The numbers below the glacier names are area-averaged elevation changes.
The equilibrium-line altitude (ELA) of the Harding Icefield was estimated by Reference Meier and PostMeier and Post (1962) to be around 600 m, and the accumulation area ratio (AAR) to be 0,68.
The more recent tidewater glacier study of Viens (1095) indicated a higher ELA of about 930-1190 m on Northwestern and McCarty, and 610-730 m on Holgate and Aialik Glaciers, respectively. Each of these glaciers has a significant part of its accumulation area well above the stated ELAs, This indicates that Harding Icefield is relatively insensitive to small climate changes (Reference BodvarssonBodvarsson, 1955).
Late-Holocene glacier fluctuations- glacial chronologies and the vegetation chronosequence in the Kenai Fjords have been studied by Post (1980a,b,c), Wiles and Calkin (1990, 1994), Reference Helm and AllenHelm and Allen (1995) and Reference Wiles, Calkin and PostWiles and others (1995). Reference RiceRice (1987) compared photographs of the icefield to USGS maps based on 1950-51 aerial photography and found that the areal extent of the icefield decreased by about 5% over a period of 34 years. Greatest loss was near sea level along the Gulf of Alaska and at 300-6OOm elevation along the north and west sides of the icefield. Also, many small glaciers below 1000 m disappeared.
The characteristics of the 13 glaciers that we investigated are given in Table 1. Figure 1 shows the outline of the icefield and the boundaries that we defined for the glaciers, the ground tracks of the profiles, snow-depth measurement sites and crossing-point locations.
The Profile Data
The airborne profiling system uses dual-frequency kinematic global positioning system (GPS) methods to determine the absolute position of the aircraft, an infrared laser ranger to measure the distance from the aircraft to the point on the surface, and a vertical-axis gyro to measure the orientation of the laser beam. Post-processing of the data gives the ground track of the profile and the absolute elevation of the surface at points spaced approximately 1.5 m along that track. The nominal accuracy of the system is 0.3 m (Reference EchelmeyerEchelmeyer and others, 1996).
GPS data quality
The quality of kinematic GPS data depends on the number and geometry of the observed satellites. All the profiles obtained in 1991 had very good GPS solution because GPS satellites were observed during the Mights, with good satellite geometry. In 1996 poor weather conditions made it difficult to fly during limes of optimal satellite configuration, and often only 1 5 satellites were observed. This caused the data for several 1996 flights to be of poorer quality (Aoalgeirsdottir, 1997). However, as shown by the crossing-point analysis following, the vertical accuracy of the 1996 data is still within about 0.35 m. During good conditions the accuracy is 0.20m.These accuracies are better than the vertical accuracy of the maps, which is 15 m al best.
Crossing points
The repeatability of the system has been closely examined by analysis of crossing points. Where one profile crosses another, the elevations of the points from the two profiles that are closest to each other can be compared. An alternative method is to compare the elevations of all points that are within 1 m of each other from the two profiles. The latter method gives a statistically more robust dataset. We have applied both methods, with the results presented in Table 2. The data in this table are divided into two groups depending on whether the aircraft landed between measurements (different flights) or not. If the trussing was from different flights, the second profile was usually measured within 24 hours of the first one. Also shown are the weighted mean differences. In the latter case, each crossing was weighted by a factor equal to one over the combination of two independent errors: the horizontal distance between the two points, transformed to a vertical difference using the local surface slope, and the estimated error in elevation as determined by the height above the surface and the till of the beam relative to the normal of the glacier surface, Points from the same flight in 1994 are divided into those from outlet glaciers and those from the upper icefield.
Table. 1. Profiled glaciers, their characteristics, dates of profile and map photography
Comparison of all the points that are less than 1 m apart horizontally, without weighting, shows that the mean elevation difference was -0.17 m in 1994 and 0.34 m in 1996. The difference between 1994 and 1996 can be explained, in part, by the difference in the quality of the GPS data for the two years, as stated above. The crossing points in the same flight had a mean difference of-0.19 and 0.20 m for 1994 and 1996, respectively. When landing between measurements, the différence for the 1994 data was small, while in 1996 the mean difference for different flights was 0.42 m.
Table. 2. Crossing points statics representing the repeated of the sysytem
When only the closest points between the two passes are examined and they are weighted by the two errors, all mean differences are 0.21 m or less. Thus, under good satellite conditions (number and geometry), we believe that the system is capable of repeating a measurement to within about 0.20 m on moderately sloped terrain, while under poor conditions the accuracy is somewhat less.
Snow Depths and Seasonal Corrections
The comparison of our elevation data with those from the topographic maps is further complicated because the two datasets were collected at different times of the year. To correct for this différence, several snow depths were measured in 1996 (Fig. 1 and Table 3). Where the snow was not measured, we extrapolated using the altitudinal gradients defined by the Aialik data for tidewater glaciers and the Ski-lak data for land-terminating glaciers.
Table. 3. Snoiv depths measured in Mar 1996'. See Figure I for location of rrteosuretnentpoints
Another complication arises because the photos that were used to make die maps were themselves taken at different times of the year. Some glaciers were photographed in June and others in August, so the snow-depth correction needed to be adjusted for different glaciers, as discussed by Adalgeirsdottir (1997). We subtracted an appropriate amount of the measured snow thickness from the profile data to correct it to the date of the map photographs (June or August), assuming a linear decrease to zero in snow depth over the melt season (no ablation measurements were made). The photography dates for each glacier are listed in Table 1; the time interval for our thickness changes can be obtained by subtracting the year of the photos from tin-profile year.
Comparison with the Maps
The profile data are given with respect to different horizontal and vertical datums than the USGS maps. The profile coordinates are obtained in the World Geodetic System of 1984 (WGS84) and height above ellipsoid, and we transformed these to the map datums (North American Datum 1927 (NAD27) and National Geodetic Vertical Datum of 1929 (NGVD29), elevations relative to mean sea level (MSL) (Reference EchelmeyerEchelmeyer and others, 1996). The U.S. National Geodetic Survey model GE01D94 (Alaska) was used to transform the vertical coordinates. All horizontal coordinates were projected to Northing and Easting of the planar Universal Transverse Mercator system (UTM). The zone boundary between UTM zones 5 and 6 bisects the Harding Icefield, so the data were transformed to either zone as required.
The elevation difference between the maps and the profiles was obtained by first digitizing the glacier contours on the maps. Then the closest point where a profile crosses a contenir was selected and the elevation difference between the profile at that point and the contour was determined. Software developed by B. Rabus, J. Sapiano, J. Gorda and L. Bombardier (personal communication, 1994-96) was used for this purpose.
As described in a later section, an additional estimate of volume change was made by comparing the profiles to the Alaska digital elevation model (DEAD, which is given in the World Geodetic System 1972 (WGS72) and NGVD29. For this comparison we used horizontal profile coordinates in WGS84 because they are sufficiently close to those in WGS72 (about 3 m in Northing and 7 m in Easting on the Harding Icefield).
Map quality
The USGS maps have a contour interval of 100 ft (30.5 m). Their stated accuracy is half a contour interval in the vertical and about 50 m in the horizontal, both of which are almost two orders of magnitude less precise than the profile data. We have found that this map accuracy can vary from region in region on a given map, depending upon the quality of the aerial photographs used in die map construction, the quality of the ground control, the steepness of the terrain and other factors, To evaluate the map quality for each glacier, we obtained most of the photographs used to construct die maps, along with the cartographers’ reports for each map and an index map showing the coverage of each photograph (USGS, Rocky Mountain Mapping Center). From the index map we found that some areas were covered by more than one flight-line, flown in different years and at a different time of the year. Also, the cartographer did not always record which photos were actually used to make the map. In most of die ambiguous cases, we examined the photos and were able to identify those that were actually used by comparing the area extent of nunataks and bedrock in the photos and on the maps. However, this was not always possible. For example, in the accumulation area above Exit, Bear and Skilak glaciers, we found that photos from different dates were used in the mapping. The contours in these areas were mismatched because of real changes that occurred between the different photo dates, and these contours were smoothed by the cartographer. The dates listed in Table 1 are our best estimates for the actual mapping photographs used for each glacier.
On the higher areas of the icefield, where the surface is relatively flat. it is often difficult to see any contrast in die photos. This lack of definition made it unfeasible to map some of the snow-covered areas stereoscopically, and the cartographer stated that contours in these regions were simply “sketched” (USGS Quadrangle reports for Seward A-8 and Kenai A-l, and personal communication from J. Sadlik. 1996). In those areas with poor contrast the elevation changes from the maps to the différent profiles lacked consistency. In the upper regions of Skilak, Exit and Aialik glaciers there is a considerable difference between the shape of the profiled and mapped surfaces. We found that this lack of surface definition in the photographs was characteristic of the areas above the snowline and crevassed areas. Using the mapping photographs, we estimated the maximum elevation for good surface definition to be about 1000 m for the main icefield and 1200 m for the southern part. We then calculated the scatter in the elevation changes above these elevations. The mean change from the profiles to the maps was o.1m, with a standard deviation of 45.5m. This large scatter is an indication of the poorly defined contours, and we assume that the standard deviation is an indicator of the random error in the contours in this upper region. Because a large fraction (0.70) of the total icefield area is above 1000 or 1200 m (Fig. lb), the volume-change calculation is strongly affected by the errors there.
Proglacial bedrock points and nunataks
Profile elevations were compared to mapped bedrock contours in proglacial areas and to nunataks within the icefield boundaries in order to obtain another estimate of map accuracy, assuming that the bedrock has not changed. The profiles crossed 19 bedrock contours in front of a total of six glaciers. The results of this comparison are shown in Table 4, where the mean difference is the apparent elevation change averaged over all the points where the profiles cross bedrock contours. The standard deviation about this mean (often called the standard deviation) indicates the typical scatter, and, as such, it is a measure of the random error in the contours. The standard deviation of the mean (often called the standard error) shows that the mean difference is not significantly different from zero. Thus, based on this analysis, no systematic errors in the maps are indicated, and the random errors are about 12 m Of less in these proglacial regions, consistent with the published accuracy.
Table. 4. Profile minus map elevation in proglacial regions and on nunataks within the icefield
Fifteen nunataks were profiled in 1996. The results were similar to those in the proglacial areas, in that no systematic error in the maps was apparent. However, the larger scatter indicates that the random error in these nunatak contours is about 24 m. This may be a consequence of sleeper slopes on the nunataks, which would cause horizontal positioning errors on the maps to translate into larger vertical errors.
Table. 5. Elevation diffrence between profiles measured in 1994 and 1996
As a result of the 1964 earthquake (magnitude: Mw = 9.2), bedrock elevations around the icefield have not been strictly constant, as we assumed. There was coseismic subsidence of about lm in the Kenai Mountains (Reference Holdahl and SauberHoldahl and Sauber, 1994), followed by post-seismic rebound of about 0.2-0.6m (Reference Cohen and FreymuellerCohen and Freymueller, 1997). However, these tectonic elevation changes are too small to be resolved from the map comparisons.
Errors due to registration of intersection points
The horizontal coordinates ol the points where the profiles intersect the contour lines were determined using a digitizing table and suitable numerical algorithms (personal communications from B. Rabus, 1994-96). The horizontal errors in map registration and in the contour-profile intersections are both about 10 m according to Reference EchelmeyerEchelmeyer and others (1996) and Reference Sapiano, Harrison and EchelmeyerSapiano and others (1998). Repeat tests were done to confirm these error estimates for our data. We found that the combined error for registration and intersection-point selection was about 13 m, which is nearly the same as that quoted by those authors, assuming the errors are independent. The vertical error is the product of the horizontal error and the tangent of the surface slope; in steep areas it can be 10 m or more. A detailed test on one glacier indicated that the mean vertical error was 3 m for 27 intersections of varying slope.
As already noted, the Harding Icefield is located within two UTM zones. This makes our comparison with the maps more complicated and introduces additional errors because of larger map distortion near the zone boundaries. Tests indicate that there could be an additional horizontal error of about 30 m near the boundary, or 1-6 m in the vertical, depending on slope.
Summary of elevation-change errors
The elevation changes determined by comparison of the profiles with the maps are subject to each of the random errors mentioned and an additional error of a few meters in the snow-depth corrections. Combining all these errors, assuming that they are independent, gives an estimate for the total error in the seasonally corrected elevation change (ice thickness) at each contour intersection. For the upper regions, above 1000 m (generally above the snowline), this random error is 52m. In the lower regions (generally in the ablation areas) it is 21 m. There may possibly be unquantified systematic errors as well, especially in the upper regions.
Long-Term Elevation Changes
The glaciers profiled in this study span a range of glacier type, aspect, size and slope. There are three land-terminating glaciers, four that now terminate in lakes, four tidewater glaciers and two that have recently retreated from tidewater. Here we choose three representative glaciers for detailed discussion; later we compare the elevation changes for all the glaciers. (The detailed results for all the glaciers are further described by Adalgcirsdottir, 1997.) These three glaciers are Exit on the northeastern side of the icefield, Kachemak on the southwest and McCarty, a tidewater glacier on the south (Fig. lb). The profile and contour elevations are shown in Figures 2-4, along with the elevation change between them. The black diamonds show the mean value (0.1 m) of all elevation differences at profile-contour crossings above 1000 m on the main icefield and above 1200 m to the south, an approximation that we use in the volume-change calculations discussed below.
Exit Glacier is a small, relatively steep glacier that terminates on land on the northeastern side of the icefield, and it has a northeastern aspect (Fig. 1). This glacier has been the subject of study by Kenai Fjords National Park (Reference RiceRice, 1987; personal communication from M. Tetrau, 1996) because of its easy access. The terminus has retreated about 500 m since 1950/1951, and it thinned 80-90 m in the lower regions (Fig. 2). At the highest elevations, the surface appears to have changed shape since the maps were made, but the contrast in the mapping photographs was poor there and the contours may be poorly drawn.
Kachemak Glacier is a land-terminating glacier with western aspect on the southwestern side of the icefield (Fig. 1). Reference Bredthauer, Harrison and BredthauerBredthauer and Harrison (1984) measured thickening on the upper glacier and a retreat of the terminus between 1952 and 1979. The photographic contrast for its upper region was better than for most glaciers on the northern icefield, and consequently there is considerably less scatter in the elevation changes there. The glacier thinned by about 20 m except at the highest elevations, where the thinning was smaller, and at the terminus, where the thinning was large (100 m; Fig. 3). The terminus retreated about 900 m.
Fig. 2. Exit Glacier, (a) Profile (1994) and contour (1950/ 1951) elevations (MSL) along the profile, (b) Elevation change. 1950/1951 to 1994. The diamonds show the mean elevation change of ell icefield contours above 1000/1200 m used to calculate the volume change in one scenario.
McCarty Glacier is a tidewater glacier with southern aspect on the south side of the icefield. It retreated about 690 m from 1950 to 1996, but it has actually been advancing since 1960 (Reference PostPost, 1980b). Its elevation-change profile is different from that of most of the other glaciers, in that the higher elevations of the glacier thickened by 20-40 m, while it thinned at lower elevations (Fig. 4). At the present terminus the thinning was about 80-100 m. Where the glacier retreated, the elevation change was estimated from bathymétrie maps (Reference PostPost, 1980b).
Elevation changes on the entire icefield
The icefield was divided into four regions, and the elevation changes for the glaciers in each region were compared (Figs 5-8). All contour crossings above 1000/1200m are show n in Figure 9. Open symbols represent contour crossings where the profile was interpolated because of a lack of profile data caused by a lake, ocean or extremely rough surface. Where there was more than one profile on a glacier, the average elevation change at a contour is shown.
The western region of the icefield consists of land-terminating glaciers with western aspect, and they are all in the precipitation shadow of the icefield (Fig. 1). The elevation changes at the lowest elevations (Fig. 5) differ among the glaciers, as expected because their termini arc at different elevations and they retreated different amounts. At higher elevations, the two southernmost glaciers, Kachemak and Dinglestadt, have similar elevation changes, while Tustu mena and Chernof show larger changes. There is considerable scatter in the elevation change in the upper regions of any one glacier. For example, the ~1300 m contour on Tustu-mena Glacier was crossed three times during one flight, with elevation changes ranging from -120 m to +45 m (!), as shown m figures 5 and 6.
Fig. 3. Kachemak Glacier, (a) Profile (1996) and contour (1952) elevations (MSL) along, the profile, (b) Elevation change. 1952-96. Diamonds as in Figure 2.
Fig. 4. McCarty Glacier, (a) Profile (1996) and contour (1950) elevations (MSL) along the profile, (b) Elevation change. 1950-96. Diamonds us in Figure 2.
Fig. 5. Elevatian change for the west side of the icefield; open symbols represent estimated elevation change where the profile was interpolated.
Two land-terminating glaciers on the northern pan of the icefield, Skilak and Exit, show similar elevation changes at higher elevations (Fig. 6), but they have much smaller elevation changes than Tustumena Glacier to the south. Skilak Glacier retreated more than 3 km, while Exit retreated only about 500m. This difference can be seen in the elevation changes at lower elevations.
The eastern region of the icefield has two tidewater glaciers. Aialik and Holgate, and one, Bear, that was a tidewater glacier but now terminates in a lagoon. These glaciers receive more precipitation than those to the north and west. The two tidewater glaciers thinned less than Bear Glacier (Fig. 7), and the thinning on Bear was similar to that on Tustumena Glacier (Fig. 6). The elevation change for Aialik Glacier, which advanced slightly over this time period, was near zero, except at higher elevations. There was an apparent strong thinning at the head of Aialik and Bear Glaciers, but these data may be in error due to poorly defined contours there.
Fig. 6. Elevation change for the north side of the icefield; open symbols represent estimated elevation change where the profile was interpolated.
Fig. 7. Elevation change for the east side of the icefield.
Four glaciers wore profiled in the southern region. Two are tidewater. McCarty and Northwestern, and the other two, Northeastern and Little Dinglestadt, have retreated from tidewater and now terminate on land. Figure 8 shows the difference between McCarty and the other glaciers: it has thinned at lower elevations but thickened above 600m. Northeastern and Northwestern Glaciers thinned at higher elevations, but at lower elevations the difference in retreat cart be seen: Northwestern retreated 4.2 km, while Northeastern retreated only 1.3 km onto land. The open symbols shown for McCarty and Northwestern Glaciers are elevation changes estimated using the bathymétric maps of Post (1980b, c); this thinning is equal to the elevation of the contour pins the Fjord depth.
Fig. 8. Elevation change for the south side of the icefield; open symbols represent estimated elevation change (from Reference PostPost, 1980b, c) and where the profile was interpolated.
Fig. 9. Elevation change for the upper regions of the icefield.
Figure 9 shows the elevation changes along the upper, flatter regions of the icefield. All contour crossings above 1000 m elevation for the main icefield and above 1200 m for the southern area are shown. These elevations arc our estimate of where surface definition was lost on the mapping photographs. The mean elevation change for all these contour crossings was 0.1 m. with a standard deviation about the mean of 45.5 m. As mentioned earlier, this large scatter is a measure of the map accuracy in the upper regions; it is likely due to poorly drawn contours.
Transects across the icefield can be drawn using our profile data (Adalgeirsdottir, 1997). East-west transects indicate large (~100m) thinning on the western, lee side of the icefield, while the eastern, more maritime side has only a slight thinning. However, these differences may be due, in part, to different glacier sizes and surface slopes. North south transects show little difference in elevation change between 100 and 800 m elevation.
Short-Term Elevation Changes
The accuracy of the profiling system allows us to obtain elevation changes over shorter periods of time by repeat profiling. We have done this for Skilak and Exit Glaciers. They were profiled both in 1991 and in 1996, Skilak on the same calendar day, and Exit two calendar days later in 1996. A second profile was also flown on Holgate Glacier in 1996, but no attempt was made to repeat the 1994 ground track, so the two profiles do not overlap as much as they did on the other two glaciers. We have also compared the profile elevation at 11 points on the upper icefield where 1996 and 1994 profiles crossed (shown as circles in Figure la). All points that were within 1m of each other were compared on Skilak, Exit and Holgate glaciers, while only the two closest points (usually within 1 m) were compared on the upper icefield.
Elevation differences and their errors from 1994 to 1996 are listed in Table 5 these are the averages among the different crossing points. The average rates of elevation change over the 2 year period for Exit, Skilak and Holgate glaciers were- 1.8, -1.6 and -2.5 m a −1. respectively. For the upper icefield, the average rate of elevation change was -2.1 m a−1. The climatic significance of these short-term changes is discussed following.
Terminus Changes
Terminus changes from the 1950s to the 1990s are listed in Table 6. To identify the terminus in our profiles, we used the position of a “kink” in the profile (where the slope changes) or the time when the terminus was reported in the aircraft. All the land-terminating glaciers are retreating: Skilak and Dinglesladt retreated the most. All of the tidewater glaciers except Aialik and McCarty are also retreating. Bear and Northeastern Glaciers have retreated onto land.
Reference ViensViens (1995) found that most of the tidewater calving glaciers in Alaska have retreated over the last 200years, likely in response to a global rise in ELA. McCarty and Northwestern Glaciers started retreating in 1900 after a long, slow advance, while the terminus of Aialik Glacier has been in the same position since 1900 (Reference PostPost, 1980a, b, c; Reference Wiles, Calkin and PostWiles and others, 1995).
We determined recent variations in terminus position for four of the tidewater glaciers (Holgate, Holgate's neighbor, Aialik and McCarty) using aerial photographs (Adalgeirsdottir, 1997). These glaciers retreated between 1950 and 1978, advanced between 1978 and 1984 and. except for Holgate, continued to advance from 1984 to 1993. McCarty Glacier retreated 1400 m between 1950 and 1978. It has advanced 700 m since 1978, bin has not yet reached the 1950s position. The other three glaciers had smaller terminus variations of about 100-300 m.
The retreat rate is clearly related to water depth. North-western Glacier retreated at 36ma−1 from 1900 to 1927 in 20-100 m deep water on the terminal shoal, but aller it reached deeper water (150-300m) the rate increased to 460 ma −l (Reference PostPost, 1980c). McCarty retreated al a rate of 25 m a−1 in shallow water, then accelerated to 800m a−1 in deeper water (Reference PostPost, 1980b).
Table. 6. Changes of each glacier from 1950s to 1990s, along with glacier area at time of mapping. Also included (in italics) are regions not profiled where elevation changes were extrapola led from neighboring glaciers
Volume-Change Calculations
To calculate the volume change of individual glaciers and of the entire icefield, we must first determine the margins and ice divides for each glacier. This is a difficult task, in part because map errors in the upper icefield make it difficult to define the ice divides. Once the boundaries are defined, the large map errors on the upper icefield will be amplified in their effect on the volume-change calculations because a significant fraction of each glacier lies in this zone of large errors.
Reference EchelmeyerEchelmeyer and others (1996) and Reference Sapiano, Harrison and EchelmeyerSapiano and others (1998) describe a method for volume-change calculation in which new contour lines are constructed from the profiles and the glacier's area is allowed to change as it thins or thickens. The original map is compared to the newly constructed one lo determine the volume change. For Harding Icefield we used a simpler method because there are large areas of the icefield where no area change occurred and because there are large areas that were not profiled, where We felt their method could not be applied with any degree of accuracy. However, we did estimate the change in area near the terminus of each glacier. This area change was then subtracted from the original area to give an estimate of recent glacier area (table 6). We then assumed that the elevation change within a given elevation band was constant across a particular glacier. This elevation change was multiplied by the map area within that band, and these volume changes were summed to give the total volume change for a glacier. The area-average thickness change was calculated by dividing the total volume change by the average of the old and the new areas.
Fig. 10. Area-averaged elevation changes for each glacier calculated with the two scenarios for change on the upper icefield: actual contour crossing and with 0 m elevation change above 1000/1200 in (MSL). The dashed line is the area -average elevation change for the whole icefield.
Some areas listed in Table 6 (Lowell and Pederson Glaciers and the areas to the east and west of Skilak Glacier) were not actually profiled In these areas we estimated the volume change by using data from adjacent glaciers in combination with actual area distributions. This allowed us to calculate a more complete volume change for the entire icefield.
The upper icefield requires special consideration in the volume-change calculation because of the large scatter in the elevation changes there. We cannot quantify the errors involved. Instead we have estimated the sensitivity of the volume change to any such errors by computing the volume changes under two scenarios: (i) using the actual elevation changes determined from the profile-to-map comparisons (Fig. 9) with their large scatter, and (ii) assuming no elevation change above 1000 m on the main icefield and above 1200 m to the south. The average elevation change of 0.1m for all the contour crossings in this area supports the second scenario (shown as diamonds in Figures 2-4). The results of these two scenarios are given in Table 6, and the area-average elevation changes are shown in Figure 10. The means of these two area-average elevation changes are shown by each glacier in Figure lb.
It should be noted that area-average elevation changes in Table 6 and Figures 1b and 10 were calculated using our measured elevation changes from the date of the maps to that of the profiles. If Sorge's law applied to the snow and ice that was lost by ablation then the long-terra average mass balance (in water equivalent), (6), would be obtained by multiplying the area-average elevation change by a density of 900 kg m−3 and dividing by the appropriate time interval. However, this law does not strictly hold when old firn is ablated (Reference KrimmelKrimmel, 1989), and to correct for this we use a density of ~850kgm−3, following Reference Sapiano, Harrison and EchelmeyerSapiano and others (1998). This value of is given in the last column of Table 6.
Comparison of the results from the different scenarios shows that scenario (i), with the actual elevation changes, gives a more negative volume change on all but McCarty Glacier. In some cases the differences are small, either because the area at higher elevations is smaller or possibly because the contours are reasonably accurate in some of these regions. In other eases, such as Aialik, Bear and McCarty glaciers, the differences are quite large, implying either that the measured elevation changes in the upper areas were significantly different from the mean, or that these contours were incorrect. It is useful to note, however, that all the resulting volume changes are negative, except that of McCarty Glacier, and the magnitudes of the changes are significantly larger than the differences between the scenarios.
The volume changes of Northwestern and Northeastern Glaciers appear to be abnormally large. For Northwestern Glacier this can be explained by its retreat in a deep Fjord, where the thickness of the ice was great. The area distribution is different for Northeastern Glacier than for the other glaciers. It has a large fractional area at low elevations (200-500m), and the elevation changes at these elevations are large.
We estimated the total volume change for the entire icefield by summing the average of the volume changes of the two scenarios for each glacier, along with the estimated changes for these areas that were not profiled. The area-average elevation change for the entire icefield is then given by this total volume change divided by the average of the old and new total areas. The total volume change is about -34km3, and the area-average elevation change is -21m. Because of unknown errors in the maps, we cannot specify an actual error Tor this elevation change. However, we can obtain an idea of the error by looking at the difference in area-average elevation changes calculated using the two scenarios. The mean difference between the two scenarios was 5 m. We use this value to estimate an error of 5 m in the area-average elevation change. This estimated error is significantly smaller than the magnitude of the area-average elevation change, so we conclude that Harding Icefield has been losing volume during the ~43 year interval between the time of the map photography and our profiling.
Use of digital elevation model to calculate volume change
The area distribution for each glacier was determined using the digital elevation model (DEM) of Alaska, which gives elevation on a 90 m grid. It is important to note that the DEM was calculated from the original maps (USGS, 1990), so it has all the errors inherent in these maps. The boundary of each glacier was digitized from the topographic maps, and the DEM was masked using this boundary. The area distribution of each glacier was then determined by counting pixels within each elevation band.
An alternative method for determining the elevation changes using this DEM was developed by H. Li and C. Lingle (personal communication, 1997). Our profiles were directly compared to a 30 m interpolation of the DEM, giving the elevation change at each DEM pixel along the profile. This method could possibly lead to substantial savings in time and effort, as it would make it unnecessary to digitize the map contours. To evaluate the quality of the DEM relative to the maps, and of Li and Lingle's method, we compared our profile-to-map changes with those calculated using their method. In most cases, the local elevation changes compared quite well, but in a few cases there were large discrepancies. Figure 11 shows an example of the elevation changes calculated using the two methods. The comparison on Northwestern Glacier is quite good, while on Northeastern Glacier the DEM is more than 100 in different than the map for a significant length of the profile. The average difference between the two methods for all the glaciers (map minus DEM) was -4.4 + 1.0 m. This is a significant error, and therefore the DEM must be used with caution. In our analyses we used the DEM only for determining the area distribution for each glacier.
Fig. 11. Comparison of elevation changes relative to the topographic maps (circles) and relative to the DEM (line): (a) Northwestern Glacier, (b)Northeastem Glacier.
Discussion
The stability of an icefield is dependent on the distribution of area in the accumulation area with respect to the ELA. The ratio of the glacier area above 1000/1200 m (our estimate of the ELA in the 1950s; Reference ViensViens (1995) estimated lower values for some of the glaciers) to the entire icefield is 0.70. Most of this accumulation area is a few hundred meters above the elevation of the equilibrium line. A 200 m rise in ELA would decrease the AAR to 0.48. This indicates that the icefield is relatively stable to small shifts in the ELA (Reference BodvarssonBodvarsson, 1955), This ratio also defines the area where the elevation changes are ambiguous due to the lack of contrast in the photos.
Harding Icefield consists of a variety of glaciers flowing from a common accumulation area. The icefield is only about 80 km long and 50 km across, so we expect that the same large-scale synoptic) climatic variations would affect the entire icefield. However, there are differences in both temperature and precipitation between the more maritime eastern side and the more continental western side of the icefield. Our results allow us to illuminate any differences in the volume changes on these two parts of the icefield, as well as to address related questions, such as: How do different types of glaciers respond to a similar synoptic-scale climate change, and how do different characteristics of glaciers influence their behavior? To this end, we investigated the correlation between the area-averaged elevation change and different glacier parameters, with the following results:
Type of glaciers: We find no obvious difference in area-averaged elevation change between tidewater, lake and land-terminating glaciers. Most of the tidewater glaciers thinned by -16 + 7 m (excluding Northeastern and Northwestern Glaciers), while the land-terminating glaciers thinned by -17± 5m.
Location on icefield and aspect: The glaciers on t he south side of the icefield appear to have thinned more than those on the north, and, related to this, glaciers with southern aspect show more thinning than those with northern aspect, Although not statistically significant, these findings may support the suggestion of Mercer (1961) that the ELA has risen more on the south side of the icefield than on the north from 1909 to 1950. Surprisingly, there were no substantial differences between glaciers on the west (continental) and those on the east (maritime).
Area and length· There is no statistically significant correlation between elevation change and glacier area or length. However, the larger glaciers appear to have thinned slightly more than the smaller ones, probably because the longer ones have more area at lower elevations.
Surface slope: A positive correlation between average thickness change and surface slope might be expected because glaciers that have smaller slopes tend to have larger changes in AAR with a change in KLA than those with steeper slopes. We find that surface slope shows a slightly better correlation with elevation change than the other glacier parameters do. This is especially true if the two somewhat anomalous glaciers, Northeastern and Northwestern, are excluded from the regression analysis (see Fig. 12).
Terminus changes: We examined the correlation between the élévation change and the fractional change in length (△L/L). If Northeastern and Northwestern Glaciers are included it appears that the glaciers that thinned the most retreated the most, as expected. However, the correlation is strongly dependent on these two glaciers and, in any case, is poor (Fig. 12. A poor correlation between thinning and retreat was also found by Echelmèyer and others (1996) and Reference Sapiano, Harrison and EchelmeyerSapiano and others (1998). For non-tidewater glaciers, Haeberli and Hoelzle (1995) found a better correlation between fractional length change and average mass balance if a measure of glacier response time (Reference Johannesson, Raymond and WaddingtonJohannesson and others, 1989) is taken into account. However, we do not have the necessary information (specifically, the ice thickness) to test this relation.
Fig. 12. Regression analyses; (a) area-average elevation change vs slope; f(b) area-average elevation change vs fractional length change. Two best-fit lines are shown in each figure, the solid line is calculated without Northwestern and Northeastern Glaciers (squares) and the dashed line is with those two glaciers included.
In all, there were no statistically significant single-variable correlations between the area-averaged elevation change and glacier parameters. This is true even if only the more accurate map comparisons below 1000/1200 m elevation are used. It may be that the situation is more complex. For example, McCarty Glacier, the only glacier with a positive elevation change over the ~43year interval, is presently an advancing tidewater glacier with southern aspect on the south side of the icefield. In the same fjord, Little Dinglestadt Glacier, with an eastern aspect, retreated from tidewater and had an elevation change close to the mean for the entire icefield. Just cast of McCarty are Northeastern and Northwestern Glaciers, both of which have large negative elevation changes, Northeastern Glacier has retreated onto land, and Northwestern Glacier retreated 4 km but still terminates in tidewater.
Can one glacier be representative for the icefield?
The glaciers of Harding Icefield had a wide range of area-averaged elevation change, from near zero to about -90 m, with an average change of-23 ± 6 m. The results shown in Figure 10 indicate that the changes on eight of the glaciers were close to the area-average change of -21 m for the ice-field as a whole. Some of the glaciers are tidewater and others are Land-terminating. This indicates that one glacier could be chosen to be representative of the icefield, but its choice requires careful consideration. The scatter in the glacier changes about this icefield average is larger than that observed by Rabus and Echelmeyer (in press) for a glacier-tzed region of Arctic Alaska, indicating that the choice of a representative glacier for Harding Icefield is more difficult than for that region.
We can compare our results for the Harding Icefield glaciers with those for two nearby glaciers, both of which are Land-terminating and relatively small. Bear Lake Glacier, 20 km northeast of Exit Glacier, thinned 12.5 m over the last 37 years (Reference EchelmeyerEchelmeyer and others, 1996; Reference Sapiano, Harrison and EchelmeyerSapiano and others, 1998), or by about 14.5m extrapolated to our 43 year period. Wolverine Glacier, about 50 km northeast of Exit Glacier, had an average annual balance of about -0.25 m (water equivalent) from 1966 to 1995 (USGS, 1997). This extrapolates to about 13 in of ice thinning over our 43 year period. These elevation changes are somewhat smaller than the thinning of 21 m estimated for Harding Icefield. However, given the large uncertainties in both measurements, they are consistent. It is interesting to note that the elevation changes at high elevations for these two glaciers are somewhat more accurate than those on the Harding Icefield because the maps are better. Bear Lake Glacier shows thickening at these higher elevations (Reference Sapiano, Harrison and EchelmeyerSapiano and others, 1998); this leads to the smaller area-averaged thinning.
Short-term elevation changes
The short-term elevation changes between 1994 and 1996 show that the glaciers are presently thinning (table 5), and because the errors are small these elevation changes are well resolved. For Exit Glacier the area-average elevation change over these 2 years is nearly constant with elevation (i.e. area average equals local change). This is in contrast to the decrease with elevation shown by the long-term elevation change (e.g. Fig. 2). Our 1996 measurements on Skilak Glacier were limited to the elevation range 600-1100 m. Over that elevation range, the 2 year change was also nearly constant and, thus, we assumed that the change over the rest of the glacier was constant as well. Two-year elevation changes on the upper icefield decrease from -4.5 m in the south to -3.0 m in the north.
When this short-term elevation change is compared to the area-average elevation change for these glaciers over I lie last ~ 43 years (Table 6, column 6), we find that the recent thinning rate is an order of magnitude larger than that measured over the longer interval, and about three times larger than that determined for the icefield as a whole. However, tins increase m thinning rate may not have climatic significance. Comparison of the apparent increased rate of thinning on the Harding Icefield with the standard deviation of Wolverine Glacier's annual mass balance record (σ = 1.25 m a−1) implies that our measured elevation changes may only reflect the scatter in annual balance, being due to short-term fluctuations in snowfall and ablation. However, it is interesting to note that a similar accelerated thinning rate has been found on McGall Glacier in Arctic Alaska (Reference Rabus, Echelmeyer, Trabant and BensonRabus and others, 1995; Rabus and Krhel-mcycr, in press).
Conclusions
Airborne altimetry profiling is an efficient and accurate way to measure elevation changes on glaciers. We found the accuracy of the system to be about 0.2 m when GPS conditions are good, meaning sufficient satellite numbers and geometry. We determined elevation and volume changes for 13 glaciers and the upper icefield by comparing the profiles to topographic maps. Errors in these changes are subject to large inaccuracies in the maps at higher elevations; these are mainly due to a lack of contrast in the mapping photographs. The published accuracy of the maps (15m in our case) appears to apply to the lower elevations of the icefield, but above the snowline the errors are about three times larger (~50m).
Harding Icefield has been thinning and shrinking since the 1950s. We estimate that the total volume loss of the icefield is about 34 km3, which corresponds to an area-average elevation change of—21 m (ice equivalent) over the ~43 year time interval, or a long-term average annual mass balance of-0.4 m (water equivalent). Even though these long-term changes are subject to considerable errors associated with the maps (estimated to be about ±5 m in the area-averaged elevation change), they do provide an important measure of how the different glaciers have been changing, and some idea of the mass of water released by these glaciers into the oceans. This average mass balance is somewhat more negative than that generally found on mountain glaciers elsewhere in the Northern Hemisphere over the past few decades (e.g. Reference Haeberli, Hoelzleand and SuterHaeberli and others, 1996; Reference Dyurgerov and MeierDyurgerov and Meier, 1997). It is also somewhat more negative than that recorded on nearby Wolverine Glacier (USGS, 1997), a record that is often used as an index for this region. It is difficult to determine from the patterns of elevation change alone whether there has been an increase in temperature, or a decrease in precipitation, or both. Inspection of climatological data from the nearby town of Seward for 1908-95(http://climnate.gi.alaska.edu/history/SouthGentral/Seward .html) shows no clear trends, but these data are of questionable quality because the location of the observation site and the time of observation have changed over the period of observation.
The wastage has probably not been uniform in time, as the rate of surface-elevation change between 1994 and 1996 appears to be much larger than the average rate between the 1950s and the 1990s. However, this short-term value may not be climatically significant, given the annual variation in mass balance of nearby Wolverine Glacier.
There seems to be no significant relation between the type, aspect, size, slope or terminus changes and the volume change. It does appear that a given glacier can be at least cumulatively representative for this region, but it must be carefully chosen. A better approach is to examine changes of several glaciers using surface-elevation profiling.
Our measurements provide an accurate baseline against which future determinations of volume change can be made. New profiles can be flown along the 1994/96 ground tracks and elevation differences calculated using similar techniques.
Acknowledgements
We are grateful to J. Sapiano, B. Rabus, J. Gorda and L. Sombardier for their work on the profiling system and software, and to M. Tetrau for providing the digitized terminus positions. We also thank J. Sadlik. at the Rock) Mountain Mapping Center for providing the index maps and reports from the time the maps were made. J. G. Gogley,W. Haeber-li and C. Lingle provided helpful comments on an earlier version of the manuscript. This work was supported by NASA grant NAGW 3727.
