1. Introduction

The Greenland ice sheet (GrIS) is an important component of the global climate system, because changes in its mass affect sea level and, potentially, Arctic climate (Reference Fichefet, Poncin, Goosse, Huybrechts, Janssens and Le TreutFichefet and others, 2003). Recent observations (Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002; Reference Luckman and MurrayLuckman and Murray, 2005; Reference Joughin, Das, King, Smith, Howat and MoonJoughin and others, 2008; Reference Van de WalVan de Wal and others, 2008) have identified seasonal accelerations of the GrIS that are apparently linked to surface melting. The accelerations can occur after supraglacial lakes sited on the ice-sheet surface (Reference Echelmeyer, Clarke and HarrisonEchelmeyer and others, 1991; Reference Lüthje, Pedersen, Reeh and GreuellLüthje and others, 2006; Reference McMillan, Nienow, Shepherd, Benham and SoleMcMillan and others, 2007) drain (Reference DasDas and others, 2008). A proposed mechanism (Reference Van der VeenVan der Veen, 2007) to link surface melting and ice acceleration is ice fracture, which may occur at the tips of crevasses beneath supraglacial lakes once sufficient water has accumulated (Reference WeertmanWeertman, 1973). Although melting at the surface of the GrIS is expected to increase over the course of the 21st century (Reference Anisimov, M.L., Canziani, Palutikof, van der Linden and HansonAnisimov and others, 2007), the current generation of climate models (e.g. Reference Alley, Clark, Huybrechts and JoughinAlley and others, 2005) fail to account for the potential feedback associated with melt-induced acceleration of the ice sheet due to an absence of data and an inability to resolve or parameterize such processes. In consequence, future sea-level projections based on the results of such models (Reference Meehl and SolomonMeehl and others, 2007) are of limited value. An improved understanding of supraglacial lake evolution at the GrIS is therefore of considerable interest.

Supraglacial lakes on the surface of the GrIS are evident in a range of satellite imagery (Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002; Reference Lüthje, Pedersen, Reeh and GreuellLüthje and others, 2006; Reference Box and SkiBox and Ski, 2007; Reference Sneed and HamiltonSneed and Hamilton, 2007). Although calculating their area is relatively straightforward (Reference Lüthje, Pedersen, Reeh and GreuellLüthje and others, 2006; Reference McMillan, Nienow, Shepherd, Benham and SoleMcMillan and others, 2007), estimating the volume of supraglacial lakes requires, additionally, a measure of depth. Remotely sensed IKONOS satellite imagery has been used (Reference Stumpf, Holderied and SinclairStumpf and others, 2003) to retrieve depth measurements from terrestrial lakes, and various studies ( ; Reference PhilpotPhilpot, 1989; Reference Maritorena, Morel and GentiliMaritorena and others, 1994) have considered the optical properties of water at other frequencies. Models have recently been developed (Reference Box and SkiBox and Ski, 2007; Reference Sneed and HamiltonSneed and Hamilton, 2007) to estimate the depth of supraglacial lakes on the surface of the GrIS from optical satellite imagery. Here we extend a method (Reference Sneed and HamiltonSneed and Hamilton, 2007) based on data from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) to determine temporal changes in the volume of water stored in a supraglacial lake in West Greenland and hence characterize its seasonal evolution. We also compare satellite-derived lake depths to independent estimates based on airborne lidar observations to assess the accuracy of the technique.

2. Study Area and Data

This study focuses on a supraglacial lake located on Jakobshavn Isbræ at the western margin of the GrIS (Fig. 1). The lake, centred at 69˚18′ N, 48˚55′ W, is situated at approximately 1020ma.s.l., and its maximum observed area during our study period was 3.4 km2.

Fig. 1. Landsat image from 1 August 2001 of the Jakobshavn Isbræ outlet glacier in West Greenland. The lake used in this study is highlighted with a black square and shown in greater detail (inset). Also shown are the airborne lidar flight-line (solid line) and a longer transect (dashed line) which was used to estimate the maximum lake depth.

We use ten summertime images recorded by the ASTER satellite between 2002 and 2005 (Table 1) to analyse temporal changes in lake water reflectance. The images were pre-processed to the ASTER On-Demand Level 2 Surface Reflectance visible and near-infrared (VNIR) data product, which incorporates atmospheric corrections at the VNIR wavelengths. In this study, we use data recorded within two ASTER radiometric frequency bands with ground pixel resolutions of 15 m; a near-infrared band (VNIR3N, 780–860 nm) was used to identify regions of shallow water and, hence, the lake bed reflectance, and a visible band (VNIR1, 520–600 nm) was used to determine lake depth. The satellite data were received in Hierarchical Data Format (HDF) and projected into the Universal Transverse Mercator (UTM) coordinate system. Because the lake-depth retrieval algorithm also requires an estimate of reflectance from optically deep water, we used ASTER images acquired over open ocean from the same satellite tracks and dates as the lake imagery, which has similar illumination conditions.

