Assessment of the surface mass balance along the K-transect (Greenland ice sheet) from satellite-derived albedos

Abstract

This paper explores the potential of using satellite-derived albedos to estimate the surface mass balance of the Kangerlussuaq transect (K-transect; Greenland ice sheet). We first retrieved surface albedos from Advanced Very High Resolution Radar data by using, among other techniques, a new cloud detection algorithm based on the relation between brightness temperature and surface elevation. We then computed the ‘satellite-derived mass balance’ (b_sat) from the mean albedo for the transect, by taking fixed values for atmospheric transmissivity and the longwave and turbulent fluxes. We found that b_sat explains 7 1% of the variance in 13 years of stake mass-balance measurements (b_m). Our method also provides good estimates of the magnitude of the interannual variability in b_m. The performance of the method degrades considerably without correction for anisotropic reflection at the surface and recalibration of the satellite sensors with dry snow at the top of the ice sheet. Sensitivity tests indicate that the method’s performance is hardly sensitive to uncertainties in parameters. Therefore, we expect that the method could be successfully applied on other glaciers and parts of ice sheets and ice caps, especially where accumulation rates are relatively small. We show that the investigated method performs best just below the mean equilibrium-line altitude.

1. Introduction

Monitoring the surface mass balance of glaciers and ice sheets is a major task of the glaciological community. However, direct measurements are limited in space and duration. This is certainly true for the ablation area of the Greenland ice sheet, for which only eight measurement series of at least 5 years exist (Six and others, 2001). Among these eight series, the measurements along the Kangerlussuaq transect (K-transect; Fig. 1), which extends almost exactly from east to west along 67˚ N on the western side of the Greenland ice sheet, currently constitute the longest record (1990/91 to present; see Greuell and others, 2001; Van de Wal and others, 2005). In this situation, estimates based on data other than direct measurements (e.g. satellite data) are extremely valuable, provided these estimates are of sufficient quality.

Fig. 1. Location of the Kangerlussuaq transect (K-transect) on the Greenland ice sheet, and blow-up of the K-transect. Mass-balance measurements were performed at all the given sites (S), and elevations are the average for the period 1995–2000. Rectangular flags indicate sites where automatic weather stations were operated.

Abdalati and Steffen (2001) used passive microwave satellite data to determine temporal variations in spatial extent of melt on the Greenland ice sheet, but they did not arrive at an estimate of the surface mass balance. This was achieved later on by Mote (2003), who also computed melt extent from passive microwave satellite data. He then converted melt extent into annual ablation by means of a degree-day model, making an assumption about the interdaily variability in the temperature. He finally estimated the annual surface mass balance of the Greenland ice sheet by combining annual ablation with annual accumulation determined from European re-analyses (ERA-40).

To the knowledge of the authors, two methods exist that use satellite-derived albedos to estimate the surface mass balance. The first method is to determine the maximum elevation of the snow-line for each ablation season. Here, the snow-line is the line separating snow that has fallen in the actual balance year (above the snow-line) from ice and firn (below the snow-line), where firn is defined as snow that has been deposited on the glacier prior to the balance year under consideration. Demuth and Pietroniro (1999), Greuell and Knap (2000) and de Ruyter de Wildt and others (2002) attempted to infer the mass balance from the satellite-derived maximum position of the snow-line. However, this method suffers from a fundamental problem, as demonstrated by de Ruyter de Wildt and others (2002) and elucidated by Figure 2. It occurs when the annual mass balance is smaller than the previous year’s mass balance. Before the end of such a balance year, the snow-line will pass the snow-line of the previous year. When the snowline continues to move upwards, firn becomes exposed at the surface. The problem is that, so far, no satisfactory way has been found to distinguish firn from snow, no matter whether synthetic aperture radar (SAR) backscatter coefficients or albedos are used. Therefore, in many years the snow-line cannot be tracked until the end of the ablation season, so no reliable estimate of the mass balance can be made.

Fig. 2. (a) In a year with a relatively positive mass balance, the snow-line will be at a relatively low elevation at the end of the ablation season. In this case, the snow-line separates snow from ice, so it can be detected from satellite imagery. (b) In a year with a relatively negative mass balance, the snow-line will be at a relatively high elevation at the end of the ablation season. The snow-line now separates snow from firn (snow that is >1 year old). Because snow cannot be distinguished from firn, the snow-line cannot be detected from satellite imagery. (Figure adapted from de Ruyter de Wildt and others, 2002.)

De Ruyter de Wildt and others (2002) suggested an alternative method. They combined satellite-derived albedos with the incoming radiation at the top of the atmosphere to estimate the amount of shortwave radiation absorbed at the surface. By performing integrations in space and over time, they then computed the radiation absorbed during the whole ablation season and averaged over the glacier area. Finally, they tested whether this quantity could be used to estimate the mass balance by comparison with direct mass-balance measurements. By imposing linear relationships, they were able to explain 76–88% of the observed variance for six outlets of Vatnajökull, Iceland.

