A distributed surface energy-balance model for a small valley glacier. I. Development and testing for Haut Glacier d’ Arolla, Valais, Switzerland

Abstract

This paper describes the development and testing of a distributed surface energy-balance model used to calculate rates of surface melting at Haut Glacier d’Arolla, Valais, Switzerland. The model uses a digital elevation model (DEM) of the glacier surface and surrounding topography together with meterological data collected at a site in front of the glacier to determine hourly or daily totals of the energy-balance components and hence of melting over the entire surface of the glacier with a spatial resolution of 20 m. The model can also be used to determine temporal and spatial variations in snow depth, snow-line position and glacier surface albedo. Calculations from the model are compared with observations made along the glacier centre line 1990, and in general the model performs very well. The correlation coefficients between calculated and measured snow-line elevation, albedo and ablation are 0.99, 0.85 and 0.81, respectively. The main source of error between modelled and measured values of these variables is probably inadequacies in the parameterization of albedo used in the model.

Introduction

The energy balance at the surface pf an ice mass determines the amount of surface melting that takes place (e.g. Röthlisberger and Lang, 1987). Prevoius studies of the energy balance and melt of ice masses have tended to fall into two categories. The first group attempts to compare computations of the energy balance with measurements of melt at just a few locations over relatively short periods of time in order to rest and improve energy-balance theory (e.g. Hay and Fitzharris, 1988; Braithwaite and Olesen, 1990; Munro, 1990; Van de Wal and others, 1992), and the second group uses energy-balance theory to model glacier mass-balance gradients and examine how possible future scenarios of climate change may alter such gradients and therefore the specific mass balance of glaciers (e.g. Murno, 1991; Oerlemans, 1992; Oerlemans and Fortuin, 1992; Oerlemans, 1993). Thus, previous studies of energy balance and melt have tended to be either non-dimensional (at a few individual points) or one-dimensional (along glacier centre lines). There have been virtually no two-dimensional studies of energy balance and melt which have attempted to model spatial variations in these factors across entire glacierized catchments. Two notable exceptions are the work of Munro and Young (1982) who developed a model to predict net short-wave radiation variations across the basin of Peyto Glacier, Alberta, Canada, and that of Escher-Vetter (1985) who developed a model to compute the net radiation and turbulent heat fluxes across Vernagtferner, Austria. Howeve, these previous distributed models have several limitations. First, their spatial resolution was limited to 200 and 100 m, respectively.Secondly, although the Vernagtferner model considered all the main components of the energy budget, the Peyto Glacier model did not account for the turbulent heat fluxes. Finally while the models considered the effects of variations in slope angle, aspect and shading on the net shore-wave radiation budget they did not consider the effects of surface albedo variations particularly accurately. Albedo values were prescribed rather than generated internally by the model. Values of 0.24, 0.25, 0.5, 0.61 and 0.74 were used to represent the albedo of ice, bare ground, firn, old and new snow, respectively, in the Peyto Glacier model (Munro and Young, 1982), while values of 0.4, 0.6 and 0.8 were used to represent the values of ice, firn and snow in the Vernagtferner model (Escher-Vetter, 1980).

This paper describes the development and testing of a distributed two-dimensional surface energy-balance model for small valley glaciers using data collected at Haut Glacier d’ Arolla, Valais, Switzerland in 1990. The model was developed as part of a wide-ranging and ongoing study of the hydrology, water quality and dynamics of the glacier. Previous studies have dealth with the glacier’s drainage-system structure (Nicnow, 1993; Sharp and others, 1993; Hubbard and others, 1995), the nature of water-storage reservoirs within the glacier (Gurnell, 1993), and the nature of chemical weathering reactions beneath the glacier (Tranter and others, 1993; Brown and others, 1994a, b).

The development of a distributed model which can be used to determine temporal and spatial paterns of glacier melt is important for several reasons. First, the model can be used to study how the relative importance of the energy-balance components to surface melt change spatially across a glacier and temporally throughout individual melt season in response to changing climatic and surface conditions. Such studies are important for the accurate prediction of changes in the mass balance of glaciers and ice sheets under possible seenarios of climate change because, in areas of high relief, the influence of topography on radiation receipts and melt is complex, and can cause significant deviations from simple elevation/melt or mass-balance relationship, which are normally used in studies of glacier response to climatic change. The spatial and temporal patterns of energy-balance contributions to melt at Haut Glacier d’ Arolla form the basis of a second, related paper that is in preparation.

Secondly, the model can be used to calculate bulk surface-meltwater inputs to a glacier. Such calculations can be combined with measurements of rainfall to compute temporal patterns of total surface-water input to a glacier over the course of individual melt seasons. These may be compared with temporal patterns of water outputs (i.e. proglacial stream-water yields) to calculate temporal patterns of bulk water storage within a glacier. Water storage is believed to play an important role in controlling rates and mechanisms of subglacial chemical weathering (Sharp, 1991) and glacier motion (Kamb and others, 1994). Calculations of bulk water storage variation for Haut Glacier d’ Arolla will be presented in a subsequent paper.

