Introduction

Medium- to long-term variations in the mass balance of glaciers and ice sheets are important indicators of climate change and provide a basis for the estimation of changes in continental water budgets and global sea level (Reference OerlemansOerlemans and others, 1998). Mass balance has traditionally been estimated using the glaciological method, which involves interpolating labour-intensive point measurements taken at stakes and snow pits across an entire glacier over time. It is therefore generally a time-consuming, costly and difficult record to derive (Reference Østrem and BrugmanØstrem and Brugman, 1991). Thus, despite their importance, detailed mass-balance measurements are carried out on only a limited number of glaciers worldwide, and longer mass-balance records ( >40 years) are available for just a handful of smaller glaciers confined to North America, Europe and the former Soviet Union. Measurements are sparse in remaining areas, especially the Andes and Himalaya and the high latitudes in general (Reference Haeberli, Bösch, Scherler, Østrem and WallénHaeberli and others, 1989).

Geodetic methods provide an alternative to the glaciological method and have been successfully employed for mass-balance monitoring of numerous glaciers (Reference HaakensenHaakensen, 1986; Reference KnudsenKnudsen, 1986; Reference Reinhardt and RentschReinhardt and Rentsch, 1986; Reference Krimmel and OerlemansKrimmel, 1989; Reference Etzelmüller, Vatne, Ødegård and SollidEtzemüller and others, 1993; Reference Willis, Arnold, Sharp, Bonvin, Hubbard, Lane, Chandler and RichardsWillis and others, 1996; Reference Sapiano, Harrison and EchelmeyerSapiano and others, 1998). The geodetic method involves comparing surveys of the glacier surface over a given period to calculate the change in volume. Glacier topographic data have traditionally been obtained from terrestrial surveys which suffer from all the limitations relevant to field monitoring of mass-balance stakes and pits (e.g. Reference Østrem and BrugmanØstrem and Brugman, 1991). However, numerous remote-sensing methods have emerged by which appropriately precise and accurate maps and digital elevation models (DEMs) can be obtained to carry out the necessary volume-change calculations. The geodetic technique does, however, suffer from two important drawbacks: (i) the elevation differences between successive surveys need to be corrected for the density of the surface materials, and (ii) it only yields the average mass balance of the entire glacier and not the spatial pattern of mass balance. Density corrections are straightforward below the equilibrium-line altitude (ELA) where the predominant process is ablation of pure glacier ice (density ∼830–917 kg m⁻³), but are a problem above the ELA because the near-surface density is not constant over the whole area and may vary between 50 kg m⁻³ for fresh powder and 830 kg m⁻³ for very wet snow and firn (Reference PatersonPaterson, 1994). However, the density error is minimal when the ELA is at about the same position in the two DEMs. Under these conditions the density profile, even in the accumulation area, is remarkably stable from year to year, especially in the late summer. Surface changes can be thought of not as peeling off a layer at the surface firn but as slicing off at the base of the ice column, so the density correction will be spatially consistent (personal communication from J. Kohler, 1999). Even under these somewhat relaxed constraints, though, careful consideration of the timing of surface surveys coupled with ground control is a prerequisite for confident application of the geodetic method.

The second drawback is of a more fundamental nature: point mass balance cannot be determined from geodetic methods alone since spatial patterns of glacier elevation change are also influenced by the pattern of mass-flux convergence and divergence which results in vertical strain at the glacier surface. However, the geodetic method may be used to test mass-balance estimates derived by the glaciological method since the net change in mass over time relates to the net volume change (Reference Østrem and BrugmanØstrem and Brugman, 1991). For example, in Norway high-quality, 1:10 000 scale topographic maps have been successfully used to both verify and challenge traditional glacier mass-balance measurements over a number of years (Reference HaakensenHaakensen, 1986; Reference Østrem and HaakensenØstrem and Haakensen, 1999). However, the net mass balance really only constitutes an indirect climatic signal; it is the spatial distribution of mass balance that directly reflects the local variables of temperature, precipitation, radiation, etc., and is therefore the better climatological indicator.

