1. Introduction
Ablation on the Greenland ice sheet is largely controlled by net shortwave radiation (Reference Greuell and KnapGreuell and Knap, 2000), and therefore variations in the surface albedo are very important. In the ablation area of the Greenland ice sheet, supraglacial lakes form during summertime. The lakes appear as ‘dark’ spots of low albedo on the ice sheet of much higher albedo. On Arctic sea ice Reference Perovich, Tucker and LigettPerovich and others (2002) found that the presence of melt ponds can lower the overall summertime sea-ice albedo from 0.65 to 0.4, and suggested that increased surface melting of sea ice will lead to greater melt-pond coverage, thus decreasing the area average albedo, and further enhance melting (Reference MorassuttiMorassutti, 1992; Reference Perovich and TuckerPerovich and Tucker, 1997).
Supraglacial lakes have been studied on the Greenland ice sheet (e.g. Reference ThomsenThomsen and Reeh, 1986; Reference Echelmeyer, Clarke and HarrisonEchelmeyer and others, 1991), on Svalbard (Reference Liestøl, Repp and WoldLiestøl and others, 1980) and on glaciers in the Himalaya (Reference AgetaAgeta and others, 2000; Reference Benn, Wiseman and WarrenBenn and others, 2000; Reference ReynoldsReynolds, 2000). Lake locations were mapped in a limited area of the West Greenland ice margin (Reference ThomsenThomsen, 1986; Reference Thomsen, Thorning and BraithwaiteThomsen and others, 1988). Comparing the position in different years, the lakes were found in the same locations (Reference Thomsen, Thorning and BraithwaiteThomsen and others, 1988; Reference Echelmeyer, Clarke and HarrisonEchelmeyer and others, 1991; Reference Jezek, Drinkwater, Crawford, Bindschadler and KwokJezek and others, 1993; Reference LüthjeLüthje, 2005). More recently, a study has measured the albedo of different surface types of the ice sheet, including the albedo of supraglacial lakes (Reference Greuell, Reijmer and OerlemansGreuell and others, 2002). However, the contribution of supraglacial lakes to enhanced ice-sheet ablation remains unstudied. This subject is addressed by the present study, in which we present a melt model that quantifies the difference between the ablation for a bare ice surface and the ablation beneath a supraglacial lake. The model is forced with climatic data measured on the Greenland ice sheet in two different years where the total surface ablation in July differed by 40%. The output from the model is used as input to a drainage model describing the vertical and horizontal evolution of supraglacial lakes in a limited area, situated around 67°20′ N, 48°57’ W (Fig. 1). Finally, we compare the result from the drainage model to the observed area coverage of supraglacial lakes as derived from moderate-resolution imaging spectroradiometer (MODIS) and Landsat 7 satellite images.
MODIS image taken on 26 August 2000 showing most of Greenland, with a resolution of (a) 1000 and (b) 250 m. (a) The arrows mark the positions of the field camp occupied in August 2003 and the climate station. (b) The area investigated for aerial coverage of supraglacial lakes.
2. Ice-Melt Model
A convenient way of solving a system with phase changes (Stefan problems) is to use the enthalpy method (Reference Alexiades and SolomonAlexiades and Solomon, 1993). Since enthalpy is defined as the total energy of the system, there is no need for tracking the phase boundary as it is a natural part of the solution. The basic structure of the melt model is the same for the case of a bare ice surface and an ice surface beneath a lake, with only small differences (commented on later). If the enthalpy is zero for a volume of ice at its melting point, Tm = 0°C, then the change in enthalpy, E, for a volume of ice at a lower temperature, T, is:
where (pc)s is the volumetric heat capacity of ice. In this case the liquid fraction, A, is zero. For a temperature greater than T m the change in enthalpy of the melted ice volume is instead given by:
where (pc)l is the volumetric heat capacity of water and L = 3.348 × 105 J kg–1 is the latent heat of fusion. The liquid fraction, A, is then 1.
If the temperature is at the melting point then 0 < E < pL and λ = E/pL. The volume of ice down to a depth Lz is divided into N control volumes, each of height ∆z = Lz/N and with a volume of ∆z × A, where A is the area of the control volume. N is fixed throughout, whereas Lz, and therefore ∆z, can be varied. All control volumes are initialized with an enthalpy corresponding to a temperature profile with a depth gradient of 2 Km–1. This assumption is based on measurements from Pâkitsoq, West Greenland (personal communication from O. Olesen and H.H. Thomsen). Since the drainage model is developed for a snow-free surface only, the calculation begins when the snow has disappeared (see section 2.3.1).
2.1. Conductive heat transport through ice
Energy transport, q, within the ice is assumed conductive and given by (Reference Alexiades and SolomonAlexiades and Solomon, 1993):
where j = 1,…, N is the ‘position’ of the volume element and kj is the effective conductivity and is assumed to be the sum of the conductivities of the different phases:
where k 1 = 0.569 W m–1 K–1 is the conductivity in the liquid state and ks = 1.88 W m–1 K–1 is the conductivity in the solid state (Reference Incropera and DeWittIncropera and DeWitt, 2002).
Since melting takes place in the surface layer and enthalpy variations only occur in the uppermost part of the ice sheet, only the top 4m of the ice were modelled. An initial grid size of 1 cm was used. At the bottom boundary, Neumann conditions were applied, i.e. no heat flux across the boundary.
2.2. Turbulent heat transfer through lake water
Since fresh water has its maximum density at 4°C, water that is heated from 0°C at the surface will be denser than the water at the bottom of the lake. Reference Taylor and FelthamTaylor and Feltham (2004) found that even for small melt ponds on Arctic sea ice, very small temperature differences between top and bottom resulted in turbulence. The heat flux between the turbulent lake and the ice bottom can be calculated using the ‘four thirds’ law (Reference Linden, Bachelor, Moffatt and WorsterLinden, 2000):
r is a dimensionless constant set to be 0.1 (Reference Huppert, Bachelor, Moffatt and WorsterHuppert, 2000; Reference Taylor and FelthamTaylor and Feltham, 2004), pwater — 1000 kgm–3 is the density of water, cl — 4217J kg–1 K–1 is the specific heat capacity (Reference Incropera and DeWittIncropera and DeWitt, 2002), β = 6.8 × 10–5 K–1 is the coefficient of thermal expansion, tl — kl/pwatercl is the thermal diffusivity, v — 1.75 × 10–6m2s–1 is the kinematic viscosity (Reference Incropera and DeWittIncropera and DeWitt, 2002), and AT is the temperature difference between the core temperature of the lake and the boundary temperature. A turbulent mixed water column is expected to have an interior (core) of a constant temperature, and the entire temperature difference is found at the boundaries (Reference Linden, Bachelor, Moffatt and WorsterLinden, 2000). Reference ThomsenThomsen and Reeh (1986) reported uniform temperature depth profiles in West Greenland supraglacial lakes with depths varying from 1.2 to 4.2 m, indicating turbulent heat transfer. In August 2003, this finding was confirmed by temperature measurements in two supraglacial lakes within the ice margin sector described in the present study; no difference in temperature with depth was observed. To achieve this in the model, the mean temperature for the lake is found by taking the average of the temperature in the lake, and all interior gridcells with a liquid fraction of 1 are then assigned this temperature.
2.3. Surface boundary condition
The melt model is forced by the net incoming energy flux at the interface between ice and atmosphere given by
Fsw and Flw are the measured shortwave and longwave incoming radiant fluxes. Fsen and FIat are sensible and latent heat fluxes, respectively, and FIw_out is outgoing longwave radiant flux, α is surface albedo. The parameter I0 represents the fraction of shortwave radiation that is not absorbed at the surface but heats deeper layers of the ice-water column.
The sensible and latent heat fluxes were calculated using a bulk method where the fluxes are proportional to the difference in temperature and vapour pressure (Reference OerlemansOerlemans, 1992; Reference Zuo and OerlemansZuo and Oerlemans, 1996).
where C (Wm–2K–1) is an exchange coefficient that, on average, accounts for the effect of variable wind speed and atmospheric stability, Tair is air temperature (K) at some reference level, T s is ice surface temperature (K), Lv — 2.501 × 106Jkg–1 is the latent heat of vaporisation, cair is the specific heat capacity of air (1.005 × 103J kg–1 K–1), p is the atmospheric pressure (Pa), esat is the saturation vapour pressure and eair is the vapour pressure (Pa) at some reference level.
The outgoing longwave radiation is calculated as 
, where ∈ is surface emissivity and σ = 5.67 × 10–8 Wm–2 K–4 is the Stefan-Boltzmann constant.
2.3.1. Albedo and transmitted fraction I0 of shortwave radiation
For the bare ice surface, the albedo was calculated from observed values of incoming shortwave radiation Fsw_in and reflected shortwave radiation Fsw_out (see section 3):
As the model is developed for a snow-free surface, the calculated surface albedo was also used to estimate the date when the snow cover had melted. When the albedo dropped below 0.6 we assumed that the snow layer had disappeared. In 1999 and 2001, this happened on 3 July at the latest. We ran the model for 30 days starting at this date. No refreezing of meltwater occurred in the model periods.
Unfortunately, few measurements of albedo vs depth of supraglacial lakes are reported in the literature. We are only aware of the measurements of Reference BoxBox (1997) from a West Greenland ice margin location at 650m elevation. Measured albedo ranged from 0.23, with pond depth less than 1 cm, to 0.15 when pond depth was 106 cm.
On the other hand, measurements of the albedo–depth relationship for melt ponds on Arctic sea ice are abundant (Reference Morassutti and LeDrewMorassutti and LeDrew, 1996; Reference Hanesiak, Barber, De Abreu and YackelHanesiak and others, 2001). Reference Ebert and CurryEbert and Curry (1993) used data from the Arctic Ocean sea ice published by Reference Grenfell and MaykutGrenfell and Maykut (1977) and Reference Grenfell and PerovichGrenfell and Perovich (1984) to parameterize pond albedo vs depth for four spectral bands. Reference Morassutti and LeDrewMorassutti and LeDrew (1996) also estimated a relationship between albedo and depth of melt ponds, based on a much larger dataset. The two parameterizations differ most at visible wavelengths. Reference Ebert and CurryEbert and Curry (1993) found a generally lower albedo that was a function of depth until around 70 cm. Reference Morassutti and LeDrewMorassutti and LeDrew (1996) found that the albedo was independent of depth for depths greater than approximately 20 cm. Reference Taylor and FelthamTaylor and Feltham (2004) took a different approach by using a two-stream model to get an approximation for the albedo vs depth. Using the method of Reference Taylor and FelthamTaylor and Feltham (2004), although with a lower albedo for bare ice (0.55 instead of 0.65) to account for the generally lower albedo of glacier ice compared with that of sea ice, the following expression for the albedo of a supraglacial lake as a function of the depth of the lake h is obtained:
This albedo-depth variation is shown in Figure 2. Compared with the few measurements by Reference BoxBox (1997), Equation (10) gives higher albedo values for h < 0.35 m, and lower values for h > 0.35 m. The maximum difference between the observations of Reference BoxBox (1997) and Equation (10) is 0.07 for h = 0.9 m. The difference decreases rapidly with increasing depth. Although one may argue that melt ponds on sea ice typically have lower albedo than supraglacial lakes, we shall use the albedo–depth relationship given by Equation (10), realizing that the magnitude of our modelled melt rates beneath supraglacial lakes may then be overestimated. This point is commented on further in section 6.
Albedo vs lake depth based on the method presented by Reference Taylor and FelthamTaylor and Feltham (2004), using a bare-ice albedo of 0.55.
The parameter I0 represents the fraction of shortwave radiation that is not absorbed at the surface but heats deeper layers of the ice-water column. In accordance with common practice, we assume that all solar energy is absorbed at a bare ice surface, rather than penetrating into the subsurface. Actually, choosing I0 = 0 is necessary to make calculated melt rates match observations.
In the case of the supraglacial lake, a value of I 0 ≠ 0 is used. According to Reference PerovichPerovich and others (1999), nearly all shortwave radiation above 700 nm is absorbed effectively at the water surface, whereas Reference GrenfellGrenfell (1979) found the fraction of incident shortwave radiation above 700 nm to be 40%. In accordance with the latter result, we use I 0 = 0.6. The part of the radiation not reflected or absorbed at the lake surface is attenuated through the lake according to Beer’s law (Reference Ebert and CurryEbert and Curry, 1993):
where z is the depth below the lake surface, and t = 0.025 m–1 (Reference Taylor and FelthamTaylor and Feltham, 2004).
2.3.2. Elevation variation of radiant fluxes
The data used as input to the surface boundary condition of the melt model were collected on the West Greenland ice margin (67°5′ N, 49°23’ W; Fig. 1) at an elevation of 1000 m. In order to estimate the melt rate at other elevations and locations, we need to assess the elevation gradient of various parameters. The shortwave incoming radiant flux is assumed to be independent of elevation (within the range 1000–1400 m). To find the change with altitude an approximation of the longwave incoming radiation is used (Reference Konzelmann, Van de Wal, Greuell, Bintanja, Henneken and be-OuchiKonzelmann and others, 1994):
where ∈ cs is the clear-sky emissivity, n the cloudiness (0 → 1), p — 4 is a constant and ∈oc = 0.952 is the emissivity of an overcast sky. ecs was approximated by Reference Konzelmann, Van de Wal, Greuell, Bintanja, Henneken and be-OuchiKonzelmann and others (1994) as:
with eair in Pa and Tair in K. b = 0.484 and m = 8 are constants. For the cloudiness we used n = 0.6. The choice of this value is arbitrary, but lies within the range of values (0.51–0.67) suggested for nearby ice-margin sites by Reference Konzelmann, Van de Wal, Greuell, Bintanja, Henneken and be-OuchiKonzelmann and others (1994), Van de Reference Van de Wal and OerlemansWal and Oerlemans (1994) and Reference Zuo and OerlemansZuo and Oerlemans (1996).
Using this approximation (Equation (12)), the longwave incoming radiation at 1000, 1200 and 1400m is calculated. The decrease with altitude is found and the result subtracted from the measured longwave incoming radiation to give an estimate of the incoming longwave radiation at other elevations.
The outgoing longwave radiation was calculated assuming an emissivity, ∈, for bare ice of 0.99 and an emissivity for supraglacial lake surfaces of 0.97, i.e. similar to the values for sea ice and melt ponds, respectively (Reference Ebert and CurryEbert and Curry, 1993).
2.3.3. Sensible and latent heat fluxes
For the exchange coefficient, C, used to calculate the sensible and latent heat fluxes (see Equations (7) and (8)), Reference OerlemansOerlemans (1992) used a constant value of 8 W m–2 K–1. Reference Zuo and OerlemansZuo and Oerlemans (1996) used a value for C that decreased exponentially with elevation. For elevations of 1400 and 1000 m, C was calculated as 3.3 and 5 Wm–2 K–1, respectively. Reference AhlstrømAhlstrøm (2003) used somewhat higher values of 8 and 10 W m–2 K–1. For the present work, both 8 and 10 W m–2 K–1 were tested, but 8 Wm–2 K–1 was chosen for all calculations because it gave the best fit to observed melt rates.
The reference vapour pressure, eair, appearing in Equations (8) and (13) is calculated as eair = esat(Tair) × qrh, where qrh is the relative air humidity assumed independent of elevation. Lacking measurements of qrh, we use the mean value qrh = 0.86 measured for a 30 day period in July 2002.
The saturation vapour pressure e sat over ice and water needed to estimate the latent heat flux (see Equation (8)) is calculated from approximations suggested by Reference KrausKraus (1972):
and
Temperature and pressure are parameters that have a strong variation with elevation. We use a summer temperature lapse rate of 0.5 K(100m)–1 (Reference Steffen and BoxSteffen and Box, 2001), and a pressure gradient of 10 Pam–1 following an approximation suggested by Reference AhlstrømAhlstrøm (2003).
3. Drainage Model
The drainage model is a modified version of the twodimensional model presented by LüReference Lüthje, Feltham, Taylor and Worsterthje and others (2006) to calculate the area covered by supraglacial lakes, using the calculated melt rates from section 2 as forcing. The elevation of the ice surface at position i,j at time t is denoted by H. The depth of a possible layer of meltwater is given by h so the surface topography of the water-covered ice surface is:
The initial surface topography 
comes from a digital elevation model (DEM) and is explained in more detail in section 4. The principle of mass conservation applied to the meltwater layer yields (Reference Lüthje, Feltham, Taylor and WorsterLüthje and others, 2006)
He(h) is the Heaviside function preventing h from becoming negative. For h ≥ 0, He(h) = 1, otherwise He(h) = 0. m is melt rate calculated as explained in section 2. pice = 900 kgm–3 is ice density. u = (u,v) is defined as the horizontal flux vector (the depth-averaged flow velocity in the layer of meltwater) with components u and v along the × and y axis, respectively. ∇ = (∂/∂x, ∂/∂y) is the horizontal gradient vector.
Water is transported to areas of lower gravitational potential by percolation through porous ice or snow, by supraglacial run-off in small streams and large rivers up to several metres in width (Fig. 3), and by subglacial run-off via moulins. We refrain from attempting to model this extremely complicated run-off pattern, the more so as the surface elevation DEM cannot resolve the streams and rivers. Instead, we make the simplest possible assumption, that u is proportional to the gradient of the surface topography:
Photographs from the field site on the Greenland ice sheet occupied in August 2003. (a) A supraglacial lake at 67°20′ N, 48°57’ W at an altitude of 1240 m. The approximate locations of sites shown in (b-d) are marked. The distance between locations shown in (b) and (d) is around 1000 m. (b) Water entered the lake by the river with an estimated velocity of 3–5m s–1 and width >2 m. (c) A small river (< 0.5m wide) of the type typically found all around the lake. (d) The outlet river that was much deeper and also broader (approximate width ≥ 5 m) than the inlet river. Water here flowed at a much lower velocity.
where Ω (ms–1) is a tuning parameter representing the change in run-off with change in slope. Equation (17) can now be written as:
The evolution of ice surface elevation is given by (Reference Lüthje, Feltham, Taylor and WorsterLüthje and others, 2006):
Equations (19) and (20) were discretized using finite differences, and the resulting algebraic equations were solved numerically. The maximum time-step used was 5 s. Meltwater that reaches the downslope boundary disappears (run-off), whereas closed boundaries are applied to the other three sides of the model.
4. Data
4.1. Input data for melt model
The data used as forcing for the melt model are shortwave incoming radiation, shortwave outgoing radiation, longwave incoming radiation, air temperature and atmospheric pressure. The data were collected on the West Greenland ice margin (67°5′ N, 49°23’ W; Fig. 1) at an altitude of 1000 m. Temperature, pressure, shortwave and longwave incoming radiation for 1999 and 2001 are displayed in Figure 4a and b, respectively.
Air temperature, atmospheric pressure, shortwave incoming radiation and longwave incoming radiation observed at 67°5′ N, 49°23’ W at an elevation of 1000m, for a 30 day period beginning on 3 July (a) 1999 and (b) 2001.
4.2. Input data for the drainage model
The surface DEM used as initial condition for the drainage model was produced from synthetic aperture radar (SAR) data from 1995 (personal communication from J. Mohr). The horizontal resolution is 100m × 100m and the relative vertical accuracy is of the order of 0.5 m. One can therefore expect that the DEM will display artificial surface irregularities with amplitudes of the order of 0.5 m. The model will accumulate water in the resulting artificially created depressions. Because of the relatively large surface slope (0.009) of the area studied, real lakes with depths less than 0.5m do not occur. To minimize the effect of artificially created surface irregularities we therefore assume that lakes are not formed as long as the water depth is less than 0.5 m.
The SAR images used to produce the surface DEM are from the fall or winter, when the supraglacial lakes are frozen over. We assume that the lakes are frozen solid during the winter and, therefore, that large subsurface water pockets do not exist at the lake formation sites at the start of the subsequent melt season. This means that the frozen lakes are assumed to be included in the surface ice layer when the melting begins.
The drainage model was applied to a 20 km × 20 km area located around a lake at 67°20′ N, 48°57’ W (Figs 1 and 3) where fieldwork was carried out in 2003. The area was divided into gridcells of 100m × 100m. Using the melt model, daily melt rates at 1000, 1200 and 1400m elevation were calculated for a bare ice surface and beneath a supraglacial lake (with an initial depth of 50 cm). Based on these results, daily melt rates were calculated by linear interpolation for each gridcell and used as input to the drainage model. The tuning parameter in Equation (18) was set to fi = 500 ms–1, but other values of fi (1, 10, 100, 250 ms–1) were also tested.
5. Results and Discussion
5.1. Melt model
Table 1 shows the model-simulated total surface ablation after 30 days of melting in 1999 and 2001 for a bare, well-drained ice surface at three different elevations. Much higher surface ablation was found in 1999 than in 2001. In Figure 5a and b the modelled surface ablation at 1000m elevation in 1999 and 2001, respectively, is compared with the observed surface ablation. The observed melt rate is more variable than the modelled melt rate, and is even negative (accumulation) for a few days. Most likely the discrepancies are due to the fact that the model does not consider precipitation either as snow or as rain. However, the overall correspondence between model and observation for the 30 day model period is excellent. As the model also captures the large observed difference between the amount of melting in 1999 and 2001, we are fairly confident of the validity of the model.
Accumulated surface ablation after 30 days of melting for a bare ice surface at three different elevations in 1999 and 2001
Accumulated surface ablation observed (dashed line) and modelled (solid line) for a bare ice surface in (a) 1999 and (b) 2001.
Table 2 shows the model-simulated total ablation after 30 days of melting for an ice surface at 1000m elevation if meltwater is not allowed to drain away. In this case, the total melting is almost doubled compared with the well-drained situation. With increasing values of the initial water depth, the total ablation increases even further (Table 2; Fig. 6a and b). Up to an initial water depth of 50 cm, the surface ablation increases rapidly, after which the effect is much reduced. The total ablation for the 30 day period is about two times greater compared with the bare-ice situation in 1999, and almost three times greater in 2001. This finding agrees very well with measurements on Arctic sea ice where the melt rate beneath melt ponds is quoted to be two to three times greater than for bare ice (Reference Fetterer and UntersteinerFetterer and Untersteiner, 1998).
Accumulated ablation after 30days of melting beneath a lake initialized with five different depths of water in 1999 and 2001
Modelled accumulated surface ablation in (a) 1999 and (b) 2001 when water cannot drain, for five different initial water layer depths. The surface ablation increases up to an initial water layer depth of 50 cm, after which the effect is much reduced.
Table 3 shows for three different elevations the model-simulated total surface ablation after 30 days of melting in 1999 and 2001 for an ice surface where meltwater is not allowed to drain away.
Accumulated ablation after 30 days of melting beneath a lake at three different elevations in 1999 and 2001. The model was initialized with 50cm of water
5.2. Drainage model
Figure 7 shows for the 20 km × 20 km model area the calculated evolution of the percentage area covered by supraglacial lakes in 1999 and 2001. The area covered by supraglacial lakes peaked in 1999 at 7.7% and in 2001 at 7.2% of the total area, reaching approximately constant values after 10–20 days. However, the process took longer in 2001 than in 1999. Figure 8a and b illustrate the modelled evolution in 2001 of the supraglacial lakes after 1 and 4 weeks of melting, respectively. After 1 week the lakes are not yet fully developed (Fig. 8a), but after 2 weeks only little change takes place. A comparison of the areas covered by the lakes after 4 weeks in 2001 and 1999 (Fig. 8b and c) shows only minor differences. The modelled maximum lake depths in 2001 and 1999 (10.16 and 10.28 m) were also very similar. For comparison, Figure 8d shows the location and magnitude of supraglacial lakes in the model area derived from analysis of a Landsat 7 image from 23 August 2000. There is good agreement between the modelled and observed positions of the lakes. However, the area covered by lakes in the Landsat 7 image is much less (2.2%) than the lake areas calculated by the model (about 7%). Comparing the location of the supraglacial lakes to a contour map of the area reveals that the lakes, as expected, form in depressions on the ice-sheet surface. The similarity of modelled area extent and lake depth in 1999 and 2001 in spite of the much greater total melt in 1999 indicates that the depressions in the ice surface fill up to a maximum water level, after which all added water drains away.
Area fraction covered by supraglacial lakes, vs time, as modelled for 1999 (solid line) and 2001 (dashed line).
(a, b) Area coverage and depth of supraglacial lakes as modelled in 2001 after (a) 1 week and (b) 4 weeks. The modelled area is 20km × 20km, located between 67°17′ and 67°28′ N, 48°43’ and 49°11’W. (c) The area covered by supraglacial lakes after 4weeks in 1999. The maximum area covered was 7.7% in 1999 and 7.2% in 2001. Maximum depths were 10.28 and 10.16m for the two years, respectively. (d) An analyzed Landsat 7 image from 23 August 2000 showing the location of the supraglacial lakes. The area covered by lakes constituted 2.2% of the surface area.
Figure 9 shows the area covered by supraglacial lakes vs time calculated for different values of the tuning parameter fi. The magnitude of fi has a large influence on the modelled extent of the supraglacial lakes, particularly for low values of fi. The average surface gradient in the area investigated is around 0.009. Hence, according to Equation (18), the mean velocity of the meltwater flowing from high to low areas is u = 0.009 fi. As mentioned in section 3, the meltwater drainage from the ice margin involves different processes with order-of-magnitude different flow speeds such as slow percolation through ice or snow, supraglacial run-off in small streams and big rivers and subglacial run-off via moulins. A reasonable assumption is that run-off is dominated by stream and river run-off with flow velocities that may reach several metres per second (Fig. 3). Even in lakes, the current can be strong. In a lake 600m in length close to 67°09′ N, 49°38’ W, for example, the central current was about 3 ms–1 (Reference Maurette, Hammer, Brownlee, Reeh and ThomsenMaurette and others, 1986). This justifies our choice of fi = u/0.009 equal to 500 ms–1. However, as long as fi is greater than 100 ms–1 the influence on the area extent of supraglacial lakes is limited (Fig. 9).
Modelled area fraction covered by supraglacial lakes In 2001 for five different values of drainage parameter Ω. When Ω> 100ms–1 the influence on the area fraction is limited.
5.3. Analysis of satellite images
To investigate the large-scale area coverage of supraglacial lakes in Greenland, MODIS satellite images from the summers of 2000, 2001 and 2002 were analyzed. We used channels 1 and 2 (0.62–0.67 and 0.841–0.876 μm; 250m nadir resolution). Supraglacial lakes occur within much of the melt zone of the Greenland ice sheet, but are especially abundant on the western margin south of 70°N. An area between 64°6′ N and 71°12′ N from the ice margin and well into the dry-snow zone was chosen for our investigation (Fig. 1b). The area was divided into four sub-areas. Cloud-free images for each of the four sub-areas from July 2000, 2001 and 2002 were geocoded and analyzed. The mean July coverage of supraglacial lakes for the total area was 0.69% in 2000, 0.93% in 2001 and 0.87% in 2002. In order to study to what extent these results were influenced by the 250m nadir resolution of the MODIS sensor, one of the subareas was also analyzed using a Landsat 7 image from 23 August 2000 using the panchromatic channel (0.52–0.9mm;15m nadir resolution). Only 6% difference in the area covered by supraglacial lakes was found. For the small area studied with the drainage model, analysis of the Landsat 7 image showed that supraglacial lakes actually covered only 2.2% of the area (Fig. 8d), i.e. much less than the approximately 7% obtained with the drainage model. One would expect the model-calculated area covered by lakes to be larger than the area observed in the Landsat 7 image, since this image was taken quite late in the melt season, when refreezing of the lake surfaces had begun (Fig. 10). However, it appears that the lakes in the model are formed at the same places as observed, but they are much larger and some of them are not even resolved as individual lakes. The explanation for the discrepancy is probably that the horizontal resolution of the DEM is not good enough to capture the drainage channels, and that vertical drainage through moulins is not included in the model.
Landsat 7 scene from 23 August 2000 showing the study area. The arrows indicate some of the locations where refreezing of supraglacial lake surfaces has started.
Near the ice margin, a significant fraction of the meltwater drains through crevasses and moulins to the glacier bed (Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002). In an area close to Swiss Camp on the West Greenland ice margin (69.57° N, 49.31° W) the area density of moulins is around 0.2 moulins km–2 (Reference Thomsen, Thorning and BraithwaiteThomsen and others, 1988; Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002). To investigate the importance of moulins, the DEM was modified slightly by introducing drainage ‘holes’ acting as sinks placed randomly in depressions of the ice surface with an area density similar to that observed near Swiss Camp. This change reduced the maximum area covered by supraglacial lakes from 7.7% to 4.1% when the model was run with the 1999 data.
6. Conclusion
Using a one-dimensional high-resolution numerical model forced by observed climatic data, we simulated summer ablation for a bare ice surface in a 30 day period. The agreement with observed surface ablation in two different years (1999 and 2001) with highly different ablation rates was excellent. Next, we simulated the melt rate beneath a supraglacial lake, using a parameterization of albedo as a function of water depth based on observations from melt ponds on Arctic sea ice. The simulation also accounted for the convective heat transfer through the water layer. The surface ablation beneath the lake was enhanced by 110% in 1999 and by 170% in 2001 compared with the ablation for bare ice.
Using the calculated melt rates as input to a drainage model, we calculated the surface evolution of supraglacial lakes within a 400 km2 area of the West Greenland ice-sheet margin. Despite much higher melt rates in 1999 than in 2001, the lake area was only 7% larger in 1999 than in 2001. Compared with satellite-derived data, the model-calculated area coverage of supraglacial lakes was three to four times too large. We suggest that this discrepancy is due to the fact that the resolution of the applied ice-surface DEM was too low and that vertical drainage through moulins was not accounted for in the model. The importance of moulins was tested by including randomly placed drainage holes in surface depressions acting as sinks. This nearly halved the area covered by supraglacial lakes.
Obviously supraglacial lakes are formed in depressions of the ice-sheet surface and, when a given water level is attained, further added meltwater will run off. Within the region where supraglacial lakes presently occur, a significant growth of the lake area in a warmer climate is therefore not to be expected unless the surface topography is changed dramatically. Furthermore, analysis of satellite images shows that less than 1% of the ablation region in southwest Greenland is presently covered by supraglacial lakes. Therefore, despite an enhanced melt rate of up to 170% beneath the lakes, the extra surface ablation caused by the supraglacial lakes within this area is minimal. However, in a warmer climate, surface melting will occur higher up on the ice sheet where small surface slopes will favour formation of large supraglacial lakes. Enhanced surface melting beneath such lakes is a hitherto overlooked feedback mechanism related to climate warming.
Acknowledgements
The elevation map was produced by J. Mohr using European Remote-sensing Satellite (ERS) data: ERS-1, 20 October 1995, orbit 22300;ERS-2, 21 October 1995, orbit 02627; ERS-1, 29 December 1995, orbit 23302;and ERS-2, 30 December 1995, orbit 03629. The fieldwork on the Greenland ice margin was performed by C. Lüthje and M. Lüthje, supported by the Arctic Technology Centre of the Technical University of Denmark (DTU), COWIfonden, the Geological Survery of Denmark and Greenland (GEUS), Hilleberg, Ketner Outdoor A/S and Fujifilm Denmark. The work on aerial classification of the MODIS images was carried out with help from N. Elmqvist and P. Ravn. We also thank P. Miller and P. Taylor for language revision and J. Walder, I. Joughin, C. Lüthje and two anonymous referees for thorough reviews that led to significant improvements in the paper.
