Linking surface hydrology to flow regimes and patterns of velocity variability on Devon Ice Cap, Nunavut

Abstract

Supraglacial meltwater reaching a glacier bed can increase ice surface velocities via hydraulic jacking and enhanced basal sliding. However, the relationships between the structure of supraglacial drainage systems, sink-point distributions, glacier flow processes and the magnitude of interannual velocity variability are poorly understood. To explore the hypothesis that spatial variations in the rate and mechanisms of glacier flow are linked to variations in supraglacial drainage system structure and sink-point distribution across an ice cap, we mapped supraglacial drainage systems on Devon Ice Cap from Landsat-7 ETM+ imagery. Spatial patterns of surface velocity and interannual velocity variability were determined using gradient correlation applied to Landsat-7 ETM+ images. Velocity variability is greater in areas close to sink-point locations, presumably because hydrologically forced basal sliding and/or bed deformation are enhanced in such areas. The distribution and characteristics of supraglacial drainage systems may play an important role in determining the distribution of regions of basal sliding, highlighting the need for knowledge of the supraglacial drainage system structure and sink-point distribution to inform efforts to model the dynamic response of Arctic ice caps to future climate warming.

Introduction

Since 2005, mean summer air temperatures in the Canadian High Arctic have been higher than at any previous time in the past 60 years. This has resulted in lengthening melt seasons and a decrease in the climatic mass balance of ice caps and glaciers in the Queen Elizabeth Islands (Sharp and others, 2011). Increasing air temperatures may also affect the dynamics of Arctic ice masses, as higher air temperatures lead to increased production of surface meltwater. If this water reaches the glacier bed it may initiate, or increase the rate of, basal sliding and/or subglacial sediment deformation, thereby increasing the ice surface velocity (Iken, 1981; Arnold and Sharp, 1992; Zwally and others, 2002; Palmer and others, 2011). Higher ice surface velocities may, in turn, lead to increased rates of mass loss, as ice is transferred more rapidly from higher elevations to lower elevations, where melting is more rapid, or discharged into the ocean via iceberg calving and/or marine melting (Van der Veen and others, 2011). Furthermore, a feedback may exist whereby changes in the extent and distribution of crevasses resulting from spatially variable changes in surface velocity (and hence surface strain rates) can alter the rate and distribution of surface meltwater travel to the bed. This, in turn, may alter basal temperatures through cryohydrologic warming (Phillips and others, 2010), and both the rate and relative importance of the processes by which the glacier flows (Burgess and others, 2005). While increasing summer melting has been shown to increase summer mean velocity in the Arctic (Bingham and others, 2003; Bartholomew and others, 2010), an inverse relationship between summer melt and both annual and wintertime velocities has been observed in Alaska, likely due to changes in the subglacial hydrological system (Truffer and others, 2005; Burgess and others, 2013). Understanding the relationships between the distributions of basal motion and of the locations where surface meltwater reaches the bed (and how these vary over time), and how surface meltwater access to the bed affects annual ice surface velocities, is therefore critical in explaining and predicting the response of ice caps to hydrological forcing.

The surface velocity field of Devon Ice Cap, Nunavut, has been mapped by Burgess and others (2005) and, more recently, by Van Wychen and others (2012). Widespread slow-flowing ice in the interior plateau region of the ice cap is mostly discharged to the ocean by fast-flowing outlet glaciers, in which velocities generally increase towards the glacier terminus. Within such outlet glaciers, along-flow variations in the direction and slope of the relationship between the ratio of surface velocity to ice thickness and the driving stress have been interpreted in terms of a downstream evolution of ice-flow processes. Burgess and others (2005) recognized four distinct ‘flow regimes’ within Devon Ice Cap, and inferred a tendency for a systematic pattern of change in both surface velocities and flow processes with increasing distance along flowlines from the ice divide to the glacier terminus. The flow regimes were categorized as follows: (1) areas where ice is frozen to the bed and flow appears to be by ice deformation alone; (2) areas around the heads of outlet glaciers, where ice likely reaches the pressure-melting point at the bed and basal sliding begins; (3) areas where further reductions in basal friction occur and the rate of basal sliding increases, likely because surface meltwater reaches the bed more widely in these areas; and (4) areas near the termini of outlet glaciers, where subglacial sediment deformation may contribute to a further increase in velocity (Burgess and others, 2005). Thus, Burgess and others (2005) proposed that downstream increases in velocity and the relative contribution of basal motion to the flow of outlet glaciers are linked to changes in basal thermal conditions, increased penetration of surface meltwater to the glacier bed, and, in some cases, a change in glacier substrate type. However, this hypothesis has not been tested by detailed comparison of the relationships between flow regimes/surface velocity, supraglacial drainage patterns and the distribution of sink points through which meltwater can enter the glacier.

Here we aim to answer the following questions: (1) Are changes in supraglacial drainage style and the distribution of locations where supraglacial meltwater enters the glacier closely linked to changes in the rate of ice flow and the inferred flow regime? (2) Are changes in flow regime related to changes in either the number of locations at which surface meltwater drains into the glacier or the mode of meltwater delivery to the glacier bed (e.g. via a sinking stream of a rapidly draining supraglacial lake)? (3) Is interannual variability in glacier velocity higher in areas close to locations where meltwater drains into the glacier than it is in areas more remote from such locations?

