A ten-year record of supraglacial lake evolution and rapid drainage in West Greenland using an automated processing algorithm for multispectral imagery

The rapid drainage of supraglacial lakes introduces large pulses of meltwater to the subglacial environment and creates moulins, surface-to-bed conduits for future melt. Introduction of water to the subglacial system has been shown to affect ice flow, and modeling suggests that variability in water supply and delivery to the subsurface play an important role in the development of the subglacial hydrologic system and its ability to enhance or mitigate ice flow. We developed a fully automated method for tracking meltwater and rapid drainages in large (> 0.125 km 2) perennial lakes and applied it to a 10 yr time series of ETM + and MODIS imagery of an outlet glacier flow band in West Greenland. Results indicate interannual variability in maximum coverage and spatial evolution of total lake area. We identify 238 rapid drainage events, occurring most often at low (< 900 m) and middle (900–1200 m) elevations during periods of net filling or peak lake coverage. We observe a general progression of both lake filling and draining from lower to higher elevations but note that the timing of filling onset, peak coverage, and dissipation are also variable. Lake coverage is sensitive to air temperature, and warm years exhibit greater variability in both coverage evolution and rapid drainage. Mid-elevation drainages in 2011 coincide with large surface velocity increases at nearby GPS sites, though the relationships between ice-shed-scale dynamics and meltwater input are still unclear.


Introduction
The Greenland ice sheet's (GIS) contribution to sea level rise between 1993-2003 was estimated at ∼ 0.21 ± 0.035 mm yr −1 by the Intergovernmental Panel on Climate Change's (IPCC) fourth assessment report (Lemke et al., 2007) ter velocity (Zwally et al., 2002;Joughin et al., 2008;Hoffman et al., 2011), and recent data suggest that variability in meltwater input timing and mechanism exaggerate flow effects (Sundal et al., 2011;Hoffman et al., 2011).During the melt season, the lower (< 1000 m a.s.l., 65-70 • N) Greenland ablation zone is rife with moulins that drain individual watersheds ranging from a few m 2 to km 2 (field observation) and provide nearly continuous inputs to the subsurface.While many lakes feed the subsurface via streams and moulins, large lakes, those with diameters upwards of ∼ 250 m (Krawczynski et al., 2009), have the potential to drain rapidly to the ice sheet bed through hydrofracture.These drainages cause localized acceleration over short (∼ 24 h) time scales (Das et al., 2008;Hoffman et al., 2011) and may also have dramatic effects on the subglacial hydrologic system, speeding the transition between cavity and channel configurations as described by Schoof (2010).In areas of thicker ice (> 1000 m ice thickness), moulins produced by lake hydrofracture may be particularly important for meltwater transport because smaller moulins and cracks can no longer reach the bed (Van der Veen, 2007; Krawczynski et al., 2009).Because of their ability to create conduits, provide large meltwater pulses, and potentially alter changes in the subglacial hydrologic system and ice flow, this study focuses on large (> 0.125 km 2 ) supraglacial lakes and their drainages.
Past studies have used visible and near-infrared satellite imagery to identify and investigate the properties of supraglacial lakes using a variety of methods; however, many employ either manual lake identification (McMillan et al., 2007;Lampkin and Vanderberg, 2011) or automated lake identification using manually adjusted identification parameters for each image (Box and Ski, 2007;Sundal et al., 2009).These studies also require manual selection of cloud-free scenes and removal of false positives due to cloud shadows or other illumination anomalies.Two more recent studies (Selmes et al., 2011;Liang et al., 2012) take advantage of the contrast between lake water and the surrounding ice in MODIS band 1 (620-670 nm, red) to identify lakes and use automated methods for image selection and processing.
A fully automated processing algorithm is vital to this study because we aim to produce a ten-year surface coverage and drainage record, including analysis of over 100 Introduction

Conclusions References
Tables Figures

Back Close
Full Landsat 7 ETM+ and over 8500 Aqua and Terra MODIS scenes over West Greenland.We automatically remove low resolution and cloudy images and filter individual lake results for localized cloud cover.We present the results from 802 lake maps of the study area-output density varied from 50 useable scenes in 2002 to 109 in 2009, with a mean of 73 MODIS and 7 ETM+ scenes per year.Our results include both watered surface area and depth results for the entire study area, but we choose to focus on 78 large (> 0.13 km 2 ), deep (> 2 m) lakes.This lake coverage and rapid drainage record aims to better describe the major surface meltwater reservoirs within our study area and the rapid transfer of this meltwater to the ice sheet bed using a fully-automated image processing algorithm for a variety of multispectral imagery sources.This processing includes image selection and filtering based on cloud cover and viewing geometry as well as surface water detection and removal of false positives caused by cryoconite and other ice surface anomalies.