Although there are 24 hours of daylight during summer at the latitude of the study area, illumination conditions alter with the changing solar zenith angle, so we refer to the satellite imagery as day- and night-time to reflect this variability. The 12 August 2005 image (Fig. 1, inset) is the only daytime image in which the lake appears to be full. Although several night-time images displaying lakes of similar area were also obtained, we were unable to estimate lake depths from these data due to unfavourable illumination conditions and due to the presence of streaks in the satellite imagery of unidentified origin. Instead, we estimated only the locations of maximum depth from the night-time images.

To study the evolution of water in the lake, it was useful to sequence the satellite imagery according to the degree of melting that had occurred prior to the date of each acquisition. Because the images were acquired during years with varying amounts of runoff, we used the quantity of positive degree-days (PDDs) accumulated to the date of each image as an alternative chronology. PDDs were calculated using temperature data recorded at four nearby meteorological stations: Swiss Camp, JAR1, JAR2 and JAR3 at elevations of 1149, 962, 568 and 323 m, respectively (Reference Steffen and BoxSteffen and Box, 2001). There is an approximately linear relationship between PDDs and elevation at these stations (Reference McMillan, Nienow, Shepherd, Benham and SoleMcMillan and others, 2007), so we were able to estimate PDDs at the elevation of the lake (1020 m) using data recorded at the Swiss Camp and JAR1 stations.

To evaluate the accuracy of our depth-retrieval model, we computed an independent estimate of the lake depth using aircraft lidar measurements acquired during NASA’s Program for Arctic Regional Climate Assessment (PARCA) and the span of the lake at maximum extent. The lidar data were recorded in May 1995 before the lake had formed in that year. We averaged ice-surface elevations from the near and far range of the lidar swath, which were approximately 90 m apart, to estimate the elevation along a transect of the lake bed (Fig. 1, inset). We then located the lake edge position along the transect using the 12 August 2005 ASTER image, and these points were assigned a depth of 0m. Inferred depths between these two positions were estimated as the difference in elevation relative to that of the lake edge. An uncertainty of ±1m is ascribed to the lidar depths to account for interannual variations in the lake bed height resulting from accumulation and ablation (Reference BoxBox and others, 2006).

3. Methodology

A method has been developed (Reference Sneed and HamiltonSneed and Hamilton, 2007) to estimate the depth of lakes on the surface of the GrIS based on their optical reflectance. The model was applied to atmospherically corrected, multispectral ASTER data to estimate the depth and volume of two supraglacial lakes in northwest Greenland (Reference Sneed and HamiltonSneed and Hamilton, 2007). We apply the same model to a sequence of ASTER data recorded in West Greenland, and we use independent lidar measurements to validate the model performance. The lake-depth retrieval model is derived from the Bouguet–Lambert–Beer law (Equation (1)), which describes how incoming radiance is attenuated as it passes through a medium. The spectral radiance, L(z,ʎ), measured at a wavelength ʎ and distance z from a target is given as

L(z, ʎ) = L(0, ʎ) exp (–K_ʎz), (1)

where L(0, ʎ) is the spectral radiance of the target and K_ʎ is the spectral attenuation of the intervening medium. In this study, the intervening medium through which the radiation passes is water and the distance, z, is the water depth.

Because the satellite imagery has been corrected for the effects of atmospheric attenuation, we assume that the spectral radiance at the lake surface is equivalent to the radiance measured by the satellite. Converting from radiance to reflectance, and solving Equation (1) for z, an expression (Equation (2)) can be derived (Reference PhilpotPhilpot, 1989) to determine the depth of the water column, given assumptions regarding lake-bed and water reflectance:

where A _d is the lake bed reflectance (or albedo), R∞ is the reflectance from optically deep water (>40 m), R_z is the reflectance measured by the satellite when the lake has a depth z, and g is the specific spectral attenuation. Once appropriate values for the parameters A _d, R∞ and g have been determined, Equation (2) provides a method for estimating lake depth from measured surface reflectance. Because this model ignores direct irradiance, its performance is most accurate when clouds are present distal from the target, least accurate under clear-sky conditions, and of intermediate accuracy during periods of low solar inclination (e.g. night). Our data are predominantly recorded during favourable conditions (night-time or in the presence of clouds). For example, we determine the maximum lake depth from an image containing 38% cloud cover. For a more detailed discussion of the lake-depth retrieval model, we refer the reader to the study in which it was developed (Reference Sneed and HamiltonSneed and Hamilton, 2007).

