Ice tectonic deformation during the rapid in situ drainage of a supraglacial lake on the Greenland Ice Sheet

. We present detailed records of lake discharge, ice motion and passive seismicity capturing the behaviour and processes preceding, during and following the rapid drainage of a ∼ 4 km 2 supraglacial lake through 1.1-km-thick ice on the western margin of the Greenland Ice Sheet. Peak discharge of 3300 m 3 s − 1 coincident with maximal rates of vertical uplift indicates that surface water accessed the ice–bed interface causing widespread hydraulic separation and en-hanced basal motion. The differential motion of four global positioning system (GPS) receivers located around the lake record the opening and closure of the fractures through which the lake drained. We hypothesise that the majority of discharge occurred through a ∼ 3-km-long fracture with a peak width averaged across its wetted length of ∼ 0.4 m. We argue that the fracture’s kilometre-scale length allowed rapid discharge to be achieved by combining reasonable water velocities with sub-metre fracture widths. These observations add to the currently limited knowledge of in situ supraglacial lake drainage events, which rapidly deliver large volumes of water to the ice–bed interface.


Introduction
Variations in the delivery of surface water to the base of the Greenland Ice Sheet induce fluctuations in ice sheet velocity on inter-annual, seasonal, diurnal and sub-diurnal time scales (Zwally et al., 2002;van de Wal et al., 2008;Bartholomew et al., 2010Bartholomew et al., , 2011Shepherd et al., 2009;Das et al., 2008;Hoffman et al., 2011). For water to access the bed of the Greenland Ice Sheet, a hydraulic pathway must first be established through often kilometre-thick ice. Given sufficient water supply, a water-filled fracture will propagate to the base of an ice mass when the overburden stress at the fracture tip is offset by the density contrast between ice and water (Weertman, 1973;Alley et al., 2005;van der Veen, 2007). Krawczynski et al. (2009) calculated that supraglacial lakes 250 to 800 m across and 2 to 5 m deep contain sufficient water to drive a fracture to the base of kilometre-thick ice. Many lakes on the Greenland Ice Sheet attain this size or larger (Box and Ski, 2007;Echelmeyer et al., 1991;Georgiou et al., 2009), of which a small proportion (13 % between 2005-2009) drain in less than 2 days (Selmes et al., 2011). Surface lakes can drain rapidly into moulins via supraglacial rivers; however, many drain by the in situ propagation of hydraulically driven fractures (e.g. Das et al., 2008), a process hereafter termed lake tapping.
Lake tapping events provide water fluxes that exceed the capacity of the subglacial drainage system, leading to transient high subglacial water pressures, hydraulic jacking and ice sheet acceleration (Bartholomaus et al., 2008;Das et al., 2008;Pimental and Flowers, 2010). The integrated effect of multiple lake tapping events, and continued water flow into the hydraulic pathways they create, has the capacity to  impact the annual ice flux in future years, especially as, in a warming climate, lakes are expected to form and drain earlier in the season (Liang et al., 2012) and at higher elevations (Howat et al., 2012). It is, however, uncertain whether this increase in water delivery will enhance the annual ice flux through a net increase in basal lubrication (e.g. Zwally et al., 2002), or reduce it due to an earlier seasonal transition to an efficient subglacial drainage system (e.g. Sundal et al., 2011). Whilst supraglacial lake-tapping events represent major perturbations of the subglacial hydrological system and provide natural experiments for process-based investigations of glacier hydromechanics, few studies have succeeded in capturing them. Boon and Sharp (2003) measured premonitory drainage events preceding the complete and rapid (< 1 h) drainage of a 6.9 m deep supraglacial pond through a 150 m thick Arctic glacier on Ellesmere Island. Das et al. (2008) observed horizontal and vertical ice motion and seismicity during the rapid (∼ 1.4 h) tapping of a large supraglacial lake through 980 m of ice to the bed of the Greenland Ice Sheet. Our study presents detailed measurements of ice motion, lake volume change and seismicity capturing the rapid tapping of a large supraglacial lake through 1.1-km-thick ice on the western margin of the Greenland Ice Sheet.