Study area
The study area is a 20 km wide swath, centered on the Sermeq Avangnardleq flow band, just north of Jakobshavn Isbrae, around 70 • N in west Greenland, and extends 80 km inland from the terminus (Fig. 1).Our lake coverage and drainage record focuses on 78 large lakes along the flow band, selected for their large surface area (> 0.125 km 2 ) and maximum depth (> 2 m, see Sect.tures any drainage events that impact ice dynamics within the ROGUE field area (Price et al., 2008;Hoffman et al., 2011).

Imagery
Our study employs imagery from the Enhanced Thematic Mapper Plus (ETM+) aboard Landsat 7 and Moderate Resolution Spectroradiometers (MODIS) on Aqua and Terra.
The ETM+'s higher spatial resolution (30 m) provides a spatially detailed albeit temporally coarse record (16 day revisit period).While MODIS's coarser spatial resolution (250 m) limits reliable classification to lakes with areas of 0.125 km 2 or larger, subdaily satellite overpasses allow for much more sensitive analyses of temporal change.
The identification algorithm uses red and near-infrared wavelengths, corresponding to ETM+ bands 3 (630-690 nm) and 4 (750-900 nm) and MODIS bands 1 (620-670 nm) and 2 (841-876 nm), respectively.We present results from imagery captured from 1 May to 30 September in 2002-2011.We exclude images with greater than 10 % study area cloud cover, determined from the Automatic Cloud Cover Assessment algorithm (Irish et al., 2006) for ETM+ and the MOD/MYD35 Cloud Mask product for MODIS.We also limit analysis to MODIS scenes captured during midday (11:00-16:00 LT) where the study area is within 30 • of nadir according to the MOD/MYD03 geolocation product to avoid issues with viewing geometry and illumination.Because of their 2330 km-wide swath, MODIS 250 m resolution images exhibit effective resolutions as coarse as 750 m near the edges of each image.When analyzing relatively small, high contrast features such as lakes, this resolution degradation can result in misclassifications.While small lakes (those nearer nominal resolution) may not be detected at all, the effective averaging that results from this resolution loss can result in wide-scale water classification where none exists.Filtering those data that fall more than 30 • from nadir limits our analyses to those pixels with dimensions no greater than 333 m in the across-track direction and reduces the likelihood of misclassification.Introduction

Conclusions References
Tables Figures

Back Close
Full

Surface water identification
We convert ETM+ and MODIS level 1 calibrated radiances to top-of-atmosphere (TOA) reflectances and reproject to the WGS-84 datum, UTM Zone 22. Surface water is identified based on a normalized difference lake index (NDLI), where R r and R nir are the red (ETM+ band 3, MODIS band 1) and near-infrared (ETM+ band 4, MODIS band 2) reflectances of any pixel, respectively.Liquid water absorbs very highly in the near infrared and can be used to delineate pooled water from the surrounding ice, which reflects highly in the infrared; however, degradation of the ice surface due to ablation and cryoconite concentration throughout the melt season drastically decreases infrared reflectivity resulting in overclassification of surface water.Because water does not absorb as highly in the visible, the red bands are not as susceptible to this phenomenon, though they also exhibit less spectral separability between ice and water.To accurately classify surface water with either band alone requires a different threshold for each image.This normalized difference index better captures the differences between liquid water and the surrounding ice and allows us to use a single threshold value for all images analyzed.To determine a threshold, we examine the distribution of NDLI at different points during the melt season (Fig. 2).The three dates in −10 % or +10 % in all ETM+ 2011 images resulted in +9.3 ± 3.7 % and −6.2 ± 1.1 % change in coverage area, respectively.

Lake depth analyses
To determine lake depth and volume, we use a modified version of the method described by Sneed and Hamilton (2007).Lake depth z can be determined by where A d is the lakebed albedo or reflectance in the red band, R ∞ is the red reflectance of optically deep water, R r is the red reflectance of the lake pixel being analyzed, and g is a constant describing the optical properties of the lake water.As R ∞ requires optically deep (> 40 m, Sneed and Hamilton, 2007) water, our estimation is limited to water in nearby fjords and Disko Bay, which may exhibit different scattering and absorption coefficients than those described above.We take R ∞ to be the reflectance of a mid-Bay pixel, free of ice, sunglint, and outwash plumes and consistent in both ETM+ and MODIS images from 3 July 2011.While g represents the characteristics of both upwelling and downwelling beam attenuation, Philpot (1989) estimates that g can be related to the diffuse attenuation coefficient for downwelling light, K d , as 1.5 K d < g < 3K d , with lower values coinciding with absorption-dominated water bodies, such as supraglacial lakes.K d is calculated from the absorption, a w , and scattering, b m , coefficients for pure freshwater (Smith and Baker, 1981; Table 1), at the appropriate sensor band wavelengths.We calibrate our estimate of g using field data from the lake near GULL (Fig. 1) by comparing the ETM+ derived depth from 3 July 2011 during the lake's filling phase to the measured depth from 28 July 2011

TCD Introduction Conclusions References
Tables Figures

Back Close
Full when the lake shore was at the same location during its gradual draining.We assume that the depths are equivalent within the resolution of the imagery.We use g = 0.86 or 2.19 K d , appropriate for the relatively clear, absorption-dominated glacial meltwater.
To select the focus lake basins, we determined the maximum depth for each pixel in the study area from ETM+ imagery between 2002-2011.Filtering out those areas with < 2 m maximum depths and, subsequently, those lakes with < 0.125 km 2 extent produced 78 distinct lakes (Fig. 1).To define our test basins, we include a 500 m radius buffer around each lake's maximum extent.Because of this buffer, the areal extents of the focus basins are significantly larger than the ten-year maximum extent of the lakes that occupy them.We filter basin-full area measurements from our surface water inventory.Such measurements represent the coverage of each basin (and of entire regions of the ice sheet) by shallow water and wet ice and, because this occurs at low elevations and late in the season, generally do not coincide with the presence of lakes.

Drainage event identification
The rapid drainage of supraglacial lakes is a catastrophic occurrence, resulting in the disappearance of the lake over a short time period.It is important to note that the time step between each image is not uniform, and because of this, assessing rapid drainage within any timing window is reliant on sufficient data density.While rapid drainage events can occur over a matter of hours (Das et al., 2008), we extend the identification window to six days to better accommodate the temporal resolution of the coverage record and assume that slower drainages (overland to moulin) occur on weekly to monthly time scales.Despite this concession, image availability prevents reliable drainage detection during some time periods.To identify rapid drainage events, we compare local maximums in lake area, A p , to subsequent lake areas, A n , within a six-day timing window.If A n is more than six days after A p , the window shifts until two points within a six day period can be compared.To qualify as a rapid drainage, a lake must lose critical area within a 6 day period; we define critical loss as 90 % of the year's maximum area, (A p − A n ) > 0.9A max .To prevent false positives, we require that the lake

Conclusions References
Tables Figures

Back Close
Full not refill within six days after the drainage.The drainage magnitude, A d , is the lake area at the local peak, A p , and the drainage time, t d , is the midpoint between the time at local maximum, t p , and time of the critically low area, t n .We define the uncertainty in drainage time as the range between these values.This is summarized in Eq. ( 4): Rapid drainages that coincide with rapid filling in a nearby basin are automatically filtered from the drainage record because they do not introduce water to the bed.A list of rapid drainage events can be found in Supplement.

Surface velocity from GPS
The ROGUE project maintains nine GPS (Trimble NetR7 with Trimble Zephyr Geodetic antennas) stations within the study area ranging from 650-1200 m.We present data from two of these stations, those nearest detected rapid drainage events in 2011.
GPS positions were determined by carrier-phase differential processing relative to a bedrock-mounted base station using Track v.1.22and final International GNSS Service satellite orbits (complete methodology described in Hoffman et al., 2011).

Volume vs. area (ETM+)
Due to the difficulty of remotely sensing lake depth and the uncertainties in this estimation, it has become commonplace to use lake surface area as a proxy for volume.
If each lake is modeled as an inverted cone, lake volume varies linearly with the ratio Introduction

Conclusions References
Tables Figures

Back Close
Full of diameter (D) to depth (z) and exponentially with surface area (Krawczynski et al., 2009;Liang et al., 2012).Because the ice sheet surface is relatively smooth at 250 m resolution and lakes tend to form in persistent depressions, low-amplitude expressions of bed topography (Echelmeyer et al., 1991), there is likely little variation in D/z among lakes, especially compared to variations in surface area (0.13 to 4.61 km 2 in the study area).When assessing all water within the study area, collectively, we find volume to be highly correlated with surface area within our predefined lake basins (r 2 = 0.896, z = 1.06 m) and take lake area to be an appropriate proxy for volume.This correlation suggests an average depth of just over a meter within the lake basins, consistent with large, relatively shallow lakes.We find no significant correlation between volume and surface area when including out-of-basin waters.We attribute both the decreased average depth and increased variance to widespread shallow standing water (r 2 = 0.443, z = 0.23).This shallow water is particularly noticeable during the second half of the melt season at lower elevations (< 900 m) and can result in anomalously high surface area estimates that are not indicative of large volumes.

Surface water area inventory
Our water inventory results include surface area data from 730 MODIS and 72 ETM+ scenes, track lake evolution through 10 yr, and identify 238 rapid drainage events.Of the 78 lakes we tracked, 73 drain rapidly at least once during the study period, though many drain catastrophically in fewer than half the years they filled (Fig. 3).Both lake occurrence and rapid drainage decrease with elevation; while fewer full lakes afford fewer opportunities for rapid drainage, we attribute the decrease of rapid drainage with elevation to increased ice thickness-ice at our GULL boreholes (888 m a.s.l.) is over 700 m thick.Despite this, several of the highest elevation lakes drain catastrophically in the majority of years that they fill.While greater ice thickness poses a barrier to hydrofracture, it also produces larger surface basins and greater lake areas and volumes.
Above 1300 m a.s.l., increased lake size appears to overcome increased ice thickness, allowing lakes to create conduits to the ice sheet bed.Introduction

Conclusions References
Tables Figures

Back Close
Full Results show high temporal variation in lake filling onset (Fig. 4), ranging from mid May (2010) to early July (2006).Lake filling occurs first at low elevations and proceeds upward, though in some years wide elevation ranges experience filling onset almost simultaneously.Lower elevations tend to lose lake area quickly, often before peak area has been reached in the other elevation bands.The middle elevation band (900-1200 m, blue in Fig. 1) hosts the most lakes but has similarly rapid turnover in all years.While there are fewer lake basins at high elevation and those lakes fill during fewer years, high elevation lakes can dominate coverage because of their large size, particularly later in the season when low and mid elevation lakes have disappeared.Widespread high elevation lake coverage is most prominent in 2004, 2007, and 2011, though notably few high elevation rapid drainage events occur in these years (4,4,0, respectively).Maximum lake coverage ranges from 14.9 (2005) to 37.6 km 2 (2007), though this accounts for only 20-50 % of the maximum potential coverage (72.8 km 2 ), calculated from the greatest extent of each lake in any given year.The large disparity between coverage and the maximum potential coverage is due not only to temporal and spatial variations in lake filling but also to variability in the mode and timing of lake drainage.
We observe large temporal variability in the timing of peak lake coverage and in coverage evolution.Rapid drainages occur most often during periods of net filling or near peak coverage, though neither the number of drainages nor their magnitude is closely related to coverage either lake-by-lake or regionally.Years with temporally narrow coverage peaks (2002,2003,2004,2006,2009, and, within the middle elevation band, 2011) exhibit temporal clustering of rapid drainage events at a wide range of elevations around the coverage peak.While these rapid drainages contribute to the subsequent loss of lake coverage, they are accompanied by significant lake area lost to non-rapid, overland lake drainage through moulins.In years without this temporal drainage clustering, lake coverage is lower in magnitude and more evenly distributed throughout the season (2005,2008,2010).In these years, rapid drainage events begin soon after filling onset in each elevation band and may prevent further lake filling.

Conclusions References
Tables Figures

Back Close
Full

Air temperature analyses
To explore relationships between local weather and lake coverage, we compare our record to air temperature measurements from the Greenland Climate Network's automated weather station (AWS) at Swiss Camp (Steffen et al., 1996).We choose Swiss Camp as a representative site because of the ready availability and reliability of AWS data and its centralized location (near the study area's center-line at 1176 m elevation).While it does not implicitly represent either low or high elevation conditions, the behavior, if not the magnitude, of cumulative positive degree-days (PDDs) at Swiss camp is representative of all of our elevation bands.Total lake area is sensitive to air temperature in the early season (Fig. 5); low elevation lakes begin to form just prior to the onset of above-zero temperatures at Swiss Camp in each year and prolonged periods of positive temperatures result in gains in total lake area.After the initiation of lake filling, the relationship between lake area and air temperature becomes less clear.In many of the cooler years, those years with fewer than ∼ 80 accumulated positive degree days, peak lake coverage occurs during the first prolonged period of above-zero temperatures.
In 2002, 2003, and 2009, total lake area fades after temperatures return to zero and is not regained despite high temperatures at the end of the season.This results in a relatively narrow coverage peak in most cases, and rapid drainage events tend to cluster around that peak.Warm years exhibit more variable behavior.In 2009 when abovezero temperatures began in mid-May and persisted through September resulting in the highest observed total PDDs (167), total lake area remained at consistently low levels at all elevations throughout the season.This combination of high PDDs and small total lake area also occurred in 2005, while the other two high PDD years (2007,2011) exhibit extensive, prolonged lake coverage.In all but one case, warm temperatures after 1 August do not contribute to total lake coverage (GC-Net data was not available to investigate the dual coverage peak in 2004), likely due to more efficient transfer of surface melt to the subsurface by existing streams and moulins.Introduction