The specific spectral attenuation, g, accounts for absorption and scattering processes that result in the attenuation of light as it passes through the water column. A number of studies (e.g. Reference Maritorena, Morel and GentiliMaritorena and others, 1994) have related g to the total attenuation in fresh water which, in turn, may be related to frequency-dependent absorption and scattering coefficients (Reference Smith and BakerSmith and Baker, 1981). Assuming there is little suspended or dissolved particulate material in supraglacial lakes (Reference Sneed and HamiltonSneed and Hamilton, 2007), the degree of attenuation is small and, in consequence, close to that of fresh water. To estimate g, we average the freshwater absorption and scattering coefficients calculated at 10nm intervals (Reference Smith and BakerSmith and Baker, 1981) across the ASTER VNIR1 waveband (520– 600 nm), and, allowing for a modest increase in attenuation relative to fresh water (Reference PhilpotPhilpot, 1989; Reference Maritorena, Morel and GentiliMaritorena and others, 1994), we obtain a value of g=0.206±0.010m–1. On average, this range of specific spectral attenuation leads to a variance of 9.5% in estimated lake depth.

To determine a value for the reflectance from optically deep water, R∞, on each date, we used adjacent images along track including sea water at a depth greater than 40 m, which is considered to be optically deep (Sneed and Hamilton, 2007). In each image, we identified an area of dark open water that was free from clouds, breaking waves and sun glint, and we calculated R∞ as the spatially averaged reflectance from this area with an associated uncertainty given by the standard deviation. Within the frequency range used in this study, differences in absorption and scattering between fresh and sea water are non-existent and small, respectively (Reference Smith and BakerSmith and Baker, 1981). By choosing pixels located far from land, the amount of particulate matter and the effects of sun glint are also minimized (Reference Smith and BakerSmith and Baker, 1981; Reference Sneed and HamiltonSneed and Hamilton, 2007). The significance of open-ocean roughness in determining R∞ is poorly known. To minimize potential uncertainties in R∞ introduced by variations in ocean roughness caused, for example, by winds, we located the darkest ocean pixels (exhibiting the lowest reflectance) in each satellite image, and we used these values to estimate R∞. We also ensured that such pixels were distal from the shore to ensure water depths greater than 40m were sampled and to avoid the occurrence of breaking waves. In consequence, our estimate of R∞ is comparable to that of the fresh water found in supraglacial lakes. During the five daytime images used in this study, R∞ varied between 0.07 and 0.14 due to variations in solar illumination and surface roughness (Table 1). However, it proved difficult to obtain images containing optically deep water on some dates, due to increased cloud cover over the ocean (Table 1), a condition that was particularly prevalent at night. There is also a clear differentiation between values of R∞ obtained at day and night (Table 1). The uncertainty in lake depths resulting from the observed variability in R∞ was <9% when observations from nearby ocean were available, and 19% at other times.

A value of the lake bed reflectance, A _d, may be estimated as the reflectance of ice that is barely covered by water (Reference Sneed and HamiltonSneed and Hamilton, 2007), and it is possible to identify such locations by considering the ratio of the observed reflectance of the lake and adjacent ice in the VNIR3N frequency range. For example, according to Equation (1), the reflectance of ice covered by 10cm of water is about 40% less than if the water were absent. First, we used measurements of the lake bed profile determined from the lidar survey to refine the estimated cut-off at which shallow water may be identified. Along a 200m cloud-free section of the 12 August 2005 ASTER image (the lake did not encompass the lidar profile on any other date), the average lake bed reflectance was 0.76±0.17, and the variance between modelled and observed lake depths was minimized when only regions exhibiting a VNIR3N reflectance greater than 46% of the adjacent bare ice were considered. We employed this threshold to identify shallow lake sections in all of the ASTER images, and computed a mean value of A _d on each date (Table 1). The lake-depth retrieval model is particularly sensitive to the choice of A _d, and the observed variance in A _d leads to average uncertainties in lake depth of 85% for lakes shallower than 3.5 m, reducing to 30% for lakes 6–10m deep. The principal component of this uncertainty is associated with our estimated lidar depth error of 1 m, which accounts for interannual variations in the lake bed depth associated with fluctuations in accumulation and ablation (Reference BoxBox and others, 2006).

