Abstract

This paper demonstrates some of the possibilities offered by GIS in arctic geomorphology. CIS is used for simulating the topography and volume of Pjetursson Glacier, Disko, Central West Greenland at different times in history. The volumes found shows good correspondance with volumes calculated from empiric formulas. An energy balance model is established in the frame of a rasterbased GIS. The model is used for demonstrating the significance of surrounding topography, on components in the energy balance.

Keywords

CIS, small arctic glacier, glacier area and volume, terrain analysis

Carsten A. Bøcker: Institute of Geography, University of Copenhagen,
Øster Voldgade 10, 1350 København K., Denmark.

Danish Journal of Geography 96: 11-20, 1996.

In the last 60 years, much work has been done to determine energy balances for melting glacier surfaces (see eg. Sverdrup 1935, Ahlmann 1953, Ambach 1974 & 85, Braithwaite 1981, Braithwaite & Olesen 1985, 88 & 90 and Wal et al. 1992), but not many with a GIS (Geographic Information System) approach. Escher-Wetter (1985) used a GIS-like data structure, calculating energy balances on Hintereisferner, Austria, by using a grid on the glacier surface. Humlum (1993) developed his own software package and used it for calculating the shortwave radiation balance of the Mitluaqkat Glacier, East Greenland. Later he used the model for calculating the energy balance on a reconstruction of the ice shield that covered Denmark during the last glaciation (Humlum & Houmark-Nielsen 1994). Etzelmüller et al. (1993) used GIS on Eriksbreen, northern Spitsbergen, in order to determine canges in mass balance, flowpattern and other parameters from 1970 to 1990 by comparing photogrammetrically derived models. They found that even if the glacier surface subsided on most of the glaciers area, crevasse pattern and surface slope did not change significantly.

The primary advantage of the GIS approach is the faciliation of performing calculations involving data with a fixed position in space, as data easily can be referenced through the data structure of the GIS. This advantage especially applies for calculations involving surface parameters. meters.In the case of a glacier surface these are e.g. slope, orientation and surface albedo.

Two types of GIS is generally available:

1) vectorbased GIS, where data are represented as lines, points and polygons. This type of GIS is very efficient for calculating distances, areas and all sorts of spatial relationships and 2) rasterbased GIS, where data are organized in a grid structure. Rasterbased GIS is efficient for spatial analysis involving several map layers especially when digital elevation models (DEM) are involved. Rasterbased GIS is also suitable for remote sensing data as the data structure for sattelite images are similar to that of a raster GIS.

In the present study two GIS was used; the vectorbased
ARC/INFO for digitizing map data and the rasterbased
IDRISI as a frame for the calculations.

Purpose

As many well tested methods exist for calculating the single elements of the energy balance for a glacier surface, the main attention of the present paper is focused on how to implement the methods in GIS and not on the methods themselves, and to demonstrate some of the possibilities of the GIS in volume calculations and terrain analysis.

The Study Area

The study area, around Pjeturssons Glacier, Disko, Central
West Greenland (69° 18'N, 53°25'W), is located app. 5 km
from Arctic Station, Godhavn (figure 1).

The study area consists of the glacier itself and its nearest surroundings. This means that the entire valley containing the glacier is included (figure 2). The study area is rectangular and spans an area of 12 km² with side lengths of 4000 m and 3000 meters orientated E-W and N- S. The lower left corner of the study area is defined as UTM zone 22 (402600 m, 7688250 m).

The valley which Pjeturssons Glacier occupies, is roughly orientated northwest-southeast and the glacier covered in 1962 app. 1.4 km². The lowest point on the glacier surface was in app. 450 m.a.s. and the highest in app. 900 m.a.s. Most of the glacier surface is oriented in northerly directions, almost 70 % of the glacier surface is oriented in the interval NW-N-NE. In addition parts of the

glacier are situated directly below a very steep northerly oriented rock wall, theoretically producing optimal conditions for maintaining a glacier on the northern hemisphere.