More recently, attempts have been made to estimate a glacier’s mass-balance distribution through combining the remote-sensing techniques employed in the geodetic method with information on the three-dimensional flow field of a glacier. Significantly, these data can now be obtained at sufficient precision and spatial resolution, and numerous studies have implemented either the mass-continuity or kinematic boundary condition at the glacier surface (Reference HutterHutter, 1983) to determine surface mass balance (e.g. Reference RasmussenRasmussen, 1988; Reference Hulbe and WhillansHulbe and Whillans, 1994; Reference Bindschadler, Vornberger, Blankenship, Scambos and JacobelBindschadler and others, 1996; Reference Hamilton, Whillans and MorganHamilton and others, 1998; Reference Gudmundsson and BauderGudmundsson and Bauder, 1999; Reference Kääb and FunkKääb and Funk, 1999). Reference RasmussenRasmussen (1988) used an iterative procedure to successfully derive both bed topography and the mass-balance distribution using measurements of surface elevation and displacement of Columbia Glacier, Alaska, derived from repeat photogrammetry. Reference Bindschadler, Vornberger, Blankenship, Scambos and JacobelBindschadler and others (1996) successfully used surface-elevation change and area-wide velocity measurements to determine mass balance at Ice Streams D and E, Antarctica. Reference Hamilton, Whillans and MorganHamilton and others (1998) measured local vertical strain rates and accumulation rates and used the method proposed by Reference Hulbe and WhillansHulbe and Whillans (1994) to give ice-sheet thickening or thinning rates across Antarctica. They determined the vertical component of velocity some 10–20 m beneath the firn surface using fixed anchors in conjunction with global positioning system survey (the “coffee-can” method), and they determined the accumulation rate from the depth of measured nuclear-bomb fallout detected from gross beta activity. Correcting these terms for firn settling and downslope motion, Reference Hamilton, Whillans and MorganHamilton and others (1998) inverted the kinematic boundary equation to yield the change in ice-sheet thickness.

At valley-glacier scale, Reference Kääb and FunkKääb and Funk (1999) used the kinematic boundary condition and high-density remotely sensed measurements of area-wide surface elevation, its change with time and surface velocity to determine the mass-balance distribution of Griesgletscher, Switzerland. By invoking a simple parameterization for the variation of strain rate with depth assuming a fixed 75% ratio of basal sliding to total surface velocity, they met with some success (±0.7 m a⁻¹) in calculating the mass-balance distribution, although their results do not agree quite so well with actual observations. Similarly, Reference Gudmundsson and BauderGudmundsson and Bauder (1999), using an identical method but applied to Unteraargletscher, also found calculated values of mass balance lay within a “reasonable range” of actual measurements. However, systematic analysis of each of the terms in the kinematic equation revealed that the simple depth-related vertical strain parameterization employed was ultimately at fault, and Reference Gudmundsson and BauderGudmundsson and Bauder (1999) conclude with a call for three-dimensional flow modelling to give an improved estimate of the vertical strain field. The intention in this paper is to do just that; to combine photogrammetric methods with three-dimensional ice-flow modelling in an attempt to derive high-resolution spatial patterns of net mass balance using Haut Glacier d’Arolla, Valais, Switzerland, as a case-study.

Approach