Thirdly, the model could be coupled with a water-routing model of glacier hydrology to predict temporal patterns of runoff in proglacial streams. Such a model would be of considerable practical value of hydroelectric companies concerned with predicting both short-term variations in meltwater runoff in response to changing meteorological conditions, and with evaluating likely runoff patterns under different climatic-change scenarios (Willis and Bonvin, in press). The coupling of the surface-melt model with a watter-routing model for Haut Glacier d’Arolla forms the basis of ongoing work.

The Field Site

The model was developed and tested using measurements made during the 1990 summer season in the catchment of Haut Glacier d’Arolla. The catchment has an area of 11.7 km² of which 6.3 km² is glacierized (Fig 1). The glacier is about 4 km long and from about 2600 m a.s.l. at the snout to about 3100 m a.s.l. has quite low slope angles (10°). However, the upper accumulation area contains a series of steep icefalls on the north face of Mont Brulé up to an elevation of over 3500 m (Fig 1). The glacier is bounded by the high mountains of Mont Collon (3637 m) and L’Evêque (3716 m) to the west, Mont Brulé (3585 m) to the south and Bouquetins Ridge (3838 m) to the east (Fig 1). These mountains cast shadows across the glacier surface and therefore play an important role in the glacier’s radiation balance.

Fig. 1. Location map for Haut Glacier d’Arolla, showing weather-station and ablation-stake locations.

Methods

The model described there is used to calculate the surface energy balance and hence the surface melt over the whole glacier surface. The model requires four main inputs: (i) a detailed digital elevation model (DEM) of the glacier surface and the surrounmding area; (ii) knowledge of solar elevation and azimuth which is used in conjunction with (i) to compute patterns shading across the glacier; (iii) the initial distribution of snow depth across the glacier which is used in conjunction with (i) to compute the intial distribution of surface albedo; and (iv) meterological data collected at a site in front of the glacier that isused together with (i), (ii) and (iii) to compute the components of the energy balance and hence melt in each grid cell of the DEM. In addition to the effect of shading, the model considers the effects of surface slope and aspect on radiation receipt at the glacier surface. It also includes a parameterization of the albedo change associated with the removal of the surface snow cover. The model uses standard lapse rate and elevation-pressure relationships to compute the turbulent energy fluxes in each DEM grid cell from the meterological variables measured near the glacier snout. The main data used for testing the model are measured daily ablation rates over the glacier surface.

Elevation data

Elevation data are needed for the formulation of both radiative and tubulent energy fluxes. For the surface of the glacier itself, elevation data were collected during the summer of 1989 and 1990 using a combination of wild theodolite and Kern electro-optical distance-measuring equipment, and a Geodimeter 400 total station. Sharp and others (1993) give full details of the field surveying. A total of 823 points was surveyed across the accessible parts of the glacier.In order to generate a DEM for the whole catchment, these data were supplemented by contour data taken from Swiss National Survey 1:25 000 topographic maps of the area. Contours from these maps were traced and then scanned, using a Datacopy scanner with a resolution of 300 dots per inch. These raster data were then converted into irregular x,y,z data using line-following software developed in the Geography Department of Cambridge University (Mayo, 1993). These two irregular data sets were then interpolated on to a regular grid with a horizontal resolution of 20 m using the “Bilinear” interpolation routine in the UNIRAS graphics package (UNIRAS, 1990). Vertical resolution is better than 0.1 m for the surface of the glacier itself, and is 10 m for the surrounding topography.

Solar-Altitude and Azimuth Data

Solar-altitude and azimuth are needed for shading and aspect calculations, and were determined using standard astronomical theory (e.g. Walraven, 1978). The methods used to calculate shading and aspect in this study are discussed in the section below dealing with the energy-balance model.

Initial snow-cover data

A snow-depth survey was carried out on 13 June 1990 at 14 2.4 mm long wooden stakes distributed along the glacier centre line (Fig 1). As relatively few points could be surveyed (because of the need to get as close as possible to an instantaneous depth distribution), snow depth for each grid cell of the DEM was calculated from a liner-regression relationship between measured snow depth and elevation rather than from a spatial-interpolation routine. However, because this survey was carried out some 2 weeks after the weather station was set up, and we wanted to use the melt model for as long a period as possible (for seasonal water-balance calculations) an initial model run was done for the start of the season from 30 May to 13 June, but with the snow depth for 13 June. The total melt calculated by the model for this period at each stake was then added to the 13 June snow depths to give an inferred snow depth on 30 May. The resulting snow-depth/elevation relationships are given in Table 1. The measured snow depth for 13 June is denoted* and the inferred snow depth for 30 May is denoted† in Table 1 ; this second relationship was used in the model runs. Because the model requires snow-depth estimates in water-equivalent units, surface snow-density meausurements were also made at each stake. Three snow pits were dug in early June 1990 to assess the vertical variation in density within the snowpack. These showed relatively little vertical variation in density (except for a thin 5-10 cm saturated zone at the base) and relatively few ice lenses. Thus, the density was assumed to be uniform through the snowpack. The mean measured snow density was 0.49 g cm ³.

Table 1. Empirical relationships used in the model — linear regression. R², coefficient of determination; N, number of data points; F, value of F statistics; p, confidence level at which the F value is significant. Figures in brackets are standard errors for linear regression slopes and intercepts