Study Site

Devon Ice Cap is located at the eastern end of Devon Island, Nunavut, in the Canadian High Arctic (Fig. 1). The surface velocity field of Devon Ice Cap was first mapped by Burgess and others (2005) using synthetic aperture radar interferometry (InSAR) applied to imagery from European Remote-sensing Satellites 1 and 2 (ERS-1/-2), and speckle tracking applied to RADARSAT-1 imagery. This mapping has recently been updated by Van Wychen and others (2012, 2014). There are significant differences in flow characteristics between the eastern and western sectors of the ice cap. The western sector overlies a plateau, and westward-flowing ice terminates on land in an 82 km wide lobe at 400–600 m a.s.l. (Burgess and others, 2005). Flow in this region appears to occur primarily by ice deformation, with surface velocities typically <15 m a⁻¹. To the south of the ice cap’s main east–west divide, this flow regime is interrupted by the south-flowing North and South Croker Glaciers, which terminate in Croker Bay and have maximum velocities of ∼210 and ∼240 m a⁻¹, respectively. To the north of the east–west divide, the northwards-flowing Sverdrup Glacier terminates in Jones Sound and flows faster than ice in the plateau region (Burgess and others, 2005) (Fig. 1). East of the main north–south divide, the topography becomes more mountainous, and several large tidewater outlet glaciers drain to the ice-cap margin from up to 60 km inland (Burgess and others, 2005). In the northeast sector, Belcher Glacier draws ice from the northern slope of the main east–west divide. Velocities approach ∼75 m a⁻¹ in its upper reaches, and increase down-glacier to ∼290 m a⁻¹ at the terminus (Burgess and others, 2005). The southeastern portion of Devon Ice Cap is drained by Southeast 1 and Southeast 2 (SE1 and SE2) Glaciers. Here fast flow extends further inland than in any other sector of the ice cap (Burgess and others, 2005). SE1 and SE2 Glaciers converge ∼40 km down-glacier from their heads, in a zone of near-stagnant ice with velocities <10 m a⁻¹. At some point in the past, SE1 appears to have flowed faster than at present and has very likely surged (Burgess and others, 2005). Recent velocity measurements, derived by applying speckle tracking to 2009 RADARSAT-2 imagery, show a deceleration of velocity near the head of SE1 Glacier since 1996, and a down-glacier extension of the region of fast ice flow was reported for the middle reaches of the glacier by Burgess and others (2005), suggesting that this glacier may be experiencing a ‘slow surge’ (Van Wychen and others, 2012). There has also been acceleration of flow in the vicinity of the terminus of SE1 Glacier, but much of the terminus regions of these glaciers remains stagnant.

Fig. 1. Devon Ice Cap. The ice-cap summit is marked with a star, and other features are labeled. Inset shows its location in the Canadian High Arctic.

The melt season on Devon Ice Cap usually begins in June or early July. It is often characterized by a series of short melt episodes interrupted by freeze-up events during periods of inclement weather (Boon and others, 2010). Total annual accumulation over Devon Ice Cap shows little year-to-year variability, and averages ∼10 cm w.e. in the northwest, increasing to ∼50 cm w.e. at the southeast margin (Koerner, 1970; Boon and others, 2010). Interannual variability in the climatic mass balance of Devon Ice Cap is controlled largely by summer melt conditions (Koerner, 2005; Gardner and Sharp, 2007). The mass balance has become increasingly negative since 1985 (Boon and others, 2010), and the 2005–09 pentad had the most negative mass balance since 1960 (Sharp and others, 2011). The climatic mass balance varies across the ice cap, with the eastern side of the ice cap losing more mass to surface melting and iceberg calving than the western sector. This is related to the large areas at low elevations in the eastern sector of the ice cap (Burgess and Sharp, 2004; Burgess and others, 2005; Boon and others, 2010). Larger drainage basins (covering 86% of the ice cap) contribute 68% of net mass loss, while the smaller drainage basins (covering 14% of the ice cap) contribute 32% to surface mass loss. Thus, the mass balance of the smaller drainage basins (which tend to be located around the periphery of the ice cap, where elevations are lower) may be more sensitive to changes in climate than that of the larger basins (Mair and others, 2005).

Background

Burgess and others (2005) identified four distinct ‘flow regimes’ on Devon Ice Cap, on the basis of along-flowline trends in the relationship between the ratio of surface velocity to ice thickness, v/h, and the driving stress, τ _d. Where there is no basal component to motion, v/h (a⁻¹) represents the mean shear strain rate through a vertical ice column (Burgess and others, 2005). This is no longer true when basal motion contributes to surface velocities, but the inverse of the ratio v/h to τ _d has units of Pa s, and is thus a measure of the ‘effective viscosity’ of the glacier system (Burgess and others, 2005). The driving stress, τ _d, is derived as

(1)