Methods

In order to calculate the energy balance, a model was developed, calculating the elements of the energy balance for a melting glacier surface, under the influence of shadowing and reflection from surrounding terrain, cloud cover, and atmospheric parameters as a sum of short- and longwave radiation balances as well as turbulent heat fluxes (Bøcker 1995). The model calculates the radiation or energy balance for any number of points in a DEM, using supplementary overlays containing slope, aspect, snow cover and surface type (albedo) as well as climate data files containing information of air temperature, humidity, precipitable water in the atmosphere, ozone layer thickness, cloud cover and wind speed. The model were programmed in Turbo Pascal as a stand-alone IDRISI application.

The extent and topography of Pjetursson Glacier and the surrounding terrain in 1962 is given by a topographic map (Blæsedalen ved Godhavn, Danish Geodetic Institute, Copenhagen 1967), which was digitized forming a base DEM.

The precision of the digitizing is estimated to be half a millimeter in the map. In the transformation of the map to UTM coordinates, a mean error was calculated to app. 4 meters in the model corresponding to 1/5 millimeter on the map, less than the digitizing precision. Further imprecision emerges as the vector data is forced into a raster grid and interpolated. The total error of the DEM after interpolation is estimated to be app. 20 meters or one pixel in the DEM.

The location of equilibrium line altitude (ELA) in 1962
was estimated from aerial photographs and digitized to-

gether with the outline of the glacier according to the topographic
map (fig. 4).

A DEM of the valley as it would appear without Pjetursson Glacier was produced (fig. 3). A series of cross profiles of the valley was drawn from the 1962 DEM, and the expected valley bottom extrapolated on these curves. This process was based on an assumption that the valley performed a high degree of symmetry, so that the glacier was expected to appear as a bulge on the south side of the valley. The north, unglacierized valley side could now be mirrored in a long axis defined as the middle of valley. Where this assumption could not be used, the cross profile of the valley bottom was based on intuition and personal experience of the appearance of former glaciarized valleys. The topographic map was updated with the new profiles and the base DEM was modified.

Another DEM, representing the appearance of Pjetursson Glacier at the beginning of this century was produced as well (fig. 4). During this period, it is expected that glaciers in the area reached their maximum extent during the Little Ice Age (LIA) (O. Humlum, pers. comm.).

The extent of the glacier was determined with a
photograph recordeded in 1906 and supported in the field
by identification of the LIA terminal moraine.

The 1906 ELA was identified in order to draw the isohypses on the 1906 glacier. This was done by following the height of the snowline as seen on the photograph. At the north edge, the ELA was determined as the height where the LIA terminal moraine vanishes. This height can be used as an indication of the ELA (Gross et al. 1977). When the glacier is expected to have a convex form in the ablation area and a concave form in the accumulation area, the isohypses can be drawn.

The extent and ELA of the 1993 glacier was determined in the field by using a handheld GPS-reciever and an altimeter, walking along the glacier where accesible and recording points with appropriate intervals. Where not accessible the extent was estimated from photographs recorded in the field.

The volume of Pjetusson Glacier was estimated for three years. This is done in principle by subtracting the empty valley DEM from the glaciarized DEM for each year and integrate ice thickness over total glacier area. However, due to digitizing imprecision and peculiarities of the interpolation method, it was necessary to compensate for errors along the edge of the glacier by using the outline of the glacier as the zero ice thickness isoline (fig. 5).

For comparison, the glacier volume was also calculated
using an empiric formula developed by Chen & Ohmura
(1990):

(1)

where V is glacier volume in 106 m3m3 and S is glacier
surface area in 106 m2.

Discussion

Fig. 4 clearly demonstrates how Pjetursson Glacier has retreated since LIA with a decrease in total glacier area of app. 40 %. It can be seen that the glacier primarily has retreated against south, hereby demonstrating the shielding effect of the northerly exponated rock wall. Despite the evidential frontal retreat from 1962 to 1993 the total