Field site and methods
The field site, Lake F (67.01 • N, 48.74 • W), is located 70 km from the terminus of Russell Glacier, West Greenland ( Fig. 1). Russell Glacier catchment represents a typical landterminating sector of the Greenland Ice Sheet which, during the melt season, accumulates supraglacial lakes (McMillan et al., 2007) and shows diurnal (e.g. Shepherd et al., 2009) and seasonal (e.g. Bartholomew et al., 2010) covariations in ice velocity and surface uplift.
Remote sensing observations indicate annually repeating growth and rapid drainage of Lake F in the same location, with a mean date of drainage between 2002 and 2009 of 14 July. During the abnormally warm melt seasons of 2003(see Cappelen et al., 2001Cappelen, 2011) Lake F tapped 3 to 4 weeks earlier compared to other years in the 2002 to 2010 period.
On 26 June 2010, Lake F was instrumented with four global positioning system (GPS) receivers, six seismometers and two water-level sensors. Five days later at 01:40 on the 30 June, the lake drained completely in ∼ 2 h. Following drainage, the bathymetry of the lake bed and the locations of fractures and moulins were surveyed. Lake volume (V ) and discharge (Q) are estimated by combining water-level measurements with the lake bathymetry. All times are expressed in Coordinated Universal Time (UTC), which is three hours ahead of West Greenland Time.

Measurements of ice motion
Four dual-frequency GPS receivers (Trimble R7 and Leica SR520) were installed surrounding Lake F (Fig. 2). GPS antennae were installed on 6 m poles drilled 5 m into the ice, which subsequently froze in and thereafter provided a record of 3-D ice surface motion. Data, sampled continuously at a 10 s interval, were processed against a bedrock-mounted reference station using the differential phase kinematic positioning software, Track v. 1.24 (Chen, 1998;Herring et al., 2010), and final precise ephemeris from the International GNSS Service (Dow et al., 2009). The reference station was located 1 km from the terminus of Russell Glacier giving baseline lengths ≤ 70 km. To improve solutions at the day boundary, 36-h files were processed and the first and last six hours of the output position time series were discarded. Assuming steady ice motion, uncertainties in the positions were estimated at < 0.02 m in the horizontal and < 0.05 m in the vertical by examining the detrended position time series for GPS NW over a 2-day period in May 2011. The output position data are characterised by high-frequency noise caused by receiver and data processing errors. To suppress this highfrequency noise without causing a shift in phase, positions were filtered with a 1-h centred moving average. To quantify differential motion between GPS receivers, relative inter-GPS separation and the rate of separation were calculated from the filtered positions at a 10-min interval.

Measuring seismic activity
The passive seismic array consisted of six GS-11-D geophones with a natural frequency of 4.5 Hz and a bandwidth of 5 to 1000 Hz, continuously recording micro-seismic velocity at a 1 kHz sampling rate on to a RefTek-130 digitiser. To improve coupling with the ice, the geophones were mounted on 15 kg concrete slabs (dimensions = 0.5 × 0.25 × 0.05 m), buried to a depth of ∼ 0.5 m and routinely reset every 3 to 5 days before they melted out. The limited number of seismic stations together with the large array aperture (1 to 2 km) and the high rate of seismicity make it difficult to correlate particular onset times with individual events. The inability to identify individual events prevents the location of seismicity using standard methods (e.g. Walter et al., 2008;Roux et al., 2010). Whilst Jones et al. (2013) demonstrate that a technique based on amplitudes may be used to locate such seismicity, the application of this method is beyond the scope of this study.
In this study, the normalised root mean square (RMS) amplitude was calculated for each seismometer using an envelope function with 1-min time windows after applying a 2-pass, 4-pole Butterworth filter. The 5 to 50 Hz passband of the Butterworth filter was selected to reduce both highfrequency noise (> 50 Hz) associated with surface crevassing (e.g. Neave and Savage, 1970) and any low frequencies (< 5 Hz) associated with the instrument response. To identify step changes in seismicity, we calculated the normalised cumulative (seismic) energy from the RMS amplitudes.