There are several reasons why the surface mass balance can be estimated from the satellite-derived albedo and therefore from the satellite-derived absorption of shortwave radiation. First, the net shortwave radiative flux generally constitutes the largest energy flux towards the surface during the ablation season (Ambach, 1979; Hock and Holmgren, 1996; Van den Broeke, 1996). The melt rate therefore increases strongly with the absorption of shortwave radiation. Secondly, and the other way around, the albedo often decreases when the melt rate increases. Due to larger melt rates, ice is exposed earlier at the surface, metamorphosis of snow occurs faster and more meltwater accumulates on the ice (Zuo and Oerlemans, 1996; Greuell, 2000). In each of these three cases, the albedo is lowered. The facts that a lower albedo causes more melt and more melt causes a lower albedo constitute a positive feedback loop, which enhances initial perturbations. Note that such perturbations can be due to temperature variations. Hence, the method somehow accounts for the effect of interannual variations in atmospheric temperature on the mass balance. The final reason for the success of the method is that accumulation affects the albedo. Less accumulation results in a lower time-averaged albedo because the transition from snow (higher albedo) to ice (lower albedo) occurs earlier in the season.

In this paper, we test the method introduced by de Ruyter de Wildt and others (2002) for the K-transect. We introduce some modifications relative to their study. One modification is the conversion of the absorbed radiation as calculated from the satellite data into the satellite-derived mass balance (b _sat). In the following sections, we first introduce and discuss the equations used to calculate the satellite-derived mass balance. Then we take a step backwards and discuss the method used to retrieve the surface albedo from the satellite data. Next, we show how the satellite-derived albedo varies in space and time and give an example of its conversion into the satellite-derived mass balance. After that, we describe the mass-balance measurements and relate the two types of mass balance. The final section contains conclusions and discussion.

2. Calculation of the Satellite-Derived Mass Balance

The easiest way to understand how we calculated the satellite-derived mass balance (b_sat) is to consider first the calculation of the absorbed shortwave radiation averaged over the transect and summed over the ablation season, 〈Qabs〉:

(1)

In this equation, a_sat is the satellite-derived albedo and I ₀ the incoming shortwave radiative flux at the top of the atmosphere. The factor τ_atm accounts for attenuation of the incoming radiation by the atmosphere. Therefore, the term in square brackets is an estimate of the shortwave radiation absorbed at any location of the surface and at the moment of a satellite pass. The inner integration then averages the absorbed radiation along the transect (coordinate x; length L; see section 5). The outer integration results in the sum of the radiation absorbed between two dates (j1 and j2). For later comparison with the mass-balance measurements, we integrated over a balance year, that is, between dates of two subsequent mass-balance measurements. We assumed no contribution to Q _abs (and therefore to b _sat) for the period 11 September to 19 May of the following year. Measurements with sonic rangers along the transect at site 6 (1015 m a.s.l.) suggest that, on average over 5 years, melt stops and starts, respectively, on these dates.

We used a time-step of 1 day for the calculations. Therefore, I ₀ is constant along the east-west-extending transect and can be put between the two integration symbols. We assumed that τ_atm can be replaced by its mean value for the ablation season and that it does not vary from year to year. This assumption implies that the influence of clouds on the incoming shortwave radiation is invariant from year to year. By taking τ_atm constant, it could also be put between the two integration symbols. Hence, we arrived at:

(2)

with

(3)

