A decade (2002–2012) of supraglacial lake volume estimates across Russell Glacier, West Greenland

. Supraglacial lakes represent an ephemeral storage buffer for meltwater runoff and lead to signiﬁcant, yet short-lived, episodes of ice-ﬂow acceleration by decanting large meltwater and energy ﬂuxes into the ice sheet’s hydrological system. Here, a methodology for calculating lake volume is used to quantify storage and drainage across Russell Glacier, West Greenland, between 2002 and 2012. Using 502 MODIS scenes, water volume at ∼ 200 seasonally occurring lakes was derived using a depth–reﬂectance relationship, which was independently calibrated and ﬁeld validated against lake bathymetry. The inland expansion of lakes is strongly correlated with air temperature: during the record melt years of 2010 and 2012, lakes formed and drained earlier,


Introduction
Meltwater runoff from the Greenland Ice Sheet (GrIS) has increased since the 1960s (Hanna et al., 2005) and has been linked to seasonal ice-flow variability through basal hydrological forcing (Zwally et al., 2002;Joughin et al., 2008;van de Wal et al., 2008;Bartholomew et al., 2010). Horizontal and vertical ice surface motion lags peak daily melt and responds to rapid supraglacial lake (SGL) tapping Das et al., 2008;Doyle et al., 2013). Furthermore, the spatial pattern of seasonal flow acceleration coincides with the development of ice sheet surface hydrology Fitzpatrick et al., 2013) where SGLs play a vital role via temporary storage and routing of surface meltwater. SGLs, which are prevalent across the entire GrIS (Selmes et al., 2011), also drain rapidly via hydraulic fracture propagation in ∼ 2 h Doyle et al., 2013), forming moulins that provide a direct hydraulic connection to the ice sheet bed. Previous studies on SGLs have quantified spatial and temporal variations in SGL area across different sectors of the GrIS (e.g. Georgiou et al., 2009;Sundal et al., 2009;Selmes et al., 2011;Liang et al., 2012), and have estimated their volume (e.g. Box and Ski, 2007;Sneed and Hamilton, 2007), whilst others have attempted to model SGL evolution (e.g. Luthje et al., 2006;Banwell et al., 2012) and measure lake-bed melt rate . The extent to which temperatures correlate with melt volume and SGL storage remains unclear. Additionally, modelling studies (e.g.  (lakes F and Z) are shown in addition to the supraglacial hydrological network, based on manual digitisation of LANDSAT and ASTER imagery between 2009-2012. The boundary of the defined Russell Glacier catchment is also shown; (b) a comparison of the catchment boundaries used in this study with those of previous studies (van Palmer et al., 2011;Mernild and Hasholt, 2009;Bartholomew et al., 2011;Chandler et al., 2013;Cowton et al., 2013), overlaid onto a Landsat 7 image from July 2012. Clason et al., 2012;Banwell et al., 2012) infer the existence of a volumetric threshold beyond which rapid lake drainage is inevitable, however this has yet to be validated by field observations.
We investigate the seasonal evolution of SGLs over Russell Glacier using 11 yr (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) of quasi-daily lake volume estimates. This period includes extremely warm years when melt records were set (Nghiem et al., 2012;Tedesco et al., 2013), and therefore provides an opportunity to determine how the evolution of SGLs are modulated by regional climatic forcing. The year 2010 was the warmest on record (Cappelen et al., 2012;Box et al., 2010), 2.5 standard deviations above the 1973-2010 average in West Greenland (van As et al., 2012). Temperatures for April, June and July in 2012 also rank amongst the warmest on record (Cappelen et al., 2012), resulting in the second warmest summer (after 2010) and the largest spatial extent of melt as yet on record (Nghiem et al., 2012). We present a field calibrated algorithm for calculating changes in SGL volume during these years and compare them with longer term patterns derived over the period 2002-2012.

Study area
Russell Glacier is an ideal region to investigate SGL dynamics as this sector of the GrIS accounts for the highest areal extent of SGLs across the ice sheet (Selmes et al., 2011). The study area covers ∼ 6500 km 2 of the western GrIS between 66.8 and 67.3 • N, centring on Russell Glacier and extends from the ice margin to ∼ 2000 m elevation (∼ 150 km inland, Fig. 1), and is hereafter referred to as the "Kangerlussuaq sector" of the GrIS. For hydrological impact assessment, this region benefits from a long-term meteorological record obtained from a transect of six automated weather stations extending well above the local equilibrium line altitude  (van As et al., 2012;van de Wal et al., 2012) and bulk proglacial discharge gauged near a bridge across Watson River, Kangerlussuaq, since 2007 .