where ρ _i is the density of ice, g is gravitational acceleration, h is ice thickness and α is the ice surface slope averaged over a distance ∼10 times the ice thickness.

Flow regime 1 (FR1) is defined by very low values of v/h, which increase with increasing τ _d. The low v/h and insensitivity to changes in τ _d in FR1 suggest that ice is cold-based and frozen to the bed (Burgess and others, 2005). FR1 is characteristic of the upper sections of flowlines and the interior plateau of Devon Ice Cap >1000 m a.s.l., where low velocities and driving stresses are common. The equilibrium-line altitude of Devon Ice Cap averaged 1166 ± 240 m a.s.l. from 1961 to 2003 (Boon and others, 2010), suggesting that much of the plateau region of FR1 experiences limited surface melt during the summer months. FR1 covers an area of ∼6300 km² and accounts for >50% of the classified flow regime area (Burgess and others, 2005). FR1 is also found in the terminus region of SE1 and SE2 Glaciers, where ice is nearly stagnant and probably in the quiescent phase of a surge cycle (Copland and others, 2003; Burgess and others, 2005; Van Wychen and others, 2012).

In flow regime 2 (FR2), which typically occurs immediately downstream of regions of FR1, v/h is higher for a given τ _d and is more sensitive to changes in τ _d than in FR1 (Burgess and others, 2005). Burgess and others (2005) interpreted this as evidence for a down-flow reduction in ‘effective viscosity’ of the glacier. This could result from warming of the ice so that it deforms more easily in response to a given stress, and/or basal ice reaching the pressure-melting point, so that basal motion begins to contribute to glacier flow. The transition from FR1 to FR2 usually occurs near the head of major outlet glaciers, coincident with the transition from a surface profile that is convex upwards to one that is concave upwards, and with the first appearance of flow stripes on the glacier surface (Burgess and others, 2005). The presence of flow stripes below the onset of FR2 suggests that basal motion contributes to surface velocities, as Gudmundsson and others (1998) showed that flow stripes develop where the ratio of basal velocity to velocity due to shear deformation of ice is large (Burgess and others, 2005). FR2 extends along significant portions of all the outlet glaciers considered in this study.

Flow regime 3 (FR3) is defined by high values of v/h relative to τ _d (Burgess and others, 2005) and by an inverse relationship between the variables that suggests that resistance to flow is controlled by friction at the glacier bed (Burgess and others, 2005). This indicates increased lubrication of the glacier bed by meltwater. Belcher, Fitzroy and North and South Croker glaciers show a progression from FR2 to FR3 as velocities increase towards their termini. Driving stress decreases down-glacier on these glaciers, which suggests that increasing surface velocities towards the termini are associated with a decrease in basal friction. The transition from FR2 to FR3 on Belcher and North Croker Glaciers coincides with areas of intense crevassing. Surface meltwater streams terminate in these crevasse fields, suggesting that penetration of surface meltwater to the glacier bed in these regions may enhance basal sliding, at least in the summer months (Burgess and others, 2005).

Flow regime 4 (FR4) is characterized by consistently low values of τ _d and large, but variable, values of v/h. This suggests very low basal friction, perhaps due to the presence of weak subglacial sediments, the deformation of which may contribute to basal motion (Burgess and others, 2005). The low mean value of τ _d, and the high sensitivity of v/h to small changes in τ _d are suggestive of flow by plastic deformation of subglacial sediments with a low yield strength (Kamb, 1991; Iverson and others, 1998; Cuffey and Paterson, 2010). Observations from Trapridge Glacier, Yukon, Canada, suggest that the yield strength of till is inversely related to the subglacial water pressure (Kavanaugh and Clarke, 2006; Cuffey and Paterson, 2010). Flow regime 4 is typically found at the terminus of outlet glaciers, including Fitzroy Glacier. However only two glaciers considered in this study (South Croker and North Croker) have large sections of FR4 near their termini (Fig. 2).

Fig. 2. Mapped flow regime distribution across Devon Ice Cap, adapted from Burgess and others (2005).

Thus, the distribution of flow regimes 2–4 is interpreted as reflecting the distribution of basal sliding and subglacial sediment deformation across Devon Ice Cap (Burgess and others, 2005). The sequence in which the four flow regimes occur along flowlines generally suggests a down-glacier progression from ice frozen to the bed in the interior of the ice cap (FR1), ice softening with transport and warming by shear heating to the onset of basal sliding (FR2), and a further increase in basal sliding as surface meltwater reaches the bed (FR3). In some regions, flow by subglacial sediment deformation may occur (FR4). However, neither this interpretation of flow regimes, nor the proposed link between flow regimes and surface hydrology, has been tested using independent evidence.

Methods

Surface drainage mapping

