1. Introduction
Tidewater glaciers – glaciers terminating at a marine interface – undergo both melting and iceberg calving, and both mechanisms transfer significant land-ice mass to the oceans. Ice loss from a tidewater glacier generally can be expressed as the water equivalent sum of surface runoff, submarine melt, volume loss due to retreat of the glacier terminus and volume loss due to iceberg calving flux. The iceberg calving flux is strongly dependent on the glacier flow rate, which is often assumed to be constant for estimates of regional total mass balance (Reference DowdeswellDowdeswell and others, 2002, Reference Dowdeswell, Benham, Strozzi and Hagen2008; Reference Burgess, Sharp, Mair, Dowdeswell and BenhamBurgess and others, 2005; Reference GardnerGardner and others, 2011). However, tidewater glaciers may experience significant flow variability at seasonal and annual timescales and this may introduce uncertainty in estimates of ice mass discharge to the oceans. Reducing this uncertainty requires an understanding of the causes of flow variability on tidewater glaciers.
Some tidewater glaciers exhibit seasonal flow variability (Reference Vieli, Jania, Blatter and FunkVieli and others, 2004) similar to that exhibited by High Arctic and alpine valley glaciers (Reference IkenIken, 1972; Reference Copland, Sharp and NienowCopland and others, 2003; Reference Sugiyama and GudmundssonSugiyama and Gudmundsson, 2004; Reference Bartholomaus, Anderson and AndersonBartholomaus and others, 2008). This variability may be missed by some methods of measuring glacier flow due to seasonal biases in their temporal sampling. For example, while feature-tracking methods applied to optical satellite imagery can measure average annual ice displacements (Reference Williamson, Sharp, Dowdeswell and BenhamWilliamson and others, 2008), interferometric synthetic aperture radar (InSAR) and speckle tracking (Reference Burgess, Sharp, Mair, Dowdeswell and BenhamBurgess and others, 2005; Reference Van Wychen, Copland, Gray, Burgess, Danielson and SharpVan Wychen and others, 2012) can measure approximately monthly displacements during fall, winter and spring. InSAR methods have also been used successfully in summer (Reference Joughin, Das, King, Smith, Howat and MoonJoughin and others, 2008), but in some cases are less effective during the melt period because of loss of coherence between radar images (Reference Van Wychen, Copland, Gray, Burgess, Danielson and SharpVan Wychen and others, 2012). Therefore it is important to have an estimate of the magnitude of seasonal flow variability when evaluating errors associated with ice flux estimates from sub-annual measurements.
Glacier flow variations may include a ‘spring event’: a sudden increase in horizontal ice velocity and simultaneous vertical ice motion that coincide with meltwater flow into moulins early in the melt season, initiating a supra- to subglacial drainage connection (Reference Iken, Röthlisberger, Flotron and HaeberliIken and others, 1983; Reference Röthlisberger, Lang, Gurnell and ClarkRöthlisberger and Lang, 1987; Reference Bingham, Nienow, Sharp and CoplandBingham and others, 2006). Additional ice uplift and fast-flow events later in the melt season have been associated with periods of strong weather-induced meltwater production and increased runoff into moulins. The enhanced rates of glacier flow in these events are usually interpreted as a result of hydrologically enhanced basal sliding, which occurs when water inputs exceed the capacity of the existing subglacial drainage network and increase subglacial water pressure to the point where water-filled cavities grow and ice–bed separation increases (Reference Iken and BindschadlerIken and Bindschadler, 1986).
Velocity variations have also been observed on High Arctic glaciers and parts of the Greenland ice sheet in association with the rapid drainage of supraglacial lakes via crevasses opened by hydrofracture (Reference Boon and SharpBoon and Sharp, 2003; Reference Van der VeenVan der Veen, 2007; Reference DasDas and others, 2008; Reference SoleSole and others, 2011). Lakes accumulate meltwater over a period of days to weeks and may deliver some or all of that water to the subglacial drainage system in an event lasting a few hours. This mechanism has the potential to rapidly perturb a slowly evolving subglacial drainage system. Both rates of meltwater production and supraglacial lake drainage events must therefore be considered when attempting to assess the causes of glacier flow variability.
Supraglacial lake drainage events have been observed directly and by remote-sensing methods. Individual lakes have been instrumented with pressure transducers and monitored in detail to determine fill and drainage rates (Reference Boon and SharpBoon and Sharp, 2003; Reference Bartholomaus, Anderson and AndersonBartholomaus and others, 2008; Reference DasDas and others, 2008). However, on larger glaciers or regions of the Greenland ice sheet where many lakes exist over a large area, it is difficult to monitor multiple lakes in this way. Satellite imagery has been used to document the spatial distribution of lakes on the Greenland ice sheet, and methods of estimating lake depth and volume have been proposed (Reference Box and SkiBox and Ski, 2007; Reference Sneed and HamiltonSneed and Hamilton, 2007; Reference Tedesco and SteinerTedesco and Steiner, 2011). Only very recently have Moderate Resolution Imaging Spectroradiometer (MODIS) images been used to determine the fill and drainage rates of these lakes at sub-daily temporal resolution (Reference Selmes, Murray and JamesSelmes and others, 2011). The temporal and spatial resolution is still insufficient to reveal details of the drainage mechanism.
Terrestrial time-lapse photography may present a way to observe lakes on large outlet glaciers at higher spatial and temporal resolution. While it may not be ideally suited for some locations where elevated vantage points are unavailable (e.g. flatter parts of the Greenland ice sheet), time-lapse photography has been used with success at many alpine and outlet glacier settings. An early prototype of an automatic time-lapse camera was deployed to observe the advance of an alpine glacier on Mount Rainier, Washington, USA (Reference Miller and CrandellMiller and Crandell, 1959). Other glacial features or processes that have been monitored with time-lapse photography include proglacial streams (Reference Humphrey and RaymondHumphrey and Raymond, 1994), snow surface albedo changes (Reference CorripioCorripio, 2004), ice melange interactions at the front of Jakobshavn Isbræ, Greenland (Reference Amundson, Fahnestock, Truffer, Brown, Lüthi and MotykaAmundson and others, 2010), iceberg plume discharge events (Reference Herdes, Copland, Danielson and SharpHerdes and others, 2012) and meltwater plume discharge (Reference MilneMilne, 2011). The combination of time-lapse photography, digital image matching techniques and photogrammetric methods has been used to determine glacier terminus velocity (Reference Krimmel and RasmussenKrimmel and Rasmussen, 1986), detect calving-line position and area change at Alaskan and Svalbard marine-terminating glaciers (Reference Motyka, Hunter, Echelmeyer and ConnorMotyka and others, 2003; Reference O’Neel, Echelmeyer and MotykaO’Neel and others, 2003; Reference Chapuis, Rolstad and NorlandChapuis and others, 2010) and calculate surface velocity fields on Greenland outlet glaciers (Reference Ahn and BoxAhn and Box, 2010). Reference Boon and SharpBoon and Sharp (2003) used a time-lapse camera to monitor lake development on John Evans Glacier, Ellesmere Island, Canadian Arctic, but to our knowledge photogrammetric methods have not been applied to time-lapse imagery of supraglacial lakes to quantify the size or rate of change of these reservoirs.
In this study, we combine photogrammetric methods with automated digital image processing and time-lapse photography to study the impact of supraglacial lake drainage events on the flow variability of Belcher Glacier, a tidewater outlet glacier of Devon Ice Cap, Nunavut, Canada. Time-lapse cameras were used to observe several of the most significant supraglacial lakes on the glacier over an entire melt season. We describe how the sequence of images from each camera was used to produce a time series of supraglacial lake area coverage change, which we used to estimate rates of water delivery to the subglacial drainage system. Comparing this dataset with multiple time series of glacier flow velocity derived from GPS measurements allows the identification of ice acceleration events that correspond in time with lake drainages. These short-duration events are part of a broader pattern of seasonal flow acceleration that occurs on this glacier. After isolating these events and measuring the corresponding ice displacement at the glacier terminus, we discuss the significance of the flow enhancement driven by lake drainage events for the total seasonal and annual ice displacement.
2. Study Site
Devon Ice Cap is among the largest ice caps in the Canadian Arctic; it has an area of 14 000 km2 and contains ∼4110 km3 of ice (Reference Burgess and SharpBurgess and Sharp, 2004). Numerous outlet glaciers, many of which have marine termini, punctuate the north, east and southern boundaries of the ice cap (Fig. 1).
Landsat-7 image (August 2000) of Belcher Glacier showing locations of time-lapse cameras, GPS stations and lakes described in Section 3. Grid coordinates are in University Transverse Mercator (UTM) zone 17X. Upper inset map shows Devon Island, which is part of the Canadian Arctic Archipelago. Devon Ice Cap (red box) is located at 75° N, 80–90° W. Map selected from the International Bathymetric Chart of the Arctic Ocean (Reference JakobssonJakobsson and others, 2008). Lower inset image (also Landsat-7, August 2000) shows the location of Belcher Glacier in the norteast quadrant of Devon Ice Cap.
Belcher Glacier is a ∼35 km long tidewater outlet glacier that flows through a steep-walled fjord and is fed by multiple tributaries. Flow rates are fastest in the terminus region, up to 336 ± 0.1 m a−1 as measured by repeat differential static GPS surveys. Airborne radar data show the grounded terminus ice is ∼250 m thick and the glacier remains grounded below sea level 15 km up-glacier from the terminus (Reference Dowdeswell, Benham, Gorman, Burgess and SharpDowdeswell and others, 2004). The glacier is a large source of icebergs, accounting for ∼17% of the net mass loss from Devon Ice Cap between 1960 and 1999 (Reference Burgess, Sharp, Mair, Dowdeswell and BenhamBurgess and others, 2005). The majority of this calving loss is due to the glacier’s high rate of ice flux rather than to retreat of its terminus, which has remained relatively stable for the past 50 years (Reference Burgess and SharpBurgess and Sharp, 2004; Reference Burgess, Sharp, Mair, Dowdeswell and BenhamBurgess and others, 2005).
3. Data
3.1. Lake observations
We have examined July and August imagery from Landsat-7 (1999, 2000), SPOT5 (2007) and Government of Canada aerial photographs (1959, 1960) to identify locations where lakes form on the surface of Belcher Glacier. Comparison of 1959–60 and contemporary imagery confirms that many of these lakes have formed repeatedly over decades.
Photographs taken on the glacier and from light aircraft during May and August fieldwork campaigns in 2006–09 provided further evidence that many lakes form annually in the same basins. During spring fieldwork, we have explored the topographic basins in which the lakes form and looked for signs of moulins or incised drainage channels leading into and out of lake basins. These can be difficult to find and map accurately due to snow cover. We also noted that some basins were ringed with stranded lake ice while others were not, which suggested that some lakes persist until after freeze-up while others drained before then. Return trips in mid- to late August confirmed that some lake basins were empty by the end of the melt season while others still contained water.
3.2. Time-lapse photography
The aforementioned imagery was used to identify potential cliff-top locations that offered viewing positions for several of the largest supraglacial lakes. In May 2007, we reconnoitered these sites by snowmobile and helicopter and began installing time-lapse camera systems. We used Harbortronics time-lapse kits, Pentax K110D 6M pixel DSLR cameras, Pentax 18–55 mm lenses and 4GB SD memory cards. Each camera recorded high-quality 3008 × 2000 JPEG files (using the in-camera Pentax RAW to JPEG convertor). More details on camera settings and time-lapse installation are provided in Appendix A.
In 2007 and 2008 we had mixed success with the time-lapse systems and our chosen viewing locations. After replacing some of the intervalometers and moving the cameras to more suitable locations, we successfully operated four time-lapse systems at a 1 hour image-capture rate over the entire spring and summer of 2009. These cameras observed the four lakes and crevasse field labeled in Figure 1.
3.3. Ice velocity time series
We deployed five Trimble dual-frequency GPS receivers that acquired continuous GPS observations at a 15 s sampling rate from May through August 2009. One receiver was used as a base station and was installed on a steel tripod anchored to bedrock at a cliff-top site. The other four were installed on the glacier surface along the center line (referred to as GPS1–4 in the rest of this paper). These four stations were used to continuously measure the flow speed of the ice; the GPS antennas were fixed to the tops of steel poles drilled 2–3 m into the glacier ice (Table 1).
On-ice GPS site details
We used TRACK, the kinematic processing module associated with the GAMIT software (Reference Herring, King and McCluskyHerring and others, 2006; http://www-gpsg.mit.edu/∼simon/gtgk/), to post-process our GPS observations differentially in kinematic mode. For each set of raw GPS observations, TRACK produces a kinematic position estimate, with 1-σ uncertainties stated for each coordinate dimension. These datasets were filtered to remove any positions with two-dimensional 1σ distance root-mean-square errors >0.05 m or vertical 1σ uncertainties >0.1 m. All remaining position estimates produced from each GPS station were concatenated into individual long time series. These time series were then sampled uniformly to produce position estimates every hour. Horizontal ice velocity time series were produced for each site by estimating the horizontal displacement between position estimates, then dividing by the intervening time-step. A 4 hour moving average was used to remove the effects of high-frequency noise.
3.4. Air temperature time series
Automatic air temperature loggers, installed in radiation shrouds, were co-located with the GPS stations to measure the local air temperature every hour over the duration of the summer. The GPS station nearest the glacier terminus (GPS4) did not have a co-located temperature logger, so we used the air temperature time series from the nearest available logger (located 5 km southwest and at +300 m elevation) adjusted by an altitude-dependent lapse rate of 4.9°C km−1 (Reference Gardner and SharpGardner and Sharp, 2009) to estimate the air temperature in the terminus area.
3.5. Ice surface elevation
A digital elevation model (DEM) of Devon Ice Cap was produced using SPOT5 satellite imagery collected in summer 2007, as part of the Spot 5 stereoscopic survey of Polar Ice: Reference Images and Topographies (SPIRIT) project. The SPIRIT DEM shows good agreement with Ice, Cloud and land Elevation Satellite (ICESat) laser altimetry in the ablation zone (Reference Korona, Berthier, Bernard, Rémy and ThouvenotKorona and others, 2009), which includes the Belcher Glacier region. We found that this DEM has sufficient spatial detail to resolve supraglacial lake basin geometry, which we used to assess the accuracy of the position and size of lakes identified via our time-lapse photography analysis methods.
4. Methods
We deployed four time-lapse cameras between 27 May and 15 August 2009 and each acquired ∼2000 images. With such a high volume of imagery, it was desirable to automate the process of finding lakes and measuring the area covered by water. A strict method is needed to quantify the size of the lake in each image and produce a time series of lake sizes from which filling and drainage rates can be derived. We developed a set of image-analysis procedures that involved the following steps:
-
1. Image pre-filtering and masking: data cleaning
-
2. Automated image classification: find the lakes
-
3. Image-space to real-space conversion: measure the lake areas
-
4. Time-series production: compare the lake areas over time
4.1. Image pre-filtering and masking
Fog cover was a persistent problem at all camera locations due to the close proximity to the open water of Baffin Bay. As our image classification procedure was unable to produce usable results from images partially or wholly obscured by fog, we required a pre-filtering algorithm to identify and remove fog-filled images from each set of time-lapse images. As a result of JPEG compression, uniformly grey, fog-filled images require much less memory storage than images showing a variety of colors. Any image with a file size below a tunable threshold was omitted from further processing.
Some of our time-lapse cameras observed multiple lakes or water-covered areas, but we were interested in monitoring and measuring individual lakes. For example, TLCam3 observed both Lake 3 and Lake 4, but we wanted to observe them individually. Likewise, TLCam4 observed a field of water-filled crevasses that we were interested in, but also a small lake far across the glacier and meltwater or marine water beyond the glacier terminus that we had no interest in. For this reason we applied a digital mask or cropped the images (preserving the original image coordinates) before applying our procedures so that individual regions of interest could be monitored.
4.2. Automated image classification
A systematic way to identify a lake in an image is to employ an image pattern recognition function that examines each pixel in an image relative to a set of criteria to determine whether the pixel fits in the target class (i.e. water) or in one or more non-target classes (i.e. snow, ice, rock, other). We used color information to provide the criteria for determining pixel class membership. Human vision also uses characteristics such as roughness or edge detection to identify transitions between snow/ice and water, but such criteria are highly variable and more challenging to implement on a pixel-by-pixel basis.
A digital image can be worked with as an X-by-Y-by-3 data array, X and Y being the pixel dimensions of the image, with three dimensions or layers of color information. The default RGB color pallet divides red, green and blue reflectance values into the three color layers. Unfortunately, lake water is not simply ‘blue’ and the reflectance values of water change as depth changes, or as lighting conditions change with weather and time of day. The significant limitation of the RGB color representation is that lighting and saturation information are mixed into each color layer, making it a poor choice for pixel classification criteria in a natural scene. Instead, we converted each time-lapse photograph into two alternative color representation systems (called HSV and L*a*b* color spaces) that define color along continuous spectra of hue and chromaticity and separate brightness and saturation information into different layers. These conversions allowed us to ignore most of the brightness variation and to focus on the hue and chromaticity layers as pixel classification criteria (see Appendix B for more detail).
We performed pixel classification using a software routine developed around the K-means (Reference MacQueen, Le Cam and NeymanMacQueen, 1967) algorithm in MATLAB. K-means is technically not a classification technique, but a data-clustering method. The user defines the number of clusters (defined by the parameter ‘k’) that the data should be sorted into. The choice of ‘k’ depends partly on how many groups within the data the user desires to find, but also on the nature and modality of the dataset itself. Some amount of reiterative testing may be required in order to find an appropriate choice of ‘k’ (Reference MacQueen, Le Cam and NeymanMacQueen, 1967).
The K-means algorithm begins by randomly selecting ‘k’ pixels from an image and uses these pixel attributes (the hue and chromaticity layer values) as the initial ‘means’ or the definition of the cluster centers. Next, all the remaining pixels from that image are assigned to the cluster with the closest mean attribute set. The cluster means are updated (becoming the average of all pixels in the cluster) and then the process is repeated until the algorithm converges on a solution where the cluster means do not change after the update step (Reference MacKay and MacKayMacKay, 2003).
The MATLAB implementation of the K-means algorithm includes numerous refinements and options for optimizing clustering results. One optimization we chose to employ involves running multiple replicates of the clustering process on each image, using new random starting pixels. Variations in the tightness of the clusters (minimum sum of distances between all pixels and their cluster mean) can occur by using different starting points; we kept the tightest result of five replicates.
After the K-means algorithm was run on an image, all pixels in the image were labeled with their cluster identification. Using the cluster associated with ‘blue’ band values, a binary mask was created to isolate the ‘lake’ feature in the image. A subsequent algorithm then drew an outline around the lake perimeter and created a list of the image-space coordinates of the lake perimeter. Finally, a rapid quality control procedure allowed us to reject any erroneous lake outlines that resulted from fog, shadow, reflections or floating lake ice obscuring parts of the lake.
4.3. Image-space to real-space conversion
The points composing the lake perimeter were then converted from pixel coordinates to Universal Transverse Mercator (UTM) coordinates. The real-space area of the lake (m2) in each image could then be calculated using the georectified perimeter.
The UTM coordinates of the lake perimeter were determined using a monophotogrammetric (MPG) technique implemented in MATLAB by M. Truffer at the University of Alaska Fairbanks. This technique adapts the methods of Reference Krimmel and RasmussenKrimmel and Rasmussen (1986) for use with digital photography. A detailed description of the technique and associated error analysis is given by Reference O’NeelO’Neel (2000).
This MPG technique has been utilized at several study sites in Alaska (Reference O’Neel, Echelmeyer and MotykaO’Neel and others, 2003, Reference O’Neel, Marshall, McNamara and Pfeffer2007) and Svalbard (Reference Chapuis, Rolstad and NorlandChapuis and others, 2010) to measure temporal changes in the calving-line position of tidewater glaciers. In these cases, the points where the glacier terminus ice cliff meets the marine water surface form an intersection line along a plane of known elevation: 0 m a.s.l. (± tidal fluctuation). Utilizing this constraint on elevation, the horizontal coordinates of any point along this intersection line may be estimated if the position and orientation of the camera and the geographic locations of several ground-control points (GCPs) in the scene are known. For our scenario, we assumed that all points along the surface of a lake were at the same elevation, thus the lake shoreline intersects a plane. We were therefore able to apply the same method to solve for the horizontal UTM coordinates of the lake perimeter with minor modifications to the code provided.
Our time-lapse cameras were initially installed only to perform a monitoring role and no photogrammetric campaign was planned; therefore, we had not measured GCPs in any of the scenes except for the terminus scene imaged by TLCam4. GCPs were therefore selected from high-resolution orthorectified imagery. An example time-lapse scene and a SPOT5 high-resolution stereo (HRS) image (captured on 20 August 2007) of the same area were examined together to locate sharp rock outcrops, nunatak edges or incised drainage channels that appeared to be persistent features that could be seen in both images. These were selected as potential GCPs. The pixel coordinates from the time-lapse image and the geographic coordinates from the SPOT image for each GCP were co-registered. Elevation values for the GCPs, lake surface and camera position were extracted from the SPIRIT DEM of the region.
The camera parameters required included the three-dimensional camera rotation angles (which were measured in the field, estimated from photographs of the camera installation or calculated based on scene geometry) and the image enlargement factor (a function of image size, sensor size and focal length).
Testing the quality of the GCPs and camera orientation parameters involved an iterative trial-and-error process: all the parameters were entered into the input settings file used by the MPG program and the program was run using a training image of the scene that included the lake at or near maximum size, with the outline drawn around it. The output of the MPG program (UTM coordinate outline of the lake perimeter and relative position of GCPs) was examined after each run to determine whether the relative positions of GCPs and lake appeared correct and the shape of the lake outline appeared true. This was verified by plotting the GCPs and lake outline on the SPOT image with the original GCPs. If the fit was not good then the program was run again after one parameter in the MPG input settings file was modified. Subtle variations in any of these parameters could produce large changes in the correspondence between the lake outline and GCP points from the time-lapse image and their projection on the SPOT image. A set of GCPs and other input parameters was settled upon when it was determined that the overall fit of the solution was as good as possible.
4.4. Time-series production
The three steps described in Sections 4.1, 4.2 and 4.3 were applied to all the images from the time-lapse cameras. The date and time were stored with the lake area calculated from each photograph to produce a time series of lake area change, beginning when the lake was first detected and ending sometime after the lake had drained or when the camera was retrieved. A time series was produced for each of the four major lakes and the water-filled crevasse area in the glacier terminus region. In order to produce an approximation of the net drainage rates of each lake, we calculated the first derivative of the lake area time series and plotted the negative values, which represent reduction in lake area (m2 h−1).
These time series of lake areas and net drainage rates were plotted on the same time axis as our GPS-based ice velocity time series and our measurements of air temperature. We made detailed observations of each individual time series, noting significant changes in surface melt conditions (inferred from air temperature), instances when lakes switched from filling to draining, and rapid variations in ice velocity. We then combined these observations to evaluate the potential relationships between melt events, lake drainage events and specific glacier acceleration events.
4.5. Lake area validation
We compared the maximum area of each lake found in the time-lapse images with the area of the same lake as seen in orthorectified aerial or satellite imagery, where available. We also used the SPIRIT DEM, which has sufficient accuracy and spatial detail to resolve the lake basins, to validate the area of the lakes. We used the Quick Terrain Modeler (QTM) software (www.appliedimagery.com) to draw a contour (using QTM’s Flood Analysis tool) at the elevations of the maximum lake surface levels on the DEM (see Fig. 3 below). We found the maximum lake surface level by adjusting a hypothetical ‘flood level’ such that water filled the basin of interest without connecting to adjacent basins or flooded down-glacier areas. The Flood Analysis tool uses a simple ‘Bathtub’ model of flooding and reports both an area and volume estimate for inundated regions of the DEM. By progressively lowering the flood level until the basin of interest was ‘dry’, we determined the basin floor elevation and thus the lake depth. Unfortunately, the DEM did not have sufficient spatial detail to resolve the terminus area crevasses, so we were unable to determine the area or volume of water-filled crevasses in the same manner.
5. Results
5.1. Image classification results
Using our automated lake detection method, we were able to identify and delineate lakes faster than was possible in trials using manual digitization methods. The K-means clustering technique enabled the identification of water-covered pixels in a variety of lighting conditions and challenging viewing situations. The segmented groups of water-filled crevasses provided the most challenging scenario attempted (a sample of the results is shown in Fig. 2). Ambiguity between water-filled crevasses and strongly shadowed crevasses often occurred, which required many of the results from the TLCam4 time series to be discarded in the quality control procedure.
Single image acquired by TLCam4 on 30 June 2009 showing the results of our automated image classification technique. In this case our technique has successfully picked out water-filled crevasses (outlined in red).
5.2. Classification challenges
One challenge we encountered was due to the changing appearance of water on the glacier surface as the melt season progressed and the surface changed from snow to slush to ice to dirty ice. The number of surface types and therefore colors recognized increased with time. K-means clustering could not converge on a stable solution if the choice of ‘k’ was inappropriate for the data. For instance, when the image was dominated by only two features, pure white snow and blue water, and we attempted to find five clusters (‘k’ = 5) then K-means either failed to converge on a solution and reported an error or attempted to divide the two logical clusters into inappropriate fractions. Conversely, if ‘k’ was set too low for an image with high variability (i.e. ‘k’ = 2 for a scene composed of snow, ice, dirt, rocks, water, shadow, etc.) then the algorithm was forced to merge clusters and produced unusable results.
To mitigate the effect of changing surface appearance over time, we divided the time-lapse images into sub-seasons (i.e. early, mid- and late summer) so that the number of clusters could be changed to optimize the clustering results for the conditions. We found that ‘k’ = 2 effectively clustered snow and water early in the summer and that clustering was more effective when ‘k’ was increased to 3 or 4 later in the summer.
5.3. Coordinate conversion results
Owing to the non-ideal methods used for GCP registration, absolute errors in the coordinate projections produced by the MPG program were quite high. The magnitude of error differed for each time-lapse camera scene, but we noted mis-registration of up to 150 m for some projected lake perimeters. However, the relative position quality was satisfactory: some lake outlines were slightly distorted, but the general shape and area of each lake outline produced by the MPG program corresponded to what can be seen in orthorectified imagery or the SPIRIT DEM.
Figure 3e–h shows a comparison of lake outlines produced by our automated method (red outline) with manually drawn outlines of the same lakes visible in a variety of orthorectified imagery (blue outline). Figure 3e shows Lake 1 at ∼25% filled state in a 1999 Landsat-7 image; the red outline represents Lake 1 at a similar fill-state in 2009. Figure 3f shows Lake 2 at or near its maximum fill level in a 1959 aerial photograph. Noticeable changes to shoreline geometry have occurred in the intervening 50 years, but the 2009 maximum lake area is very similar. Figure 3g shows the red outline of Lake 3 at maximum fill-state in 2009 and a blue outline of the lake at ∼70% fill-state, as seen in a georeferenced DigitalGlobe (WorldView2) quicklook image captured on 24 June 2011. Figure 3h shows Lake 4 at maximum size in a 2007 SPOT5 image compared with the red outline of Lake 4 at maximum size in 2009. Comparisons of lake areas based on these red and blue lake outlines are presented in Table 2.
Lake area validation. (a–d) Each lake at its maximum fill level in the original oblique-angle time-lapse photograph. The red outline is the lake outline detected by our image classification procedure and corresponds to the set of points converted from image-space to real-space (UTM) coordinates. (e–h) Comparison of the lake outlines produced by our method (red) with outlines of the same lakes visible in orthorectified imagery (blue), which were drawn manually. Since these images are not coincident in time, some differences are expected. The red arrows represent the view angle of the time-lapse camera. (i–l) Comparison of the lake outlines from the time-lapse images (red) with the lake basins found in the SPIRIT DEM.
Lake area comparisons
Since we could not always find the lakes at fill-states useful for comparison in the available orthorectified imagery, we used a second method to validate the area of the four lakes. As described in Section 4, we used QTM’s Flood Analysis tool and the SPIRIT DEM to estimate the area, volume and depth of our selected lake basins.
The lake basins found via this method are shown as blue areas outlined in black in Figure 3i–l. Areas for each lake found via this method are shown in Table 2. Lake volume and depth estimates are given in Table 3. Comparison of the maximum lake areas found via the DEM method and our photogrammetry method shows good agreement; the maximum area difference is 13%. As we are only interested in areal changes of the lakes, and not in absolute positions, we consider our lake outline coordinate conversion results acceptable.
Lake volume estimates. Each volume and depth estimate is reported with an uncertainty found by raising and lowering the lake surface elevation by 0.25 m, which is the limit of resolution of the DEM
6. Observations
6.1. Ice velocity and air temperature events
The air temperature time series in Figure 4e can be divided into two periods: spring from day 170 to 190 (19 June–9 July) and summer from day 190 to 217 (9 July–5 August). This division corresponds closely to the pattern observed in the horizontal ice velocity time series: the glacier exhibited different flow characteristics before and after day 190. Below we describe several of the distinct events observed during these two time periods.
In the spring period, air temperature rose above 0°C on most days and fell below 0°C on most nights. Time-lapse imagery from this period showed snow melting and water beginning to accumulate in supraglacial storage. The horizontal ice velocity at stations GPS1–3 remained consistently at or near the mean winter flow rates recorded at these locations and there was little to no variability in vertical ice motion. However, ice near the terminus began to accelerate during this period, as seen in the GPS4 horizontal velocity time series.
Acceleration event A1, days 172.5–174.3 (21–23 June), marks the first departure from the relatively stable winter flow rate of 0.6 m d−1 in the terminus region. This event was composed of two pulses of high flow speed and simultaneous abrupt jumps in vertical ice motion (up to 0.2 m). The end of this event was followed by a gradual increase in ice velocity leading up to event A2.
Acceleration event A2, days 181.0–185.5 (30 June–4 July), began with a rapid acceleration of ice flow to ∼1 m d−1, which was sustained for 19 hours followed by a 24 hour deceleration to 0.85 m d−1. Event A2 culminated with a rapid return to ice flow above 1 m d−1, after which flow speeds return to rates seen prior to day 181. A slight ‘bulge’ can be seen in the vertical motion of GPS4 during this event, though it is less distinct than the rapid uplift seen during event A1. The dual peaks in ice velocity seen during event A2 correspond closely in time to twin peaks in air temperature recorded in this region. Similar air temperature fluctuations were recorded at the measurement sites up-glacier, but the ice velocities measured at stations GPS1–3 did not vary in the same way as at GPS4.
In the summer period, air temperature remained consistently above 0°C, apart from two short negative excursions in the higher-elevation regions of the glacier. The air temperature generally ranged between 0°C and +6°C, with the exception of one particularly warm day during which a peak temperature of +13°C was reached in the lowest-elevation regions of the glacier. Horizontal ice velocity became much more variable during this period, and several specific events are described below. In addition, between days 200 and 204 (19–23 July) ice flow at GPS2 and GPS3 exhibited distinct diurnal variations that were in phase with diurnal temperature variations.
Acceleration event A3 was a wave of increased ice velocity that was recorded along most of the glacier between days 191 and 197 (10–16 July). It initiated on day 191 in the mid-upper glacier, when ice velocity at GPS2 rapidly increased in 24 hours from a background rate of ∼0.18 m d−1 to a rate of 0.64 m d−1 . Rapid vertical motion of the ice lagged the onset of horizontal acceleration by ∼16 hours. The ice surface at GPS2 then rose 0.15 m up to 0.2 m and remained near this elevation for the duration of event A3. Acceleration at GPS3 began as the first peak in velocity was reached at GPS2, at the beginning of day 192. Velocity at GPS3 peaked and began to slow down at the beginning of day 193, at which point GPS4 began to accelerate. Near the end of day 193, GPS2 and GPS3 simultaneously accelerated again and velocity at GPS4 continued to climb. During this second episode of acceleration, GPS3 recorded its highest velocities in 2009 (0.87 m d−1), and GPS4 its second highest (1.52 m d−1). Velocity peaked at GPS2 (0.82 m d−1) on day 195. Event A3 concluded at the beginning of day 197, when the velocities at GPS2–4 returned to pre-A3 velocity. GPS3 and GPS4 did not show significant variation in vertical ice motion during this event.
Near the end of event A3 (day 195.5, 14 July) when GPS2 velocity fell rapidly, horizontal and vertical velocities at GPS1 began to vary. Vertical uplift of ice at GPS1 increased gradually from day 195.5 to peak at day 199.2 (18 July). During this 88 hour period, GPS1 horizontal velocity accelerated above its pre-summer mean and fluctuated several times, often in phase with velocity at GPS2. Between days 201 and 204 (20–23 July) both GPS1 and GPS2 exhibited distinct diurnal velocity fluctuations accompanied by vertical motion undulations, which occurred in phase with diurnal air temperature variations in the same area. The ice remained elevated until the surface at GPS2 began to fall on day 205.5. Lowering followed at GPS1 on day 206.8. The ice surface elevation at both stations returned to the pre-uplift elevation by day 210. These events can be seen in Figure 4 and are enlarged in Figure 5.
Acceleration event A4 (days 197.3–198.5, 16–17 July) was rapid and led to the fastest ice flow rates measured in the 2009 GPS4 time series (1.59 m d−1). Acceleration also occurred at GPS3 at this time, though ice velocity did not exceed the peak reached during event A3. This short-duration event concluded when the velocities at GPS3 and GPS4 returned to pre-A4 values. Air temperature was decreasing slightly at the initiation of event A4.
Acceleration event A5 (days 204.3–205.3, 23–24 July) involved rapid and near-simultaneous acceleration of ice at GPS4 and GPS3, followed by acceleration at GPS2 ∼1 hour later. Velocities at GPS4 and GPS2 peaked rapidly and then began to fall, while velocity at GPS3 remained elevated for nearly the full duration of the event. This event began ∼12 hours before the onset of an extreme positive temperature anomaly (up to 13°C) on day 205 (24 July).
Acceleration event A6 (days 205.3–206.1, 25–26 July) immediately followed event A5. At GPS4, the peak velocities recorded during events A5 and A6 were approximately the same, while GPS3 and GPS2 recorded lower velocities during event A6.
6.2. Meltwater storage and drainage events
The collection of water-filled crevasses near the glacier terminus was the first supraglacial reservoir to fill and drain. This area of the glacier (the lower ∼5 km) is very heavily crevassed, with crevasses running approximately perpendicular to the glacier flow direction. However, they are not long straight fissures; these crevasses have been curved and contorted by ice flow, such that the crevasse mouths form openings of various widths (up to 10 m across) and lengths (tens to hundreds of meters). Available satellite and photographic imagery shows that the network of drainage channels visible on the ice surface in the upper- and mid-glacier regions is progressively destroyed as the ice is advected into the terminus region and becomes more crevassed. This makes surface drainage inefficient in this region of the glacier and we observed meltwater collecting in crevasse mouths over a wide area (Fig. 2), with no obvious routes for evacuation over the surface.
Although in Figure 4a the ‘filled crevasse area’ is plotted as a single aggregated meltwater storage reservoir, it is actually a collection of many small unconnected ponds that fill and drain independently of each other. Individually none of these crevasse ponds formed a large lake, but they covered a total area comparable with that of Lake 1. As described in Section 4, it is not possible to estimate the volume using the SPIRIT DEM.
The two most active phases of crevasse pond drainage (D1 and D2) were days 171.2–172.6 (20–21 June) and days 180.1–185.4 (29 June–4 July). D1 consisted of two pulses of drainage over ∼34 hours and D2 consisted of six drainage pulses over 5 days. (These drainage pulses are bounded within the blue boxes marking events D1 and D2 in the lake area change plot of Figure 4.) During each of these periods, the time-lapse imagery showed multiple individual crevasses fill with water and then drain rapidly. Because the total area of these crevasse ponds was aggregated, the rapid changes that occurred at individual crevasses appear as a few discrete drainage pulses in the plot of lake area change.
Drainage event D3 occurred between days 193.5 and 194.7 (12–13 July) when Lake 2 became the first of the large lakes to drain. The majority of the lake drained in the first 4 hours of the event, leaving a small remnant lake that continued to decrease in size over several more hours. It is difficult to identify a single crevasse or moulin that may have been responsible for draining the lake. The lake basin lay directly below an icefall, and numerous large crevasses were visible in the area. The only visible surface drainage channels lead into the lake basin from up-glacier and none appear to flow out, probably because there is a slight topographic rise just down-glacier of the lake basin.
Drainage event D4 (days 197.3–199.1, 16–18 July) involved the rapid complete drainage of Lake 3, the largest of the lakes observed in this study. The time series of Lake 3 area change shows that drainage of this lake began very rapidly, slowed briefly and then continued rapidly, emptying most of the lake within 14 hours. The time-lapse photography revealed several large deep crevasses within the lake basin immediately after drainage.
Drainage event D5 occurred between days 201.4 and 203.2 (20–22 July) when Lake 4 drained rapidly, leaving water in only the deepest part of its basin. The time-lapse imagery shows that this lake did not drain across the glacier surface. This lake is very close to the lateral margin of the glacier, and its down-glacier end lies in a deep trench formed along a medial moraine. A large crevasse extending from the edge of the glacier intersects this trench and provides the only indication visible on the surface of a potential englacial drainage route.
Drainage event D6 occurred after Lake 2 partially refilled at approximately day 202 (21 July) and then drained rapidly between days 203.9 and 204.6 (22–23 July). This second drainage of Lake 2 left no remnant lake behind.
The drainage of Lake 1 was not identified as a specific event as it is very different from the crevasse pond drainages (D1 and D2) and the rapid lake drainages (D3–D6) described above. This lake filled in ∼18 days and drained slowly over a similar time period. Maximum size was reached on day 198 (17 July). We examined the time-lapse imagery from this date and noticed that water in the lake drained through a snow dam which had filled a pre-existing incised drainage channel that intersected the southern margin of the lake basin. The breach of this snow dam resulted in the first significant drainage pulse that can be seen in the Lake 1 area change time series on day 198 (Fig. 5). The next day, a second, parallel, channel developed and the flow of water leaving the lake increased slowly as water flow progressively enlarged and deepened the supraglacial channels. However, these channels did not become deep enough to rapidly or completely drain the entire basin. The Lake 1 area change time series demonstrates that the lake decreased in size during the coldest parts of the day and grew during the warmest parts of the day, when the rate of meltwater influx matched or exceeded the rate of outflow from the drainage channels. For this reason, the drainage of Lake 1 appears as a series of small drainage pulses, the first eight of which occur with ∼24 hour periodicity.
Observations of the glacier surface drainage system from satellite imagery, field photography (from aircraft and from cliffs overlooking the glacier) and ground survey during the 2008 melt season (Reference DuncanDuncan, 2011) show that the outflow from Lake 1 drains into a large moulin ∼5 km from the lake and 1.5 km west of GPS1 (moulin location shown in Fig. 1).
7. Discussion
Examination of the summer 2009 horizontal ice velocity time series and air temperature measurements supports the hypothesis that variations in meltwater production drive ice flow variability. During the ‘summer’ period, air temperature demonstrates general covariance with ice velocity, and in some instances short-term variations in ice velocity occur in phase with diurnal temperature variations. However, some of the ice velocity variations (including the most rapid acceleration events observed) cannot be explained as a simple function of variable surface meltwater production. Instead, we suggest that several of these ice acceleration events are driven by supraglacial lake drainage.
7.1. Lake drainage typology
Based on the results presented, we identify three mechanisms for the drainage of supraglacially stored water:
-
Type 1 (crevasse pond drainage): Large crevasses fill with water forming a region of small discontinuous ponds. Lack of visible surface drainage routes and rapid asynchronous drainage of individual crevasse ponds suggest that the crevasses empty via their bottoms.
-
Type 2 (slow lake drainage): Supraglacial lake overfills its basin. Net storage ceases and net drainage occurs as water spilling out of the lake progressively erodes larger surface channels. Drainage rate is similar to filling rate (on the order of 102–103 m3 h−1 over multiple days), and incomplete drainage may leave a remnant lake in the deepest part of the basin.
-
Type 3 (fast lake drainage): Supraglacial lake drains via a crevasse or moulin opening within the lake basin. Drainage is complete and rapid (on the order of 104 m3 h−1 for <2 days).
To clarify, our use here of the term ‘slow drainage’ does not indicate an insignificant rate of meltwater drainage. The slow and fast drainage rates may be one or two orders of magnitude apart, but either drainage rate would be sufficient to sustain crevasse fracture propagation to the glacier bed and/or overcome rates of refreezing in cold-ice conduits, according to the findings of Reference Van der VeenVan der Veen (2007).
In the subsections below, we describe how the observed drainage events and associated ice acceleration events fit within this typology.
7.2. Type 1: crevasse pond drainage
Drainage of the crevasse ponds (events D1 and D2) was associated with localized acceleration of the terminus area of the glacier (events A1 and A2) well before any of the ice up-glacier had accelerated beyond the stable flow velocity of late winter/early spring.
7.2.1. Events D1 and A1
Approximately 26 hours after initiation of crevasse pond drainage event D1, ice acceleration and uplift (event A1) are seen in the GPS4 time series. The ice velocity during this event has a multi-peak pattern, similar to that of pond area change rate during event D1. The delay between the onset of drainage and ice acceleration could suggest that these are unrelated events. However, since this is the first drainage event of the year, the delay may instead indicate that water draining from crevasse ponds required some time to establish drainage routes through the thick ice in this region.
A previous investigation of the Belcher Glacier terminus region (Reference MilneMilne, 2011) utilized bathymetric mapping of the bay in front of the glacier (Ocean Mapping Group, 2008) and an ice thickness profile parallel to the terminus measured by airborne radio-echo sounding to assess whether the glacier front was floating. The bathymetry data revealed that the glacier encountered deep water (∼220 m) along a ∼1000 m wide medial segment of the terminus. Reference MilneMilne (2011) showed that the hydrostatic pressure exerted by the sea water across this section of the terminus would likely generate negative subglacial effective pressure in an area extending <500 m up-glacier from the margin. As GPS4 was placed along the glacier center line ∼500 m from the ice margin, it would have been located on the edge of this narrow band of ice that is potentially at or near the floating point. Because of the very low (close to negative) effective pressure at the glacier bed in this region, relatively small inputs of meltwater drainage could potentially increase ice–bed separation and cause measurable changes in basal sliding.
The ice uplift and speed-up at GPS4 following event D1 is consistent with the assertion that a hydraulic connection was made to the glacier bed and a brief spike in water pressure resulted in temporary ice–bed separation. Immediately after this acceleration event, both the ice velocity and ice surface elevation briefly decreased and then returned to pre-event velocity and elevation. This rapid reorganization of flow speed and surface elevation suggests that ice encountered momentary friction at the bed as the water supplied from the crevasse ponds was evacuated from the region and subglacial water pressure renormalized. If evacuated water formed a meltwater plume at the front of the glacier at this time, we were unable to see it because of remaining sea-ice cover.
Vertical ice surface motion could also have been caused by ice flowing over an obstacle at the glacier bed. We do not have sufficiently detailed information on bed topography to rule out the existence of bedrock undulations; however, it seems unlikely that roughness at the glacier bed could be expressed on the glacier surface over the horizontal distance of ∼1.18 m traversed by GPS4 during event A1. Alternatively, the observed uplift and subsidence could have been due in part to varying rates of longitudinal strain and associated vertical extension of the ice caused by the rapid localized velocity anomalies at GPS4 in relation to the surrounding ice.
7.2.2. Events D2 and A2
The second sequence of crevasse pond drainage, event D2, began on day 180, and a clear acceleration of ice flow (event A2) is apparent in the GPS4 time series within ∼18 hours. Ice velocity dropped to near the pre-acceleration velocity within hours of the cessation of the last drainage pulse in D2 (Fig. 4).
Lake area, GPS and temperature time series plotted on the same time axis. (a) Area of four lakes and the water-filled crevasse region measured from time-lapse photography. (b) Lake area change rate, our proxy for net drainage rate. Vertical shaded boxes delimit drainage events (D1–D6) and are color-coded to corresponding lakes. (c) Horizontal (xy) ice velocity measurements from four GPS stations. See Figure 1 for location of GPS stations. Boxes outlined in black delimit acceleration events (A1–A6), which are referred to in the text. (d) Cumulative change in ice surface elevation, relative to day 170, after correction for downslope motion. The GPS2, GPS3 and GPS4 time series have been offset by 0.2, 0.4 and 0.6 m, respectively, to improve viewing. (e) Air temperature measurements co-located with two of the above GPS stations.
D2 involved many crevasse ponds distributed over a ∼1.5 km2 area. Some of these crevasses drained for the first time during D2, though some may have drained during D1 and refilled if their exit conduits resealed. Several periods of freezing temperatures occurred between events D1 and D2, which would have reduced meltwater production and flux into established drainage conduits, perhaps allowing them to close via creep closure and refreezing. At the beginning of D2, the shorter delay (relative to event D1) between the initiation of drainage and ice acceleration suggests that the englacial drainage routes between the surface and sub-glacial drainage systems in this region were quickly reestablished. Alternatively, the shorter delay may have been the result of the greater number of water-filled crevasses (and thus higher water volume) that drained during event D2.
Compared with event A1, very little vertical motion was observed at GPS4 during event A2, which is superimposed on a longer-term gradual acceleration at GPS4 that began after A1 and continued through A3. This gradual speed-up of already fast-flowing ice in the lower glacier region would have increased longitudinal extension and the associated ice thinning and may have obscured transient vertical uplift. In general, neither GPS4 nor GPS3 demonstrates significant vertical surface motion during any of the subsequent horizontal acceleration events, while GPS2 and GPS1 show clear vertical uplift during the period of peak meltwater production and runoff (as indicated by air temperature). Both GPS3 and GPS4 are located on areas of the glacier with very steep surface gradients compared with the nearly flat surfaces at GPS1 and GPS2. It is possible that subtle vertical ice motion superimposed on the steep downslope flow of ice is harder to detect at GPS3 and GPS4 than at the other GPS sites.
7.3. Type 2: slow lake drainage
Type 2 drainage is exemplified by the drainage of Lake 1, which transitioned from a state of net storage to net drainage when drainage channels developed the capacity to remove water from the lake faster than influx from local melt production. The drainage channels emptying the lake connected to the existing supraglacial drainage network, and the released lake water augmented the daily runoff into the nearest sink point.
During the initial drainage period of this lake (days 198– 199, 17–18 July), the ice in the vicinity of the moulin exhibited its maximum seasonal uplift and flow velocity (see GPS1 in Fig. 5). Over the next 72 hours there appear to be coherent variations in lake area change, diurnal temperature, ice velocity and glacier uplift rates at GPS1 and GPS2. These events are marked by arrows in Figure 5. However, although these events seem to be linked, it is not possible to separate the relative contributions to enhanced ice flow made by lake outflow and runoff of daily meltwater production. The direct impact on ice velocity is not known.
Enlarged view of events from days 190 to 209. (a) Horizontal ice velocity at GPS1 and GPS2. (b) Vertical ice surface displacement at GPS1 and GPS2 corrected for downslope motion. (c) Air temperature at GPS2. (d) Lake 1 area change. The red arrows highlight the apparent synchronization of events in all four time series.
7.4. Type 3: fast lake drainage
Type 3 drainage is exemplified by Lakes 2–4, which fill gradually over a period of many days to weeks and then drain in a few hours when a drainage conduit opens, leaving an empty basin or tiny remnant of the original lake. The drainage of Lake 3 and the secondary drainage of Lake 2 (events D4 and D6) had the most pronounced impact on ice velocity, most notably in the terminus area but also at GPS3. The drainage of these lakes appears to accelerate the ice to short-lived velocity peaks. However, the fast drainage of Lake 4 (event D5) made no apparent impact on ice flow, and during event D3 there was a complex relationship between the drainage of Lake 2 and the ongoing evolution of supra- and subglacial drainage connections.
7.4.1. Events D4 and A4
Event D4, the drainage of Lake 3, coincided with rapid ice acceleration at GPS3, followed by acceleration at GPS4. The fastest flow rates at GPS4 in 2009 (1.59 m d–1) were recorded during this event, but peak velocity at GPS3 did not exceed the value reached during A3. Air temperature was decreasing at the beginning of event A4, so this acceleration event was probably not coincident with a sudden increase in meltwater production.
7.4.2. Events D6, A5 and A6
The extreme air temperature anomaly (+13°C) on day 205 caused high rates of meltwater production. We expect this meltwater was quickly routed to moulins via the surface drainage system that existed at this late stage in the melt season. Rapid delivery of runoff to the subglacial drainage system was the likely cause of acceleration event A6, but cannot explain event A5, which began 12 hours before the rise in air temperature. It is more likely that event A5 was triggered by D6, the final drainage of Lake 2. The peak rate of D6 lake area change coincided with the near-simultaneous acceleration of ice at GPS4 and GPS3 on day 204.
7.4.3. Event D5
Though event D5 was the second largest drainage event observed, drainage of Lake 4 did not cause noticeable flow acceleration at either of the nearest GPS stations (GPS3 and GPS4). Indeed, ice velocity was decreasing gradually during event D5. We have been unable to determine how water drains from this lake. However, a large crevasse that connects to the down-glacier end of the lake basin presents the possibility that this lake may drain towards the glacier margin and may not deliver water to the central region of the glacier.
7.4.4. Events D3 and A3
Acceleration event A3 began before drainage event D3, but velocities recorded at GPS2–4 peaked after D3. The following evidence suggests that the lake drainage did not trigger this acceleration event, but may have enhanced it. Day 190 marked the beginning of a 1 week warm period when air temperatures remained consistently above freezing. Snowpack diminished progressively, revealing bare glacier ice, and the resulting reduction in surface albedo will have contributed to higher melt rates. Time-lapse imagery shows the occurrence of a slush flow in the vicinity of GPS3, and the development of surface drainage channels along the entire mid- and upper glacier. All lakes observed on the glacier, except the crevasse ponds, grew during this period. Event A3 began at GPS2 on day 191.1 and then spread progressively down-glacier; GPS3 began to accelerate at day 191.6 and GPS4 at day 193.2. We interpret this period of flow acceleration as a ‘spring event’ (Reference Röthlisberger, Lang, Gurnell and ClarkRöthlisberger and Lang, 1987) brought on by high rates of meltwater production and the connection of many surface drainage channels to moulins and crevasses.
The ice surface uplift recorded at GPS2 early on day 192 may have been caused partly by longitudinal compression and vertical extension since, at that time, velocity at GPS2 was temporarily higher than velocity down-glacier at GPS3 (Fig. 4c and d). However, this cannot explain why the surface remained elevated and continued to rise when the downstream ice began flowing faster later in the day. Thus the ice uplift associated with acceleration at GPS2 more likely indicates ice–bed separation in response to the injection of surface meltwater into the subglacial drainage system.
At day 193.3, GPS3 decelerated to ∼0.6 m d−1 while GPS4 continued to accelerate (Fig. 4). The ice between GPS3 and GPS4, already undergoing extensional strain, will have temporarily experienced increased rates of longitudinal strain. Lake 2 is located in the zone between GPS3 and GPS4 and it is plausible that increased longitudinal stretching of the ice in this region caused the opening of crevasses, leading to the drainage of Lake 2. Lake 2 began draining on day 193.5, and peak lake area change rate (our proxy for net drainage rate) occurred at day 193.9. Peak velocity was recorded at GPS3 2 hours later (at day 194.0) and at GPS4 on day 195.0, 7 hours after Lake 2 ceased draining at day 194.7. Following these peaks, velocity fell gradually at both GPS3 and GPS4, returning to pre-A3 values by day 197. Air temperature also began decreasing after day 195, briefly returning to 0°C in the mid-glacier region on day 197, which will have decreased meltwater production rates.
The enhanced ice velocity recorded during event A3 seems to be a response to the combined influence of increased surface runoff and stored meltwater gaining access to the subglacial drainage system. While the event began as a ‘spring event’, the rapid injection of water from Lake 2 (∼366 400 ± 32 000 m3) into the subglacial drainage system likely contributed to the peak velocities and/or duration of acceleration event A3.
7.5. Impact on seasonal and annual ice displacement
In order to understand how the acceleration events induced by lake drainage affect the total flow of the glacier, we focus further attention on the terminus region. Ice displacement measured in the terminus region provides a good indication of the dynamic flux of a tidewater glacier. Our GPS4 measured ice displacement along the center line and within ∼500m of the glacier margin, the fastest-flowing region of the glacier.
The locations of the lakes and the relative timing of their drainage likely influenced the impact that drainage events had on the subglacial drainage system and in stimulating accelerated ice flow in the terminus region. The crevasse pond drainage events that we observed occurred very close to the glacier terminus, a region potentially close to flotation, and occurred very early in the melt season, before highly developed drainage routes had evolved. This combination of factors allowed a relatively small volume of meltwater drainage to produce a perturbation to a poorly developed subglacial drainage system, resulting in a measurable response in ice flow localized to the terminus region. In extreme contrast, the slow drainage of Lake 1 occurred far up-glacier late in the melt season when a highly developed subglacial drainage system was efficiently evacuating high volumes of ongoing melt. While there may have been some localized effect, there was negligible impact on ice flow acceleration in the terminus region. By comparison, the fast lake drainages that we observed occurred in the mid- or late melt season, when some form of subglacial drainage system was already in place, and occurred in the lower-mid-glacier. Only some of these events (i.e. D4 and D6) were sufficiently rapid and appropriately located to generate short-term hydrologic perturbations and stimulate measurable ice acceleration near the terminus.
The total ice displacement at GPS4 during the high-velocity events (A1, A2, A4 and A5) that we believe were driven by drainage events (D1, D2, D4 and D6) was 7.98 ± 0.14 m. However, the displacement resulting from ice velocity in excess of the background rate (i.e. the mean pre- and post-event velocity) during these events was only 1.34 ± 0.14 m. For comparison, the excess displacement resulting from the ‘spring event’ (event A3) was 1.08 ± 0.14 m. The summer displacement (in excess of the mean winter flow rate of ∼0.6 m d−1 at GPS4) that occurred between days 172 and 217 was 14.70 ± 0.14 m. The total annual ice displacement was 336.8 ± 0.1 m (May 2009–May 2010).
In summary, the additional ice displacements produced by the lake events, the spring event and the total summer accelerated flow represent 0.4%, 0.3%, and 4.4%, respectively, of the annual ice displacement near the glacier margin. However, it is difficult to meaningfully isolate these events, which are all part of the larger pattern of seasonal melt-enhanced ice flow. The drainage of crevasse ponds early in the melt season and of larger lakes in the mid-melt season likely served to initiate and enhance the overall seasonal meltwater-induced acceleration by opening connections between the surface and subglacial drainage networks. The prime example was event D3: while it was difficult to determine the immediate ice response to this single perturbation, it almost certainly played a role in enhancing the spring event that was observed along the entire glacier.
8. Conclusion
We have presented a new method for determining and monitoring the morphology of a supraglacial lake from time-lapse photography of a glacier surface. This has allowed us to estimate the area change rates (a proxy for net filling and drainage rates) of several lakes for which we have no other hydrological measurements. The photogrammetric methods we have employed are somewhat coarse and the results could be improved if GCPs were collected within the area covered by each camera view during the measurement period. Lake areas derived using our lake identification and photogrammetry procedures are within 14% of those found by manually digitizing the outline of the same lake in remotely sensed imagery. These results are acceptable for the purpose of determining the timing and relative magnitude of supraglacial lake drainage events and demonstrate what can easily be accomplished using oblique terrestrial photography.
A typology of lake-filling/drainage styles has been defined and can be linked to modes of ice velocity variability. Type 3 rapid lake drainages contribute to rapid, short-duration increases in glacier velocity. The observed link between these events is consistent with the hypothesis that drainage of some supraglacial lakes can quickly inject large volumes of water into the subglacial drainage system, raising subglacial water pressure and promoting ice flow acceleration.
The lake-driven rapid ice flow events observed here constitute ∼10% of the summer enhanced ice displacement, or 0.4% of the total annual ice displacement, as measured near the glacier terminus. This value may understate the significant impact that lake drainages, particularly those early in the melt season, may have on driving open conduits, priming the subglacial drainage system and initiating seasonally enhanced velocity.
Acknowledgements
This work was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery and International Polar Year (IPY) Special Research Opportunity (SRO) grants to M. Sharp. B. Danielson has been supported by the Dr S.M. Blair Research Graduate Scholarship, the Steve and Elaine Antoniuk Graduate Scholarship in Arctic Research, the Queen Elizabeth II Graduate Scholarship (Faculty of Graduate Studies and Research, University of Alberta) and by grants from the Canadian Circumpolar Institute and Northern Scientific Training Program (Department of Indian and Northern Affairs Canada). Logistical support for fieldwork was provided by the Polar Continental Shelf Project, Natural Resources Canada. We thank the Nunavut Research Institute and the communities of Grise Fjord and Resolute Bay for permitting us to work on Devon Ice Cap. We thank Shad O’Neel and Martin Truffer for sharing their monophotogrammetry code, and Alex Gardner, James Davis and Hannah Milne for their invaluable assistance with field data collection.
Appendix A
At each camera site, a 2 m steel tripod and guy-wires were anchored to exposed rock using expansion bolts to provide a solid foundation for the camera system.
The camera systems were produced by Harbortronics and include a weatherproof enclosure with glass window and mounting brackets, lithium battery with solar panel and charging regulator, an adjustable shoe-type camera mount, and a DigiSnap2000 intervalometer/shutter controller. The camera used was a Pentax K110D 6M pixel DSLR with an 18–55 mm lens and 4 GB SD memory card. At each installation, the camera lens was set to the desired focal length and manually focused on the scene. Auto-exposure and shutter speed were selected, so that the camera could adjust to changing lighting conditions, but a pre-programmed ‘daylight’ white balance was selected in order to provide consistent color cast through the resulting image series.
Appendix B
We converted all time-lapse images from RGB color space into the following alternative color spaces:
HSV color space
-
Hue (H): contains the ‘color’ information. Hue values (0.0–1.0) map to a color wheel ranging through ROYGBIV.
-
Saturation (S): represents the intensity or purity of a color. 0 would produce a grey pixel, while 1.0 would produce pure unmixed colors.
-
Value (V): contains the brightness information of each pixel. 0 is no brightness (black) and 1.0 is full brightness.
L*a*b* color space
-
L* dimension: luminance on a black-to-diffuse white axis.
-
a* dimension: chromaticity, indicating where color falls along the red–green axis.
-
b* dimension: chromaticity, indicating where color falls along the blue–yellow axis.
These conversions allowed brightness information to be isolated in the V and L* layers while H, a* and b* (and optionally S) layers were used as criteria for assessing class membership.
Documentation on the above color-space definitions is available in articles by G. Hoffmann, which are available online: CIELab Color Space – http://www.fho-emden.de/∼hoffmann/cielab03022003.pdf; Color Models: RGB / HLS / HSB – http://www.fho-emden.de/∼hoffmann/hlscone03052001.pdf. Lower inset map is a Landsat-7 image of Devon Ice Cap, with the location of Belcher Glacier outlined in the orange box. Both Landsat-7 images are from August 2000.