Meteorological data

Meteorological data needed to drive the model were collected using a Delta-T automatic weather station situated 100 m in front of the glacier at an elevation of 2547 m (Fig 1). Incoming short-wave radiation, air temperature, wind speed, wind direction, relative humidity and precipitation were measured at 10 min intervals. Hourlt totals of precipitation and hourly averages of the other variables were logged. These hourly values constituted the data used by the model. A second weather station was established from 9 July to 23 August at an elevation of 2846 m on the glacier surface (Fig 1). Data from this station were used to test some of the assumptions regarding atmospheric lapse rates used in the model, and to test how well short-wave radiation variations over the glacier surface were represented by the model.

Albedo data

The albedo of a glacier surface plays a fundamental role in the absorption of short-wave radiation. As one of the aims of this study was to parameterize and then model the change in albedo over the course of a melot season, albedo was measured at approximately weeklt intervals at each of the centre-line ablation stakes during the 1990 ablation season. These data were used to test and adjust the parameterization of albedo used in the model. Measurements of reflected radiation in the visible-near infrared range (400-1200 m) were made using the photoelectric sensor of a Milton multi-band radiometer. This sensor had a field of view of 15°, and the radiometer was held perpendicular to the surface at a heitht of 30 cm. At this height, the area sampled was approximately 50 cm², and the recorded values were insensitive to the precise angle of the sensor.

The reflectance of the glacier surface was compared with a reference surface, for which a Kodak grey card was used (Milton, 1989). Detailed calibration showed that the reflectance of this card in the 400-1200 nm range was 21.5%. and that this reflectance did not deteriorate during use. The reflectance of the glacier surface (R) was calculated from the photoelectric voltage measured from above the glacier surface and above the grey card as follows:

(1)

where V _s is the voltage over the glacier surface, V _r is the voltage over the grey card, V ₀ is the offset voltage and K is the reflectance of the grey card. The offset voltage is included in this expression because, even with the sensor in complete darkness, a small current flows through the sensors. The offset voltage (i.e. the voltage recorded with the sensor shielded form light) was measured at the start and end of each weekly albedo survey, and the mean of the two readings was used.

At each stake, three sets of measurements were made at three closely spaced locations (at the stake itself, and at locations 1 m either of the stake transverse to the direction of glacier flow) in order to try to remove small-scale variations in albedo. The albedo measurements used for model calibration (see Equations (9)-(11) below) were the arithmetic mean of these three measurements. Generally, the three values lay within ±0.05 of the mean value.

Ablation data

Measured ablation data form the main test of the model. Data were collected at the 14 stakes distributed along the glacier centre line(see Fig 1). These were drilled into the glacier until approximately 10 cm remained above the surface. Ablation rates were measured by lowering a plastic disc 10 cm diameter over each stake until the outer edge of the disc just touched the snow or ice surface. The amount of surface lowering was calculated by measuring the length of stake above the disc and comparing this with the previous day’s reading. This figure was then converted into a daily ablation rate. Stakes were redrilled into the surface as required, generally when approximately 30 cm of stake was left below the surface.

When snow was present at a site, surface snow-density measurements were also made by carefully inserting a cylinder of known volume into the edge of a shallow trench dug at the time of measurements. The snow was then weighed using a spring balance and the density was calculated. This technique did have problems, however. Great care was needed when inserting the cylinder, in order to avoid compacting the snow sample, and very wet or very dry snow tended to fall out of the cylinder very easily. For this reason, three measurements were made at each site and the mean value was taken to represent snow density. As a further safeguard, if any one measurement was different from the other two by more than 10% it was rejected and another sample taken. If no snow was present, a constant density of 900 kg m⁻³ for ice was assumed. Knowledge of the density of the surface material allowed the measured surface lowering to be converted into water equivalent units, which could then be compared with the values calculated by the model.

The Energy-balance Mode

The glacier surface energy-balance model used in this study incorporates elements from several previous energy-balance studies. These generally agree that short-wave solar radiation in the dominant source of melt energy (e.g. Munro and Young, 1982; Braithwaite and Olesen, 1990; Oerlemans, 1993; Paterson, 1994), and the model described here treats this source of energy in the greatest detail. The model also includes turbulent heat fluxes and long-wave radiation, which are acknowledged to make a smaller, but nevertheless significant, contribution to surface melt energy. It is assumed that the energy available for melt (Q _M) can be obtained from:

(2)

where Q ^* is the net short-wave radiation flux, I ^* is the net long-wave radiation flux, and Q _S and Q _L are sensible and latent turbulent heat fluxes, respectively. For convenience, energy fluxes are expressed in equivalent ablation units (mm of water per unit time) by dividing the energy flux by the latent heat of fusion of water (L) and the density of water (ρ _w) (Table 3). These four energy-balance components are calculated for each grid cell of the DEM of the glacier surface and surrounding topography, and are then used to compute hourly values Q _M for each grid cell over the course of the 1990 melt season from 30 May to 24 August. The computation of each energy-balance component is now considered in detail. Parameter values used in the model that are derived from data collected for this study are given in Table 2 ; other parameter values, typically standard physical constants or empirical values used in other studies, are given in Table 3.