Measurements of lake dynamics
Two pressure transducers (P1 and P2, Solinst M15 Levelogger) were installed in Lake F, logging pressure in metres of water head at two minute intervals with a specified accuracy of ±1 cm (Fig. 2). The records of P1 and P2 were compensated for changes in atmospheric pressure using a third Solinst Levelogger located at GPS NW.
Post-tapping, six transects were surveyed across the lake bed with a Leica SR520 GPS receiver recording at 1 Hz. The transects were differentially corrected and interpolated to form a digital elevation model (DEM) of the lake bathymetry with a grid spacing of 10 m. Time series of lake volume (V ) and discharge (Q) were calculated by combining the DEM with the water-level data. We estimate the uncertainty in the lake volume and discharge at ± 8400 m 3 and ±130 m 3 s −1 , respectively.
To extend the lake volume record, a time series of Lake F volume was estimated from daily-acquired atmosphericallycorrected Moderate-resolution Imaging Spectroradiometer (MODIS) images by applying the method of Box and Ski (2007). Uncertainty in this method was estimated at ±15 % by comparing MODIS-derived lake volumes with an independently collected lake bathymetry dataset.
To investigate the extent and timing of rapid-draining lakes within the Russell Glacier catchment, an automatic lake classification was applied to daily-acquired, cloud-free MODIS images. Lakes were classified using the Normalised Difference Water Index (NDWI) following the method described in Huggel et al. (2002). An empirically determined NDWI threshold was used to distinguish between water and other objects with a low NDWI (e.g. ice with a low albedo). The lake classification was trained using lake perimeter measurements derived from Landsat 7 images and differential GPS surveys. In combination with the NDWI threshold, thresholds for both the red and blue bands were used to further reduce misclassification of pixels with a similar spectral signal to water. Images with partial cloud cover were manually inspected. We define rapid draining as lakes that disappear within a 4-day interval, and the date of drainage as the midpoint between the date it was last seen and the date it disappeared. Typically, slow-draining lakes took much longer to drain than the 4-day threshold, and there is a clear distinction between lakes that drain rapidly and those that drain slowly. Figure 1 shows the date of drainage of rapidly draining lakes and the maximum extent of all lakes across the Russell Glacier catchment in 2010. The drainage network shown in Fig. 1 was created using hydrological modelling software (ArcGIS hydrological toolkit) from a 30 m resolution DEM derived from Système Pour l'Observation de la Terre (SPOT) data acquired on 2 July 2008.

Mapping hydraulic potential gradients
Basal and surface elevation DEMs collected by skidoo-based radio echo sounding following the method of Pettersson et al. (2011) were used to calculate the gradients of hydraulic potential, assuming basal water pressures were everywhere equal to the ice overburden pressure (Shreve, 1972). The resulting hydraulic potential gradients ( Fig. 11) can only be used to approximate the direction of subglacial water flow due to the variability in subglacial water pressure and basal hydraulic conductivity.