Calculating the volume of water stored in each lake was complicated by the presence of surface ice in many of the ASTER images. On 12 August 2005, when no surface ice was apparent, we estimated the average water depth along the deepest transect passing through the lake (dashed line, Fig. 1, inset) and multiplied the result by the surface area to estimate the lake volume. On all other dates, surface ice was present, so we developed an alternative method to estimate their volume. First, we identified the location in each lake with the lowest reflectance value in the VNIR1 band, and took this to correspond to the point of maximum depth, on the assumption that there exists little surface ice above deep water (lake ice generally forms first over shallow water, typically near the margins, and last over deep water, typically towards the centre, because the lake bed is formed of material (ice) more reflective than water, resulting in lower albedos and greater absorption of solar radiation in deeper water). We used the 12 August 2005 image, which showed no ice cover, to locate the position of maximum depth. In the other images used in the study, there appeared to be less surface ice cover or no ice cover over this area when compared with the shallower areas. We then computed an initial estimate of their volume based on the maximum depth and the lake’s semi-major and semi-minor axes (Table 2), assuming an ellipsoidal profile. To account for departure from such an idealized profile, we used the 12 August 2005 ASTER image to determine an empirical volume-scaling factor, calculated as the ratio of the volume retrieved using the product of surface area and average depth to the volume retrieved using the ellipsoidal representation. We then multiplied the ellipsoidal volume for each image by this factor scaled relative to the lake’s surface area in each image. For images where the surface area was <50% of that in the 12 August 2005 image, we kept the ellipsoidal approximation for the volume, including an additional error term to account for the uncertainty in bed geometry.

4. Results and Discussion

To assess the accuracy of the depth-retrieval model, we compared estimates of lake depth derived from the ASTER data on 12 August 2005 with those determined from the lidar transect (Fig. 2). Uncertainties in radiometric depths are estimated as a combination of the estimated errors in A _d, R∞ and g. We calculated the root-mean-square (RMS) deviation between the radiometric and lidar depths, and across the entire lidar transect the average deviation in height was 0.7 m. The smallest deviations were found in the region

Fig. 2. Scatter plot of lake depth estimates based on a radiative transfer model and airborne lidar surveys. Data in regions of partial cloud cover and cloud absence are marked in black and grey, respectively.

where there was no cloud cover, and within the cloud-free sections the RMS height deviation was 0.3 m. We also observe that as depth increases, the RMS deviation reduces, indicating that the model performs better in deep water. Depths across the lidar transect varied between 0 and 4 m. For depths less than 2 m, the RMS deviation is 0.9 m; for depths of 2–4m it is 0.6 m. This agrees with the pattern of variance for A _d, the primary source of uncertainty, which we observe decreases with an increase in radiometric depth (Table 1).

To estimate the maximum lake depth during the melt season, we identified a transect that passed through the deepest points of the lake (shown as the longer dashed line in Fig. 1, inset) on 12 August 2005, the date on which our imagery resolved the maximum lake area. We retrieved reflectance values in the VNIR1 band at 15 m intervals across this transect and used the depth-retrieval model, combined with our estimates of A _d and R∞, to convert these reflectances to depths. We estimate that the lake had a maximum depth of 9.6m and a mean depth of 5.5m on this date. From the average radiometric depth and lake area we calculated a lake volume of (18.6±3.7)×106m3 on 12 August 2005 (Table 2). Given the temporal resolution of our imagery, we were unable to determine whether this value represents the maximum lake volume during the summer melt season.

Because the satellite data were acquired during different years with different degrees of summer melting, we used an alternative chronology to investigate the seasonal evolution of supraglacial lakes. Lakes are fed by runoff from melting snow and ice as summer progresses, and positive degree-days (PDDs), the daily averaged sum of hourly temperatures above 0˚C, have been shown to be approximately proportional to the quantity of runoff produced as a result of snow or ice melting (Reference BraithwaiteBraithwaite, 1995). We computed the quantity of PDDs accumulated each year to the date of each satellite image (Table 2) from nearby meteorological observations (Reference McMillan, Nienow, Shepherd, Benham and SoleMcMillan and others, 2007), and from these data we were able to sequence estimates of lake depth derived for different years according to the degree of melting (Fig. 3) to describe their seasonal evolution. PDDs were estimated at the elevation of the lake (1020 m) by fitting a linear model to temperatures recorded at the Swiss Camp and JAR1 meteorological stations sited above (1149 m) and below (962 m) the lake, respectively. To estimate the associated error, we fitted a linear model to PDDs recorded at all four meteorological stations and computed the average difference between observed and modelled PDDs to be 12, equivalent to 5% of the average annual total.