Table 2. Empirical relationships used in the model — non-linear regression. R², coeffiecient of determination

Table 3. Parameter values used in the model

Short-wave radiation

The short-wave radiation flux (Q ^*) is given by:

(3)

where α is the surface albedo, Q ^′ _i is the instantaneous flux of direct and diffuse solar radiation received on the surface at a particular location (W m⁻²), L is the latent heat of fusion of water (J kg⁻¹) and ρ _w is the density of water (kg m⁻³). Q ^′ _i for each grid cell of the DEM was calculated from global radiation (Q ^′) measured at the weather station near the glacier snout. This value was then adjusted to take account of the effects of slope, aspect and shading by the surrounding topography as detailed below.

Slope angles and aspects for a given grid cell are obtained by examining four of the grid cells surrounding it:

(4)

(5)

where z^′ is the angle of the slope from the horizontal (always positive), A^′ is the aspect of the slope, measured as degrees away from due south, positive to the west, and negative to the east, h is the interpolated elevation of a given prid cell, subscripts “i” and “j” are the x and y grid positions of the grid cell in question (where grid cell (i = 1, j = 1) is the southwest corner of the DEM grid, which is oriented north/south), and s is the length of the side of the gride cells in the DEM (20 m).

Topographic shading is determined as follow. All grid cells in the DEM are intially designated as being “in-sum”. For each grid cell of the DEM, in order of nearness to the direction of the Sun, the computational algorithm “walks” away from the start cell alon the path of the solar beam, in steps equal to the DEM grid size. At each step, two tests are made. If the elevation of the new grid cell is above the height of the solar beam at that point, the start grid cell of the current walk is designated as “in-shade”, the current walk stops and the next start grid cell is selected. If the new grid cell is lower than the height of the solar beam and has been designated on a previous walk as “in-sum” othen the start grid cell stays designated as “in-sum”, the current walk stops and the next grid cell of the current walk. If neither of these conditions has been met by the time the walk reaches the edge of the DEM, the original grid cell remains designated as “in-sum”. This procedure is quite efficient, as generally one of the conditions is met quite quickly, so most walks do not need to go all the way to the edge of the DEM.

If a given grid cell is designated as “in-sum”, the radiation received by the radiometer (Q ^′)is converted to the equivalent radiation received by a surface normal to the Sun’s rays (Q ^′ _n) using the relationship

(6)

where z is the angle of the Sun above the horizon.

