Introduction
Compared with traditional approaches, satellites enable us to observe more glaciers, especially those that are difficult to reach, more frequently. Different remote-sensing methods can be used to measure glacier extent, height, velocity and surface mass balance (Konig and others, 2001; Reference De Ruyter de Wildt, Oerlemans and Björnssonde Ruyter de Wildt and others, 2002).
In this paper, we discuss surface mass balance. Traditional approaches require many stakes to obtain a good spatial distribution of accumulation and melt (Reference PeltoPelto, 2000). Satellites, however, have the potential to deliver mass balance at a higher spatial resolution because they can be used to provide representative images of surface albedo (Reference Reijmer, Knap and OerlemansReijmer and others, 1999). How mass balance and surface albedo are related is discussed below.
Before mass balance can be inferred purely from satellite images, the procedures need to be tested thoroughly. Many authors have already worked on this subject. Some have tried to relate the position of the snow-line at the end of the ablation season to the mass balance because the snow-line is a surrogate for the equilibrium line, which can be used to infer the mass balance. For example, Reference Hall, Williams, Barton, Sigurðsson, Smith and GarvinHall and others (2000) and Reference De Ruyter de Wildt and Oerlemansde Ruyter de Wildt and Oerlemans (2003) used synthetic aperture radar (SAR) to look at the snow-line and equilibrium-line positions in Iceland. They found that the equilibrium line could be obscured by the firn line of the previous year if that firn line lies below the snow-line. Reference Hall, Williams, Barton, Sigurðsson, Smith and GarvinHall and others (2000) also argued that Landsat images were preferable to SAR images, because they better distinguish the glacier from the surrounding moraine. Another disadvantage of the SAR method is that it requires a good image at a specific moment, something that can be quite difficult in Iceland because of its frequent cloudiness.
An alternative method to using the satellite-derived snowline position is to use satellite-derived surface albedo or potential absorbed radiation at the surface, which is derived from albedo. Among the first authors to use Advanced Very High Resolution Radiometer (AVHRR) images to derive the albedo of a glacier were Reference Koelemeijer, Oerlemans and TjemkesKoelemeijer and others (1993). Subsequent authors improved the technique by correcting for surface slope (Reference Knap, Brock, Oerlemans and WillisKnap and others, 1999; Reference De Ruyter de Wildt, Oerlemans and Björnssonde Ruyter de Wildt and others, 2002). Recent papers have improved the technique further by adding an atmospheric correction, and a correction for the anisotropic reflection of radiation by either ice (Reference Stroeve, Nolin and SteffenStroeve and others, 1997) or snow (Reference De Ruyter de Wildt, Oerlemans and Björnssonde Ruyter de Wildt and others 2002) or both (Reference Klok, Greuell and OerlemansKlok and others, 2003).
In this paper, we present a dataset of calculated radiation from 1991 to 2002, which we compare with in situ measurements of surface mass balance. The images are from Vatnajökull, Iceland (Fig. 1), one of the largest ice caps in Europe. They were taken by the AVHRR instruments aboard US National Oceanic and Atmospheric Administration (NOAA) satellites, specifically NOAA-11, NOAA-14 and NOAA-16. The AVHRR instrument has a rather coarse resolution (1.1 km at nadir) compared to some other satellite instruments (e.g. those on the Landsat satellites (30 m at nadir)), but it is adequate for an ice cap as large as Vatnajökull (8300 km2). An advantage of AVHRR is that the satellite passes over the area several times a day. A negative aspect of Vatnajökull is that it lies in an active volcanic area, so its surface can be affected by factors unrelated to incoming radiation, like volcanic ash, which can change the albedo.
Fig. 1. Location of Iceland (inset) and Vatnajökull. Several drainage basins are marked: T, Tungnaarjökull; K, Köldukvíslarjökull; D, Dyngjujökull; BU, Bruarjokull; BE, Breiaamerkurjokull; E, Eyjabakkajokull.
Images from 1991 to 1999 have already been presented by Reference De Ruyter de Wildt, Oerlemans and Björnssonde Ruyter De Wildt and others (2002). Because many steps in the procedure have been changed in our study, these images were reprocessed to produce a homogeneous AVHRR radiance time series. The method described below requires visual inspection of images at several steps, especially for the detection of clouds, but it is computationally efficient.
Albedo and Mass Balance
Reference De Ruyter de Wildt, Oerlemans and BjörnssonDe Ruyter De Wildt and others (2002) showed that the mean surface mass balance (b) correlates well with the potential absorbed radiation (Q pot,net) averaged over a glacier surface over the melt season, and to a lesser extent with the average surface albedo (a). This is because the amount of shortwave radiation absorbed by the surface greatly influences the amount of summer melting. The winter balance is also taken into account by this method: if more snow has fallen, it will take longer to melt and the surface albedo will remain high for a longer period. To some extent, this method is sensitive to summer snow, which can have a great effect on albedo for a few days, while not being very important for the overall mass balance (Reference Greuell and OerlemansGreuell and Oerlemans, 1986).
Net potential radiation is defined as
Q pot is the amount of radiation that would reach the surface if the atmosphere were completely transparent. It can be calculated from standard astronomical theory (Reference WalravenWalraven, 1978). The mass balance (B) is related to the summer Q pot,net averaged over the whole ice cap (〈Q pot,net〉):
where
where D is the duration of the summer and A is the area of the observed surface.
It is possible that an important part of the ice cap is not visible during part of the year. This will lead to a misinterpretation of the mass balance. Therefore we used a curve to describe 〈Q pot, net〉 with a priori knowledge about its course during the summer. A weight factor is introduced which gives the images with the least clouds the largest weight. Reference De Ruyter de Wildt, Oerlemans and BjörnssonDe Ruyter de Wildt and others (2002) found that the following Gaussian curve describes the changes in 〈Q pot, net〉 well:
where a, b and c are tunable constants and day is the day of the year. Figure 2 shows 〈Q pot, net〉 vs the day of the year for two extreme years on the drainage basin Bruarjokull. Every point in the graph depicts a single image. Before applying Equations (2) and (3) to estimate the mass balance, we describe the method to retrieve the surface albedo.
Fig. 2. Examples of the fit of Equation (3). This figure shows the most extreme years for Brúarjökull, namely 1993 (triangles and dashed line) and 1995 (circles and solid line).
Satellite Image Processing
The AVHRR images were purchased from Dundee Satellite Receiving Station, Scotland. We attempted to acquire a cloud-free image every 2 weeks for April-October 1991-2003. However, it was not possible to obtain completely cloud-free images for every date, so a cloud-masking procedure was necessary. The images for the years 1991-94 were obtained from NOAA-11, those from 1995 to April 2001 were derived from NOAA-14, and all later images were acquired by NOAA-16 (Table 1). Bands 1 and 2 were used to determine the albedo of the surface, while bands 1, 3 and 4 were used for cloud masking.
Table 1. Dates of the satellite images used in the analysis
Cloud masking
It is often difficult to detect clouds over an ice cap, because clouds, ice and snow frequently have similar reflectance characteristics in the visible and infrared wavelengths. Thus, cloud masking becomes a matter of pattern recognition. Reference De Ruyter de Wildt, Oerlemans and BjörnssonDe Ruyter De Wildt and others (2002) found that the differences between bands 3 and 4 could be used to aid in the detection of clouds. A parameter R is defined as:
where c 3 and c 4 are the raw count values per pixel in the AVHRR bands 3 (3a for NOAA-16) and 4. If R is higher than a chosen threshold, the pixel is classified as cloud. Threshold values are not equal for every image, and anomalies like cloud shadows are not detected. However, in addition to this procedure, we also examined a composite image of bands 1 (blue), 3 (red) and 4 (green). By including the visible band it was easier to see the extent of the ice cap and it was also possible to distinguish clouds from patterns in the ice itself and this enabled us to better tune R. Portions of the ice cap that were still suspect after the best possible R had been chosen were removed manually.
Geolocation
We used a digital elevation model to produce a mask of the ice-cap margin. This mask was moved over the image until the fit with the margin was optimized. These fits were checked manually afterwards. We estimate that the accuracy of the geolocation is within one to two pixels.
Calibration
The raw count numbers of the AVHRR images in bands 1 and 2 were converted into top-of-the-atmosphere radiances. Reference Rao and ChenRao and Chen (1995, Reference Rao and Chen1999) and http://noaasis.noaa.gov/NOAASIS/ml/n16calup.html provide calibration coefficients for bands 1 and 2 for the AVHRR instruments aboard NOAA-11, NOAA-14 and NOAA-16 respectively. Instrumental drift is taken into account for NOAA-11 and NOAA-14, while the instruments aboard NOAA-16 have shown no systematic instrumental drift so far. However, the calibration coefficients have been updated several times, so the sensor is not completely stable. The zenith of the sun and the variable distance between the Earth and the Sun are also taken into account.
Atmospheric correction and BRDFs
The incoming radiation at the top of the atmosphere is not equal to the radiation received at the surface, because of the influence of the atmosphere. We used the 6S (Second Simulation of the Satellite Signal in the Solar Spectrum) model to determine this influence (Reference Vermote, Tanré, Deuzé, Herman and MorcetteVermote and others, 1997). The model incorporates descriptions of both the anisotropic reflection of the atmosphere and the surface. The anisotropy of the latter is described using bidirectional reflection distribution functions (BRDFs). The BRDF for ice was derived by Reference Greuell and de Ruyter de WildtGreuell and de Ruyter de Wildt (1999) using measurements made on Morteratschgletscher, Switzerland, and the BRDF for snow was determined by Reference KoksKoks (2001) on the basis of measurements made on Glacier du Géant, Italy. We present the results with and without correction for anisotropy of the surface.
To distinguish between ice and snow we used the broadband albedo threshold of 0.43 determined by Reference Greuell, Reijmer and OerlemansGreuell and others (2002). We selected this value because it was determined on Vatnajökull itself. If the broadband surface albedo after all corrections lay below 0.43, the surface was classified as ice and otherwise as snow. If the images contradicted each other, i.e. the image with the snow BRDF classified the pixel as ice, while the ice BRDF image classified it as snow, the mean of the two broadband albedos was taken. With this method it is not a priori possible to say if a surface is ice or snow, because the broadband albedo is used to determine that. Therefore we had to calculate the albedo for all images with both BRDFs and select the correct surface type of each pixel afterwards.
The relation between the reflectance at the top of the atmosphere and the surface albedo depends on the sun– target–satellite geometry, water vapour, ozone and visibility. The latter is a measure for the aerosol content of the atmosphere. We generated a look-up table (LUT) that describes the relation between all possible zenith and azimuth angles of the sun and the satellite. The calculations were done for a standard sub-arctic summer profile with constant values for water vapour (10.6 kg m–2), total ozone (341 Dobson units) and visibility (109 km), all at 1500 m a.s.l.
To justify the use of constant values, we repeated all processing steps with atmospheric properties specific to each date for the years 2001 and 2002. Water vapour was taken from the ERA-40 re-analysis by the European Centre for Medium-Range Weather Forecasts, ozone content from NASA (http://www.toms.gfc.nasa.gov/ozone/ozone.html) and visibility from Fagurhólmsári station (63˚53’ N, 16˚39’ W) just south of Vatnajökull. This justified the use of the LUT values because the difference between the albedo calculated with both methods was only a few per cent at most.
Slope correction
All calculations above assume a horizontal surface, but this is seldom the case. We used the method of Reference Knap, Brock, Oerlemans and WillisKnap and others (1999, equation 4) to correct for the slope of the surface. The correction factor is defined as:
where 9 s is the solar zenith angle, y is the angle between the normal vector of the surface and the vector pointing to the sun, and a is the ratio of diffuse to direct solar radiation under clear-sky conditions. Its value depends on the wavelength (0.138 and 0.079 for bands 1 and 2, respectively).
Narrow-to-broadband conversion
The calculations described above compute the narrowband albedos in AVHRR bands 1 (α1) and 2 (α2). These narrowband albedos have to be converted to a broadband albedo (a). We used several narrow-to-broadband (NTB) conversion equations and compared the results, namely:
Reference De Ruyter de Wildt, Oerlemans and Björnssonde Ruyter de Wildt and others (2002):
Reference Greuell, Reijmer and OerlemansGreuell and others (2002):
Reference Greuell and OerlemansGreuell and Oerlemans (2004):
Results
This section is divided into two parts. The first deals with the uncertainties in the satellite-derived albedo, and the second deals with the mass balance determined using the calculated albedo. The results are presented in comparison with a standard run. This standard run uses an isotropically reflecting surface and NTB conversion equation (6c). Equation (6c) was chosen because it was the most recent NTB conversion equation available. The average albedo of the whole ice cap is shown for the years 1991–2002 in Figure 3 for the standard run.
Fig. 3. Calculated albedo for the whole ice cap, 1991–2002. The standard run uses an isotropically reflecting surface and NTB conversion equation (6c). This figure clearly shows the interannual variability.
Satellite-derived albedo
The differences between the isotropic and anisotropic cases are typically on the order of a few per cent, as shown in Figure 4, and the correlation between both cases is large (0.98). The anisotropic albedo tends to be higher than the isotropic albedo, but there are exceptions. Not accounting for the anisotropy of snow and ice will, however, lead to errors in the calculated surface albedo. Reference Knap, Brock, Oerlemans and WillisKnap and others (1999) did not correct for anisotropy and they considered this the main source of errors in their method. Reference De Ruyter de Wildt, Oerlemans and BjörnssonDe Ruyter de Wildt and others (2002) corrected for snow, but not for ice, while Reference Klok, Greuell and OerlemansKlok and others (2003) used BRDFs for both snow and ice. However, Reference Greuell and OerlemansGreuell and Oerlemans (2004) argue that it is better to use no BRDF at all than one made for a glacier with completely different surface properties. Reference Reijmer, Bintanja and GreuellReijmer and others (2001) also mention that the BRDFs they derived for Antarctica cannot be used for other places. Albedo variations on Morteratschgletscher are mainly caused by the amount of water at the surface (Reference Klok, Greuell and OerlemansKlok and others, 2003), while the dust content is more important for Vatnajökull (Reference Greuell, Reijmer and OerlemansGreuell and others, 2002). Until these calculations are repeated with BRDFs specifically determined for Iceland, it is impossible for us to choose between the isotropic and the anisotropic case. Unfortunately, such BRDFs do not yet exist.
Fig. 4. Radiation calculated with an isotropically reflecting surface vs radiation calculated with an anisotropically reflecting surface. The correlation between the two cases is high (R = 0.985).
The NTB equations by Reference Greuell and OerlemansGreuell and Oerlemans (2004) (Equation (6c)) and by Reference LiangLiang (2001) (Equation (6d)) are based on modelling and were subsequently compared with measurements, while Equations (6a) and (6b) are solely based on measurements. Equations (6c) and (6d) have been determined for a wide range of possible albedos. This is important, because extrapolation of NTB conversion equations beyond the range of albedos that were used as input for the equations is hazardous. However, by using different NTB equations we could estimate the effect of the NTB step on the calculated albedo. A scatter plot of the standard run albedo vs the other albedos is shown in Figure 5. Equation (6a) is a preliminary form of Equation (6b), so we expected that this version would perform worst. It was only included because it had been used in a previously published paper (Reference De Ruyter de Wildt, Oerlemans and Björnssonde Ruyter de Wildt and others, 2002). However, the standard run albedo is more strongly correlated with Equations (6a) and (6d) than with (6b). Although Equation (6b) is more recent than (6a), it is based mostly on measurements of surfaces with low albedo. Indeed, the largest differences between the Equation (6b) and (6c) results are found among the higher albedos, with differences as large as 13%. The correlation between 〈Q pot,net〉 over the summer for these NTB conversion equations is strong in most cases (around 0.9). This comparison has shown that selection of a NTB conversion equation is no trivial matter, because it can lead to significant variations in the calculated albedo and also in 〈Q pot,net〉, as shown in the next subsection.
Mass balance determined by satellite compared with in situ measurements
The mass balance was measured on the various drainage basins using a stratigraphic method. This method has been described in detail in Reference Björnsson, Pálsson, Guðmundsson and HaraldssonBjörnsson and others (1998c). Briefly, it consists of measuring changes in thickness and density along a number of selected flowlines which span the elevation ranges of the drainage basins. The estimated error of the mass balance per outlet is <20%. The measured mass-balance data were obtained from various reports from the Science Institute, University of Iceland (Reference Björnsson, Pálsson, Guðmundsson and HaraldssonBjörnsson and others, 1997, Reference Björnsson, Pálsson, Guðmundsson and Haraldsson1998a, Reference Björnsson, Pálsson, Guðmundsson and Haraldssonb, Reference Björnsson, Pálsson, Guðmundsson and Haraldssonc, Reference Björnsson, Pálsson, Guðmundsson and Haraldsson1999, Reference Björnsson, Pálsson and Haraldsson2002).
Figure 6 shows the potential absorbed radiation for the standard run vs the mass balance for several years for six drainage basins combined (Tungnaarjökull, Köldukvíslarjökull, Dyngjujökull, Brúarjökull, Breiðamerkurjökull and Eyjabakkajokull). The correlation between measured mass balance and 〈Q pot,net〉 varies strongly between drainage basins (Table 2). Correlation coefficients are low if there are not enough mass-balance measurements and if not enough good images through the year are available. In particular, it was difficult to find good images at the beginning and at the end of the melt season, because clouds are more prevalent in those periods. The residual standard deviation (m w.e.) is low in all cases (Table 2), showing that the data points are strongly related to each other even if the correlation is low. The choice of the NTB conversion equation has a large effect on the correlation coefficient and there is no single NTB conversion equation that performs best for all drainage basins (Table 2).
Fig. 6. Mean summer 〈Q pot,net〉 for the standard run vs mean measured mass balance (B) over the whole area where mass-balance measurements were taken. The largest outlets were given the largest weight in the calculation of the mean. Every point shows a year, and only those years for which mass-balance measurements were available for all outlets are included.
Table 2. Correlation coefficients (r) and residual standard deviations (RSD) in m w.e. for measured mass balance vs mean summer 〈Q pot,net〉 for six drainage basins of Vatnajökull. B is the mean measured mass balance per outlet in m w.e. and n is the number of years for which mass-balance measurements were available. (6a), (6b), (6c) and (6d) refer to the different NTB conversion equations discussed in the text
The correlations we determined are generally lower than those found by Reference De Ruyter de Wildt, Oerlemans and Björnssonde Ruyter de Wildt and others (2002). The addition of the years 2000-02 is partly responsible for this, as 2001 and 2002 decrease the correlation substantially. The images from these years come from the NOAA-16 satellite, while all other images were acquired by NOAA-11 and NOAA-14. Köldukvìıslarjökull gives by far the worst correlation coefficients for all three NTB equations, contrary to Reference De Ruyter de Wildt, Oerlemans and Björnssonde Ruyter de Wildt and others (2002) who found a correlation of 0.92. However, if the years 2001 and 2002 are removed, the correlation for the mean summer 〈Q pot, net〉 and the measured mass balance becomes 0.8˚ for NTB conversion equation (6a) at Köldukvìslarjökull. Using data from different satellites may not be as straightforward as it seemed at first and this should be dealt with in future studies.
Assuming that the combination of the six drainage basins is representative of the whole ice cap both in space and time, we can estimate the total mass balance of the whole ice cap using the linear relation shown in Figure 6. This is shown for the years 1991-2002 in Figure 7. We intend to validate this result in future work by comparing our findings with modelling studies and by adding more satellite images. The derived relation is unique for Vatnajökull, but this method can be applied to other ice masses which are large enough compared to the AVHRR pixel size; it has already been applied to part of Greenland by Reference Greuell and OerlemansGreuell and Oerlemans (2005).
Fig. 7. Estimated mass balance for Vatnajökull based on the relation illustrated in Figure 6.
Conclusions
We tested a method to infer the mean surface mass balance of Vatnajökull from NOAA AVHRR satellite images for the years 1991-2002. The measured mass balance correlates well with mean summer 〈Q pot,net〉 if enough cloud-free images can be found in an individual year and if the temporal variation in measured mass balances is high enough. The correlation between measured mass balance and mean summer 〈Q pot,net〉 for all outlets combined was 0.92, with a residual standard deviation of 0.18 m w.e.
The calculated albedo is greatly affected by the selected NTB conversion equation and, to a lesser extent, by the choice of BRDF. It is affected slightly by variations in atmospheric water vapour but not ozone or visibility, although these may be important for other ice masses. The correction for anisotropy of the surface is an important step in the satellite-image processing procedure, but until BRDFs are made for ice and snow in Iceland, this step will remain an important source of error.
Finally, switching between AVHRR instruments may cause a serious error in retrieved albedos. A recalibration like that presented by Reference Greuell and OerlemansGreuell and Oerlemans (2005) may improve our results.
Acknowledgements
We thank M. de Ruyter de Wildt for his help with starting this project, A. Brooks and N. Lonie of Dundee Satellite Station for providing the AVHRR images, G. Gísladόttir of the Icelandic Meteorological Office for the visibility data, and B. Brock and A. Klein for critical and helpful reviews. Our research was sponsored by SPICE (Space borne measurements of Arctic glaciers and implications for sea level), European Union grant EVK2-2001-00262 SPICE.