Exploring the relationship between supraglacial drainage characteristics and the distribution of areas of inferred basal sliding requires mapping of supraglacial channel networks, supraglacial lakes, drainage sink points and crevassed areas. These features were digitized manually from imagery acquired by NASA/USGS (US Geological Survey)’s Landsat-7 Enhanced Thematic Mapper Plus (ETM+) sensor on 7 July 1999. Surface drainage channels were delineated using both the 15 m panchromatic band and a 30 m red, green, blue (RGB) false-color composite (bands 5,4,2). Sink points were identified as points where surface drainage channels ended abruptly without re-establishment further down-glacier. Not all streams have associated sink points, as some streams appear to re-establish down-glacier, likely vanishing for short sections under snow bridges, and other supraglacial streams are routed to the glacier termini or margins. Whilst the 15 m panchromatic band of Landsat-7 was of sufficient resolution to identify larger meltwater channels, supraglacial meltwater features that were <15 m were not detected using this imagery. Lakes larger than 30 m × 30 m were classified automatically using the supervised spectral angle mapping (SAM) classification in ENVI 4.3, based on a spectral library that was created using training classes and a pixel purity index created from the 7 July 1999 Landsat-7 ETM+ scene. The classification output was checked manually for consistency, and corrections were made in ArcGIS for any non-lake areas that were classified as lakes or visible lakes that were not detected by the classification. Although the location of supraglacial lakes is likely constant from year to year, as their locations are controlled by bedrock topography (Lampkin and Vander-Berg, 2011; Sergienko, 2013), Landsat-7 ETM+ imagery from 2000–11 was inspected over Devon Ice Cap to look for changes in the distribution, density or structure of supra-glacial drainage across the ice cap.

It was not possible to calculate lake area change throughout the 1999 melt season, due to the low temporal resolution of the available Landsat-7 ETM+ imagery (which was partly due to cloud cover over Devon Ice Cap). Unfortunately it is not possible to use Moderate Resolution Imaging Spectroradiometer (MODIS) imagery to extend this analysis, due to the pixel size of MODIS imagery, which is too large to capture most of the relatively small lakes on Devon Ice Cap.

Gradient correlation

In order to explore the year-to-year variability in annual mean ice velocity and how this relates to the ice-flow regimes defined by Burgess and others (2005), an 8 year record of annual mean ice velocity was produced from sequential clear-sky panchromatic band (15 m resolution, channel 8 (520–900 nm)) images from Landsat-7 ETM+. Summer images were used where possible, to ensure the minimum amount of snow cover in each scene. Images were chosen to be approximately one calendar year apart, as long as cloud cover allowed. Georectified sequential image pairs were processed using an automated gradient correlation algorithm developed by Haug and others (2010) to find displacement distances of ice surface features (e.g. crevasse fields and rock debris). This method is insensitive to the scan line striping found in Landsat-7 ETM+ imagery acquired after 2003, making it highly suitable for this study.

Orientation images were first created from the original georectified image pair, based on image intensity differences in both the x- and y-directions (Haug and others, 2010). The complex orientation images, f _o and g _o, are created from

(2)

where

where sgn is the signum function, i is the complex imaginary unit, f the image at t = 1 and g the image at t = 2 (Haug and others, 2010). The new complex orientation images, f _o and g _o, consist of a real and imaginary part, where the real matrix represents the image intensity differences in the x-direction and the imaginary matrix represents the image intensity differences in the y-direction.

The orientation images were then divided into 64 × 64 pixel windows (960 m × 960 m). This window size was suitable for all outlet glaciers studied, as the maximum scene-to-scene displacements (4–500 m) were captured by this window size. The spacing between matching windows was 15 pixels, giving a densely populated grid. The correlation surface, P(x, y), was computed using

(3)

where F ₀(u, v) is the fast Fourier transform (FFT) of the reference window from f _o(x,y), is the complex conjugate from the FFT of the search window from g _o(x, y) and IFFT is the inverse fast Fourier transform (Haug and others, 2010). Sub-pixel accuracy was achieved by fitting a parabolic function to the maximum displacement and two surrounding points using the methods of Haug and others (2010).

Displacements determined by the orientation correlation were converted to annual mean velocities, calculated as the measured displacement divided by the number of days between the image pair and multiplied by 365. Results were checked for reliability by ensuring that the pattern of ice motion was realistic, and that vectors of ice motion were aligned with the ice-flow direction expected from the surface topography. As a further check, the patterns and magnitudes of ice velocities were compared with velocity maps derived by Burgess and others (2005) and Van Wychen and others (2012). While the vast majority of matches produced realistic ice motion patterns and velocities, some erroneous matches were generated, due to small-scale differences in shadow and illumination between scenes. Matches over slow-flowing areas often produced inaccurate displacement directions. Any matches that were not aligned (±45°) with the local aspect of the glacier or had displacement magnitudes less than the uncertainty (17 m) were automatically removed. Any remaining erroneous matches were removed manually.

In order to test the hypothesis that ice-flow variability increases in areas close to locations where meltwater enters the glacier, velocity standard deviation maps were also created in ArcGIS 10 from the calculated velocity rasters, using data for each pixel for which more than two years of velocity data were available.

The terminus area change for each outlet glacier in this study was calculated by manually digitizing the position of the terminus of each outlet glacier in georeferenced sequential image pairs. The total area of terminus change was calculated in ArcGIS 10 by summing the total area of terminus retreat and advance across the length of the calving face between sequential image pairs.

Statistical significance was determined at the 95% confidence interval.

Error analysis

The errors associated with gradient correlation come from two different sources: the co-registration of the image pair (the systematic error; σ_{ν syst})_, and the feature-tracking process itself (the random error; σ_{ν rand})_. The systematic error associated with the co-registration was kept to <1 pixel, so with a pixel size, P _s, of 15 m and time separation, Δt, of 1 year, the systematic error for the velocity, σ_{v syst,} becomes P _s/Δt or ±15 m a⁻¹ (Berthier and others, 2003). The random error associated with the feature-tracking algorithm is <1 pixel, with accuracy ranging from 1/8 to 1/4 pixel (Argyriou and Vlachos, 2007). Taking a conservative estimate of the random error to be ±0.5 pixels, which corresponds to a ∼8 m a⁻¹ error, the total error associated with the procedure used is

(4)

(Berthier and others, 2003). This uncertainty is a significant source of error for ice with a velocity of ∼30 m a⁻¹, but less significant for outlet glaciers with much higher velocities (Berthier and others, 2003). Terminus region velocities of outlet glaciers in the study area vary significantly between years, exceeding 200 m a⁻¹ in some years while dropping below 100 m a⁻¹ in others. When the terminus region velocity is 200 m a⁻¹ the error is approximately ±8.5%. Velocities were checked for accuracy by including matches over stable ground, for which the orientation was inconsistent, and the calculated velocity was always <4 m a^−1.

Results and Discussion

Interannual velocity variability

Annual velocities of all outlet glaciers studied show year-to-year variability. Annual mean center-line velocities at the terminus of each glacier are shown in Table 1. For all outlet glaciers (except the surge-type glaciers, SE1 and SE2), the maximum standard deviation of annual velocity occurs near the glacier terminus. Velocity variability generally decreases up-glacier (Fig. 3), although there are localized regions with high standard deviations of velocity in the mid- to upper reaches of all the outlet glaciers studied (Fig. 3). There is no significant relationship between annual mean center-line velocity and annual terminus area change for any glacier in this study, except South Croker Glacier. However, significant differences in the mean and standard deviation of velocity are found between the four ice-flow regimes on Devon Ice Cap (Table 2).

Fig. 3. Standard deviation of mean annual velocity on Belcher Glacier, determined from all available annual velocity fields. Supraglacial drainage channels and the locations of their sink points are also shown.

Table 1. Image dates (dd/mm/yy) used for each image pair, with annual mean velocity at the center line at the terminus

Table 2. Mean and standard deviation (s.d.) of all annual mean surface velocities measured in each flow regime on Devon Ice Cap

Year-to-year velocity variability and supraglacial hydrology

Maps of supraglacial drainage patterns on Devon Ice Cap were created from 1999 Landsat-7 ETM+ imagery. They show surface drainage channels ending upstream of zones of high interannual velocity variability, and visible sink points are common in the mid- and upper reaches of most outlet glaciers (Fig. 3). A t-test shows a statistically significant difference between the standard deviations of annual velocity in areas around sink-point locations and similarly sized areas around randomly generated control points with the same elevation distribution. A higher standard deviation of velocity is found near sink-point locations (Table 3). This suggests that the magnitude of the interannual velocity variability on Devon Ice Cap is linked to proximity to water-input locations. Thus, year-to-year differences in water input to these regions may be affecting the subglacial water pressure and basal friction, leading to year-to-year differences in the amount of basal motion. This could be linked to annual variability in meltwater production, leading to interannual changes in the amount, rate or timing of surface meltwater delivery to the glacier bed. Changes in the hydraulic efficiency of the subglacial drainage system may also play a part in controlling the velocity response of glaciers to meltwater inputs (Kamb, 1987; Kavanaugh and Clarke, 2001; Bingham and others, 2003; Bartholomew and others, 2010). Meltwater that reaches an inefficient subglacial drainage system early in the melt season produces a larger velocity perturbation than meltwater reaching an efficient subglacial drainage system later in the melt season. This suggests that the timing of meltwater inputs in relation to the stage of subglacial drainage system development may be an important control on the magnitude of velocity variability (Kavanaugh and Clarke, 2001; Bartholomew and others, 2010). In addition, variations in the number of high-magnitude meltwater input events during a summer melt season may influence the annual velocity variability, as large volumes of meltwater can overwhelm the capacity of the subglacial drainage system, leading to velocity perturbations (Schoof, 2010; Danielson and Sharp, 2013).

Table 3. Results of Student’s t-tests for the significance of difference in the standard deviations of mean annual velocities in regions of Devon Ice Cap that are located close to sink points of supraglacial streams, and a randomly selected collection of points with the same distribution with respect to ice-cap surface elevation

There is a difference in the relationship between velocity variability and the proximity to sink points between regions above and below 600 m elevation. At higher elevations, velocity variability increases near sink-point locations and decreases in the vicinity of the randomly distributed points. Below 600 m elevation (typically closer to the termini of outlet glaciers, that are characterized by FR3 and FR4 and are extensively crevassed), there is little difference in velocity variability between sink-point locations and randomly selected points. This may be because surface meltwater can access the glacier bed more easily (and widely) in these crevassed areas, resulting in more uniformly high velocity variability. This suggests that velocity variability becomes more spatially homogeneous in regions with more water-input points and greater water transfer to the glacier bed. Alternatively, velocity variability near the termini may arise from local forcing at the terminus, such as tidal flotation, removal of sea ice and loss of buttressing, or dynamic readjustment to iceberg-calving events.