Regional scale lake dynamics
Within Russell Glacier catchment in 2010, 45 % of the lakes were classified as rapid (< 4 days) draining. The earliest rapid drainage occurred in late May and the latest in mid-July and in general lower-elevation lakes drained earlier than those located at higher elevations. Rapid lake drainage events occur in clusters; multiple lakes within the same elevation band drain simultaneously (Fig. 1). The rapid tapping of Lake F coincided with the rapid drainage of several adjacent lakes and an isolated peak (30 June to 2 July) in the discharge of catchment-wide proglacial Watson River (see

Formation and drainage of Lake F
In 2010 Lake F began to form on 5 June, attaining its maximum extent on 24 June with an area of 4.5 km 2 and a vol- ume of 1.8 × 10 7 m 3 (Fig. 3). Following the installation of the pressure transducers on 26 June the lake volume steadily decreased at a mean rate of 13.8 m 3 s −1 from 1.1 × 10 7 m 3 to 7.4 × 10 6 m 3 immediately prior to rapid tapping. This period of low discharge amounts to a volume of 3.6 × 10 6 m 3 and could be entirely accounted for by supraglacial discharge into Lake Z and moulin M4 (see Figs. 1 and 2). On 28 June 2010, water from Lake F travelled through a slow-flowing series of elongate ponds before a < 1 m wide supraglacial stream fed the water into moulin M4. During the slow period of drainage, the majority of water left Lake F via a 5-m-wide supraglacial river feeding into Lake Z.
Rapid discharge (here defined as Q > 50 m 3 s −1 ), associated with the tapping of the lake via in situ fracture propagation, occurred between 01:40 and 03:15 on the 30 June 2010 with the discharge peaking at Q max = 3300 m 3 s −1 at 02:47. Lake F had completely drained by 03:50. Immediately prior to rapid tapping, Lake F had a maximum depth of 9.9 m, an area of 3.8 km 2 and a volume of 7.4 × 10 6 m 3 .

Observations
On 29 June 2010, a healed crevasse, similar to Fig. S5 of Krawczynski et al. (2009), was observed from the western margin of the lake running through the lake in a W-E direction. Closed moulins were observed in the approximate positions of M1 and M5. It is likely that these relic features were formed by lake tapping events in previous years.
At 04:50 on 30 June 2010, less than 2 h after the end of rapid discharge, ∼ 0.3 m deep standing water was observed in the centre of the lake, overflowing across the clean-cut On 1 July the location, dip and strike of fractures were surveyed. The main fracture, F1, was mapped for 3 km but extended beyond this as a thin (< 1 cm wide) crack. F1 and F2 were sub-vertical, dipping towards the north and west, respectively (Fig. 2). Differential vertical displacement was observed along F1 with the northern hanging wall displaced typically 0.1 to 0.3 m (but up to 1 m near the centre of the lake) above the southern foot wall. This structure can be interpreted as a reverse dip-slip fault and evidence for transverse (cross-flow) compression (Fig 5a). The largest vertical displacement was measured in the deepest region of the lake, ∼10 m east of M2.
Along F1 a number of ice blocks, detached from the ice surface, had subsided into the fracture or been uplifted by floatation (see Fig. 4a). The structure of the subsided blocks is that of a high-angle normal fault with a dropped 2 to 5 m wide graben (Fig. 5b). Normal faults are evidence for transverse (cross-flow) extension across F1 (Cloos, 1955;Price and Cosgrove, 1994). Similar supraglacial fracture structures associated with hydraulic fracturing were observed following the Skeiðárjökull and Sólheimajökull jökulhaups in 1996 and 1999, respectively (see Fig. 12 of Roberts et al., 2000).
Five moulins (M1 to 5) were also mapped on 1 July (Fig. 2). The largest moulin (M1), ∼ 10 m in diameter, was explored to a depth of 45 m below the surface (Fig. 4b). At 45 m depth M1 continued vertically downwards, but access was restricted by fallen blocks of ice. Water entering M2 could be heard to flow englacially along F1 at a shallow depth before descending vertically down M1.
The main ∼ 5-m-wide supraglacial river flowing into Lake F from the north was intercepted by a fracture, forming three moulins collectively named M3 (Fig. 2). The evolution of the M3 moulins was observed over a number of days and is consistent with Stenborg (1968). Initially, water flowed into the clean-cut fracture at three discrete points and began to cut channels due to the frictional heat of melting. Over time the channels became wider and deeper. The channel of the largest moulin incised the fastest and ultimately captured all the flow. A similar evolution was observed for M2, which also continued to receive water throughout the melt season. Moulins M1 and M5 were not connected to supraglacial streams after rapid tapping.

Ice displacement
The horizontal and vertical ice motion recorded by the GPS receivers preceding, during and following the lake tapping event are depicted on Fig. 6. The first abnormal GPS motion occurred on 29 June when GPS SE accelerated to the west, concomitant with ∼ 0.2 m of uplift. GPS NE and GPS SW also accelerate, albeit with a lower magnitude. The northern two GPS receivers (GPS NW and GPS NE) show anomalous motion in the north-south plane, including transient reverse ice flow at GPS NW. The mean daily vector for 29 June is of a lower magnitude for the western GPS stations compared to the eastern GPS stations, implying a compressive strain regime that is also evident in the decreasing separations between GPS NW-SE, NW-NE and SW-SE (Fig. 7a, b and f). A detailed time series of discharge (Q), rate of discharge (dQ/dt), uplift (Z) and rate of uplift (dZ/dt) is presented in Fig. 8. At 01:15 on 30 June 2010 the uplift rate at GPS SW suddenly increases, leading to a maximum 0.34 m of uplift at 02:09. In this period, GPS NW and GPS NE are uplifted by 0.07 m and 0.1 m, respectively, and there is no discernible uplift at GPS SE. Coincident with the start of rapid discharge at 01:40, the uplift rate at GPS NW increases. In the interval 01:40-02:00, the lake volume decreases by 1.75 × 10 5 m 3 (2.4 % of the lakes pre-tapping volume) and GPS NW and GPS SW move north-west ( Fig. 6a and b) with slight compression (Fig. 7c). At 02:00, GPS NW continues to move in the north-west direction, accelerating and reaching a maximum displacement to the north of 0.3 m. At the same time (02:00), GPS SW reverses in direction, moving to the south-west (Fig. 6b) coincident with a step increase in the discharge ( Fig. 8a and b), seismicity (Fig. 9), uplift rate at GPS NW (Fig. 8c), and separation rate between all the GPS receivers with the exception of GPS SW-SE (Fig. 7). Rapid separation continues until 03:00 when GPS SW reverses to the north (Fig. 6b) and GPS NW reverses to the south (Fig. 6a), causing the separation rates to become negative (Fig. 7). The maximum discharge Q max of 3300 m 3 s −1 occurs at 02:47 simultaneous with the peak uplift rate at GPS NW of 0.8 m h −1 (Fig. 8). The maximum relative separation between GPS NW-SW, NW-SE, NW-NE and SW-NE is attained at 03:00 (Fig. 7), simultaneous with a lull in seismicity across all six seismometers (Fig. 9) and a transient 1 m 3 s −2 increase in the discharge rate (dQ/dt) which remains negative (Fig. 8b)

Interpretation and discussion
On the basis of the observations described above, the rapid in situ tapping of Lake F on the 30 June 2010 can be decomposed into three episodes: (1) initial drainage (01:40-02:00); (2) fracture opening (02:00-03:00); and (3) fracture closure (03:00-03:15). The horizontal and vertical ice surface velocities during each episode are illustrated on Fig. 10. These episodes are bounded by the duration of rapid discharge and it should be noted that closure of fractures extended beyond 03:15 (see Fig. 7). The timings of each episode, together with the time of Q max , are indicated on Figs. 6, 7 and 8. Episode 1 begins at the onset of rapid discharge and ends when GPS SW changes direction to the south, causing north-south extension across the lake (Fig. 10a). In episode 1, a small (1.7 × 10 5 m 3 , 2.4 % of the lakes pre-tapping volume) amount of water drains coincident with uplift concentrated at GPS SW and an increase in seismicity for the western half of the seismic array (S1-3, Fig. 9). Little or no increase in seismicity is recorded by the eastern seismometers (S4-6, Fig. 9).
As the fractures open in episode 2, Q rapidly increases peaking at 02:47 (Q max = 3300 m 3 s −1 ) simultaneous to the maximum rate of uplift at GPS NW (dZ/dt = 0.8 m h −1 ), indicating that water rapidly attained the ice-bed interface causing hydraulic ice-bed separation. The divergence of GPS NW-SE, NW-SW, SE-NE and SW-NE (Fig. 7a,  e) are interpreted as the opening of the main fracture F1 (Fig. 10b). Likewise, short-term longitudinal (with-flow) extension between GPS NW and GPS NE of ∼ 0.2 m (Fig. 7b), involving the reverse motion of GPS NE commencing at 02:00 (see Fig. 6d), is interpreted as the opening of subsidiary fracture F2. The opening of F2 involves the displacement of GPS NE to the east up the bed slope (Fig. 11). As soon as discharge begins to decrease after 02:47 the force holding F2 open reduces and, aided by the bed slope, GPS NE reverses in direction to the west (Fig. 6d), closing F2. The circular path of GPS NE during lake tapping (Fig. 6d) can be interpreted as the combined effect of fractures F1 and F2 opening and closing. In episode 2, 90 % (6.7 × 10 6 m 3 ) of the lakes pre-tapping volume drained compared to just 6.8 % (5.1 × 10 5 m 3 ) in episode 3. At 03:00, the boundary between episodes 2 and 3, the fractures attain their maximum width, and the short respite in fracturing is coincident with a lull in seismicity evident in the records for all six seismometers (Fig. 9). This quiescence suggests that seismicity is predominantly generated by the opening and closure of fractures and not by water flow which was continuing at a high rate (Q = 1450 m 3 s −1 ).
Immediately after 03:00, there was a step increase in seismicity particularly evident at seismometer S4 (Fig. 9j), as the inter-GPS separation rates between GPS 1-3, 1-4, 1-2, 3-4 and 2-4 became negative (Fig. 7) Fig. 9. Normalised seismic RMS amplitude (a-f) and normalised seismic cumulative energy (g-l) for the passive seismic array ordered anticlockwise from west (S1) to east (S6). The time of rapid (Q > 50 m 3 s −1 ) discharge (grey shade) and the time of maximum discharge (dashed vertical line) are indicated. closed (Fig. 10c). Episode 3 is characterised by decreasing discharge and uplift (Fig. 8). Following the end of rapid discharge, the rate of closure reduces, but it takes several hours before the inter-GPS separation rates stabilise (Fig. 7). We attribute the elevated seismicity post-03:15 ( Fig. 9) to continuing fracture closure (Fig. 10c) and the subsidence of GPS NW and GPS NE ( Fig. 8c and f). Following lake tapping, the inter-GPS separations show a stable increase in length between GPS NW-SW/blackbox [CE]For the consistency of all pages, please change the hyphens to endashes in the directions within Fig. 7 graphs and GPS SW-NE and a stable shortening between GPS SW-SE and NW-SE (Fig. 7). This supports the observation of juxtaposed extensional and compressional supraglacial fracture structures (Fig. 5). The separation between GPS NW-NE and GPS SE-NE reduces to ∼ 0 m by 12:00 on 30 June, suggesting the total closure of F2 and the most eastern section of F1, respectively. This observation is consistent with the small (cm-scale) fracture width measured in the field.
The characteristics of the rapid tapping of Lake F are consistent with Das et al. (2008), who observed the rapid (1.4 h) drainage of a 4.4 × 10 7 m 3 supraglacial lake through 980 mthick ice on the western margin of the Greenland Ice Sheet at 68.7 • N. The reverse and circular motion of the singular GPS station of Das et al. are comparable to those of GPS NE which was similarly located on westerly flowing ice, north of a longitudinal (with-flow) fracture and west of a transverse (cross-flow) fracture. The step increases in seismicity observed during the tapping of Lake F (Fig. 9) are also comparable (see Fig. S3 of Das et al., 2008).
Prior to lake tapping, a healed crevasse consistent with Fig. S5 of Krawczynski et al. (2009), was observed running from the western margin of the lake easterly towards its centre. Closed moulins were observed in the approximate positions of M1 and M5. It is likely that the healed crevasse and the closed moulins are relic features formed during tapping events in previous years. It is possible that the rapid tapping in 2010 was the re-opening of fractures and moulins formed in previous years; however, our observations do not reveal whether tapping involved the formation of a new fracture or the re-activation of a healed crevasse.

Initiation mechanism
Drainage of water into moulin M4, located to the west of Lake F, during the slow-discharge period prior to rapid tapping could theoretically cause localised uplift and acceleration leading to longitudinal (with-flow) extension. Alley et al. (2005) assert that the tensile stress caused by the acceleration of downstream ice may be important for initiating hydrofractures. However, discharge into 1-m-wide M4 was relatively low and prior to the tapping of Lake F there is no evidence for longitudinal (with-flow) extension (Fig. 7). On the contrary, the two western GPS receivers (GPS NW and GPS SW) were uplifted and accelerated several hours later than the two eastern GPS receivers (GPS SE and GPS NE), causing longitudinal (with-flow) compression across the lake in the hours preceding lake tapping ( Fig. 7b and f).
The observation of a compressive strain regime prior to lake tapping is in agreement with Krawczynski et al. (2009), who found that water-filled crevasses can propagate without longitudinal (with-flow) tension and that a given volume of water has the propensity to propagate a water-filled crack further in regions with less tension (or even slight compression), as thinner cracks require less water to remain water filled.
Although it is possible that small drainage events were masked by variations in supraglacial discharge, the lack of notable changes in the lake volume prior to 01:40 on 30 June 2010 (Fig. 8a) suggests that premonitory drainage events as observed by Boon and Sharp (2003) did not occur.  Fig. 11. Map of the hydraulic potential gradients and subglacial topography for Lake F. The black arrows are vectors of hydraulic potential gradient scaled by the metres of hydraulic potential change per metre. The gradients of hydraulic potential were calculated assuming basal water pressures were everywhere equal to the ice overburden pressure. The lake margin immediately prior to lake tapping is shown together with the locations of moulins, fractures and GPS receivers. The contour interval for the basal topography is 10 m.

Hydraulic pathways
We assert that during lake tapping, discharge occurred along most of the fractures' lengths. Assuming first that the inter-GPS separation is entirely a result of fracturing, and second following Krawczynski et al. (2009) that the fractures are parallel-sided with constant width, the fractures' crosssectional area at the time of peak separation (03:00) can be estimated. Interpolating the maximum separation between GPS NW-SW, NW-SE, SW-NE and SE-NE gives a peak fracture width averaged over F1's 2.6 km-wetted-length of 0.4 m and an estimated maximum cross-sectional area of 842 m 2 . This fracture width agrees well with the modelling results of Krawczynski et al. (2009) that suggest an idealised conical lake of a similar size to Lake F (diameter = 2.2 km, area = 3.8 km 2 ) would drain via a 0.4-m-wide fracture across the width of the lake in ∼ 2 h. The maximum cross-sectional area of F2 is estimated at 140 m 2 . Combined, F1 and F2 have a total cross-sectional area at 03:00 of 982 m 2 . Dividing the discharge at 03:00 (Q = 1450 m 3 s −1 ) by the combined cross-sectional area gives a mean flow velocity of 1.5 m s −1 . Hence, due to the length of the fractures, rapid discharge can be achieved by combining reasonable water velocities with sub-metre fracture widths. Assuming an ice thickness of 1100 m, the total volume of the fractures at 03:00 is 1.1 × 10 6 m 3 , or 15 % of the lake volume prior to tapping. The actual volume of water that drained between 01:40 and 03:00 was much larger at 6.8 × 10 6 m 3 . During rapid discharge the frictional heat of turbulent flow would preferentially melt the sections of the fracture transporting the greatest flux, leading to the concentration of flow (Walder, 1982) and the development of moulins. There is a clear distinction between moulins that formed during rapid discharge and moulins that formed afterwards. Immediately, post-tapping moulins M1 and M5 were fully formed but dry whilst moulins M2 and M3 were clean-cut fractures accepting water. We assert that discharge was initially concentrated down the largest (∼10 m diameter) moulin M1 (Fig. 4b).
When the water level fell below the elevation of M1 surface water flow would have been concentrated through F1 in the deepest region of the lake, ∼10 m east of M2 (Fig. 2).
At 04:50 on 30 June 2010, < 2 h after rapid tapping, there were no supraglacial channels in the lake bed. A shallow (∼ 0.3 m deep) pond of water in the deepest region of the lake was overflowing the clean cut edge of F1. Over a number of days, channels formed by the frictional heat of preferential water flow, concentrating water flow into discrete moulins. Continued water flow into moulins would work to keep moulins open for the remainder of the melt season. On 1 July, sections of the fracture were water filled and water draining into moulin M2 could be heard to flow laterally along F1 at a shallow depth before draining into moulin M1, suggesting that post-tapping large sections of the fracture were closed.

Vertical ice surface motion
Ice surface uplift is typically attributed to bed-parallel motion, vertical strain and ice-bed separation due to high subglacial water pressures (Hooke et al., 1989). Fault displacement (Fig. 5) may also cause vertical ice motion (e.g. Walder et al., 2005). All four factors must be considered when interpreting the vertical motion of a GPS receiver installed on the ice surface.
The bed slope within the study area is highly variable (Fig. 11) and may be responsible for some of the differential ice motion (Fig. 6). In contrast to the smooth vertical motion of GPS NW, GPS SE and GPS NE the vertical motion of GPS SW is characterised by sudden steps coincident with the start and end of the fracture opening episode. GPS SW is located on the strongest subglacial gradient of all the GPS receivers, south of a conical subglacial peak (see Fig. 11). Horizontal motion along the inclined bed slope can satisfactorily explain the vertical motion of GPS SW. At 00:00 on 30 June 2010 the trajectory of GPS SW was perturbed to the north-west, coincident with ∼5 cm of uplift. This vertical motion can be explained by north-westerly motion up the bed slope (Fig. 11). The subsequent southerly-motion of GPS SW down the bed slope between 02:00 and 03:00 is coincident with subsidence. This lowering is simultaneous with uplift at GPS NW, GPS SE and GPS NE, suggesting that the water delivered to the bed during rapid discharge did not access the bed beneath GPS SW. Finally, at 03:00 when GPS SW moves north up the bed slope there is a second low magnitude (10 cm) period of uplift exclusive to GPS SW. Sugiyama et al. (2008) observed the greatest uplift near to the drainage centre during the subglacial drainage of icemarginal lake Gornersee in Switzerland and we therefore assert, like Das et al. (2008), that surface uplift was likely greater near the centre of the lake. The highest-magnitude uplift observed in this study of 0.9 m by GPS NW is the most consistent with the 1.2 m of uplift measured by Das et al.. The bed slope underneath GPS NW slopes down towards the north-west (Fig. 11) so north-west acceleration during the fracture opening episode is not responsible for the observed uplift at GPS NW. Vertical strain cannot account for the uplift as extension (which would cause thinning and lowering) was observed for all the inter-GPS separations involving GPS NW during this episode (Fig. 7a, b, c). Fault displacement in the form of vertical offset of the fracture walls in a reverse dip-slip fault is attributed to compressional strain (Fig. 5b) which is unlikely to have occurred during the fracture opening episode when all inter-GPS separations were extensional (Fig. 7). As neither motion along an inclined bed slope, vertical strain or fault displacement can explain the observed vertical motion of GPS NW during the fracture opening episode, the uplift at GPS NW can be attributed to icebed separation resulting from high subglacial water pressures caused by the delivery of a large quantity (6.7 × 10 6 m 3 or 90 % of the lakes pre-tapping volume) of water to the icebed interface.

Subglacial water routing
We argue subglacial water delivered to the ice-bed interface along F1 would be preferentially routed through a subglacial valley to the north-west (Fig. 11) . Field measurements of the fracture structure and the differential motion of the GPS receivers support this assertion. Based on the locations of F1 and F2, we can conceptualise the ice mass structurally into three semi-independent blocks: the southern, north-eastern and north-western. The direction of dip of sub-vertical fractures F1 and F2, to the north and west respectively (see Figs. 2 and 5), together with the permanent offset of the northern hanging wall above the southern foot wall of F1, suggest that during the fracture opening episode, the northwestern block was preferentially uplifted and ejected to the north-west. This is consistent with the observation of the greatest horizontal and vertical motion of the north-western block on which GPS NW was located ( Fig. 6a and e) and the highest rate of seismicity recorded by the most western seismometers (S1-3, Fig. 9). In contrast to the slow subsidence of GPS NW, GPS SE and GPS NE following the end of rapid discharge, the vertical motion of GPS SW rapidly returned to a steady height (Fig. 8), suggesting that water was not routed or stored at the bed in this area.

Conclusions
Detailed measurements of ice surface motion, discharge and seismicity during the rapid in situ drainage of a large annually-tapping supraglacial lake through kilometre-thick ice on the western margin of the Greenland Ice Sheet contribute to the currently limited knowledge of rapid supraglacial lake tapping events.
Horizontal ice motion during lake tapping is dominated by ice tectonic deformation involving the transient opening and closure of multiple fractures. We assert that during rapid discharge, water flowed down most of the fractures' lengths. By reconstructing the peak cross-sectional area of the fractures from the inter-GPS separations, we find that the fractures' kilometre-scale lengths allow rapid discharge to be achieved by combining reasonable water velocities with submetre fracture widths.
The maximum uplift rate of 0.8 m h −1 occurred simultaneously with the maximum discharge of 3300 m 3 s −1 providing evidence that water rapidly attained the ice-bed interface, raising subglacial water pressures above overburden over a large area of the bed. Basal topography and the gradient of hydraulic potential exerted control on water routing, horizontal ice motion and uplift. The greatest horizontal displacement and vertical uplift was observed above the preferential subglacial drainage route. Lake tapping events rapidly deliver large pulses of surface water to the bed of the Greenland Ice Sheet, causing transient ice-bed separation and acceleration; however, it remains unclear what impact this water delivery will have on the annual ice flux.