The method proposed here is slightly different to that used in the two-valley glacier studies mentioned above, insofar as we invoke the mass-continuity equation (Reference PatersonPaterson, 1994) and not the surface kinematic equation to indirectly determine the mass-balance distribution b(x, y). This difference reflects the fact that we are using three-dimensional flow modelling to compute the flux-divergence field, that is the net expression of all flow components, rather than just the vertical strain at the surface. It also reflects our use of surface velocity measurements to constrain the flow-modelling effort rather than provide a direct input field, which is required by the kinematic equation. The advantage of our method is that it requires considerably fewer data points at a lower spatial density. This permits the mass-balance calculation to be made across areas where high-density surface velocity measurements are not available. The disadvantage of this approach is that it does rather crucially depend on successful modelling of the three-dimensional flow regime, a difficult task given the need to specify area-wide boundary conditions at the glacier surface and bed. However, since the kinematic method also crucially hinges on indirect determination of the three dimensional flow for the vertical strain field (Reference Gudmundsson and BauderGudmundsson and Bauder, 1999), then given the datasets available and the efficacy of the previous flow-modelling effort, we pitch for what we consider the lesser of two evils.

Within a Cartesian co-ordinate system x, y, z, with z the vertical ordinate, positive upwards, the rearranged equation for mass continuity (Reference PatersonPaterson, 1994) gives:

(1)

where t is time, h is the ice elevation of the ice surface, H(x, y) is ice thickness, is the vertically averaged velocity vector and ∇ is the horizontal difference operator.

To calculate the distribution of mass balance b(x, y), the lefthand terms in Equation (1) require estimation. The first term, ∂h/∂t, can be readily calculated using the geodetic methods outlined above, that is differencing two DEMs of the glacier surface for successive years. The next term, the local mass-flux convergence or divergence, is derived using numerical three-dimensional ice-flow modelling optimized with actual surface velocity measurements.

To date, three ice-flow models have been documented, based on both finite-element and finite-difference numerical techniques which are capable of computing a glacier’s three-dimensional flow field at sufficient resolution (Reference BlatterBlatter, 1995; Reference HansonHanson, 1995; Reference GudmundssonGudmundsson, 1999). On the basis of previous success at modelling the dynamics of Haut Glacier d’Arolla (Reference Hubbard, Blatter, Nienow, Mair and HubbardHubbard and others, 1998), here we use the first-order finite-difference algorithm of Reference BlatterBlatter (1995) to solve the three-dimensional stress and strain field. Boundary-condition requirements for application of the first-order model are:

1. glacier surface DEM

2. subglacial bedrock DEM

3. surface velocity measurements.

In addition to the glacier geometry, all of the above models require some degree of optimization of their respective flow parameters (rate factor, flow exponent and sliding parameters, where appropriate) by recourse to actual strain measurements. Although, in ideal circumstances, full three dimensional strain measurements would provide the perfect constraint on the basal and surface boundary conditions, such extensive data are an impractical expectation, and measurements of surface velocity alone must suffice if this method is to have any real applicability. However, down to a limit generally agreed to be of the order of one ice thickness (Reference Balise and RaymondBalise and Raymond, 1985; Reference MacAyealMacAyeal, 1992; Reference Bahr, Pfeffer and MeierBahr and others, 1994), the higher the resolution of such velocity data available, the greater the degree of model optimization, and hence the greater confidence that may be placed in the resulting flow field.

Objectives

Specifically, our aims in this study at Haut Glacier d’Arolla are to:

Study Area

Haut Glacier d’Arolla covers an area of 6.3 km² and is situated at the head of the Val d’Hérens, a southern tributary of the Rhône valley (Fig. 1). The glacier has a northerly aspect, is about 4 km long and comprises two upper basins which lie at about 3000 m a.s.l. and which combine to feed a glacier tongue extending down to about 2560 m a.s.l. Slopes over much of the glacier are gentle, except above the upper southeast basin where the glacier is fed by a series of steep icefalls rising up to 3500 m a.s.l. on the north face of Mont Brûlé. The glacier is surrounded by the Bouquetins ridge to the east and the mountains of L’Évêque and Mont Collon to the west. Over the last decade the ELA has been well above 2800 m a.s.l. (Reference OerlemansOerlemans and others, 1998), leading to a strong negative mass balance with about 2.5–3.0 m a⁻¹ of surface ablation across the lower tongue and > 100 m of retreat since 1989.

Fig. 1.

Haut Glacier d’Arolla, showing locations of surface-velocity markers and the mass-balance network along the glacier centre line. The glacier surface is contoured at 100 m intervals.

Methods and results

Glacier surface-elevation change

DEMs of the glacier surface were constructed from analytical photogrammetry of vertical photographs taken about 1000 m above the glacier on 17 September 1992 and 20 September 1993. The methods employed in the production of the DEMs are standard and are described in detail by Reference Willis, Arnold, Sharp, Bonvin, Hubbard, Lane, Chandler and RichardsWillis and others (1996). In brief, a least-squares adjustment procedure was used to fix the DEMs into the Swiss Grid Co-ordinate System by optimizing 15 on-glacier and 8 off-glacier markers to three Swiss Survey Trigonometric Points. Some parts of the upper glacier are blanked out due to difficulties in obtaining high-precision elevation data from snow-covered textureless surfaces. In these zones, lack of sufficient visual contrast in the aerial photos prohibits the ray-intersection procedures on which analytical photogrammetry is based. Hence, the resulting 20 m resolution DEMs are largely for firn- and ice-covered areas of the glacier, with vertical errors typically less than a few centimetres over much of the tongue but up to ± 10 cm over parts of the upper glacier. Horizontal errors associated with the digitization of markers and subsequent least-squares adjustment procedure are about ±35 cm.

Subtracting the 1992 DEM from the 1993 DEM gives the spatial pattern of surface-elevation change, where a decrease in surface elevation is represented by negative values, and surface-elevation gain by positive values (Fig. 2). Thinning occurred over most of the glacier, with maximum values of about 4.5 m on the lower eastern glacier margin. Some areas of the glacier experienced thickening of up to 0.5 m, especially the lower west margin where the shading effects of Mont Collon combined with a thick debris mantle appear to have severely limited ablation compared to the rest of the lower tongue. However, net mass wastage is the predominant trend, with an average surface elevation loss of 0.86 m, equating to a net volume change of −5.34 × 10⁶m³.

Fig. 2.

The pattern of surface-elevation change, September 1992–September 1993, calculated from the difference between the 1992 and 1993 DEMs. The pattern is contoured at 0.0 m, with positive values (dark shading) indicating glacier thickening and negative values (light shading) indicating glacier thinning

Modelling the annual flux-divergence field

The first-order flow model of Reference BlatterBlatter (1995) is used to compute the spatial pattern of annual flux divergence across the glacier. This three-dimensional model has been applied successfully to Haut Glacier d’Arolla at a horizontal grid resolution of 70 m with 40 scaled vertical layers, verified against diagnostic field data and used successfully to elucidate the annual flow regime between 1994 and 1995 (Reference Hubbard, Blatter, Nienow, Mair and HubbardHubbard and others, 1998). On the explicit supposition (which will be tested later) that the annual flow regime at Haut Glacier d’Arolla has remained essentially unchanged year to year from 1992 to 1995, we use the identical 1994/95 boundary conditions to derive the 1992/93 flow field, and the reader is referred to the above study for specific details.

To replicate the annual flow regime, Reference Hubbard, Blatter, Nienow, Mair and HubbardHubbard and others (1998) conduct three independent numerical modelling experiments, each constrained by a contrasting basal velocity boundary scenario defined by measured surface displacements made at 34 surface markers (Fig. 1) between 1994 and 1995:

winter base flow: assuming zero basal motion and providing a means by which the value of the flow-law rate factor A of 0.063 a⁻¹ bar⁻³ could be determined by tuning to winter surface velocities observed between 3 February and 8 February 1995;

summer flow: a basal motion distribution representative of mean melt-season conditions defined by the residual of summer surface velocities (observed 30 June–3 September 1994) and the above observed winter surface velocities;

enhanced spring flow: an enhanced distribution of basal motion representing conditions at the bed during a major glacier-wide speed-up event which lasted 10 days defined by the residual of the observed spring-event surface velocities (measured 19–29 May 1994) and observed winter surface velocities.

Each of the above flow scenarios represents a characteristic snapshot of the distinctive flow regimes operating over the course of a year at Haut Glacier d’Arolla. Thus, although diurnal variations in the morphology and hydraulics of the subglacial drainage system affect ice dynamics (Reference Hubbard, Sharp, Willis, Nielsen and SmartHubbard and others, 1995; Reference Gordon, Sharp, Hubbard, Smart, Ketterling and WillisGordon and others, 1998; Reference Nienow, Sharp and WillisNienow and others, 1998), these diurnal patterns are aggregated into the three flow scenarios, which in turn are combined to yield the annual flow regime. Reference Hubbard, Blatter, Nienow, Mair and HubbardHubbard and others (1998) combine the flow fields associated with each snapshot using a time-weighted average consisting of 28 weeks’ winter base flow, 22 weeks’ summer flow and 2 weeks’ enhanced spring flow to give the overall 1994–95 annual flow regime. This composite pattern compares well with observations at Haut Glacier d’Arolla between 1994 and 1995, accurately reproducing both the observed pattern of surface velocity and the subsurface velocity distribution in cross-section derived from repeat, down-borehole inclinometry measured over the same period (Reference Harbor, Sharp, Copland, Hubbard, Nienow and MairHarbor and others, 1997; Reference Hubbard, Blatter, Nienow, Mair and HubbardHubbard and others, 1998).

What is significant is the critical effect that the basal boundary condition has on the modelled flux-divergence field. A time-weighted annual flow composite derived from representative seasonal snapshots with a significant basal motion contribution yields a considerably more sensible three-dimensional flow field than if the model were tuned to annual surface velocities alone even if the constraining measurements are from a later year. To substantiate the assumption that the 1994/95 flow regime can be used to constrain the 1992/93 modelling, a comparison is made between our modelled annual composite flow regime and actual measurements of surface velocity for 1992/93 (Fig. 3). Area-wide spatial surface velocity vectors were derived from the measured displacement of prominent surface features on the aerial photographs from which the high-precision DEMs were constructed. Corresponding moulins, boulders and crevasse edges that feature in both sets of aerial photographs were digitized with respect to fixed reference points, and displacements were calculated accordingly. Spatial precision is good, with the estimated accuracy for initial feature location as about ±20 cm and that associated with digitization as about ±35 cm, yielding an rms error of about ±40 cm for marker digitization and hence an rms error of about ±57 cm for the total displacement vector.

Fig. 3.

Modelled annual surface velocity distribution at 70 m resolution comprising the time-weighted average of three representative flow snapshots: winter base flow, normal summer flow and enhanced spring flow at Haut Glacier d’Arolla (grey arrows), compared to actual velocities derived from the displacement of surface features between September 1991 and 1992 (black arrows and inset).

Thirty surface displacement vectors were derived for the period 1992/93 and are both consistent with each other and in good agreement with the pattern of modelled annual surface velocities (Fig. 3). Although spatial coverage of observed velocity vectors is limited to the glacier tongue, a plot of observed against modelled surface velocities indicates a high level of fit (R ² = 0.84) with a mean accuracy of ±0.7 m a⁻¹ which equates to an estimated error of < 15% for flow modelling (Fig. 3). Such a nice fit is hardly surprising given the highly optimized nature of the three flow snapshots used to construct the annual flow pattern. However, such a result does lend strong support for the assumption of interannual consistency of the flow regime at Haut Glacier d’Arolla, over the short term at least, and the estimated error of 15% provides approximate bounds on the flux-divergence modelling.