Conclusions References
Tables Figures

Back Close
Full

Rapid lake drainage-induced ice flow
We observe the coincidence of rapid drainage events and surface speedups in 2011 at the two sites nearest the detected drainages (Fig. 6).Due to temporal uncertainty in rapid drainages, the specific contributions of each drainage to the surface velocity record is unclear, though drainages appear to coincide with the three of the four major, short-duration, mid-season speedups similar to those seen by Hoffman and others (2011).Most noticeable is the clustering of drainage events around the 1 July surface velocity peak.The two largest drainages around the peak occurred within 10 m of 33N0's (Fig. 6a) elevation, and that station recorded a larger surface positive velocity anomaly during the 1 July speedup than the station at GULL as well as a more pronounced slowdown following peak velocity (Fig. 6b).This slowdown is perhaps indicative of local reorganization of the subglacial hydrologic system following the nearby drainages, once the water supply begins to diminish.The final peak (11 July) may be the result of flow coupling associated with drainage outside the study area or a rain event since we observe no coincident rapid drainages; however, the number of usable images limits reliable drainage detection after 5 July.

Conclusions
We present a fully-automated system for identifying and tracking supraglacial lakes in both ETM+ and MODIS imagery and use it to track 78 large lakes along Sermeq Avangnardleq between 2002 and 2011.We observe wide variations in filling onset and lake coverage evolution.While filling generally progresses inland and upward throughout the season, we observe interannual differences in filling onset date in low, mid, and high elevation zones of the study area ranging from a few days and over a month.Rapid drainage events are observed during all phases of filling but most often near the seasonal coverage maximum or, in years with strongly elevation-staggered filling, near local peaks in lake coverage for each elevation range.Our limited surface ve-Introduction