glacier area for the two years are almost the same, as the glacier in 1993 extents further up on the plateau north of the glacier than in 1962. The area of the 1993 glacier is disputable, as it is uncertain that parts of the accumulation area are permanent or temporary features.

The accumulation area ratio (AAR), the ratio between accumulation area and total glacier area, for a glacier in balance is app. 0.65 (Gross et al. 1977). However, if the glacier recieves a substantial part of its accumulation through e.g. avalanches or snow drift, the AAR for a glacier in balance can be lower, even down to 0.5 (O.Humlum pers. comm.).

The AAR in 1906 was 0.47, implying that the glacier had begun its retreat from its maximum position during the LIA. In 1962 the AAR was only 0.19, a clue to massive retreat. Since 1962 the AAR has increased to 0.41 for the total glacier area, signifying negative mass balance, but considering the two parts of the glacier individually, the AAR for the small part of the glacier is 0.19 and for the bigger part 0.50. The increase of the accumulation area on the plateau combined with an AAR of 0.50 could suggest that the glacier recieves material from snowdrift, thus implying a near neutral mass balance.

Humlum (1987) investigated 205 glaciers and firn areas on southwestern Disko and found that high situated glaciers oriented northerly and northwesterly had experienced much larger recession since the LIA than glaciers oriented otherwise. Humlum (1987) argues that they recieve a substantial part of their accumulation by snow drift from southwesterly direction. As these winds are prevailing in the 850 mb level, Humlum (1987) suggests that extensive backwasting of the glaciers can be contributed to a decrease in wind speed in the 850 mb level since LIA.

Humlums (1987) theory could explain the frontal retreat since LIA, as the orientation of the valley makes it a trap for snow transported by south-easterly winds. If it is true that the glacier today is in a near neutral mass balance, one could postulate a recent increase of wind speed in the 850 mb level. However, a confirmation of the postulate would demand careful studies of upper air meteorological data and a repetition of Humlum's investigation (Humlum 1987) with new aerial photos.

In 1906, maximum ice thickness was app. 120 meters in the uppermost part of the glacier with another local maximum of app. 120 meters on the glacier tounge, separated by a subglacial treshold (fig. 5). These values were bigger than expected, but as the volume calculations show, apparently parentlyreasonable. The subtraction of overlays yielded a volume of 1.23 10⁸ m3m³ while the calculation gave a volume of 1.11 10⁸ m³. This result is far better than expected and the app. 10% difference can easily be contributed to the build in imprecision in both methods.

In 1962, the glacier shows the same pattern of thickness
as in 1906 but with a maximum of app. 35 meters on the

glacier lounge and app. 100 meters in the upper part. The volume calculations gives good agreement between the two methods: subtraction of overlays gives a volume of 4.03 10⁷ m3m³ and calculation gives 4.52 10⁷ m³, again a difference of app. 10%.

Since 1962 the ice thickness has decreased substantially as the the maximum ice thickness in 1993 is app. 55 meters in the upper part of the glacier, according to figure 5. As the glacier has been divided into two, the volume were calculated individually for the two parts, yielding a total volume of 4.32 10⁷ m³. This has to be compared to the 1.7110⁷ m3m³ found by subtracting overlays. This difference of app. 60% could be explained by a glacier surface set too low in the DEM. However, an alternative explanation is sought in the fact that the accumulation area of the glacier has extended towards south and west on the plateau where it is rather thin, yielding a large surface area but a small volume. This could result in a situation where formula 1 is not valid. In addition, it was previously mentioned that the area of the glacier in 1993 was somehow overestimated. All three explanations is probably influencing the result: An overestimation of the surface area and an assumption of a large thin glacier area on the plateau will make the calculated volume become smaller, whereas an elevation of the glacier surface will increase the estimated volume, thus making the results of the two methods converge.

The method for determining the isohypses on the glacier surface for the years 1906 and 1993, is far from ideal. If arial photographs are available, elevation models can be obtained by photogrammetry, as done in Krimmel (1989) and Etzelmüller & Sollid (1995). However, the method has its limitations: the cartography must be done extremely carefully (Krimmel 1989), in case of newly fallen snow, a large glacier surface lacks contrasts and modelling in that area is not possible (Etzelmüller & Sollid 1995). Krimmel (1989) obtains good results for the mass balance over a 27 year period on South Cascade Glacier, Washington, but poor results for shorter intervals within the period and concludes that surveing is by far the best method, but off course not retrospectively possible. Krimmel (1989) uses a grid resolution of 100 m, compared to 5-10 metres in Etzelmüller & Sollid (1995) and 25 metres in the present study. This could introduce errors, especially along glacier edges. However, as computers gets faster and faster, the problem with a coarse grid resolution, should be minimized.

Topographic Influence on the Radiation Balance

Figure 6 shows the components of the energy balance, calculated by the model from hourly climate data for 0803, a cloudfree day with an almost constant wind (< 3 m/s) in the ablation area. The incoming shortwave radiation (S_in) initiates at app. 0500 hours where the sun appears on the local sky hemishere. Only the diffuse radiation is calculated as the sun is still hidden behind the local horizon. S_in increases at a constant rate until app. 0700 hours where the slope of the curve increases dramaticaly when the sun passes the local horizon and beam radiation is added to S_in. The same effect is seen at sunset where the sun passes the local horizon app. two hours before disappearing from the sky hemisphere. The influence of the horizon on S_in is reflected in both the outgoing (S_out) and net shortwave radiation (S*). S_out is symmetric with S_in as the surface albedo is chosen constant during the day. Several authors (e.g. Bailey et al. 1989) argues that the albedo is varying with solar altitude. A constant surface albedo therefore represents a type of compromise often involved in modelling. The shortwave components of the energy balance (S_in, S_out and S*) are all symmetric around app. 0200 hours, displaced two hours from midday due to surface geometry. Two components of the energy balance, the incoming

longwave radiation (Lin) and the sensible heat flux (SHF) are lineary proportional to air temperature, which results in parallel courses of the two. It is assumed that the glacier has a constant surface temperature of O°C through all of the ablation season. As a result the outgoing longwave radiation is a constant flux.

(Lout) Figure 7 shows the same calculation as in figure 6, only removing the influence of the surrounding terrain, as if the point was situated on an infinite plain. It is the shortwave components of the energy balance in figure that are of major interest as the turbulent heat fluxes and the longwave radiation components are independent of the surrounding

terrain. Comparison of S_in in figure 6 and 7 respectively, shows that S_in is bigger in figure 6 around solar zenith. This might be in contrast to what could be expected; that the topographic shielding would block the scattered diffuse radiation (I_s) and thus would be bigger on a perfectly flat surface as in figure. 7, than in the bottom of the valley as in figure. 6. Instead, the energy loss by the topographic blocking effect of diffuse radiation is more than compensated for by the effect of reflection of diffuse radiation by the surrounding terrain (I_r). This off course results in a corresponding increase in S* and Q*.

In table 1, the effects of topographic shielding and terrain reflection of diffuse radiation are separated and shown together with the beam radiation (I^^) as mean daily energy flux densities. It is seen that I^^, as expected, is biggest in figure 7 whereas the total diffuse radiation is biggest in figure 6. However, the most interesting feature in table 1 is that S_in for figure 6 and 7is almost equal. When integrating over the whole day, the difference in S_in is found to be only 888 KJ between figure 6 and 7. In other words: S_in is only 3.5% lower for a point down in a valley than for a point on a perfectly flat plain. However, the difference is a function of the sky view factor, i.e. how deep in the valley the point is situated.

In figure 8, S_in was calculated for each point in the entire test area in the same manner as in figure 6 and 7. This calculation illustrates the problems involved with choosing a suitable grid resolution. The calculation of sky view angles for each of the 30000 points in the DEM took app. six days on a very fast PC and demanded 90 million calls to the calculation routine. Figure 8 (A) is S_in calculated with the proper topographic features, (B) shows S_in calculated as if each point in the DEM was situated on an infinite, perfectly flat surface. (C) is the ratio (A)/(B). By using a ratio instead of a numerical difference, other effects affecting the magnitude of S_in is left out of account (e.g. surface albedo). (D) is the DEM of the study area for comparison.

The distribution of S_in on the surface (figure 8A) follows almost what could be expected. The map shows highest and lowest radiation on southerly and northerly oriented slopes respectively, whereas the plateaus receive quite constant radiation. An interesting feature is the beam of

high radiation at the valley mouth in the north-eastern corner of the map (framed in figure 8 A) This feature shows the importance of terrain shadowing and it is likely that the higher radiation is reflected in vegetation in this portion of the valley slope. In figure 8 B, all effects of topography and surface geometry are removed and S_in is therefore only dependent on terrain altitude. The only thing that separates Figure 8 B from Figure 8 D is the unit of the scale bar. The terrain altitude is not a very important parameter as S_in only changes app. 5 W/m² in app. 750 m. Even though very easily predicted, the results shown in map B are very relieving for a model developer, as it explicitly confirms the models internal validity.

Figure 8 C is characterized by a large white areas. This is where the surface recieves more shortvawe radiation under influence of topography than without this influence. Not only on southerly oriented slopes, but also on plateaus and in the valley. This shows how important it is to understand the influence of topography on incoming shortvawe radiation in order to perform reliable model calculations.

The surface geometry also affects the global radiation. Figure 9 shows the influence of (a) slope with a constant aspect of due south and (b) aspect with a constant slope of 30° respectively, calculated without influence of topography.Figure 9 (a) shows the decrease in day length with increasing slope and an increase in midday S_in with increasingslope until somewhere near 60° from where S_in decreases. The influence of aspect is reflected in a displacementof the period with shortwave radiation for the surfaces oriented against east, south and west respectively, with a higher midday S_in value for the surface oriented

against south. At the surface oriented against north, the period with direct sun is seperated in two; app. one hour from sunrise and from app. one hour before sunset. The above is straightforward, but how do the two effects inter act?

In order to illustrate the three-dimensional relationship between S_jn, slope and aspect, the daily mean S_in was calculated for 32 combinations of aspect and slope and a surface was interpolated from these values (figure 10). The generel outline of the surface looks as expected, except for the fact that, according to figure 10, a wall (slope = 90°) oriented due west or due east recieves a higher mean S* than a wall oriented due south. The explanation is that at noon the sun is higher on the sky than in the morning or in the afternoon. This results in a lower angle of incidence of the direct solar radiation on the wall facing south, than at a wall facing east or west. As the amount of incoming radiation is proportional to the angle of incidence, the south-facing wall will recieve the least shortwave radiation.

The surface form in figure 10 is a function of several location-bound parameters (e.g. latitude and surface albedo) as well as time of year, but it is believed that if the same calculations were performed elsewhere or at another time of year, the general outline of the surface would remain the same, only differing on the magnitude of S*.

Concluding Remarks

In the present paper some applications of GIS in the field of arctic geomorphology has been presented. It was shown that GIS is a very efficient tool for calculating volume changes, alltough some limitations were detected. This especially applies for errors along borders of glaciers, where the chosen resolution of the DEM is of great importance. Despite the limitations it is the firm belief of the present author that the future within mass and volume studies of glaciers is closely connected with GIS and computer modelling.

The topographic influence of the radiation balance was discussed and it was shown how important it is to understand the influence of topography on incoming shortvawe radiation in order to perform reliable model calculations.

References

Ahlmann, H.W. (1953): Glacier variationsand climatic fluctuations.
Bowman memorial lectures, series three. The American
Geographical Society.