Flow regimes, supraglacial drainage and year-to-year velocity variability

Flow regime 1

There are distinct differences in the supraglacial drainage systems of regions with different flow regimes (Fig. 4). FR1 contains the largest number of supraglacial lakes (665), the largest lake area (5 742 900 m²) and the highest number of sink points (102) of all the flow regimes on Devon Ice Cap (Fig. 4). However, because it covers such a large area, it also has the lowest density of lakes, sink points and crevassed regions (Fig. 5). The plateau region of Devon Ice Cap accounts for the majority of the FR1 area, although localized regions of FR1 are also found in the mid- to upper regions of some outlet glaciers, and in the terminus regions of SE1 and SE2 Glaciers, where the majority of FR1 sink points are located. The high-elevation plateau of the ice cap contains 473 lakes with an average area of 9184 m², but no sink points have been found. Lakes in this region of the ice cap do not always drain (Fig. 6), and likely refreeze in situ at the end of the melt season. Lakes that do drain appear to drain across the surface of the ice cap, thus they do not appear to be a direct source of supraglacial meltwater to the subglacial drainage system.

Fig. 4. Relationships between (a) flow regime and the distribution of (b) sink points, (c) the number and total area of supraglacial lakes and (d) the number of crevassed regions and total crevassed area on Devon Ice Cap.

Fig. 5. Differences in (a) sink-point density, (b) crevasse density, (c) lake size and (d) lake density for different flow regimes on Devon Ice Cap.

Fig. 6. Panchromatic Landsat-7 ETM+ imagery from 21 July 2002 and 14 September 2002, showing supraglacial lakes at high elevation on the plateau above SE2 Glacier. These lakes are located in a region characterized as FR1. There is little change in area of two of the three lakes between the image dates, and ice forming in the center of the lakes in the second image suggests that the lakes refreeze in situ and do not drain either englacially or supraglacially.

FR1 has the lowest mean and standard deviation of annual mean velocities of all four regimes. This supports the interpretation of Burgess and others (2005) that ice is frozen to the bed in FR1, as large year-to-year differences in ice velocity are not expected where internal deformation is the predominant form of ice flow. The limited opportunity for supraglacial lake water to penetrate to the glacier bed in the plateau region of the ice cap is also consistent with this interpretation. Changes in surface meltwater production in this region are unlikely to lead to changes in ice velocity, as surface meltwater has limited access to the glacier bed.

Flow regime 2

FR2 has 0.18 lakes per km², a higher density than any other flow regime on Devon Ice Cap. However, the fraction of surface area covered by lakes in FR2 is lower than in FR3 and FR4 (Fig. 5). The large number of small lakes in FR2 may provide the water needed to establish hydrological connections to the bed. Such connections are established either when lakes drain over the glacier surface via supraglacial streams to sink points, by which water can access the bed, or when they drain suddenly through their floors. Lake-drainage pathways in regions of FR2 are detectable in Landsat-7 ETM+ imagery as large drainage channels that extend downstream from the lakes to sink-point locations (Fig. 7) in FR2 and FR3. This suggests that supraglacial lake drainage over the ice surface is more common than in situ lake drainage in FR2.

Fig. 7. Panchromatic Landsat-7 ETM+ imagery from 28 June 2002 and 31 August 2002 showing supraglacial lakes on the upper reaches of North Croker Glacier. The lakes are located in a region characterized as FR2. There is a clear decrease in lake size between the image dates, showing lake drainages at this location. Clearly defined drainage channels (indicated by arrows) extending down-glacier from the lakes suggest that lakes drain across the surface of the glacier via supraglacial streams.

The density of sink points more than triples between FR1 and FR2, suggesting that meltwater penetration to the bed via sink points is more widespread in FR2 than in FR1 (Fig. 5). However, the distribution of sink points is not uniform throughout FR2. Over 70% of sink points are found <600 m a.s.l. The low elevations of sink points and the observations of supraglacial lake drainage by over-ice flow in FR2 suggest that large volumes of surface meltwater do not reach the glacier bed in FR2. Rather, surface meltwater appears to flow supraglacially from lakes to isolated input locations at elevations close to the transition between FR2 and FR3. The limited opportunities for meltwater to access the bed in FR2 suggest that surface meltwater does not widely affect glacier flow in FR2. This supports the interpretation of Burgess and others (2005) that the reduction in ‘effective viscosity’ seen in FR2 is associated with the softening of basal ice due to transport and warming of basal ice to the pressure-melting point, which will make it more deformable and perhaps initiate basal sliding.

The mean ice velocity nearly doubles from FR1 to FR2, likely due to the onset of basal sliding as ice reaches the pressure-melting point at the bed (Burgess and others, 2005). Both the standard deviation of annual velocity and the variation of the standard deviation of annual velocity from the mean are greater for FR2 than FR1 (Table 2). This is possibly because sink points at lower elevations in FR2 allow surface meltwater to access the bed, where it can reduce friction and locally enhance basal sliding. Year-to-year variability in velocity will occur in the regions surrounding sink points, as rates of seasonal sliding are sensitive to changes in the supply of surface meltwater to the bed (Danielson and Sharp, 2013). Variations in proximity to sink points may account for the increase in the standard deviation of the standard deviation of velocity in FR2, as isolated regions of high interannual velocity variability are located close to sink points.

Flow regime 3

FR3 has the highest density of sink points, suggesting that meltwater delivery from the surface to the bed may be widespread throughout this flow regime (Fig. 5). Furthermore, the near-tenfold increase in the density of crevasses between FR2 and FR3 suggests that crevasses may provide a means of transferring surface meltwater to the bed in FR3. Lakes are less abundant in FR3 than in FR2, but both the mean lake size and fractional lake area coverage are greater than in FR2 (Figs 4 and 5).

The inferred onset of widespread basal sliding in FR3 is suggested by the inverse relationship between v/h and τ _d (Burgess and others, 2005). This is consistent with the observed differences in supraglacial hydrology between FR2 and FR3, which allow greater delivery of surface meltwater to the bed, which likely enhances basal sliding. The onset of crevassing and the increased area coverage of supraglacial lakes, which are often coincident in location with crevasses, suggest that crevasse hydrofracture may generate supraglacial/subglacial connections in FR3. Crevasse hydrofracture requires both the existence of fractures that can propagate and large volumes of water to drive crevasse propagation to the bed (Van der Veen, 2007), both of which exist in this region. In FR3, lakes that decrease in size during the melt season (implying lake drainage) do not appear to be drained by large supraglacial channels as they are in FR2. This suggests that they drain through their floors by crevasse hydrofracture (Danielson and Sharp, 2013), rather than by runoff over the glacier surface (Fig. 8). Lake-drainage events form moulins that remain open for the remainder of the melt season (Das and others, 2008), and can increase the number of locations at which surface meltwater can access the glacier bed, leading to more widespread basal sliding in FR3 than in FR2.

Fig. 8. Imagery of South Croker Glacier from 28 June 2002 and 31 August 2002. Supraglacial lakes identified by red circles are in regions of FR4, while green circles show a lake in a region of FR3. The approximate location of the transition from FR3 to FR4 is shown by the blue curve. There are clear differences in the size of lakes between the beginning and end of the melt season, suggesting lake drainages occurred. There are no visible supraglacial streams draining the lakes in the FR4 regions, suggesting that lake drainages in these regions are likely occurring englacially through the ice beneath the lakes, rather than across the surface of the ice.

FR3 has the highest mean annual velocity of all the flow regimes, and the standard deviation of annual velocity is higher than in FR1 or FR2. This supports the suggestion of Burgess and others (2005) that FR3 is coincident with widespread basal sliding driven by increasing surface water inputs to the bed (especially since many of the sink points in FR2 lie immediately upstream of the boundary between FR2 and FR3). Year-to-year differences in the production of surface meltwater may lead to differences in the amount of basal sliding in this region. Lakes that form here are larger than those in FR2, and they appear to drain suddenly as a result of crevasse hydrofracture on the lake floor. Sudden lake-drainage events on the Greenland ice sheet have been linked to regional increases in ice velocity (Das and others, 2008; Joughin and others, 2013). Year-to-year differences in the amount of surface melting will alter the volume of meltwater contained within these lakes. This may lead to year-to-year differences in the ability of lakes to drain suddenly, as the crevasse-filling rate controls whether or not fracture propagation to the bed will occur (Van der Veen, 2007; Danielson and Sharp, 2013). This may provide an additional source of year-to-year velocity variability, as multiple sudden lake-drainage events provide a large input of surface meltwater to the glacier bed, causing periods of increased velocity in some years, but not in others.

The increase in variability of annual mean surface velocity from FR2 to FR3 suggests an increase in the variability of surface meltwater input to the bed. This could be caused by year-to-year differences in summer air temperatures and melt intensity. Annual melt rates increase down-glacier in response to higher temperatures at lower elevations, and a reduction in mean surface albedo, as the length of time for which bare ice is exposed each summer increases. In low-snow years, when the ablation zone is snow-free for a larger portion of the melt season, the amount of surface meltwater being transported downstream will be high, as melt rates will be high across the catchment area, due to the lower albedo of bare ice. The converse is true in high-snow years or years with cool summers. Thus, the increase in velocity variability between FR2 and FR3 may also be due, in part, to the increasing contrast in surface meltwater availability between low- and high-snow years, which is likely to be accentuated at lower elevations because of stronger surface albedo feedback effects.

Flow regime 4