Conclusions References
Tables Figures

Back Close
Full  PDD years (2003PDD years ( , 2009)).High PDD years may exhibit similar peak coverage and temporally localized drainage events (2011) or highly variable coverage and sporadic drainages (2010).
From 2003-2007, estimates increase to an average ∼ 0.5 mm yr −1 (Cazenave and Llovel, 2010).Surface melting and calving are each responsible for about half of this mass loss (Van den Broeke et al., 2009), though they are not independent.The introduction of meltwater to the base of the ice sheet can accelerate flow through sliding by up to 200 % of the win-Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2.4).The Real-time Observations of the Greenland Under-ice Environment (ROGUE) field program targeted this area in 2010-2012, installing nine surface GPS stations in the region to monitor surface velocity.In 2011, the ROGUE drilling campaign instrumented eight boreholes to monitor ice motion, ice temperature, and subglacial hydrology, installed automated weather stations, and collected surface hydrology data at two of the GPS station sites (FOXX: 69 • 26 47 N, 49 • 52 43 W; GULL: 69 • 27 9 N, 49 • 42 52 W).Our chosen swath cap-Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 2
represent typical early, mid, and late season coverage states, respectively.NDLI gradually increases throughout the melt season as the surface transitions from snow to bare ice.Meltwater appears in the upper tail of the distribution with NDLI increasing with water depth.To avoid misclassifying shallow water and wet snow or ice as lake water, we place our threshold near the end of the distribution tail during peak coverage (Fig.2b), when, by inspection, most of the standing water is in lake basins.We attribute the late season upward shift in the distribution to widespread shallow water and wet ice in the lower ablation zone that does not necessarily indicate an increase in lake volume.A threshold of NDLI = 0.21 was used for all images.Adjusting the threshold by Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |