Introduction
The mass balance at a point of the glacier surface represents the direct climate influence on the glacier at that point. Determining the local mass balance and its distribution pattern on the glacier surface contributes to the monitoring of climate variations and to the understanding and modelling of glacier reaction to climate. Traditional methods of mass-balance determination using stakes and snow pits are time-consuming and expensive, even for a rudimentary network. Just deciding where the "representative" points should be located is in itself a major undertaking. Finally, it would be desirable to obtain mass balance for many glaciers in a geographic region, but this is rarely done because of the tremendous costs involved. Our aim is therefore to devise and test a remote-sensing method of direct mass-balance determination (cf. Reference Haeberli, Haeberli, Hoelzle and SuterHaeberli, 1998). In this work, we use photo-grammetrically derived velocity fields and a kinematic ice-flow model to determine mass balance over larger areas than is possible with traditional stake measurements.
Direct determination of mass balance by modelling the ice flow is not often done due to the need to specify area-wide boundary conditions at the glacier surface and bed. Reference RasmussenRasmussen (1988) obtained both bed topography and mass-balance distribution iteratively, using repeated aero-photogrammetric measurements of elevation and displace-ments on Columbia Glacier, Alaska, U.S.A. Reference Fastook, Brecher and HughesFastook and others (1995) derived bedrock elevations from repeated aerophotogrammetry. To help analyze an ice core drilled at Dyer Plateau, Antarctica, Reference Raymond, Weertman, Thompson, Mosley-Thompson, Peel and Mulvaney.Raymond and others (1996) calculated the mass balance of the vertical ice column at the drill site. For that purpose, not only the boundary conditions but also a profile of vertical velocity from borehole measurements was used. More often than with the above point-by-point approaches, the discharge fluxes for deriving the mass balance of glacier sections have been calculated. For valley glaciers, longitudinal sections with known average change in elevation and bounded by transverse profiles of velocity measurements are frequently used (e.g. Reynaud and others, 1986). Reference Bindschadler., Vornberger, Blaukenship, Scambos and Jacobel.Bindschadler and others (1996) derived the mass balance of ice-stream sections by means of area-wide velocity measurements using a similar approach.
Here we present a method which requires area-wide surface elevation, its changes with time and area-wide surface velocities for input. Glacier bed topography is known from radar measurements. Using these input data and some simple assumptions of ice mechanics, the three-dimensional ice-velocity field at the surface, and subsequently, the mass balance, can be calculated at every point of the glacier surface.
Owing to the size of many alpine glaciers, the resolution of satellite imagery is not adequate to accurately determine surface topography and its temporal three-dimensional variations. Analytical aerophotogrammetry, on the other hand, is a powerful and precise remote-sensing tool for monitoring three-dimensional geometrical changes of gla-cier surfaces. Nevertheless, its applicability is limited to areas with sufficient optical contrast. Thus, we applied the method described here not on accumulation areas but only on glacier tongues. The promising possibilities as well as the restrictions, the results of ice-flow and mass-balance computation, are shown for the tongue of Griesgletscheir, Valais, Swiss Alps, between autumn 1990 and autumn 1991.
Site Description
Mass-balance measurements on Griesgletscher (Fig. 1)started in 1961. The annual mass balance has been determined using the glaciological method based on a set of slakes. In addition, the photogrammetric method has been applied at intervals of some 5 years, with results that compare well with the glaciological method (Reference Siegenthaler, Kasser, Aellen and SiegenthalerSiegenthaler, 1986; Reference Aellen and FunkAellen and Funk, 1988; Reference Funk, Morelli and StahelFunk and others, 1997). Bind-schadler (1981) simulated the future evolution of Griesgletscher according to different mass-balance assumptions. Reference Hambrey, Milnes and SiegenthalerHambrey and others (1980) investigated the relation of ice structures and the strain regime. Reference Vieli, Funk and BlatterVieli and others (1997) modelled the three-dimensional ice flow and future behaviour of the whole glacier according to measured and estimated future mass balance. Since 1997, 33 stakes have been drilled over the whole glacier and surveyed several times a year to obtain the surface velocity field. High-resolution aerial photography from a height of about 2000 m above ground (image scale about 1:15 000) is repeated every autumn by the Swiss Federal Office of Cadastral Surveys, covering the entire glacier on one photo strip of about 13 images (Fig. 1 ).
The tongue of Griesgletscher between the icefall (left) and an artifical reservoir lake (right). Section of an aerial photo taken by the Swiss Federal Office of Cadastral Surveys, 10 September 1991. Griesgletscher is part of the Swiss glacier monitoring network and the object of numerous investigations. Inset: Location of Griesgletscher in the Swiss Alps.
Kinematic Boundary Condition at the Glacier Surface
Basic Equations
Our calculations of the local ice flow and mass balance follow the kinematic boundary condition at the surface. In this section we describe the basic equations and the special ice-mechanical assumptions we used for our computation.
From the mass conservation in a vertical column over the entire ice thickness, the local mass balance b at a point (x, y) can be expressed as the difference between the change in surface elevation within time, ∂z s/∂t, and the horizontal flux differences over ice thickness, ∂q x/∂x and ∂q v/∂y:
(Reference Paterson.Paterson 1994, p. 256-258). The horizontal flux is the integral of the horizontal velocity components, u and v, over ice thickness, and thus the flux differences can be written as
where zb is the glacier bed elevation, zs is the surface elevation, us and vs are the horizontal components of surface velocity, ub and vb are the horizontal components of basal velocity, ∂z s/∂x and ∂z s/∂y are the surface slope components and ∂z b/∂x and ∂z b/∂y are the bed-slope components. Inserting Equation (2) into Equation (1) gives
where ws is the vertical ice velocity at the surface
Equation (4) can be written as
Where ⋵xx and ⋵yy are the horizontal strain rates . Assuming incompressibility of ice (⋵xx + ⋵yy + ⋵zz = 0), Equation (5) can be expressed as
where ⋵zz is the vertical strain rate. Equations (3) and (6) represent the kinematic boundary condition at the glacier surface (Reference HutterHutter, 1983; Reference Paterson.Paterson, 1994). This geometric relation is valid at every point of the glacier surface. Some of the iec-flow components and the connection to geodetic stake measurements and photogrammetric measurements are depicted in Figure 2.
Photogrammetric and geodetic measurements and their relation to some components of ice flow which appear in the kinematic boundary condition at the surface.
Ice-mechanical assumptions
In our study, we try to determine all the components on the righthand side of Equations (3) and (6) and, subsequently, to calculate the local vertical ice velocity at the surface, ws, and the local mass balance, b. For that purpose, we use repeated standard photogrammetry to obtain surface slopes and surface elevation changes. A special photogrammetric technique is used to measure surface velocities. Bed topography and ice thickness are derived from ground-based radio-echo soundings. The basal velocity (sliding and sediment deformation) is estimated by comparing an existing ice-flow model for Griesgletscher (Reference Vieli, Funk and BlatterVieli and others, 1997) and measured surface velocities. We obtained a ratio of basal sliding to total surface velocity of about 75% (see below, section on calculated ice flow and mass balance).
To determine the vertical strain at the surface, ∈ zz s , it is necessary to know the variations of vertical strain rate with depth. These could be obtained by borehole measurements at single points (e.g. Reference Gudmundsson., Iken and FunkGudmundsson and others, 1997) or by ice-flow modelling. Reference Gudmundsson and BauderGudmundsson and Bauder (1999) give first results of a study on Unteraargletscher, Swiss Alps, which combines the approach using the kinematic boundary condition with a three-dimensional flow model (Reference GudmundssonGudmundsson, 1997). On Griesgletscher, there are no borehole measurements available, and Reference Vieli, Funk and BlatterVieli and others' (1997) flow model of Griesgletscher does not include sliding, which is assumed to have a marked influence on straining due to the high sliding ratio. Therefore, we decided to estimate the column-averaged vertical strain using the ice thickness, (h = zs-zb), the vertical strain rate at the surfaced ⋵zzs and a simple assumption for the variation of strain rate with depth:
In the case of no internal deformation and pure sliding, γ would be 1(i.e. vertical strain rate constant with depth). In the case of no sliding, the vertical strain rate at the glacier bed is 0 which, together with the assumption of linear variation of vertical strain rate with depth, gives γ a value of 0.5. As our estimate we use γ ≈ 0.75, considering the high sliding ratio for Griesgletscher tongue, and we introduce a large error estimate for γ in our calculations.
The vertical strain rates at the surface in Equation (7) can be derived from the horizontal ones with -⋵zzs = ⋵xxs + ⋵yys. The horizontal strain rates at the surface in turn were calculated from the surface velocity field, measured photogrammetrically in our study, using the expression
Surface Topography and Surface Kinematics
Photogrammetric methods
Aerophotogrammetric determination of digital terrain models (DTMs) and subsequent comparison of multitemporal DTMs is an effective and well-established technique for determining glacier surface elevation, surface slopes and the changes in surface elevation. The method uses monotemporal stereo models, composed of at least two overlapping photographs taken from different places (spatial baseline). The terrain point A (see Fig. 3) is computed by intersecting two spatial rays, each fixed by the known projection centres and the projections (A1'(t1)) and (A2'(t1)) of the selected terrain point. Repeating the procedure at other points gives a DTM at time t1. Repeating the procedure using photographs taken at time t2 gives the elevation of the same location (point B) and a DTM at time t2. From both DTMs, the changes in surface elevation throughout the area can be calculated.
Geometry used for determining surface elevation and its changes by monotemporal photogrammetric stereo models and for determining surface displacements by multitemporal stereo models. See text for explanation.
To determine surface displacements, a new photogrammetric method to measure surface displacements with high precision was developed and applied, based on the earlier concepts of Reference FlotronFlotron (1979), Reference ArmenakisArmenakis (1984) and Reference SaksonoSaksono (1984). The method was implemented on the computer-aided photogrammetric compilation system (analytical plotter) (Reference KääbKääb, 1996). It uses multitemporal stereo models composed from aerial photographs taken at different times [temporal baseline) and from different places (spatial baseline) (e.g. photo 1 (t2) and photo 2 (t1) in Fig. 3 ).
Between times t1 and t2 the surface point A has moved with the glacier to point C. A marker on the surface or features like closed crevasses or rough ice are suitable targets. The projection of such a point (A2'(t1) is chosen on the photo at time t1. Intersecting the spatial ray fixed by this image point and the known projection centre with the terrain surface, represented by the DTM at time t1, gives the three-dimensional ground coordinates of point A. This procedure is called monoplotting. The image point C1'(t2), which corresponds A' 2(t1), can be found using the stereoscopic overlap. The operator cancels the terrain movement between times t1 and t2 and t-2 by displacing one image while simultaneously looking at the multitemporal photographs. This simultaneous and stereoscopic observation supports the identification of corresponding points, improves the accuracy of the measurements and indicates whether a local displacement measurement is consistent with surrounding displacements. In the same way as for the spatial position of point A, the position of point C is computed by monoplotting. Thus, area-wide spatial displacements of the glacier surface can be deduced.
The glacier surface for which the velocity field is determined must show corresponding features in every photograph, and the displacements must be greater than the accuracy of the method. Accuracy increases with the time between the photographs, but the number of corresponding features may decrease. The photogrammetric technique described here works especially well for analyzing the creep of mountain permafrost (rock glaciers) and is also suitable for observing glacier flow and slope movements (e.g. Reference Gudmundsson., Iken and FunkGudmundsson and others, 1997; Reference Kääb, Haeberli and GudmundssonKääb and others, 1997,Reference Kääb, Gudmundsson and Hoelzle1998; cf. Reference Knizhnikov, Gelman, Osipova and TsvetkovKnizhnikov and others, 1998. for application of a similar technique ). It is very difficult to determine displacements on snow-covered terrain, where features are subtle and change rapidly. The accuracy of the method is assumed to be about ±30μm rms in image scale or, considering the photogrammetric parameters of the Griesgletscher missions, about ±0.4m rms on the glacier. The accuracy of the displacements was checked by repeating measurements using semi-independent multitemporal stereo models and by comparing the results with geodetic stake measurements on several glaciers and rock glaciers. The comparison for Griesgletscher indicated an accuracy (rms) of the displacement measurements of about ±0.6 m a-1 in the case of a 1 year interval between the two photo missions.
Surface elevation
For the photogrammetric measurements two stereo pairs of photographs, taken on 22 August 1990 and 10 September 1991, were used. The orientation was computed from artificial ground markers painted on stable bedrock beside the glacier for which coordinates were known from a geodetic network. DTMs at both times were composed of a regular grid of elevation points with a spacing of 50 m. Each DTM includes about 1000 points on the glacier surface. Contour lines were interpolated from these DTM points and the margin line of the glacier (also determined photogrammetrically) using the DTM interpolation software HIFI (height interpolation using finite elements; Reference Ebner and ReinhardtEbner and Reinhardt, 1984). Below the steep icefall, the surface slope of the tongue decreases to a nearly planar middle part at about 2600 m a.s.l., before increasing again towards the terminus (Fig. 4).
Surface contour lines of the Griesgletscher tongue (in m a.s.l.), derived from a 1990 DTM, and positions of the stakes for determination of mass balance and ice velocity. Coordinates (Swiss system) are in metres.
Changes in elevation for the period 1990-91 were computed as the difference between the two DTM grids. A Gaussian low-pass filter averaging nine weighted differences to obtain a new smoothed value at the centre point of the 3 x 3 cell filter mask was applied to reduce noise caused by terrain roughness and random measurement errors (Fig. 5). The accuracy of one height measurement is about ±0.4 m, so the elevation differences have a standard deviation of about ±0.6 m a -1 rms. Over the period 1990-91 the tongue thinned by up to 4ma -1 in the planar middle part (Fig. 5), a considerable thickness loss amounting to about 3% of total ice thickness. During the same year, most other glaciers in the Swiss Alps also experienced significant thinning due to the extremely warm summer. The raw elevation differences outside the glacier margin are at the error level, with some larger amounts caused by (moraine) erosion and melting of ice remains. In Figure 5, depicting the filtered data, the remaining non-zero elevation differences adjacent to the glacier are mostly an artifact of the low-pass filtering, which smooths sharp edges. This can be seen from the fact that the elevation differences decrease to small amounts at greater distances from the glacier. The filtering effect is especially evident at the lake, which despite being at different levels in 1990 and 1991, does not show a uniform change in elevation (Fig. 5).
Changes in surface elevation on the tongue, 20 August 1990 to 10 September 1991.
Surface displacements
The photogrammetric monoplotting technique for determining surface displacements, outlined earlier in this section, was used to derive surface velocities and, subsequently, to calculate surface strain rates and estimate surface strain and basal velocity. Two separate sets of displacements were measured for the period 1990-91, each consisting of about 700 points on the glacier surface. Two multitemporal stereo models were constructed for that purpose, each composed of different photographs of the photo strips of 1990 and 1991 and each newly oriented (i.e. one model of photo 1 (t2) and photo 2 (t1), and the other of photo 1 (t1) and photo 2 (t2); cf. Fig. 3). The semi-independent sets of displacements help reduce distortion effects in single photographs, errors in the orientation procedures and, most important, errors in single measurements. General error influences, such as bad ground-control points or camera errors, are detected to a lesser extent. Figure 6 shows the raw velocity field obtained from one of the displacement measurement sets. It gives an impression of the basic data quality, while the flow field itself is discussed below (Fig. 7). The velocity vectors in Figure 6 have a density similar to that of the DTMs, i.e. 50 m. Sufficient targets at both times were found in most regions, and measurement was fairly straightforward. In some regions missing targets, low surface contrast or large surface changes prevented determination of surface displacement. And at a few places, similar but not truly correspondent targets were identified, leading to obvious data blunder.
Raw surface displacement vectors measured by simultaneous monoplotting using a multitemporal photogrammetric stereo model.
Surface velocity field of the Griesgletscher tongue during 1990-91, derived from two measurement sets like Figure 6. The results can be compared with surface velocity from stake measurements ( inset).
To compare the two velocity fields, the scattered data were interpolated to a grid with a regular spacing of 50 m using a module of the IMSL mathematical software library (least-square fit to neighbours using quadratic polynomial; Reference AkimaAkima, 1978). Extrapolations greater than 50 m from original input data were deleted. Each vector in the first gridded dataset was compared with the corresponding vector in the second set, and a final value was computed as the average of both vectors. Velocities differing >15% in one or more of the two vector components (u,v) in both sets were eliminated, as were velocities measured only in one dataset. The resulting average velocity field was smoothed by the same Gaussian low-pass filter described earlier. After that filtering, data gaps of not more than one gridcell (50 m) were closed by interpolation. The velocity field resulting from this post-processing procedure (Fig. 7) was used for the following calculations. The horizontal velocity accuracy at a single point, deduced by comparison with stake measurements (inset in Fig. 7), was about ±0.6 m a-1 rms; the accuracy of the average velocity values was about ±0.4 m a-1 rms.
The highest annual surface velocities during 1990-91 on the Griesgletscher tongue occur in the icefall (up to 40 m a -1 ). The surface velocities decrease to 9ma -1 in the planar middle pari. The photogrammetric results reveal most ice in the ablation area coming from the southern (S) part of the accumulation area. Thus, the flowlines of the northwestern (NW) part of the tongue are shifted against the NW valley flanks. The large ice velocities directly at the NW margin suggests glacier sliding. Below 2550 m a.s.l., ice flow increases to >10ma -1 and then decreases again to 6—7 m a -1 near the terminus.
Surface strain rates
To model local ice flow and mass balance, we estimate the vertical strain from the vertical strain rate at the surface (cf. Equation (7)). For that purpose, we derive the horizontal strain rates from the surface velocity field using Equation (8). The resulting ⋵xy are sensitive to noise in the raw dis-placement measurements (cf. Fig. 6), so we apply the smoothed surface velocities (cf. Fig. 7). We use the surface velocity field described in the previous subsection to calculate horizontal strain rates following Reference NyeNye (1959), in which deformation at a gridpoint is computed from velocities at that point and its four nearest neighbours (see also Reference Bindschadler., Vornberger, Blaukenship, Scambos and Jacobel.Bindschadler and others, 1996). Figure 8 shows the horizontal strain rates at the glacier surface transformed to their principal axes. The average strain-rate accuracy, deduced by redundancy in Nye's algorithm, is about ±0.0002 a-1, only a few per cent of the calculated deformations. The vertical strain rate at the surface (Fig. 9) is derived from the incompressibility condition, (⋵xx + ⋵yy + ⋵zz = 0).
Horizontal principal surface strain rates during 1990-91 calculated from the velocity field. The crevasses of 1990 were mapped photogrammetrically.
Vertical surface strain rates during 1990-91 calculated from the horizontal strain rates using the incompressibility condition.
The increasing velocities in the upper part of the icefall cause marked extending flow. Strong horizontal compression, and hence vertical extension below the icefall, is due to decelerating flow. As expected, the horizontal strain-rate pattern at the margins (Fig. 8) shows a form and orientation typical for shear deformation. Also, below 2600 m a.s.l. the central glacier experiences very low strain rates as often as longitudinal compression. At 2600-2580 m a.s.l. the ice flow against the NW valley flank causes transverse horizontal compression due to converging flow and, additionally, longitudinal horizontal compression due to decreasing ice velocity. Both together result in strong compression and a rotation of the principal strain-rate axes of about 450 with respect to flow direction. Diverging flow and increasing velocity produce horizontal extension in the middle of the glacier tongue at about 2560-2580 m a.s.l. In most other parts of the tongue, horizontal compression and associated vertical extension predominate, as generally expected for ablation areas.
Figures 8 and 9 also include the 1990 crevasse pattern mapped photogrammetrically. Crevasses are expected where horizontal extension prevails. The orientation of the crevasses should theoretically be perpendicular to the direction of the largest horizontal extension (principal axis) (Reference VaughanVaughan, 1993; Reference Harper., Humphrey and PfefferHarper and others, 1998). This orientation is consistent for most mapped crevasses on Griesgletscher (Fig. 8; cf. Reference Hambrey, Milnes and SiegenthalerHambrey and others, 1980). Exceptions may occur, for example, where crevasses advect from regions of horizontal extension into regions of horizontal compression, without having yet totally closed. Small zones of horizontal extension, indicated by crevasse formation but without corresponding evidence in the strain-rate pattern, may have been missed in the strain-rate calculations, possibly due to inadequate resolution of the velocity grid spacing, or to faulty interpolation or excessive smoothing of the velocity-field. Finally, the observed pattern of crevasses is the result of some cumulative strain experience prior to 1991, and is not necessarily the result of the strain-rate pattern measured for the period 1990-91.
Determination of Glacier Bed Topography
To estimate the vertical ice velocity at the surface following Equation (6), bed slope and ice thickness have to be known. Both are derived from radio-echo soundings performed in spring 1988 and 1990 along II profiles using a monopulse-radar device (Reference Funk, Gudmundsson and HermannFunk and others, 1995). The method used to interpret the data is explained in Reference Fabri.Fabri (1991). The accuracy of the measurement method as deduced from comparisons of radar-derived ice thickness with deep ice drillings on several other glaciers is estimated to be about 5-10% of the ice thickness (Reference Haeberli and FischHaeberli and Fisch, 1984). Additional radio-echo soundings for Griesgletscher are planned, especially in the steep middle part of the glacier where no measurements have been possible so far.
A regular 50m gridded digital elevation model (DEM) and contour lines (Fig. 10) were interpolated from the ice-thickness profiles using the software HIFI (Reference Ebner and ReinhardtEbner and Reinhardt. 1984). The most notable features of the bed topography are the overdeepening between profiles R6 and R9 and the associated transverse bedrock riegel between profiles R9 and RIO. This riegel cannot be deduced from the radar measurements alone, and has been inferred from the bedrock topography outside of the glacier. The existence of the interpolated riegel seems to be confirmed by the comparison between modelled ice flow and stake measurements (see following section). Ice thickness of the tongue is up to 150 m in the over-deepening, and about 100 m over the riegel, both relative to the 1991 surface elevation.
Glacier bed topography (in m a.s.l.) of the Griesgletscher tongue interpolated from profiles of radio-echo soundings (R5-Rll) and bedrock topography of the surrounding ice-free area.
Calculated Ice Flow and Mass Balance
So far we have determined a DTM of the glacier surface, surface velocities and surface strain rates between 1990 and 1991, and a DEM of the glacier bed. To model vertical ice velocity at the surface (ws) following Equation (6), assumptions for the basal velocity and the relation between the measured surface strain rates and those at depth (cf. Equation (7)) are necessary.
Reference Vieli, Funk and BlatterVieli and others (1997) simulated the ice flow of Griesgletscher using surface and bed elevation data measured by photogrammetry and radio-echo soundings, respectively. Applying a flow model of Reference BlatterBlatter (1995), they calculated ice deformation neglecting basal sliding. Comparing the modelled ice deformation at the surface with the velocities measured photogrammetrically gives a ratio of basal sliding to total surface velocity of about 75%. At most places on the glacier tongue the sliding ratio differs by only a few per cent from this value. Deviations in the icefall and at the terminus are artifacts of the flow model's grid spacing (150 m) and of insufficient bed data. Thus, we decided to use the average sliding ratio of 75% over the whole tongue. The basal velocity is assumed to have the same direction as the surface velocity. For the above approach, we introduce an error of 10-20%, The vertical strain at the glacier surface was estimated from the ice thickness and the vertical strain rate at the glacier surface using Equation (7) and a γ of 0.75, and introducing a large error for this assumption (Table 1).
Approximate average, maximum and error of terms of the kinematic boundary condition at the glacier surface for Griesgletscher tongue (values in flow direction)
Using the measured photogrammetric and geophysical data, and the assumptions made above, it is possible to compute the vertical ice velocity at the surface and the local mass balance at every point from Equations (3), (6) and (7). Before describing the results, we discuss the accuracy of the calculations. To analyze the error propagation and the final accuracy of the calculated ice flow and mass balance, we estimate errors for the input data, the ice-mechanical assumptions and the final results (Table 1). To obtain the average accuracy of our model (rms), on the one hand, and to estimate a maximum error range, on the other hand, we take both average and maximum values for the input data, and analyze error propagation for both cases. The resulting maximum model errors are given below in parentheses. For the error analysis we neglect the sleep part of the icefall. Slopes and velocities of the estimates are transformed to the flow direction. The errors in bed slope and ice thickness are rough but maximum estimates. Consequently, the error of γ has an effect of up to ±0.2 m a -1 (max. 1.1) on the vertical ice velocity, while the error of basal velocity has an effect of only ±0.04 m a -1 (max. 0.9). These estimates show that in the special case of the Griesgletscher tongue, the rather rough assumptions for sliding velocity and total vertical strain have a minor effect on ws due to the low bed slope and low strain rates at the surface (cf. Equation (6)). All errors taken together cause an uncertainty in the calculated vertical ice velocity of about ±0.3 ma -1 (max. 1.5), which is dominated by the uncertainty in γ. Errors in elevation change, surface velocity and surface slope contribute an error in the calculated mass balance of ±0.5 m a -1 (max. 0.6). Thus, the total error of the mass balance modelled for a point is about ±0.7 m a -1 (max. 1.7).
The calculated vertical ice velocity shows some interesting features (Fig. 11). The large negative vertical velocity in the icefall is due to the steep bed slope and is therefore to be expected. The region along the NW margin where the glacier flows towards the NW valley flanks is a conspicuous feature, in which the model calculation yields positive vertical velocities of up to +3 ma -1. This may be due to the ver-tical strain (cf. Fig. 9). In addition, the flow direction in this area is rotated toward the margin, causing the along-flow bed slope to increase, and hence could result in positive vertical velocities and decreasing horizontal velocities, i.e. compressive flow. Positive vertical velocities at the surface have been calculated upstream of the transverse bedrock riegel, mainly as a result of the increasing bed slope at the riegel. In the remaining parts of the planar middle area the ice flows approximately horizontally, as expected because of the low strain rates and flat bed. At the steeper terminus, below a surface altitude of about 2550 m, it has been calculated that the ice flows downwards in two stages which are the result of bed-slope variations.
Calculated vertical ice velocity for 1990-91, including comparisons with stake measurements ( inset).
Some calculated values of vertical ice velocity can be compared with stake measurements (inset in Fig. 11). The differences show an agreement within the range of the error estimate of ±0.3 m a -1. At stake 7, above the transverse riegel, the stake measurement of the vertical velocity confirms the interpolation of this riegel.
Using the vertical surface velocity, the surface slope, the elevation changes and horizontal surface velocity, the distribution of the mass balance for 1990-91 can be calculated (Equations (3), (6) and (7)). The results, plotted on a 100m grid, are shown in Figure 12. The modelled mass-balance distribution mainly results from the distributions of vertical ice velocity and change in elevation. The glaciological mass-balance measurements at stakes 8-10 agree with calculated values within the error estimate of ±0.7 m a -1. At slake 7 no mass balance was available. At stake 6 in the lower part of the icefall, the resulting mass balance does not agree well with the stake measurement, even though the calculated vertical velocity appeared reasonable. Additional photogrammetric measurements were made in this area without detecting any gross errors. The horizontal velocity corresponds well with the stake measurement. Our favoured hypothesis is that this represents a gross error in the stake measurement; however, a difference between modelled and measured mass balance, in principle, may also he caused by the different periods of photogrammetry and stake measurement. The aerial photos used were taken on 22 August 1990 and 10 September I991, while the stakes were measured on 25 September 1990 and 24 August 1991. Different melting rates in the non-overlapping periods would influence the mass-balance measurement and therefore affect all comparisons between model and field data.
Calculated mass-balance pattern of the Griesgletscher tongue for 1990- 91, including comparisons with stake measurements (inset).
Despite the mass balance at stake 7, all differences between the calculations and the measured vertical ice velocity and mass balance are within the error estimates. Therewith, they are the expected result of the accuracy of the photogrammetric and geophysical methods used, the ice-mechanical assumptions and the error propagation through the kinematic boundary condition.
The modelled mass-balance distribution (Fig. 12) shows a high spatial variability. The values range up to nearly-6 m and vary by a factor of up to 2 at the same altitude. Even taking into account the accuracy of the applied method, these differences remain significant. Several influences may play a role, analysis of which will require further investigations. The obtained mass-balance pattern may be an expression of the solar radiation pattern, the snow distribution and a height gradient of air temperature. The spatial variations of the solar radiation, affecting the snow- and ice melt, depend on the topography of the glacier and of the surrounding terrain. Some low mass-balance values at the SE margin of the glacier may be attributed to the shading of the adjacent mountain ridge, while higher mass-balance values at the NW margin may reflect the absence of such shading. The ice mass-balance pattern is influenced by the snowmelt pattern, and therefore also by the distribution of snow depth over the glacier, which we do not know. The general decrease of mass-balance values to the SW is expected due to the height gradient of air temperature. Enlarged area of the ice surface, and thus high melting rates, may further influence mass balance at zones with crevasses or other high activity such as deformation (e.g. crevasses around slake 7, or strong deformation at the NW margin; cf. Figs 8 and 9).
Conclusions and Outlook
The kinematic boundary condition at the glacier surface provides a valuable tool for using area-wide photogrammetrically derived velocity fields and elevation data to model ice flow and mass balance. Ice-mechanical assumptions and data noise can have significant influence on calculated vertical ice velocity and mass-balance values. The effect of errors and assumptions will increase with increasing basal velocity, increasing bed slope and increasing surface slope (cf. Equations (3) and (6)).Therefore, the method is most suitable for planar glacier sections. Additionally, surface strain rates and surface strain can he expected to be low-under such conditions. Applying the model to ice flow over steeper glacier sections like icefalls will be difficult, but, on the other hand, may provide a useful alternative means of determining mass balance for such difficult-to-reach areas.
In addition to these model-inherent restrictions, there are limits imposed by photogrammetric techniques. Photogrammetric determination of elevation and displacement requires sufficient and persistent surface contrast and corresponding features, and is therefore restricted to ablation zones, or to crevassed or rock-covered areas in the accumulation zone. Contrast-giving features on the snow, such as sastrugi or dust, allow for elevation measurement, but often are not persistent enough for displacement measurement. Photogrammetry works much better on contrast-rich debris-covered ice than on the debris-free ice which predominates on Griesgletscher. Therefore, the method described here is especially suitable for determining the mass balance of debris-covered glaciers, which is difficult with other approaches.
To summarise, the modelling of vertical ice flow and mass-balance distribution using area-wide photogrammetric and geophysical data is most promising for ablation areas with low slope (cf. Bauder, 1996; Reference KääbKääb, 1996). Under the above restrictions, this approach has the potential to reduce the expenditure involved in traditional stake-based mass-balance programmes, for example by helping in siting representative stakes. Finally, this approach provides the spatial variability of the mass balance, which might contribute to improving three-dimensional models of glacier flow, energy-balance models and interpretation of mass-balance processes.
The future goal of modelling glacier mass balance using photogrammetric and geophysical data is to contribute to a remote-sensing strategy for monitoring mass-balance variations. For that purpose, it is necessary to apply repeated annual photogrammetry to obtain a time series of mass-balance distribution. In the case of such time series, the effect of some time-independent errors will be markedly reduced by analyzing only relative mass-balance variations. In fact, first results of a corresponding 20 year series (Reference KääbKääb, 1996) gave a reasonable mass-balance curve and suggest that even strong errors in glacier bed determination and in the ice-mechanical assumptions were mostly eliminated by calculating differences between annual mass balances.
Acknowledgements
We gratefully acknowledge the careful and constructive reviews of two anonymous referees. Thanks are due to H. Gudmundsson and W. Haeberli for many useful discussions and to numerous colleagues for their help with fieldwork. The study was part of the Swiss National Research Programme No. 31 on "Climate Change and Natural Disasters". The photogrammetric investigations would not have been possible without the valuable aerial photographs taken by the Swiss Federal Office of Cadastral Surveys.
MS received 18 September 1998 and accepted in revised from 26 March 1999