The most striking feature of FR4 is the very high crevasse density in this regime, where crevasses cover >50% of the area. Distinct sink points are less apparent in FR4 than in the other regimes (Fig. 4). This is likely due to the extensive crevassing, which prohibits the development of supraglacial stream networks. Streams are not readily visible on the ice surface, probably because water is draining englacially through crevasses, from which it may drain to the bed. The increase in crevassing from FR3 to FR4, in conjunction with the relatively low number of sink points, suggests that the distributed drainage of surface meltwater through crevasses is likely the dominant form of meltwater delivery to the bed in this flow regime.

FR4 has the second highest density of supraglacial lakes (after FR2), and the largest average lake size. The large area of and, probably, volume of water stored in the supraglacial lakes, combined with the extensive crevassed areas in close proximity to (or beneath) these lakes may favor hydrofracture of crevasses to the bed in this region, causing supraglacial lake-drainage events that provide large inputs of water to the subglacial system. Figure 8 shows supraglacial lake-drainage events that occurred on South Croker Glacier in 2002. The mechanisms of lake drainage appear to differ between FR2, FR3 and FR4, in that lakes in FR3 and FR4 do not drain supraglacially via stream networks to sink points, but rather drain by crevasse hydrofracture through their floors.

Widespread lake drainage via crevasse hydrofracture throughout FR4 suggests that large volumes of supraglacial water sporadically reach the glacier bed in such regions. Large and variable inputs of water delivered to the subglacial system may enhance basal sliding via hydraulic jacking (Iken and Bindschadler, 1986; Bingham and others, 2003), as they can overwhelm the capacity of the subglacial drainage system (Schoof, 2010). High subglacial water pressures in FR4 may also enhance subglacial sediment deformation, as the yield stress and deformation rate of till are both dependent on water pressure (Kamb, 1991; Iverson and others, 1998). Conversely, the high variability in annual velocity observed in FR4 may be due to a decrease in annual velocity, driven by decreasing winter velocities following summers with high melt rates. This has been observed in Alaska, where there are increases in the efficiency of the basal hydrological system, driven by high summer melt rates, and a reduction in the amount of basal water and basal sliding through the fall and winter months (Truffer and others, 2005; Burgess and others, 2013).

The highest standard deviation of annual mean velocity occurs in FR4. This is consistent with (though not proof of) the interpretation of Burgess and others (2005) that sediment deformation contributes to glacier motion in this region. Furthermore, the major outlet glaciers draining Devon Ice Cap are grounded below sea level at their termini (Dowdeswell and others, 2004). If subglacial conduits are connected to the ocean then it is likely that basal water pressures are elevated in the regions grounded below sea level. If the basal water pressure near the terminus approaches or exceeds the critical value for flotation and/or sediment deformation, then small year-to-year differences in driving stress may have a dramatic effect on glacier velocities. Alternatively, the high year-to-year variability in glacier velocity found in FR4 could be a result of changes near the glacier terminus, such as thinning and flotation, that reduce back-stress and lead to enhanced surface velocities and terminus retreat (Moon and Joughin, 2008; Nick and others, 2009). However, terminus area change rates calculated from sequential image pairs do not show any significant correlation with terminus velocity at any glacier considered in this study, except South Croker Glacier. Changes in ice thickness at the terminus, submarine melt rates and year-to-year differences in the thickness and extent of proglacial sea-ice melange are unknown, but have been shown to influence terminus velocity in other regions (Holland and others, 2008; Amundson and others, 2010; Lindsey and Dupont, 2012). These influences may contribute to the year-to-year variability in ice velocity at the termini of the glaciers in this study, but we have no evidence to support this.

Conclusions

There are distinct differences in the supraglacial drainage characteristics of regions of Devon Ice Cap with different flow regimes. These are linked to variations in the density of locations at which surface runoff can enter the glacier and potentially reach the glacier bed. The mechanisms of meltwater delivery to the bed vary between flow regimes: little opportunity for meltwater penetration in FR1; surface runoff to localized sink points at low elevations in FR2; flow into widespread sink points, or drainage by crevasse hydrofracture, in FR3; and drainage by crevasse hydro-fracture (and perhaps rapid supraglacial lake drainage) in FR4. These differences in supraglacial hydrological characteristics between flow regimes are consistent with the interpretation of flow regimes proposed by Burgess and others (2005), in that they suggest that the distribution of regimes with flow by basal sliding and/or sediment deformation is, to some extent, controlled by the spatial pattern of surface meltwater delivery to the glacier bed. Furthermore, a feedback may exist whereby changes in basal thermal regime and/or basal substrate from FR1 to FR4 lead to increases in crevasse density, which, in turn, affect the density and distribution of both supraglacial meltwater drainage and meltwater sink points, which might lead to changes in basal temperature due to cryohydrologic warming (or cooling) (Phillips and others, 2010). Different flow regimes, characterized by differences in the volume and mechanism of meltwater delivery to the glacier bed, display distinct differences in the degree of year-to-year variability in annual mean velocity. As flow regime goes from FR1 to FR4 and the inferred contribution of basal velocity to overall ice velocity increases, so too does ice interannual velocity variability.

The results presented here support the argument of Burgess and others (2005) that the distribution of regions of Devon Ice Cap where basal sliding or sediment deformation contributes to ice surface velocities is linked to the distribution of locations where supraglacial meltwater is delivered to the glacier bed.