Ambach, W. (1974): The influence of cloudiness on the net radiation
balance of a snow surface with high albedo. Journal of
glaciology, 13, 67: 73-84.

Ambach, W. (1985): Climatic shift of the equilibrium line: Kuhns
koncept applied to the Greenland Ice Sheet. Annals of
Glaciology, 6: 76-78.

Bailey, E.J. et al (1989): The radiation balance of alpine tundra,
Plateau Mountain, Alberta Canada. 21, 2: 126-134.

Braithwaite, R.J. (1981): On glacier energy balance, ablation and
air temperature. Journal of glaciology. 27, 97: 381-391.

Braithwaite, R.J. & Olesen, 0.8. (1985): Ice ablation in west
Greenland in relation to air temperature and global radiation.
Zeitschrift Gletscherkd. Glazialgeol. 20: 155-168.

Braithwaite, R.J. & Olesen, 0.8. (1988): Winter accumulation
reduces summer ablation on Nordbogletscher, S. Greenland.
Zeitsch. für Igletscherkunde und glacialgeologie. 24, 1: 21-30.

Braithwaite, R.J. & Olesen, 0.8. (1990): A simple energy balance
model to calculate ice ablation at the margin of the
Greenland Ice Sheet. Journ. of Glaciology, 36, 123: 222-228.

Bøcker, C. (1995): Simulating energy balance and ablation on a
small arctic glacier, using a GIS-approach. Unpublished
masters thesis. Institute of Geography, Univ. of Copenhagen.

Chen, J. & Ohmura, A. (1990): Estimation of alpine glacier
water ressources and their changes since the 1870s. lAHS
Publ. 193: 127-135.

Escher-Wetter, H. (1985): Energy balance calculation for the
ablation priod 1982 at Vernagtferner, Oetztal alps. Ann.
Glaciol. 6: 158-160.

Etzelmüller, 8., Vatne, G., Ødegård, R.S. & Sollid, J.L. (1993): Mass balance and changes of surface slope, crevasse and flow pattern of Eriksbreen, northern Spitsbergen: an application of a geographical information system (GIS). Polar Research, 12, 2: 131-146.

Etzelmüller, B & Sollid, J. L. (1995): Net mass balance of selected polythermal glaciers on Spitsbergen, Svalbard, measured by digital photogrammetry. Polar Research, in press.

Gross, G., Kerschner, H. & Patzelt, G. (1977): Melodische Untersuchungen die Schneegrenze im alpinen gletschergebiete. Zeitschrift für gletscherkunde und glacialgeologie, 12: 223-251.

Humlum, O. (1987): Glacier behaviour and the influence of
upper air conditions during the Little Ice Age in Disko,
Central West Greenland. Geografisk Tidskrift 87: 1-12.

Humlum, O. (1993): Towards an ablation model for the Midluagkat
Glacier, SE Greenland: Formulation of a shortvawe
radiation balance model. Submitted.

Krimmel, R.M. (1989): Mas balance and volume of South Cascade Glacier, Washington 1958-1985. In: Oerlemanns, J. (ed.) (1989): Glacier fluctuations and climatic change, pp. 193-206, Kluwer Academic Publishing.

Humlum, O. & Houmark-Nielsen, M. (1994): High deglaciation
rates in Denmark during the late Weickselian - Implications
for the paleoenvironment. Geografisk Tidsskrift 94: 26-37.

Sverdrup, H.U. (1935): The scientific results of the Norwegian- Swedish Spitzbergen expedition in 1934, part IV, The ablation on Isachsens plateau and on the Fourteenth of July Glacier in relation to radiation and meteorological condition. Geografiske Annaler 17:145-166.

Wal, R.S.Wvan der, Oerlemans, J. & Hage, Van der, C. (1992): A study of ablation variations on the tounge of Hintereisferner, Austrian Alps. Journal of glaciology 38, 130: 319-324.