We set τ_atm = 0.62, in accordance with 5 years of measurements of the incoming shortwave radiative flux at site 6. The extraterrestrial flux I ₀ can be calculated very accurately with standard equations (Garnier and Ohmura, 1968). This leaves as the only unknown. It was retrieved from Advanced Very High Resolution Radiometer (AVHRR) images, but this could only be done for clear-sky days. We converted the values for days with suitable images (N days with albedo into a continuous time series by using a Gaussian filter (Fig. 3):

(4)

with the weights w_ij computed according to:

(5)

In this equation, di and dj are day numbers and is a timescale (4days).

Finally, b_sat was computed by two operations:

(6)

Fig. 3. Satellite-derived surface albedo, averaged over the K-transect, and extraterrestrial radiation and shortwave radiation absorbed at the surface during the 1995 ablation season. Dots give all the satellite-derived albedos. On the basis of these points, we compute the temporal variation in the albedo (upper continuous line) by using the Gaussian filter of Equations (4) and (5). We combine this albedo with the incoming radiation at the top of the atmosphere to calculate the amount of shortwave radiation absorbed at the surface (lower continuous line). Finally, b _sat is computed by integrating the absorbed shortwave radiation curve minus the threshold over time, omitting periods when the absorbed shortwave radiation is less than the threshold. Triangles show the shortwave radiation absorbed at the surface on days with suitable images (and correspond therefore to the dots). De Ruyter de Wildt and others (2002) fitted a three-parameter Gaussian function (dashed line) to the triangles and then integrated this function over time to determine the total amount of absorbed radiation.

First, by adding a term representing the sum of the longwave radiative and turbulent fluxes (Q0), we converted the estimated absorbed shortwave radiation for each day into the total energy flux from the atmosphere. Melting only occurs when the total energy flux is positive, so we did not allow negative values. We set Q ₀ = −48 W m–² , in accordance with measurements and calculations (Greuell and Konzelmann, 1994) for ETH Camp (1155 m a.s.l.; about 300 km north of the K-transect). We show later that, by introducing the parameters τ_atm and Q0, the mean and the standard deviation of the satellite-derived mass balance becomes similar to the mean and the standard deviation of the measured mass balance. Secondly, we divided the total energy flux by the latent heat of fusion Lf (0.334 × 10⁶ J kg^-1). Note that Equation (6) ignores aspects of the surface mass balance such as accumulation and refreezing of meltwater. Implications of the neglect of these processes are discussed in section 7.

Our method incorporates some modifications with respect to the work of de Ruyter de Wildt and others (2002). The first modification is the conversion of energy into mass. This has the advantage that the result of the calculations can be directly compared with the mass-balance measurements. De Ruyter de Wildt and others (2002) also did not use the parameters τ_atm and Q₀. Another modification is our method for calculating the total amount of absorbed radiation from values on the days of the clear-sky satellite images (Fig. 3). Whereas we perform an interpolation and extrapolation to values of the surface albedo, then multiply by I ₀ and finally integrate over time, de Ruyter de Wildt and others (2002) first multiplied by I ₀ for the days of the clear-sky satellite images, then fitted a three-parameter Gaussian function to the values of the clear-sky days and finally integrated the Gaussian function over time. The disadvantage of the method of de Ruyter de Wildt and others (2002) is that it puts constraints on the form of the function describing absorbed radiation as a function of time (e.g. that the function is symmetric with respect to a specific date (which is found by fitting the function to the available data)). In our method, such constraints do not exist.

3. Retrieval of the Surface Albedo from the Avhrr Data

The method we used to retrieve the surface albedo from the satellite data strongly resembles the methods employed by de Ruyter de Wildt and others (2002), Klok and others (2003) and Greuell and Oerlemans (2005). We will therefore only outline the retrieval method, but we will discuss deviations from the methods described in the above-mentioned literature in more detail. We performed the following processing steps:

1. 1. Selection of potentially cloud-free images was accomplished by browsing through the archive of the Dundee Satellite Receiving Station (http://www.sat.dundee.ac.uk). We ordered 388 images for the years of the mass-balance measurements (1990–2003), which corresponds to an average of 28 images per year. With a few exceptions, a maximum of one image per day was taken. Since the accuracy of the retrieval deteriorates with increasing satellite zenith angle (Greuell and Oerlemans, 2005), we selected the daytime image with the smallest satellite zenith angle for the K-transect. All our images were acquired between 1402 and 1842 GMT (solar noon is at 1525 GMT). They had a pixel size of 1 by 1 km and covered the area 64.5–70.5˚ N, 40–56˚W.

2. 2. Correction of the geolocation was carried out with a ground-control point at the intersection of the K-transect and the ice-sheet margin (the location of the former Greenland Ice Margin Experiment (GIMEX) base camp). Precision is one to two pixels.

3. 3. Three cloud-masking procedures that are applied to assemble the AVHRR Polar Pathfinder (APP) dataset have been discussed and tested for the Greenland ice sheet by Stroeve (2001). She showed that two of these methods lead to a severe overestimation of the frequency of cloudy cases. The third method (‘cmask1’) yields almost the correct frequency of clear and cloudy cases at lower elevations, but Stroeve (2001) does not show how successful the method is for individual images. Therefore, and since cmask1 relies on continuous sequences of daily images, which we did not have, we employed a new algorithm for cloud detection. Because we were interested in the K-transect only, we first superimposed a window of 12 by 115 km over the transect, extending from the latitude of the northernmost mass-balance site (site 5 at 67.10˚N) to the latitude of the southernmost mass-balance site (site 7 at 66.99˚ N) and from the longitude of GIMEX base camp (50.25˚ W) to the longitude of the point midway between mass-balance sites 9 and 10 (47.61 ˚ W). This is the window for which we also computed the mean surface albedo and therefore b_sat. For cloud detection in this window, we calculated the brightness temperature in AVHRR band 4 (10.3-11.3 μm) of the pixels contained in the window. The brightness temperature approximates the surface temperature or, if clouds are present, the cloud top temperature. We then made scatter plots of the brightness temperature against elevation, which was taken from the digital elevation model (DEM) produced by Bamber and others (2001). A unique relation between brightness temperature and elevation and a monotonic decrease of temperature with elevation can be expected for clear skies all over the transect. Figure 4a and b correspond, respectively, to a melting and a freezing surface under clear skies along the entire transect and do indeed show such a relationship. Deviations from this behaviour, such as are visible in Figure 4c, must be caused by clouds. In order to remove images with clouds from the batch of images, we fitted a second-order polynomial to the points in the scatter plot belonging to each image (see Fig. 4) and computed the residual standard deviation. Images with a residual standard deviation >0.5˚C were excluded entirely from further analysis (25% of the images).

4. 4. Conversion of the AVHRR counts in band 1 (0.58-0.68μm) and band 2 (0.725-1.10μm) into effective narrowband reflectances was performed with calibration coefficients of Rao and Chen (1995) for NOAA-11, Rao and Chen (1999) for NOAA-14 and from the Internet (http://noaasis.noaa.gov/NOAASIS/ml/n16calup.html) for NOAA-16. The coefficients for the sensors on board NOAA-11 and NOAA-14 are time-dependent in order to account for degradation of the instruments. The coefficients for NOAA-16 are constant with time.

5. 5. Conversion of the effective narrowband reflectances in bands 1 and 2 into planetary narrowband reflectances was achieved by accounting for the effects of the solar zenith angle (with respect to a horizontal surface) and the varying distance between the sun and the Earth.

6. 6. Corrections for the atmosphere and for anisotropic reflection by the surface were made simultaneously to account for the effect of the atmosphere on the directional distribution of the signal (Rahman, 1996). The radiative transfer model we used to calculate the atmospheric correction (6S (Second Simulation of the Satellite Signal in the Solar Spectrum); see Vermote and others, 1997) couples a prescribed bidirectional reflectance distribution function (BRDF, a function describing the directional distribution of the reflected radiation) of the surface with the atmosphere.

   With 6S we produced look-up tables (LUTs) for the two AVHRR bands and three different BRDFs, namely one for melting glacier ice (Greuell and de Ruyter de Wildt, 1999), one for melting snow (Koks, 2001) and one for isotropic reflection. Each LUT contains a five-dimensional matrix of the surface narrowband albedo for various values of the planetary narrowband reflectance, for various values of the elevation and for various values of the three angles describing the sun–satellite– target geometry. To obtain values of the surface albedo for actual cases, we applied an equation of the type proposed by Klok and others (2003) for interpolation in terms of the planetary albedo and elevation, and we interpolated linearly between angles describing the sun– satellite–target geometry.

   • Atmospheric composition was as follows:

   • standard sub-arctic summer profiles of pressure, temperature, water vapour and ozone (McClatchey and others, 1972);

   • a mixture of 15% (by volume) dust-like aerosols, 70% (by volume) water-soluble aerosols and 15% (by volume) oceanic aerosols;

   • an optical thickness for aerosols equal to 0.065 at 550 nm (corresponding to a visibility of 109 km) as measured by Stroeve and others (1997) at ETH Camp.

7. 7. Correction for the effect of the sloping surface on the incoming radiation was carried out as by Klok and others (2003). Slope and orientation of the surface were taken from the DEM published by Bamber and others (2001).

8. 8. Narrowband-to-broadband (NTB) conversion was performed according to equation 13 of Greuell and Oerlemans (2004).

9. 9. Steps 1–8 resulted in three albedos for each pixel, namely one for each type of BRDF. The albedo for isotropic reflection was only used in step 10 and therefore not for the K-transect. According to helicopter measurements carried out by Greuell and others (2002), an albedo of 0.60 separates snow from ice for the K-transect. So, if both ice BRDF albedo and snow BRDF albedo were larger than this threshold, the snow BRDF albedo was taken. If both albedos were smaller than the threshold, the ice BRDF albedo was taken, and if the two values were on either side of the threshold, the two albedos were averaged.

10. 10. We recalibrated the satellite sensors. Many of the processing steps 1–9 introduce errors in the final result, the surface albedo. Though some of these errors possibly will not cancel out when averaged over the transect and summed over the season, they are not expected to vary from year to year. Examples of such errors are errors in cloud detection and errors in NTB conversion. In that case, the correlations that our analysis produces would not be affected. This is definitely not true if the assumed temporal variations in the calibration coefficients of the satellite sensors are incorrect.

To check whether such problems with the calibration coefficients occurred, we produced time series of satellite-derived narrowband albedos for an area with constant surface albedo (Fig. 5a and b). We chose a rectangle near the top of the ice sheet (250 by 75 km) with elevations ranging between 2600 and 2970 m. According to the maps by Abdalati and Steffen (1997), melt is very rare within that rectangle, so the surface albedo should be fairly constant in space and in time. In order to detect surfaces within the rectangle that were possibly melting at or before the time of image acquisition, we calculated the spatial variation of the albedo of the pixels from the rectangle. Melting is expected to enhance the spatial variation in the albedo because melt rate and duration will be a function of elevation. Images with a standard deviation of the albedo >0.02 were removed, leaving a total of 131 clear-sky, dry-snow images (on average, nine per year).

For calculation of the albedo within the dry-snow rectangle we used almost the same retrieval method as for the K-transect. The only difference was the assumption of isotropic reflection. This choice was based on unpublished measurements of dry-snow BRDFs (personal communication from E. Hendriks, 2004), which, for the relevant sun–satellite–target geometries, are closer to isotropic reflection than to the melting-snow BRDFs that we employed for the K-transect.

Figure 5a shows that the band albedos vary considerably with time. In the averages for each year, the variations are smooth for periods during which the images were collected by the same NOAA platform. Abrupt changes occur at times of changes in platform and therefore in sensors. We checked the possibility that the interannual variations displayed in Figure 5a were related to slow variations in sun–satellite–target geometry. Such a relationship could have been caused by errors in atmospheric corrections and/or by the unaccounted anisotropic reflection by the surface. However, we found no relation between any angle describing the sun–satellite–target geometry and the apparent albedo variations. We conclude that the apparent variations in surface albedo are due to unaccounted variations in the calibration coefficients of the satellite sensors.

On the basis of the band albedos (shown in Fig. 5a for band 1), we recalibrated the satellite sensors. According to the measurements by Grenfell and others (1994) at the South Pole, the clear-sky albedo of dry snow is about 0.96 in AVHRR band 1 and about 0.80 in the broadband for the solar zenith angles of most of our images. Because it is difficult to derive a value for AVHRR band 2 from the values of Grenfell and others (1994), we determined that value with the two-stream model of Reijmer and others (2001). We found that an albedo in band 2 of 0.88 corresponds to the mentioned band 1 and broadband albedo. We then multiplied the satellite calibration coefficients by factors that only varied by year and band number, so that for each year the average surface narrowband albedos in the dry-snow rectangles became equal to prescribed values of 0.96 (band 1) and 0.88 (band 2). The time series of the satellite-derived narrowband albedos for the dry-snow rectangle after recalibration are displayed in Figure 5b. Much of the residual scatter can be ascribed to the unaccounted anisotropy of the reflection and to natural variations in the surface albedo. The multiplication factors needed to achieve the corrections are listed in Table 1. They range between 0.981 and 1.097 (mean: 1.052) for band 1 and between 0.960 and 1.168 (mean: 1.032) for band 2.

Fig. 4. Brightness temperatures calculated from AVHRR band 4 (10.3–11.3 μm) counts as a function of surface elevation. All pixels from the K-transect window with an elevation >800m are plotted for three images, corresponding to a clear sky and a melting surface (a), a clear sky and a freezing surface (b) and a partly cloudy sky (c), respectively. Best fits of second-order polynomials are plotted as dashed lines, but are not visible in (a) and (b). Residual standard deviations are 0.11˚C (a), 0.33˚C (b) and 2.96˚C (c).

Fig. 5. Retrieved surface albedo in AVHRR band 1 for a ‘dry-snow rectangle’ near the top of the Greenland ice sheet before (a) and after (b) recalibration of the satellite sensors. Only albedos for cloud-free and dry-snow conditions within the rectangle (according to our algorithms) are considered. Each circle corresponds to a single image. Yearly averages are given by triangles. The dashed line shows the assumed mean clear-sky albedo for the rectangle.

Table 1. Correction factors for the calibration coefficients of AVHRR sensors, 1990–2003. We multiplied calibration coefficients from the literature and the Internet (see text) by these factors. They were obtained by retrieving the surface albedo for a rectangle near the top of the Greenland ice sheet for which albedo variations are assumed to be small

Figure 6 compares our satellite-derived albedos with in situ measurements made by a weather station at site 6 (1015m a.s.l.). The in situ measurements cover six ablation seasons and were made with Kipp & Zonen CNR 1 sensors, which have an appropriate spectral response (1 between 305 and 2800 nm; 0 for other wavelengths). The main source of error in individual measurements probably is tilt of the sensors relative to the underlying surface, but on the timescale of 6years we expect this error to be random. For validating the satellite-derived albedos, we selected mean in situ measured albedos for hours embracing the time of image acquisition. As satellite-derived albedo we took the value of the pixel that contained site 6. Taking satellite-derived albedos averaged over squares of 5 by 5 pixels surrounding site 6 instead of taking the value for a single pixel had a negligible influence on the comparison.

Fig. 6. AVHRR-derived surface albedos as a function of the measured surface albedo for site 6 (1015 m a.s.l.), situated within the K-transect window. Each point represents the value for a single clear-sky image.

A major problem of this type of validation is a scale problem as noted earlier by Hall and others (1990) and Knap and others (1999). The in situ measurements made at a height of about 4 m above the surface represent an area with a diameter of roughly 25 m. This is two orders of magnitude smaller than the resolution of AVHRR data (1.1km at nadir). Most probably, this difference in scale explains outliers like the low (0.2–0.4) measured albedos corresponding with high (0.6–0.8) satellite-derived albedos. While the snow-pack under the ground-based sensors had melted away and the exposed ice was perhaps partly covered by water, most of the surrounding terrain was at the same time still covered by snow. Despite these outliers, for most points (71%) the absolute difference between calculated and measured albedo is <0.1. Disregarding the outliers, most satellite-derived albedos <0.5 are somewhat too small, whereas satellite-derived albedos >0.5 are somewhat too high. Following the analysis by Greuell and Oerlemans (2005) we might blame such a systematic difference mainly on insufficient knowledge of BRDFs (our BRDFs are based on measurements made in the Alps) and on errors in the (revised) calibration coefficients of the satellite sensors.

4. Satellite-Derived Albedo and Satellite-Derived Mass Balance

Figure 7 shows selected albedo profiles for 1995, and Figure 3 the corresponding averages for the K-transect window. Most of the features visible in these figures are as expected. The albedo tends to increase with distance from the margin (and therefore with elevation). In addition, it tends to decrease with time at all locations until the beginning of August. This decrease becomes faster once the slush line has passed a location. With the passing of the slush line, a uniform snow cover is succeeded by a mixture of snow, slush and ice, and later by ice alone (Greuell and Knap, 2000). It is unclear whether smaller temporal variations, like the increase in the albedo along the entire transect between 21 and 31 July, are real or artifacts due to inadequacies in the retrieval method. Between 7 and 18 August, the albedo increased abruptly as a result of snowfall. After 18 August, it decreased slightly again, but never reached the low values of July. A noteworthy feature is the minimum in the albedo at about 50 km from the margin. This can be ascribed to accumulation of meltwater on the ice surface (Greuell, 2000).

Fig. 7. Computed albedo profiles along the K-transect for nine days in 1995 (dates are dd-m-yy). The arrows point towards the positions of the slush line on different days according to the algorithm of Greuell and Knap (2000). This position could not be determined for 18 August and 8 September. Only the part of the transect left of the dashed line is considered in the computation of b _sat.

Figure 3 also demonstrates how the mass balance is estimated from the satellite-derived albedos. First, albedo values for each day of the ablation season were computed from the 21 values available for 1995 by using the Gaussian filter of Equations (4) and (5). Next, we inserted the daily values of the satellite-derived albedos into Equation (6) to calculate the satellite-derived mass balance for each balance year.

5. Mass-Balance Measurements

Surface mass-balance measurements along the transect (Greuell and others, 2001; Van de Wal and others, 2005) have been carried out since 1990 by members of the Institute for Marine and Atmospheric Research Utrecht (IMAU). This record therefore currently consists of 13 balance years (1990/91–2002/03). The K-transect extends from west to east at 67˚ N. At each of seven sites, ranging in elevation between 340 and 1520 m, height changes along at least (and usually) two stakes were measured. In 1993/94, an eighth mass-balance site, located in the accumulation area at 1850 ma.s.l., was added (site 10). The measurements were performed between 18 and 26 August each year, with the exception of 1990 (21 July), 1991 (5 August) and 1992 (12 August). This makes the duration of the first three balance years considerably longer than a year, but these differences are accounted for in the satellite-derived mass balance by adapting the integration limits j ₁ and j ₂ in accordance with the measurement dates.

Per site and per balance year the mass balance was calculated by averaging over the available stakes. We estimated the error in the mass balance for individual sites on the basis of differences between stakes at the same site (see Braithwaite and others, 1998). In the case of two stakes, we took half of the difference between the two values. We then rounded upwards in order to account for unmeasured spatial variations on a scale larger than the distance between the stakes. We found that the error in the mass balance for individual sites increases from about 100 mm w.e. at site 9 (near the mean equilibrium-line altitude (ELA)) to about 300 mmw.e. at site 4 (near the ice-sheet margin). Mass-balance profiles for several years and the average for all years are plotted in Figure 8. On average, the ELA lies at about 1500m a.s.l. Therefore, the part of the K-transect considered here (up to halfway site 9 and site 10) belongs almost entirely to the ablation facies. Only the surface at the highest considered elevations might belong to the soaked facies (Benson, 1961). The interannual variation decreases with elevation and is dominated by year-to-year variations in ablation, which is caused by the dryness of the area, with accumulation rates of about 200 mm a^–1 (Ohmura and others, 1999).

Fig. 8. Some of the mass-balance profiles measured along the K-transect. The balance years 1991/92 and 2002/03 are characterized by extremes in the net mass balance.

Though we did test the satellite method for individual sites, the main test was for the net (transect mean) mass balance. Net mass balance was computed by averaging over the individual sites, weighting each site with the distance along the transect between the two points halfway towards the neighbouring sites. The highest site was excluded because at that site the measurements started 3years later than at the other sites. In view of the error in the mass balance for individual sites, the error in the annual net balance is about 80 mmw.e.

6. Relation Between Satellite-Derived and Measured Mass Balance

Figure 9 compares measured net mass balance (b _m) and satellite-derived mass balance (b _sat). The optimal linear fit is (the estimated mass balance b _m,est and b _sat in mmw.e.):

(7)

b _sat explains 7 1% of the variations in b _m (correlation coefficient = 0.84; residual standard deviation in b _m = 214 mmw.e.). The scatter plot confirms an almost linear relation between b _m and b _sat, which is not necessarily expected in view of the complex relationships between the albedo, net shortwave radiation, melt, accumulation and mass balance.

Fig. 9. Scatter plot of the transect mean measured mass balance against the transect mean satellite-derived mass balance. Individual points are labelled by corresponding balance years. The line represents the best fit according to Equation (7).

In view of the simplicity of the method, it is remarkable that the magnitude of b _sat (average for all years: –1142mm w.e.) is similar to the magnitude of b _m (average for all years: –1202mmw.e.). The same is true of the interannual variability in b _sat and b _m. Standard deviations over all 13 years are 366 and 395 mmw.e., respectively.

We also tested the method for individual sites, by taking the mean satellite-derived albedo for squares of 5 by 5 pixels surrounding the sites to calculate the satellite-derived mass balance. The lowest two sites, located within a few kilometres of the margin, were excluded from the analysis because the tundra often influences the satellite-derived albedos for these sites. Figure 10 shows the correlation coefficient between b _m and b _sat, and therefore the method’s performance as a function of elevation. Correlation coefficients are significant at the 99% level at all elevations, except for the elevation of the highest site (significance <95%). The performance is optimal just below the mean ELA.

Fig. 10. The method’s performance characterized by the correlation coefficient between measured and satellite-derived mass balance as a function of the surface elevation. Correlation coefficients are significant at the 99% level when they are higher than the values given by the dashed lines. This limit varies with elevation because the number of balance years varies with elevation.

7. Conclusions and Discussion

We have demonstrated that the surface mass balance along the K-transect can be estimated on the basis of AVHRR-derived surface albedos. It appears that the satellite-derived mass balance (b_sat) explains 7 1% of the variations in the observed annual mass balance.

Excluding two of the thirteen available balance years, b _sat would have explained a considerably higher percentage of the variance (90%). For the balance year 1991/92, b _sat is too small by 444 mmw.e. with respect to the best fit of Equation (7). Perhaps the discrepancy is partly due to high amounts of accumulation in 1991 and 1992 (Mote, 2003), but only exceptional amounts could explain the entire discrepancy, which is more than twice as large as the mean annual precipitation (200 mm). We also checked whether the deviation could have been caused by the eruption of Pinatubo, Philippines, on 15 June 1991. During this eruption, huge amounts of aerosols were ejected into the stratosphere, which reduced the incoming radiation in the subsequent years. We did not take this effect into account. According to model simulations by Stenchikov and others (1998), the reduction in the incoming radiative flux (shortwave plus longwave) at the location of the K-transect during summer 1992 was about 1Wm^–2 . If this change in the energy flux is entirely translated into a change in the amount of melt, the total amount of melt decreases by only about 30 mmw.e. We conclude that the outlier of 1991/92 is an order of magnitude too large to be ascribed to the direct radiative effect of the aerosols ejected by Pinatubo.

For the balance year 1996/97, b _sat is too large by 569 mmw.e. with respect to the best fit of Equation (7). A small part of this deviation (about 160 mmw.e.) could be due to the anomalous mass-balance measurement for that balance year at site 6 (Fig. 8), which is perhaps related to the formation of meltwater streams around stakes.

Both the images for 1992 and 1997 are close to the averages for the whole period in terms of the angles describing the sun-satellite-target geometry and in terms of the acquisition time. Extremes in the geometry could have led to extreme errors in b_sat, which are linked to errors in, for example, the BRDFs. Extremes in acquisition time could have caused deviations in b_sat linked to diurnal variations in the surface albedo, as observed for the K-transect by Greuell and others (2002). In summary, we found no satisfactory explanation for the two outlying balance years, except perhaps for a significant contribution by large amounts of precipitation during 1991/92.

By introducing r_atm and Q0 we turned the equations of de Ruyter de Wildt and others (2002) into a simple equation that sums all the energy fluxes between the atmosphere and the surface. Note at this point that r_atm and Q₀ were not used as free parameters to be determined by statistics but were derived from independent data. Another adaptation of the equations of de Ruyter de Wildt and others (2002) was division of the total energy flux by the latent heat of fusion so that energy is converted into mass. These two adaptations enabled the direct comparison of our product (b _sat) with the measured mass balance (b_m). We found a close agreement between both the mean values and between the interannual variabilities of b _sat and b_m.

We have performed a number of tests to assess the sensitivity of the method’s performance for uncertainties in parameters. The upper part of Table 2 lists the results of these tests in terms of changes in the mean and in the standard deviation of b_sat and in terms of the explained variance. Bold numbers correspond to large changes in the performance. It appeared that reasonable changes in the prescribed beginning and end of the ablation period, in the width of the transect window and in the dry-snow albedos have a small effect on all results. Moreover, the explained variance is hardly affected by any of the parameter changes. However, the mean value of b_sat is very sensitive to variations in τ_atm and Q0. Therefore, the good agreement that we found between the mean values of b_sat and b_m seems fortuitous and we conclude that the method is not suitable for assessing the absolute magnitude of the surface mass balance. This is even more so in face of the fact that Equation (6) neglects crucial aspects of the energy and mass balance, such as accumulation, removal of winter and nocturnal cold waves by conduction and refreezing of meltwater, and diurnal variations.

Table 2. Sensitivity of the method’s performance for various changes in parameters and processing steps. Changes in mean b _sat, in standard deviation of b _sat and in explained variance are relative to the values presented in section 6 (–1142mmw.e., 366 mmw.e. and 7 1 % , respectively). Bold values represent large changes in performance

Table 2 also shows that the calculated interannual variability in b_sat, quantified by the standard deviation, is rather robust against the uncertainty in the various parameters. The only exception is the sensitivity to changes in τ_atm, but perhaps this parameter can be estimated with sufficient accuracy from the percentage of cloud-free images. We conclude that our method is suitable for estimating the interannual variability of the surface mass balance. One might argue that this conclusion is dangerous because accumulation and refreezing of meltwater are not considered in Equation (6) and therefore not considered in the sensitivity analysis. However, in dry areas like the K-transect, accumulation and refreezing are small. Consequently, the interannual variations in these processes form only a small contribution to the interannual variations in the surface mass balance. Interannual variations in ablation are dominant (Greuell and others, 2001). Hence, we expect that the neglect of accumulation and refreezing hardly affects the computed interannual variability in b_sat. The same is true for the explained variance.

We expect that our method is also suitable for estimating interannual variations in the surface mass balance of other glaciers and parts of ice sheets and ice caps. In theory, application of the method in regions that are more maritime than West Greenland could be problematic. In wetter regions, fewer cloud-free satellite images are available, and annual variations in accumulation are more important than in dryer regions. However, de Ruyter de Wildt and others (2002) successfully applied the method to the maritime ice cap Vatnajökull.

We also tested the effect of several processing steps on the results (bottom half of Table 2). It appears that we would have obtained similar results if we had omitted cloud detection (on the pre-selected images!) or if we had used linear interpolation to obtain values of for days without suitable images. However, the explained variance would have been much lower without applying the BRDF correction (by 6%) and without recalibration of the satellite sensors (by 9%). These large sensitivities indirectly indicate (but do not prove) that for the K-transect the BRDFs that we used are closer to reality than the assumption of isotropic reflection and that the derived corrections improve the accuracy of the calibration coefficients of the satellite sensors. Those corrections reach maximums of 10% for band 1 and 17% for band 2. Such an order of magnitude is in line with differences between calibration coefficients proposed by various authors (see Tahnk and Coakley, 2001). In future, recalibration can be improved when more knowledge on the albedos and BRDFs of dry snow in the AVHRR bands becomes available. Alternatively albedos could be derived from data provided by Moderate Resolution Imaging Spectroradiometer (MODIS) or Along Track Scanning Radiometer (ATSR) sensors, which are both calibrated with on-board devices and consequently deliver data with small uncertainty (better than 2% for MODIS according to Justice and others, 2002; 3% for ATSR according to http://www.atsr.rl.ac.uk/documentation/docs/visiblecalibration/index.shtml). A disadvantage of those data is that collection in suitable bands (at least one visible channel is needed) started only in 2000 (MODIS) and 1995 (ATSR).

As Figure 10 demonstrates, the method works best just below the mean ELA. We explain this by the duration of the season with ice, and therefore snow, at the surface. This quantity probably has a maximum in its interannual variation at the elevation for which the method performs best. Such interannual variations have a large effect on b _sat because they cause relatively large variations in the albedo due to the difference between snow and ice albedo. Interannual variations in the duration of the season with ice (snow) at the surface are smaller at the lowest and highest elevations of the transect. At the lowest elevations, ice is exposed at the surface even at the beginning and the end of many ablation seasons, as suggested by measurements with sonic rangers at site 5 (490ma.s.l.). Apparently, strong winds blow snow off the ice. At the highest site (site 10, 1850 ma.s.l.), snow always covers the surface and melt rates are low, so interannual variations in albedo are relatively small. Therefore, and because the interannual variability in the mass balance is smaller than elsewhere along the transect, the method performs poorly for site 10. Given this variation with elevation of the method’s performance and its explanation, we expect that the method will perform better for glaciers where substantial ablation occurs even at the highest considered elevations and substantial accumulation occurs even at the lowest considered elevations. This might be why de Ruyter de Wildt and others (2002), using the method of this paper, could explain much more of the variance in the observed annual mass balance of six different outlets of Vatnajökull (76–88%) than we could explain for the K-transect (71%).

Finally, we want to comment on other methods that can be used to estimate the surface mass balance. Greuell and Knap (2000) have studied how the surface mass balance of the K-transect could be estimated by inferring the annual maximum elevation of the slush line from AVHRR images. Despite the fundamental problem of that method, discussed in the introduction of the present paper, the ‘slush line method’ could explain the same percentage in the observed variance as our ‘albedo method’. Greuell and Knap (2000) also studied the potential of the degree-day method to estimate the mass balance. As input they employed temperature data from Kangerlussuaq, situated about 25 km west of the lowest point of the K-transect. They could explain 82% of the observed variance with the degree-day method, so for the K-transect that method seems superior to the two methods using satellite data as input. We speculate that an energy-balance model would perform even better than a degree-day model. Also, by combining two or several of the mentioned approaches, results could be improved. Methods using meteorological data as input lose some of their power when they are applied at greater distances from meteorological stations, whereas the performance of the methods based on satellite data is independent of that distance. At larger distances from meteorological stations, the methods based on satellite data might therefore be more successful than methods based on meteorological data. In fact, the method dealt with in this paper allows estimation of the interannual variation in the surface mass balance along the entire margin of the Greenland ice sheet in a spatially continuous way.