Fig. 3. Evolution of lake volume according to the number of PDDs accumulated prior to the image acquisition date in each season.

Our data show that the supraglacial lake began to form between dates on which the cumulative PDDs reached 10–35. During the early part of the melting season (from ∽PDD 10 to ∽PDD 80) the lake volume increased at a rate of (0.12±0.03)×106m3d–1 (we converted from PDDs to a daily rate by matching PDDs to the date at which they occurred in 2005). After ∽PDD 80, the rate at which the lake filled accelerated to (0.37±0.06)×106m3d–1 (Fig. 3). We interpret this change in the rate of lake filling as an indication of the transition from snow- to ice melting at the ice-sheet surface, which occurs after a fraction of the winter snowpack has melted (Reference ReehReeh, 1991). Reference BraithwaiteBraithwaite (1995) used an energy-balance model to investigate variations in the rates of snow and ice melting under a range of climatic conditions. The effect was parameterized using empirically determined linear coefficients (PDD factors, measured in mmd˚C–1). The PDD factor for snow was found to be 33–44% of the corresponding value for ice surfaces, and it is suggested that the difference is due to lower energy fluxes occurring within snow. In addition, during the snowmelt phase, some of the meltwater will refreeze within the snowpack, resulting in even less water contributing to filling of the lake when compared with the ice-melting phase.

The proportional increase in the rate of filling (3.05) is comparable to the ratio of ice- to snow-melting factors obtained empirically in a range of studies elsewhere at the GrIS (Reference BraithwaiteBraithwaite, 1995). Our data resolve a peak lake volume of (18.6±3.7)×106m3 on PDD 124. Complete lake drainage had occurred by PDD 153, which suggests that during this period (corresponding to a ∽6 day interval in 2005) the lake drained at a minimum rate of (2.66±0.53)×106m3d–1. In reality, satellite- (Reference Box and SkiBox and Ski, 2007) and field-based (Reference Boon and SharpBoon and Sharp, 2003; Reference DasDas and others, 2008) studies have observed that lakes drain more rapidly, so it is likely that the actual drainage rate is significantly higher.

According to previous work (Reference Lefebre, Gallèe, van Ypersele and HuybrechtsLefebre and others, 2002), winter snowfall at the GrIS melts at a rate of 3 mmw.e. mass per PDD. Based on this rate of snowmelt and assuming that the transition between slow and rapid lake filling (PDD∽80) reflects the time at which all winter snow has melted, we estimate a winter snowpack of ∽24 cmw.e. This value is in agreement with the results of a climate model re-analysis (Reference Hanna, Huybrechts, Janssens, Cappelen, Steffen and StephensHanna and others, 2005) which determine that net accumulation in the Swiss Camp runoff area is 30±6 cmw.e. in 2001.

5. Conclusions

We employed a radiative transfer model (Reference Sneed and HamiltonSneed and Hamilton, 2007) to estimate temporal changes in the depth and volume of a supraglacial lake sited on the surface of Jakobshavn Isbræ. The model is based on data recorded by ASTER between 17 June 2002 and 12 August 2005. At its largest extent, the lake covered an area of 3.38 km2 and we estimate its maximum depth was 9.6±1.0 m. We estimate further that the lake volume at this stage was (18.6±3.7)×106m3. When compared to independent estimates of lake depth based on airborne lidar observations of the lake bed elevation, we find the RMS deviation between modelled and observed depths was 0.3m for cloud-free satellite imagery.

By considering the number of PDDs accumulated prior to each lake image date, we were able to sequence estimates of lake volume according to the degree of surface melting, which varies considerably from year to year (Reference BoxBox and others, 2006). The lake begins to pond after around ten PDDs, and lake drainage occurs between dates on which the accumulated PDDs reached 124–153. We estimate that in 2005 the lake drained between 12 and 24 August at a minimum rate of (2.66±0.53)×106m3d–1. An interesting feature of the lake evolution is that, after the accumulated PDDs exceed about 80, the rate of filling increases abruptly by a factor ∽3. We interpret this as corresponding to the period during which melting of the previous winter snowpack ceases, when melt rates are observed to increase by a similar factor (Reference ReehReeh, 1991; Reference BraithwaiteBraithwaite, 1995). A more detailed understanding of the extent to which supraglacial lakes affect flow of the GrIS will require an assessment of rates of water storage and drainage across a wider catchment.