The patterns of flux divergence derived from each modelled flow snapshot, expressed as vertical displacement per year, exhibit large spatial variations (Fig. 4a–c). Through mass conservation, the net vertical displacement across the whole glacier must, by definition, equal zero, which is indeed the case. However, large local departures do occur and account for up to ±4.5 m a⁻¹ of vertical displacement. The length scale of this variability is typically hundreds of meters, indicating that the flow regime is complex. Some general trends can be identified and are most apparent in the winter base-flow snapshot (Fig. 4a). The upper basins tend to lose mass as a result of flux divergence in a zone of down-glacier flow acceleration and glacier thickening. Conversely, the glacier tongue tends to experience mass gain through flux convergence within a zone characterized by general glacier thinning and deceleration of the flow field. A notable exception to this pattern can be observed in the confluence zone, where high-velocity flow from the southwest basin converges with the southeast arm of the glacier and results in a zone of marked positive displacement up to 3 m a⁻¹.

Fig. 4.

Modelled annual vertical displacement at the glacier surface resulting from flux divergence for (a) winter base flow, (b) mean summer flow, (c) the enhanced spring flow and (d) the time-weighted annual composite pattern. Displacement is measured in metres and is contoured at 0.0 m

The general patterns observed in the winter base-flow scenario are compounded under mean summer sliding conditions; the tendencies for flow acceleration across the upper glacier and deceleration across the tongue are enhanced, resulting in a more pronounced pattern of flux divergence and convergence in these respective zones (Fig. 4b). The enhanced spring snapshot results in a highly disturbed pattern of vertical displacement, as would be expected under such extreme flow conditions (Fig. 4c). An area of >1 km² across the lower glacier tongue experiences strong deceleration, resulting in a zone of maximum vertical uplift of about 5 m a⁻¹ which is compensated by a zone of intense surface lowering up-glacier as a result of net flow acceleration. Although the distribution of these zones of vertical uplift and deflation is strongly related to the form of the imposed basal velocity perturbation, the magnitude of vertical displacement is in general agreement with estimates of actual uplift observed during the spring event (Reference NienowNienow, 1997).

Computed pattern of mass balance

The computed pattern of mass balance, calculated from the residual of the annual flux-divergence composite (Fig. 4d) minus the 1992/93 pattern of elevation change across the glacier, exhibits strong spatial variations reflecting the datasets from which it is was derived. Longitudinal variations are particularly marked, varying from around −5.0 m at the snout to nearly +1.0 m in the upper basins (Fig. 5). Lateral variations are also apparent, indicating relatively higher ablation and hence greater wasting towards the centre of the glacier compared to the margins. This is particularly pronounced along the western margin, where both shadowing from Mont Collon and a thick debris mantle suppress ablation.

Fig. 5.

The modelled pattern of mass balance (uncorrected for density), 1992/93, derived from the residual of the pattern of surface-elevation change minus the pattern of modelled surface displacement due to flux divergence. Inset is a bivariate analysis and least-squares linear regression between measured and modelled mass-balance interpolated onto the mean position of the stakes installed at the glacier centre line.

Comparison with measured mass balance

To evaluate the accuracy of the computed mass-balance pattern and to correct for density, a comparison is made with specific mass-balance observations at stakes installed along the glacier centre line (Fig. 1). Direct comparison of computed mass balance, interpolated onto the mean position of each of the stakes between 1992 and 1993, with the actual measurements at these stakes over the same period reveals a good correspondence (R ² = 0.91; Fig. 5). Furthermore, since the modelled mass balance is at this stage uncorrected for density, the regression line gradient of 1.13 yields a mean relative density of 0.88 which is in good agreement with the values specified for glacier ice (Reference PatersonPaterson, 1994).

Although the direct comparison of computed and observed values is very encouraging, the stake network which runs up the centre line (Fig. 1) is not optimally located to confirm the lateral spatial heterogeneity predicted in the computed mass-balance distribution. Moreover, measurements from the 14 mass-balance stakes do not constitute a statistically definitive test of the method presented. However, comparison of the density-corrected modelled mass-balance–elevation relationship for all ∼12 000 points on the glacier with that of the stake network (Fig. 6) indicates that the modelled values straddle observed mass-balance measurements fairly evenly.