The radiation received by the sloping glacier surface (Q ^′ _i, Equations (3) is then

(7)

where A is the solar azimuth, defined as degrees away from due south, positive to the west of due south and negative to the east.

If a grid cell is designated as “in-shade”, however, only diffuse radiation is allowed to reachg the surface. Because of the lack of detailed cloud-cover records, this is treated more schematically than direct solar radiation. Following Oerlemans (1993), diffuse radiation from the sky is assumed to be one-fifth of the meausred global solar radiation in all cases. Reflected radiation from the surrounding topography is then added to this value. Thus, the radiation received by a shaded surface is calculated using a view-factor relationship:

(8)

where α_m is the mean albedo for the whole basin (Table 2 ). In this equation, the first term on the righthand side represents the radiation reveived from the part of the sky hemisphere visible from the grid cell in question (assuming isotropic radiation from the sky), and the second term represents the reflected radiation received from that part of the ground hemisphere visible from the grid cell, again assuming isotropy (Munro and Young, 1982).

The slopes and aspects calculated for each grid cell of the DEM are shown in Figure 2 and 3, respectively. As can beeb, slopes on the glacier are generallly below 10°, except on the north face of Mont Brulé and on parts of the western firn basin above 2900 m, and the glacier generally faces between 270° and 45°. Thus, the general effect of slope and aspect variations is to reduce the incident solar radiation received by the glacier surface. These effects are discussed in the results section below.

Fig. 2. Slope angles on Haut Glacier d’ Arolla calculated from the DEM.

Fig. 3. Surface aspect on Haut Glacier d’ Arolla calculated from the DEM.

The effect of shading of the glacier surface by the surrounding topography on short-water radiation receipts is also important, due to the enclosed nature of the basin. Figure 4 shows the effect of shading at 0700, 0900, 1800 and 1900 h on 21 June. The main glacier tongue is the most affected by shading; in the morning it is shaded by Bouquetins Ridge, and in the afternoon by Mont Collon and L’ Evêque (cf. Fig 1). By contrast, the western and eastern firn basins of the glacier are generally not shaded until very late in the day. Thus, daily radiation totals received by the lower glacier are generally lower than those for the upper parts. As the summer progresses, these patterns remain broadly similar, but the main tongue of the glacier remains shaded for longer in the mornings, and becomes shaded earlier in the afternoon due to the shorter days and lower solar heithts after the summer solstice. Again, these effects are discussed below.

Fig. 4. Three-dimensional computer images of Haut Glacier d’ Arolla, looking south-southeast towards Mont Brulé, with simulated shading patterns generated by the model for 21 June. (a) 0700 (b) 0900 (c) 1800 (d) 1900 h.

Grid-cell albedos are calculated within the model using relationships of the type proposed by Oerlemans (1992, 1993), as follows. The parameterization used is an empirical one based on observed albedo variations on many valley glaciers. Oerlemans (1992, 1993) assumes that albedo of ice surfaces is elevation-dependent, and is generally lower at low elevations, due to increased debris and dust in the ice. Thus, a map of this background ice albedo (α_b) can be defined as:

(9)

where h is the elevation of an individual grid cell (m a.s.l.), E is the elevation of the equilibrium line n(m a.s.l.) and a ₁ and a ₂ are adjustable parameters. The actual surface albedo is then affected by the presence of snow, and by “weathering” of the surface (be it ice or snow). Oerlemans assumes that the depth of any snow cover is important, as this determines how much of the background surface can be “seen” through the semitranslucent snow cover, and that “weathering” of the surface can be approximated by the total amount of melt that has occurred since the beginning of the melt season. Oerlemans uses an exponential curve for the effect of snow depth, and assumes a linear relationshi of albedo with “Weathering”, thus:

(10)

where α_s is the resulting surface albedo, α_os is the standard albedo of old snow, d is snow depth (m w.e.), M is the cumulative melt (m w.e.) and a ₃ and a ₄ are adjustable parameters. Thus, if no snow is present (d = 0), the exponential term equals 1, and therefore the surface albedo is the same as the background albedo, less the effect of M. In the model, the standard albedo of old snow has a value of 0.72, a value lower than that of fresh snow, to allow for a standard aging effect (Oerlemans, 1992, 1993). Field experience at Arolla suggested that even very small amounts of new snow falling during the ablation season had an important effect on albedo, and this was taken into account using a similar exponential relationship relating the depth of new snow and the albedo of the underlying surface:

(11)

where α_s(n) is the resulting surface albedo in the presence of new snow, α_ns is the standard albedo of fresh snow, and p is the depth of fresh snow (m w.e.). The large absolute value of the multiplier of p (5000) ensured that small snowfalls affected the albedo. The albedo of fresh snow was taken as 0.85. Precipitation in a given cell was assumed to fall as snow if the air temperature in that cell at the time 1°C or (cf. Oerlemans, 1993). Subsequent melt was first removed from the depth of new snow. If any new snow remained unmelted after 3 days, its depth was added to the depth of any old snow, the depth of new snow was reser to zero, and the albedo was then calculated using Equations (10). This procedure simulates, albeit crudely, the aging of newly fallen snow.

Measurements of ice albedo on Haut Glacier d’ Arolla were generally lower than those reported by Oerlemans (1993), and so the parameters were adjusted using measured albedo. The parameters in Equations (9) (a ₁ and a ₂) for the background albedo were calculated using the first measured albedo value after the snow line had retreated above any particularly stake; the parameters in Equations (10) (a ₃ and a ₄) were calculated using all the measurements of albedo, snow depth and cumulative melt, together with the background albedo as estimated by Equations (9). In both cases, non-linear ordinary least-squares curve-fitting techninques were used. The resulting parameter values, correlation coefficients and standard errors are given in Table 2. The final surface albedo had maximum and minimum values of 0.85 and 0.12 imposed if necessary.

Long-wave radiation

The long-wave radiation flux is given by

(12)

Where I↓ is the incoming long-wave radiation, and I↑ is the outgoing long-wave radiation, both in W m². Assuming that the glacier surface is always at 0°C, and it radiates as a black body, I↑ has a constant value of 316 W m⁻². I↓ is given by Braithwaite and Olesen (1990):

(13)

where

is the effective emissivity of the sky, σ is the Stefan-Boltzmann constant (W m⁻¹ K⁻⁴), and T _a is the absolute air temperature (K). Effective emissivity is calculate as:

(14)

where k is a constant depending on cloud type; n, from 0 to 1, is cloud amount; and

is clear-sky emissivity. Ohmura (1981) gives k values for eight different cloud types, but as data collected for this study did not include this information, a constant value of 0.26 was used, the mean of the values for altostratus, altocumulus, stratocumulus, stratus and cumulus cloud types, following Braitwaite and Olesen (1990). Twice-daily to 5 July. The daily means of these data were then regressed against various meteorological parameter to try to find a relationship that could be used in the model (as the meterological data collected by the weather stations did not include cloud cover). The most effective regression relationship was based simply on the daily temperature range, and is given in Table 1. Following Ohmura (1981) the clear-sky emissivity is given by:

(15)

Turbulent heat fluxes

Turbulent heat fluxes for each cell are calculated using the relationship derived by Ambach (1986), based on energy-balance measurements on the Greenland ice sheet. For a melting glacier surface (i.e. with a temperature of 0°C and vapour pressure equal to the saturated vapour pressure at 0°C), and assuming adiabatic stratification in a Prandtl-type boundary layer (in which the vertical fluxes of energy and momentum are constant with height, and nearly all of the total increase in wind speed from the ground to the free atmosphere takes place (Kraus, 1973)), the relationship are:

(16)

and

(17)

where K _s and K _l are coefficients, P is atmospheric pressure (Pa), T is air temperature (°C), V is wind speed (m s⁻¹), and δ_e is the difference between the vapour pressure of the air and the saturation vapour pressure at the glacier surface (Pa). The subscript “2” indicates that the measurements were made at 2 m above the surface. Numerical values of K _s and K ₁ include the adjustment for latent heat of fusion and density of water, and vary depending on whether the surface is ice or snow (due to their different roughness), and for K ₁ on whether the surface is experiencing condensation or evaporation (see Table 3 ). Temperature for each cell in the DEM is calculated using a standard mean atmospheric lapse rate (see Table 3 ) and the difference in altitude between the cell in question and the weather station. The δ_e values for each cell is calculated using measured relative humidity (Which is assumed to be constant with elevation) and the lapse-rate corrected air temperature. The measured relative humidity from both weather stations showed a diurnal cycle, with similar maximum values recorded at around dawn, when air temperature were at their lowest. The lower weather station showed much larger diurnal variations, but this station malfunctioned quite frequently, with anomalously low (or even negative) values recorded, though only during daylight hours. Thus, whether the increased diurnal range at this location, even on days when no extreme values were recorded, was a real phenomenon or an artifact of the weather station could not be determined with confidence. In view of these problems, the relative-humidity data from the upper weather station were used when available. During periods when the upper weather-station data were not available, and the lower weather station was malfunctioning (deemed to be when humidity was recorded as being less than 40%), a value of 75% (the mean value for the season from both weather stations) was used. This problem shoule not be of great importance, as latent heat fluxes form the smallest of the four energy-balance components (see below).

Model Results

Validation

The model produces three main output data sets which can be compared be compared with field data in order to evaluate the model’s performance. These are the pattern of snow-line retreat during the course of a melt season, the pattern of albedo change over the glacier surface through time, and the pattern of ablation Over the glacier surface through time. In this paper, we examine how well the model outputs match observed variations in these data sets, using calculations and measurements made along the centre-line stake network shown in Figure 1. More detailed spatial and temporal variations in these data sets will be dealt with in a subsequent paper.The accuracy of the model is assessed on the basis of the Pearson correlation coefficients and regression standard errors, summarised in Table 4.

Table 4. Regression relationships between modelled measured variables. Figures in brackets under regression slopes and intercepts are standard errors; r, Pearson correlation coefficient; N, number of data points; M, Modelled parameter; O, observed parameter

Snow-line retreat

The observed and modelled pattern of snow-line elevation through time is shown in figure 5. The model performs remarkably well; the only obvious discrepancy is the snowfall (and associated fa11 in snow-line elevation) in the first week of July, which is not captured by the model, though no retreat of the model snow line occurs during this period. This mismatch would seem to be caused by the model underestimating snowfall events. This could be due either to the use of only one weather station to determine conditions over the whole basin, or to the weather station itself under-recording snowfall events. The correlation coefficient between observed and calculated snow-line elevation is 0.99.

Fig. 5. Modelled and observed snow-line elecation for 1990. Solid line, observed; dashed line, modelled.

Albedo

Field albedo measurements were used to tune the model equations used to calculate albedo, and so cannot really be used to calculate albedo, and so cannot really be used as an objective test of the model. Nevertheless it is instructive to see how well the equations used can simulate the changing albedo of the glacier surface. Measured and modelled albedo for four of the stake locations through the course of the melt season are shown in figure 6. For the two lower stkes (B and F), the model slightly underestimates the albedo of the glacier surface for the first half of the summer. The general downward trend in albedo is captured, and the lowest model and measured albedo values agree quite well. The model does not capture the slight end-of season rise in albedo; this effect is difficult to explain, as no snow fall during this period. For the two higher stakes (I and N), the model performs rather better. Again, there is some discrepancy stakes. After the end of June, the modelled albedo decreases in a way very simila to the observed values. The relationship between all measured and modelled albedos is shown in figure 7. For 106 data points, the correlation coefficient is 0.85. The slope of the regression relationship is slightly greater than 1; together with the small negative intercept, this suggests that the model marginally overestimates high albedo values, and underestimates low values. The good relationship between measured and modelled albedo values is not surprising, given that the parameters in the albedo relationship used by the model were adjusted using the measured albedo values, but it does show that the basic form of the relationship, and the main controlling variables (elevation, snow depth and cumulative melt) chosen by Oerlemans (1993) and used in this study can explain much of the albedo variation observed on alpine glaciers. The main problem with the relationship is that the debris concentration, which is one of the main controls on the background ice albedo, is only very approximately dependent on elevation at the scale of an individual glacier, where local flow patterns and debris availability strongly influence ice debris content. Thus, Equations (9) performs rather worse than Equations (10) (see Table 2 ); this inaccuracy affects all the model albedo estimates, however.

Fig. 6. Modelled and observed albedo variations for four stake locations: (a) stake B, (b) stake F, (c) stake I (d) stake N. solid line, observed; dashed line, modelled.

Fig. 7. Scatter diagram of modelled vs observed albedo values. Dashed line, fitted relationship (see Table 4 ).

Ablation data

Time series of measured ablation rates from the main test of model performance. Measured and modelled ablation for stakes B, F, J and N shown in figure 8. The magnitude of the modelled ablation agrees well with the measured values, and the general shape of the curves through the melt season also matches quite well. However, the modelled ablation is generally less variable at shorter time-scales than the measured ablation, which shows quite large short-term oscillations. Some of this discrepancy may be associated with difficulties in measuring ablation rates over short time intervals (Müller and Keeler, 1969; Braithwaite and Olesen, (1990).

Fig. 8. Modelled and observed ablation variation for four stake locations: (a) stake B, (b) stake F, (c) stake I, (d) stake N. Solid line, observed; dashed line, modelled.

The relationship between all the measured and modelled ablation values is shown in figure 9. There is obviously some scatter, but the overall correlation for the 985 data points is 0.81. The slope of the regression equation is less than I which suggests that the model under-predicts high and overestimates low ablation values. The small positive value of the intercept (0.97) also means that the model overestimates low ablation. This may help explain the smoother nature of the modelled ablation curves compared to the measured values. For both coefficients thestandard errors are very small (0.02 and 0.08, respectively), which suggests that the predictive power of the model is good. The standard deviation of the difference between the measured and modelled daily ablation values for all the stakes is 12.9 mm w.e.d⁻¹, confirming the explanatory power of the model. This is somewhat more accurate than the standard deviations of 13.6 and 18.9 mm w.e.d⁻¹ found by Braithwaite and Olesen (1990). The mean difference between the measured and modelled daily ablation values for all the stakes over the summer is - 0.3 mm w.e.d⁻¹.

Fig. 9. Scatter diagram of modelled vs observed ablation values. Dashed line, fitted relationship (see Table 4 ).

In order to further identify the origins of the discrepancies between the measured and modelled ablation values, the differences between the observed and predicted melt values (hereafter called the model error) were correlated with the observed melt values to investigate whether the model was systematically predicting melt less accurately at high or low ablation values (typifying high-pressure or low-pressure climatic conditions). However, there was no significant correlation (r = 0.006 ), suggesting that the model was performing equally well at both high and low melt values. This suggests that the errors result either from the measured ablation values (cf. Müller and Keeler, 1969; Braithwaite and Olesen, 1990) or from a simplificationin the model which is not directly weather-related. This latter possibility was investigated by correlating the model error against the individual energy-balance components. These gave very low, statistically insignificant correlation coefficients (see Table 4 ) except in the case of the short-wave radiative flux (r = 0.24). This regression of the model error against the short-over wave radiative flux gave a negative slope coefficient which implies that the model is over-predicting melt rates during periods of high short-wave radiation inputs. There are two possible reasons for this. The first is that the modeld does not accurately predict short-wave radiation arriving at the glacier surface; the second is that the model does not accurately predict the absorption or solar radiation was tested using data from the second weather station on the glacier surface (Fig 1.) The coefficient between the modelled incoming short-wave radiation at this site and the measured values was 0.95, suggesting that the model predicts incoming radiation very well. Thus, the more likely source of error is in the absorption of short-wave radiation by the glacier surface. This confirms the suggestion above that glacier-surface albedo is still not represented as accurately as desired. This is especially true on the lower glacier where the model generally shows lower albedo values than those measured, particularly during the early part of the summer. Under-prediction of the albedo would then lead to higher than expected melt rates, especially during periods of high radiation receipts.

Energy partitioning

The model, having been tested, can also be used to evaluate the importance of the different contributors to melt energy, and the role that topography plays in melt energy. For the model run discussed above, the mean daily contribution to the surface energy budget are shown in Table 5. As can be seen, the dominant contributor of melt energy is short-wave radiation, which supplies a daily mean for the whole glacier surface for the whole season of 24.5 out of a total energy budget of 17.3 mm w.e.d⁻¹. Long-wave radiation is the main cause of heat loss from the glacier surface (−10.4 mm w.e.d⁻¹), and is largely responsible for the total energy budget being smaller than the short-wave inputs. The turbulent fluxes are less significant; sensible heat provides a small net input of energy (4.4 mm w.e.d⁻¹), and latent heat transfer is a small source of heat loss (−1.2 mm w.e.d⁻¹). This is in broad agreement with most previous studies of glacier energy balance, which are dominated by radiative heat transfer (e.g. Braithwaite and Olesen, 1990; Munro, 1990), with small positive contributions from sensible heat transfers, and smaller positive (e.g. Munro, 1990) or negative (e.g. Braithwaite and Olesen, 1990) contributions from latent heat.

Table 5. Model energy-balance components. Units are mm w.e.d⁻¹

However, for this study, the short-wave radiative heat flux is rather higher than in most other studies. This would seem to have two possible causes. First, the albedo of Haut Glacier d’ Arolla (once snow has melted) is rather low; the mean of all measured ice-surface albedo values is 0.19, which is at the lower end of Oerlemans’ (1993) “dirty-ice” category, and is lower than the values of 0.24 and 0.4 used by Munro and Young (1982) and Escher-Vetter (1985), respectively, in their energy-balance models. This would increase the absorption of solar radiation. Secondly, 1990 seems to have been a warm and sunny year. The snow line retreated very rapidly,to very high elevations (over 3000 m). The average July and August temperature measured at a nearby Swiss meteorological station at Bricola was approximately 0.5°C higher than the long-term average from 1968 to 1994. Similarly, solar radiation received during these two months was nearly 7% greater than the long-term average.

This effect of low ice albedo is confirmed if the results are partitioned into ablation values recorded on snow or ice surface (see Table 5 ). For ice, short-wave radiation supplies 33.6 out of a total gain of 26.0 mm w.e.d⁻¹, but for snow (which has a higher albedo) it only supplies 8.3 out of total loss of −11.9 mm w.e.d⁻¹. Sensible heat flux is higher for ice, and is slightly negative on snow; latent heat fluxes are negative on both surfaces, and, again, long-wave radiation is the dominant form of energy loss on both surfaces.

This overall loss of energy from snow surfaces seems contradictory, as it implies snow should not melt. However, it should be borne in mind that melt occurs only when the total energy balance is positive. This typically occurs only during daylight hours, when high short-wave radiation inputs counteract any possible energy losses. During darkness, the surface energy balance is often negative, due larglely to the emission of long-wave radiation. For snow surfaces, this effect dominates the seasonal total energy balance. If the various energy contribution are summed for periods when melt is occurring (i.e. when the total energy flux is positive) rather than for the whole season, a slightly different pattern emerges (Table 5 ). Sensible heat flux is positive for both surfaces, (though lower for snow, due to the colder temperatures at the higher elevations at which snow is found). Latent heat flux is small for both surfaces, but is positive for ice and negative for snow. This again seems to be due to the occurence of snow at higher elevations. Short-wave radiation is the main source of energy inputs, but ∼4 times smaller on snow-covered surfaces than for ice surfaces. This difference is largely due to albedo differences, though the snow on the north face of Mont Brulé (where aspect greatly reduces short-wave radiation receipts) also accounts for some of this reduction. The smaller negative values for long-wave radiation show that even during daylight the long-wave flux remains negative. The smaller values are due simply to only about half of the loss occuring during daylight.

The role of topography

The model was also used to evaluate the role that topography plays in melt-energy receipts given that short-wave solar energy forms the dominant component of melt energy. Three additional model runs were completed, in which, first, aspect was not taken into account; secondly, shading was not taken into account; and thirdly, neither was accounted for. The effects of these on energy receipts are shown in Table 5. Excluding the effects of aspect increases short-wave energoy receipts by 14.7%; excluding shading increases receipts by 5.2%; and excluding both increases receipts by 20.8% (equivalent to 5.1 mm w.e.d⁻¹).

Conclusions

We have developed a three-dimensional distributed energy-balance model to calculate spatial and temporal variations in the energy-balance components and hence melt across small valley glaciers. We have developed and tested the model using data collected from Haut Glacier d’Arolla, Valais, Switzerland. The model is based on simplifications of energy-balance theory derived by Ambach (1986) and successfully used by Braithwaite and Olesen (1990) to model the ablation at the margin of the Greenland ice sheet. It treats the short-wave radiation balance in most detail as this energy-balance component is generally acknowledged to make the biggest contribution to the melt of alpine glaciers (Munro and Young, 1982). The effects on the short-wave radiation budget of surface slope, aspect and albedo and of shading from surrounding topography are accounted for in the model. Surface albedo is calculated internally by the model and is based on parameterizations developed by Oerlemans (1992, 1993). We believe that the model represents an improvement on previous distributed models (e.g. Munro and Young, 1982; Escher-Vetter, 1985), which have a lower spatial resolution and have no parameterization for albedo. However, though generally quite effective, the albedo parameterization used here cannot explain all the variance in albedo. In particular, local spatial variations in debris concentrations in the ice on individual glaciers cause large variations in the surface albedo, which are not accounted for. This problem will be difficult to solve, as such variations depend strongly on the flow patterns and debris availability of the glacier concerned, rather than simply on elevation, snow cover and cumulative melt.

Over the 1990 summer melt season, the mean difference between the modelled and observed ablation at 14 stakes distributed along the glacier centre line was just −0.3 mm w.e.d⁻¹. with small adjustments to some of the parameter values this error could be reduced to zero if desired. The rather small standard deviation for the daily ablation of 12.9 mm w.e.d⁻¹, which is somewhat smaller than those derived from the Greenland study of Braithwaite and Olesen (1990), indicates that the model in this study generally performs well.

The model results confirm that solar-radiation inputs are the largest source of melt energy to alpine glaciers. The model also shows that topographic effects play a very important role in determining receipts of solar radiation. Given the importance of solar radiation to glacier senergy budgets, these effects should not be ignored in energy-balance studies, especially for glaciers in areas of high relief. Neglecting these effects would lead to an over prediction of melt, which could be very significant in mass-balance or water-balance studies of glacierized catchments. The model developed here can also be used to study the spatial and temporal variations in the relative importance of the individual energy-balance components to surface melt of valley glaciers and to compute surface meltwater inputs as part of a study of the water balance of valley glaciers. Such studies for Haut Glacier d’ Arolla will from the basis of subsequent papers.