Methods
To investigate supraglacial lake dynamics, we extracted lake area, depth and volumes from 11 yr of daily-acquired MODIS imagery. The following sections detail the methods, limitations and uncertainties of this technique. The methods for deriving ancillary data (surface runoff data and proglacial discharge) and their spatial coverage are discussed. Previous studies classified supraglacial lakes in MODIS imagery using manual digitisation (McMillan et al., 2007;Lampkin and Vanderberg, 2011), band thresholds (Box and Ski, 2007;Liang et al., 2012), and object-orientated approaches (Sundal et al., 2009;Johansson et al., 2013). Here, we opted for a semi-automatic strategy to classify SGL extents, using a modified normalized difference water index (NDWI) (Huggel et al., 2002;Yang and Smith, 2013, Eq. 1) and freely available RSGISlib software (http://www.rsgislib. org). Before classification, MODIS bands within the visible spectrum (0.46 to 0.57 µm; bands 3 (blue) and 4 (green)) were sharpened from 500 m to 250 m resolution, using the ratio of bands 1 (red) at 250 m and 500 m resolution (L. Gumley et al., http://earthdata.nasa.gov/sites/default/files/field/ document/MODIS_True_Color.pdf). The approach for identifying water cover was to isolate regions with a high reflectance in the blue spectral band (band 3), utilising their striking colour contrast against the surrounding ice. Thresholds for both red and blue bands were used in combination with the NDWI threshold to reduce misclassification of pixels with spectrally similar signals to SGLs. A cloud filter was also included in the classification using MODIS band 6 (1628-1652 nm) which is sensitive to solid and liquid phase water, and final SGL classification was subsequently carried out on pixels which had band 6 reflectance of less than 0.15 consistent with Box and Ski (2007).

Depth and volume estimation
Methods for calculating the volume of water contained in SGLs from satellite imagery have been established based on the principles of light transmission through water. Lake reflectance in MODIS band 1 reduces exponentially with depth allowing a depth-reflectance function to be extracted (Box and Ski, 2007). However, other studies have used a physically based retrieval to derive depth based on lake bottom albedo using the Bouger-Lambert-Beer law and applying the technique to ASTER scenes providing estimates of lake volume (Sneed and Hamilton, 2007). Numerically modelled depths using this technique were subsequently validated with in situ measurements (Sneed and Hamilton, 2011;Tedesco and Steiner, 2011).
In this study we opt to calculate lake volume using an empirically derived depth-reflectance relationship. This relationship was calibrated using bathymetry data of two SGLs surveyed in July 2010, referred to here as lake F (67.00 • N, 48.71 • W) and lake Z (66.96 • N, 48.78 • W) (Fig. 1a). On 26 June 2010, lake F was instrumented with water level sensors and on 29 June the lake perimeter was surveyed using a Leica System 500 Global Positioning System (GPS) receiver differentially corrected against a bedrock-mounted reference station. Following lake drainage on 30 June , six elevation transects were surveyed over the lakebed using the GPS recording a location point every second. The combination of lake bathymetry and depth time series allowed for shoreline reconstruction . The second, undrained lake was surveyed by mounting a depth sounder (Garmin X10) and GPS receiver (Leica System 500) on an inflatable raft, whilst the perimeter was surveyed using a similar GPS set-up. The depth sounder has a stated accuracy of ±0.1 m. The GPS transect data and the extrapolated margins for both lakes were interpolated onto a 10 m × 10 m grid using kriging interpolation to yield lake bathymetry maps (Fig. 2). Lake perimeter measurements were overlain on MODIS imagery from corresponding days. MODIS pixels that contained more than 50 % water were accepted as lake pixels, constraining the minimum and maximum band thresholds for water. SGL bathymetries (Fig. 2) were interpolated to 250 m resolution to match the MODIS band 1 resolution, so that data could be extracted and plotted to reveal the relationship between depth and surface reflectance (Fig. 3a). Following Box and Ski (2007), we approximated the relationship between depth and surface reflectance using a least-squares fit, where D is depth and R is MODIS band 1 reflectance: The equation used in this study differs slightly to the one used by Box and Ski (2007) in that lake depth has a shallower upper and lower limit in relation to values of low and high reflectance respectively. This refinement is a result of our lake bathymetry measurements. Equation (2) was applied to all pixels classified as water across the study region to yield depth estimates, which were subsequently integrated over the SGL area to calculate water volume.

Surface melt, pro-glacial discharge and catchment delineation
Surface runoff data were derived from a surface energy balance model using AWS measurements and MODIS (MOD10A1 product) derived surface albedo as input (van As et al., 2012). The model iteratively searches the surface temperature for which the energy fluxes are in balance. When this fails due to a limiting melting-point temperature, the surplus energy yields snow or ice melt. Model input data were interpolated to 100 m elevation bins to determine the distributed melt patterns in the region. Model uncertainty for runoff totals was calculated to be between 6-14 % (van As et al., 2012). Discharge measurements were recorded at a gauging station near the bridge crossing the Watson river at Kangerlussuaq (Fig. 1a). Discharge was derived using pressure transducers and a stage-discharge relationship, and at high stages with an acoustic doppler current profiler . The catchment described in Hasholt et al. (2013) has been refined in this study. The new catchment outline, hereafter referred to as the Russell Glacier catchment ( Fig. 1a and b) was derived by manually digitising the locations of all major supraglacial streams using panchromatic (15 m resolution) and multispectral (30 m resolution) ASTER and LANDSAT images between 2009-2012. Multiple images acquired at different times in the melt season were used to map streams in each year to maximise coverage of the network. Using this map of the supraglacial hydrological network the Fig. 3. (a) Lake F and Z MODIS band 1 reflectance vs. depth with best fit curve; (b) the relationship between lake area estimates derived from manually digitising 15 m resolution panchromatic LANDSAT 7 images and the automated MODIS classification; and (c) the relationship between independently measured lake depth for 3 SGLs and modelled lake depth using the MODIS depthreflectance relationship. Hasholt et al. (2013) catchment was edited, ensuring that the catchment boundary did not cross streams, and instead encompassed the stream network for each lake ( Fig. 1a

Area
We evaluated the accuracy of the classification by comparing the area of individual SGLs extracted from MODIS data with a SGL area data set derived from finer spatial resolution (15 m) panchromatic LANDSAT ETM+ images. In total, 62 lakes (free of surface ice) were manually digitized from four panchromatic LANDSAT ETM+ images from different years and at different stages of the melt season (26 June 2008, 17 July 2004, 20 July 2005, 8 August 2006. The comparison shows a strong agreement between the two data sets ( Fig. 3b) with strong correlation (r 2 = 0.97) and a root mean square deviation (RMSD) of 0.14 km 2 , comparable with previous studies (e.g. 0.11 km 2 , Selmes et al., 2011, and0.22 km 2 , Sundal et al., 2009). This RMSD value was multiplied by the number of SGLs occurring each day to provide error margins for total SGL area estimates.
The minimum area of any individual SGL identified from MODIS images is restricted by pixel size (0.0625 km 2 ) thereby introducing potential bias in total lake area and its size distribution. Previous work using higher resolution products calculated that small lakes (< 0.1 km 2 ) account for 12 % of the total lake area based on analysis of a single mid-season ASTER image (Sundal et al., 2009). In our study SGLs smaller than 0.0625 km 2 were identified from a LANDSAT ETM+ image (26 June 2008), indicating that they account for 7.5 % of the total SGL area at this time. From the same image, 11.1 % of SGLs were found to be smaller than 0.1 km 2 in line with the Sundal et al. (2009) study.

Volume
Depth surveys were undertaken on two SGLs located to the north of the study area on 14 and 15 August, 2005, referred to as Lakes A and B in Box and Ski (2007), and on Lake Disco on 22 July and 8-9 August, 2007. The results of these surveys were overlain on the corresponding MODIS derived depth maps enabling comparison (Fig. 3c). There is good agreement between the MODIS-derived depth and the raftmeasured depth (r 2 = 0.79, RMSD = 1.47 m). To calculate percent volumetric error, the point measurements from Lake Disco surveys were gridded and interpolated onto a 10 m grid, with the lake margin delimited from a Landsat 7 image acquired on the 8 August 2007. The 9.5 × 10 −3 km 3 calculated SGL volume agrees within 15 % of the volume from MODIS from 22 July (11.3 × 10 −3 km 3 ) and 8-9 August (11.0 × 10 −3 km 3 ) 2007. A further comparison between the recorded volumes of Lakes A (7.1 × 10 −3 km 3 ) and B (7.8 × 10 −3 km 3 ) from Box and Ski (2007), and the estimated volumes using the method described in this study (5.9 × 10 −3 km 3 and 6.7 × 10 −3 km 3 ) reveals a volumetric difference of 15.6 % and 13.7 %, respectively.

Cloud Cover
A key limitation to monitoring SGL evolution is the temporal and spatial continuity of available imagery. Cloud cover prevents the daily monitoring of SGLs and therefore the formation and drainage of some SGLs could have been missed. To limit the impact of clouds, SGL classification was only carried out on pixels that had a MODIS band 6 reflectance value of less than 0.15, consistent with Box and Ski (2007). Following automatic classification, images were manually inspected for cloud cover, and images which contained cloud that obscured the view of lakes were removed from the data set. This removal of cloudy images reduced the false positive identification of lakes due to cloud shadows, which in some cases have spectrally similar signatures to optically deep lakes. Although images included in the classification were selected when the Terra satellite was at or near nadir, the changing solar zenith throughout the course of the melt season meant that cloud shadows were present on some images. To eliminate this potential error each image was manually checked for the presence of clouds and their shadows. Furthermore, given that optically thin clouds reduce reflectance values in the red band -on which the depth reflectance function is applied -the removal of cloud-covered images is essential to accurately and consistently derive lake volumes. Cloud cover therefore limited the number of images used to evaluate SGL evolution each year (denoted as n in Fig. 4) and this constraint may introduce bias when identifying rapidly draining lakes. To identify such lakes we selected a four-day threshold as it provided the optimal interval for identifying fast draining lakes on a data set characterised by intermittent gaps due to cloud cover.

Prolonged blackout
Extended periods of cloud clover, termed blackouts, represent a problem for the analysis of lake drainage patterns as we are unable to quantify the amount of water stored in lakes which form and drain during these periods, therefore our results will likely underestimate lake storage and loss. During a prolonged blackout period with consistent cloud cover the degree of surface melting and its contribution to lake formation and expansion is unknown, an inherent limitation of this study. However, given that the average life cycle of a lake is ∼ 10 days, we note only 13 blackout periods which extend beyond this time frame between the start of May and the end of August, over the entire 11 yr analysis. Even given a perfect data set, an uncertainty exists regarding lakes loss/gain during the time interval between satellite passes and hence results of this analysis must be considered minimum estimates of area and volume. Analysis of supraglacial lake evolution using optical satellite imagery will always face the problem of blackout periods due to cloud cover and is therefore a key limitation of the method. Although some studies have demonstrated the use of SAR data for monitoring lake evolution due to its ability to penetrate clouds, it is not on its own a replacement, due to its poor temporal coverage and can only be used to complement other visible-near infrared data sets. Despite the restrictions imposed by cloud cover and relatively low spatial resolution MODIS imagery remains the best source of data available for monitoring SGL evolution.

Classification
The classification we applied is sensitive to the presence of water and does not identify lakes covered by floating ice or snow. Therefore using our technique, lake area and volume will be underestimated in regions or seasons with a high proportion of ice/snow-covered lakes.

Depth-reflectance technique
The accuracy of the volume calculation relies on the calibrated lake depths being representative of the whole SGL population. Lakes deeper than our calibration lakes may exist elsewhere on the ice sheet: the maximum depths of lakes F, Z, A, B and Disco were 16 m, 16 m, 11.5 m, 10 m and 15.2 m, respectively. Box and Ski (2007) note that there is a cut-off to the depth-reflectance relationship, beyond which increased depth fails to reduce pixel reflectance. The depth-reflection function itself limits the accuracy of the volume results as it relies on accurate calculations of atmospherically corrected surface-reflection and bathymetry measurements. Previous remote sensing studies have presented lake area estimates (e.g. Selmes et al., 2011), but Liang et al. (2012) go a significant step further by using a first-order, conical lake area to volume approximation to derive lake volume estimates. Tedesco and Steiner (2011) compared MODIS (band 3 and band 4) and LANDSAT (band 1 and 2) derived lake depths using the method of Sneed and Hamilton (2007) with in situ measurements. Tedesco and Steiner (2011) revealed that MODIS band 3 and 4 tends to overestimate and underestimate lake depth respectively when compared to in situ data. Here we used MODIS band 1 to derive lake depth because it is more sensitive to water depth than bands 3 and 4 (Box and Ski, 2007). To calibrate the relationship between reflectance and depth, we used bathymetry data from deeper lakes (up to 16 m, Fig. 2) than Tedesco and Steiner (2011), whose calibration lake had a maximum depth of 4.55 m. Tedesco and Steiner (2011) also found differences between the depth values derived from MODIS on board TERRA and AQUA, with values estimated by MODIS on TERRA being closest to those estimated by Landsat and measured in situ. On this basis we only used images from the TERRA satellite to derive depth estimates. We also acknowledge the limitations of deriving volumes from mixed pixels, for example at the lake margins, where the presence of ice lowers the overall pixel reflectance leading to an underestimation of lake depth and hence volume.

Results
Our aggregate data set reveals that an average of 200 SGLs (> 0.0625 km 2 ) are observed each year within the 6500 km 2 Kangerlussuaq sector of the GrIS (Fig. 1a). However, the total number of lakes varies from year to year (one standard deviation, ±20 lakes) and 21 % reoccurred annually (Fig. 1a). Mean area of an individual lake was 0.68 km 2 but they can grow as large as 8.2 km 2 (66.92 • N, 48.15 • W). SGLs formed 21 days earlier in 2010 compared to the 11 yr mean, with small lakes forming at low elevations by 9 May (Fig. 4 and Table 1). The maximum total volume of SGLs at any one time between 2002 and 2012 occurred on 23 June 2012 (DOY 175) when over 0.25 km 3 of water was stored across the Kangerlussuaq sector of the GrIS coincident with the warmest month in four decades (Cappelen et al., 2012).
The largest individual SGL volume was 0.038 km 3 on 14 July 2011 (DOY 195), during the second warmest July on record. Intra-annually, the timing of lake formation, periods of peak storage (Fig. 4) and lake-drainage events (Fig. 5) are highly variable, but it is clear that during high melt years (2010 and 2012), formation and drainage occurs much earlier in the season (regional maximum volume, Tv max , was attained 38 and 20 days earlier than the 11 yr average) than in comparatively low melt years (e.g. 2002 and2006;Figs. 4, 5 and 6a). The average maximum extent of lakes between 2002 and 2012 at each 200 m elevation interval illustrates that smaller lakes tend to form at lower elevations (< 1000 m, Fig. 6b). There is also large inter-annual variability in the fraction of rapidly draining lakes (defined here as those tapping within a 4 day period), with 28 % of lakes classified as such since 2002 (Fig. 7e). There is no obvious correlation with surface runoff, although the peak number of rapidly tapping lakes did occur in the peak melt year (2012). However, in record warm years, characterised by a greater number of cloud free MODIS images, it is possible that a sampling bias is introduced; that the greater number of rapid draining SGLs identified is a consequence of more complete temporal   lan et al., 2007;Liang et al., 2012;Doyle et al., 2013). Using 2010 as an example (Fig. 8a), lakes can be seen storing and draining meltwater within five elevation bands. During 2010, 83 % of SGL volume was lost in lakes located between 1000 and 1600 m elevation (Fig. 8b) with only 9 % and 7 % of water lost from lakes below 1000 m and above 1600 m elevation, respectively. Although most lakes form in the same locations each year, there is significant inter-annual variability in the maximum areal extent of SGLs within each 200 m elevation band (Fig. 9a). The maximum areal extent of lakes during relatively high temperature summers was 158 km 2 in 2003, 154 km 2 in 2007, 136 km 2 in 2010, 160 km 2 in 2011 and 160 km 2 in 2012, significantly larger than in other years. For example, SGL maximum areal extent was 40 % greater in the record melt year of 2012 compared with the low melt year of 2006. Furthermore, during 2012 lakes occupied a larger area (49 %) above 1400 m (96.6 km 2 ) than in previous years (11 yr average of 64.7 km 2 ) (Fig. 9a). The years 2003 and 2007 gave the second and third highest elevation lake extents above 1400 m a.s.l., 33 % and 25 % greater than the 11 yr mean, respectively. The maximum areal extent of lakes above 1400 m for each year in the period 2002-2012 was compared with mean monthly temperatures of June and July (Cappelen et al., 2012), revealing a positive correlation (r 2 = 0.78) supporting findings of SGLs forming at higher elevations in the highest temperature years (Fig. 9b). With the exception of the lower temperature 2002 and 2006 melt seasons, there is a significant temporal trend in both lake volume loss and the maximum elevation of lakes over the last decade (Fig. 7d, f).
Seasonal patterns of SGL meltwater storage closely correlate with melt rates derived from an energy balance model (van As et al., 2012) for 2009/10 (Fig. 4). Peaks in surface meltwater runoff and storage in lakes generally precede increases in discharge recorded from the Watson River (van As et al., 2012;Hasholt et al., 2013). Volume lost from lakes was compared to peaks in the river discharge record to determine if specific SGL-drainage events can be identified (Fig. 5).
Our results indicate that SGL drainage does coincide with short-term perturbations in Watson River discharge, providing tentative evidence for the impact of changing supraglacial water storage on catchment-wide hydrology.
Lake-drainage maps (Fig. 10) illustrate the seasonal lake evolution, highlighting the spatial and temporal variations of rapidly (< 4 days) draining and relatively slow (> 4 days) draining lakes. In general, SGLs within the study area fit consistent drainage patterns each year, remaining rapid or slow draining and not switching between years. There is no obvious relationship between lake size (volume/area) and its drainage mode; fast or slow. In some parts of the study area, SGLs were observed to drain in clusters, i.e. lakes within close (< 6 km) proximity were observed to drain on the same day or within a couple of days of each other (Fig. 11), possibly impacting on the proglacial discharge hydrograph (Fig. 5).

Inland expansion of supraglacial lakes
Using a field-calibrated and independently validated depthreflectance function, we tracked the area and volume of 200 ± 20 individual SGLs in 11 melt seasons (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) across the Kangerlussuaq sector of the GrIS (Figs. 1a, 7a and The Cryosphere, 8, 107-121, 2014 www 10). The recurrence of SGLs in the same positions each year (Fig. 1a) is in agreement with previous work (Thomsen et al., 1988;Echelmeyer et al., 1991) and indicates that the lakes do not advect down-glacier with ice flow. Lampkin and Vanderberg (2011) and Sergienko (2013) demonstrated that ice surface topography is coupled to underlying basal topography, confirming that most SGLs form above the same subglacial basins each year. Within our study area, the persistent SGLs (21 % of the population) are preferentially located away from the fast-flow units as identified by Palmer et al. (2011). The largest lake area observed within the study area was 8.25 km 2 , coinciding with a maximum volume of 0.038 km 3 , comparable to findings by Box and Ski (2007), although larger lakes have been observed in other studies (16.8 km 2 by Selmes et al. (2011)). On average, the maximum volume of any single lake (Lv max ) occurs on 25 July (DOY 206), ∼ 12 days after the regional maximum volume (Tv max ) (Table 1; Fig. 4) because larger (area) lakes tend to form at higher elevations, later in the season and also take longer to fill (Fig. 6b). The total number of SGLs within any year shows no relationship with Tv max and Lv max (Fig. 7b, c). For example in 2012, 200 lakes were observed, however Tv max was 80 % greater than the 11 yr mean, and SGLs formed at higher elevation than in all previous years (Fig. 7d). Lake area varies with elevation, with smaller lakes forming at low (< 1000 m a.s.l.) and high elevations (> 1600 m a.s.l.), and the largest lakes forming in between (Fig. 6b). The majority (84 % in 2010) of water stored in lakes occurs within the 1000-1600 m elevation band (Fig. 8b), reflecting the effect of ice sheet hypsometry on lake morphology. We speculate that steeper slopes and crevassing associated with fast flow at lower elevations within this region Fitzpatrick et al., 2013) prevents lakes from growing as large before drainage occurs. At higher elevations (above 1600 m), lower temperatures and the shorter melt season limits meltwater availability and lake growth. However, SGL expansion in this upper melt zone is sensitive to positive temperature and melt anomalies and may be used to indicate longer term climate trends (Howat et al., 2013). Availability of surface runoff, governed in part by surface temperatures and winter snow pack, determines the initiation of lake formation (Fig. 4, Johansson et al., 2013). The timing and duration of individual lake-drainage events varies interannually (Fig. 10), explaining the large differences between the total volume of water stored in lakes in each year (Fig. 4). As previously stated, during record warm years (2010 and 2012) SGLs formed earlier and attained their Tv max on DOY 156 and DOY 174, 38 and 20 days earlier than the 11 yr average, respectively ( Fig. 4 and Table 1). The record melt of 2010 was exceeded in 2012 during which the largest melt extent within the satellite era (∼ 97 % of the ice sheet) was observed (Nghiem et al., 2012). During 2012, the melt season lasted up to two months longer than the 1979-2011 mean (Tedesco et al., 2013), and although we recorded no increase in the number of lakes (Fig. 7a), the total volume of SGL (Tv max ) was ∼ 80 % greater than the 11 yr mean and lakes were observed at higher elevations (> 1800 m) and occupied a greater area > 1400 m (96.6 km 2 ) than in previous years (Figs. 7b, d and 9a Lake-drainage events have contributed directly to hydrograph anomalies in the discharge record from Watson River. For example, the rapid drainage of a number of large SGLs between DOY 175-184 (24 July-3 August) during 2010 (Fig. 5, Doyle et al., 2013) coincided with a peak in discharge without a corresponding peak in melt (Figs. 4 and 5).
Our data reveal considerable inter-annual variation in water storage and loss (Figs. 4 and 6a), with no two years following the same pattern. The rapid formation and drainage of lakes during the record melt year of 2012 strongly contrasts with the slow formation of lakes during the relatively cool year of 2006 (Fig. 6a). Assuming that all lakes within the study area encompassing the Kangerlussuaq sector of the GrIS region drain via the Watson River, a comparison between cumulative volume loss and discharge reveals that lakes temporarily store and release on average 13 % of total discharge (Table 2). However, if only lakes within the boundaries of the defined Russell Glacier catchment are considered to drain via the Watson River, the mean lake contribution to discharge is 7 % (max 11.5 %). Given the uncertainties associated with catchment delineation and deriving lake volumes, the estimates of lake storage are in line with modelling results over this region (Leeson et al., 2012), which have shown that supraglacial lakes stored a maximum of 12 % of total runoff during 2003. Using 2011 as an example, lake storage within the Palmer et al. (2011), van As et al. (2012 and Russell Glacier catchment was 3.5 %, 8.7 % and 7.6 %, respectively (Table 2). Previous studies note the uncertainties associated with defining a catchment for a single outlet from the GrIS (van As et al., 2012;Clason et al., 2012)  is probable that water piracy occurs between catchments as catchment boundaries change depending upon the state of the supraglacial and subglacial systems. Under a future warming climate (Meehl et al., 2007), SGLs are expected to expand further into the ice sheet interior as melt-season extent and duration increase (Howat et al., 2013). We present evidence supporting previous work (Liang et al., 2012;Howat et al., 2013), illustrating that in higher temperature years SGLs formed at higher elevations (Figs. 7d and 9b). Lakes occupied a significantly larger area at higher elevations (> 1400 m) during 2012 (Fig. 9a) than in previous years (49 % greater than the 11 yr average). This coin-cides with the highest mean June/July temperatures recorded in Kangerlussuaq (Fig. 9b) and the greatest Greenland melt extent on record (Nghiem et al., 2012). Furthermore, we find a positive correlation between the maximum areal extent of lakes above 1400 m a.s.l. and mean monthly temperatures of June and July (Fig. 9b). Preferential absorption of solar radiation by SGLs, causes expansion of their existing basins (Box and Ski, 2007;Tedesco et al., 2012), creating a positive albedo feedback and encouraging SGLs to reform in the same place in successive years. Despite the high melt season of 2010, the maximum extent of SGLs above 1400 m only matched the 11 yr mean (Fig. 9a). However, during 2010, The Cryosphere, 8, 107-121, 2014 www.the-cryosphere.net/8/107/2014/ SGLs in each of the lower elevation bands were larger than all other years (Fig. 6b), indicating that the large volume of melt in 2010 primarily occurred across lower elevations. A similar pattern occurred in the warm summer of 2003, where SGLs occupied a below average extent above 1600 m (Fig. 9a), but the average maximum lake size between 1400-1600 m was 60 % greater than 11 yr mean (Fig. 6b).

Rapid lake-drainage events
Rapid SGL-drainage events provide a mechanism for surface water to penetrate through kilometre-thick ice, establishing hydraulic pathways which meltwater continues to drain through the remainder of the melt season Doyle et al., 2013). Over the 11 annual melt periods of this study, on average 28 % of the SGL population rapidly drained (Figs. 7e and 10). As the number of SGLs occurring higher up-glacier increases (Figs. 7d and 9a), more SGLs will rapidly drain, potentially expanding the region affected by melt-induced accelerated flow. We also note that at higher elevations, a number of SGLs appear to drain twice in higher melt years. For example, in 2012 a 5.2 km 2 SGL located at 1400 m drained completely on 23/24 June, subsequently refilled and by 5 July had completely drained again. Lake refilling has been observed in the field in automatic camera imagery (Brizgys and Box, 2005). Such observations indicate that the subglacial hydrological system may not readily adapt to surface meltwater inputs as has recently been postulated (e.g. Sundal et al., 2011). This suggests that basal hydrostatic creep-closure rates across the upper ablation zone where ice thickness is greater than 1 km exceed those processes which help to maintain an efficient hydrological system for the remainder of the melt season (Schoof, 2010;Sole et al., 2013). Spatial and temporal clustering of SGL-drainage events (Figs. 10 and 11) provides evidence that neighboring SGLs dynamically trigger synoptic-scale tapping. Box and Ski (2007) calculate that a lake containing 0.0983 km 3 of water, upon drainage and without being laterally constrained will spread over a subglacial area of 97 km 2 , assuming a void thickness of 1 m. We suggest that the drainage of neighbouring SGLs acts to mechanically trigger one another via short-term perturbations in the regional stress/strain regime transmitted over length scales on the order of tens of kilometres. Large surface rifts and fractures have been observed in SGLs that have drained rapidly Doyle et al., 2013), and are the primary mode of drainage identified in these studies. Many SGLs across the GrIS appear to drain in this manner each year and we find no common critical threshold of lake volume or depth which appears to control drainage initiation. Based on our results, we hypothesise that the rapid drainage trigger mechanism is site specific, governed by the interplay of seasonal routing and delivery of meltwater moderated by the regional topographic setting and stress/strain regimen together with the inherent mechanical resilience of the ice surface in the locality.

Conclusions
Using 502 cloud-free atmospherically corrected MODIS images, we derived water volumes from ∼ 200 seasonally occurring SGLs between 2002 and 2012 using a field-calibrated and validated depth-reflectance relationship. We find smaller lakes form at lower elevations (< 1000 m) and larger area lakes form at higher elevations later in the season, taking longer to fill. Although SGLs only occupied a relatively small region of our study area (2 %), they store a disproportionate volume of bulk runoff (7-13 %) and have important implications for ice dynamics through the release of surface meltwater into the subglacial hydrological system via rapid in situ drainage or through overflow into moulins. We find evidence of SGLs draining in clusters, causing short-term perturbations in discharge gauged at the Watson River, and infer that one drainage event dynamically triggers rapid tapping in neighbouring lakes. We find that lake size does not influence its drainage mechanism, and that there is no evidence for a critical lake depth or volume threshold to initiate rapid drainage.
In years of high summer temperatures, SGLs form and drain earlier in the season (e.g. 38 and 20 days earlier than the 11 yr mean during 2010 and 2012, respectively), and cover a larger surface area (e.g. 40 % greater in the record melt year of 2012 compared with the cooler year of 2006). Furthermore, inland expansion of SGLs is strongly correlated with air temperature (r 2 = 0.78), with lakes occupying a greater area within the upper ablation zone during warmer years. For example, in 2012, lake area was 49 % greater above 1400 m compared with the 11 yr mean and extended further inland (> 1800 m) than previously recorded. In a warming climate, spatial and temporal expansion of SGLs with concomitant surface to bed coupling may impact strongly on inland ice sheet flow dynamics and interior mass drawdown.