Fig. 6.

Bivariate plot of specific annual mass balance, 1992/93, measured at the 14 stakes along the glacier centre line and the ∼12 000 modelled mass-balance values (density corrected by a factor of 0.88 derived from Figure 5) vs elevation.

Discussion

The good correspondence between measured and computed mass balance indicates that the approach presented here has the potential to provide an alternative method for determining spatial patterns of glacier mass balance. However, the method requires substantial qualification insofar as our chosen case-study has: (i) a self-contained and consistent geometry with a low aspect ratio, rendering it amenable to a modelling approach; (ii) good surface velocity measurements to optimize the modelling effort; and (iii) a temperate, isothermal regime. Furthermore, Haut Glacier d’Arolla is predominantly undergoing net ablation, and thus difficult snow- and firn-density corrections have been neglected.

The first qualification raised implies numerous cases where flow modelling is inapplicable: glaciers characterized by a high aspect ratio which have icefalls or other steep and irregular geometry; those experiencing strong decoupled behaviour (i.e. calving or actively surging glaciers); and those where basal accretion or melting processes contribute a significant proportion to the total mass balance.

The second point is less problematic but still of primary concern. For this method to be effective as a realistic means of deriving high-density patterns of mass balance, the required velocity data must be determined remotely. Currently in use or in development are numerous remote-sensing technologies, techniques and platforms capable of deriving annual, seasonal and intra-seasonal surface velocity fields at the necessary spatial and temporal resolutions for many glaciers worldwide. Surface velocities such as those determined in this study can be obtained by tracing the displacement of crevasses, boulders, ogives and other prominent surface features observed in pairs of aerial photographs or digital images (Reference RasmussenRasmussen, 1988; Reference Rentsch, Welsch, Heipke and MillerRentsch and others, 1990; Reference Etzelmüller, Vatne, Ødegård and SollidEtzelmüller and others, 1993), or the cross-correlation of sequential satellite imagery (Reference Bindschadler and ScambosBindschadler and Scambos, 1991; Reference Whillans, Jackson and TsengWhillans and others, 1993) and interferometric synthetic-aperture radar (InSAR) (Reference Frolich and DoakeFrolich and Doake, 1998). For valley glaciers, Reference Kääb and FunkKääb and Funk (1999) have developed a new photo-grammetric method to provide high-density surface displacements at high precision using multitemporal stereo models composed by aerial photographs taken with different temporal and spatial baselines. However, the method requires corresponding features in each photograph and does not work on snow-covered terrain or where features are subtle and change rapidly.

The third qualification, that Haut Glacier d’Arolla is temperate, is the biggest stumbling-block precluding almost all high-latitude ice masses with cold or polythermal regimes. Coupling a thermal component to the numerical flow modelling is one possible direction, but introduces yet another significant level of complexity and further degrees of freedom. For example, the surface temperature field, the geothermal heat flux and the thermal history of larger ice masses are required to provide boundary conditions to the problem. However, for glaciers and smaller icefields with relatively short response times of < 1000 years, coupled thermomechanical modelling may be a possible solution given good bounds on the geothermal heat flux and the surface temperature field. This latter field may be constrained by integration of multiple infrared images at the appropriate resolution or by ultra-high-frequency ground- or airborne-based radar systems which have had some success at locating englacial thermal boundaries (Reference BjörnssonBjornsson and others, 1996).

In conclusion, the method presented performs well at Haut Glacier d’Arolla, where abundant velocity data are available to optimize the flow modelling, where the thermal regime is isothermal and the predominant process is ablation. The question remains as to the general applicability of the approach given the amount and diversity of data required and the range of skills needed to collect them. Whether it could provide a practical means of deriving mass-balance information on more than a handful of glaciers, for which extensive datasets are available, remains to be seen.