Rapid development and persistence of efficient subglacial drainage under 900 m-thick ice in Greenland

Intensive study of the Greenland Ice Sheet’s (GrIS) subglacial drainage has been motivated by its importance for ice dynamics and for nutrient/sediment export to coastal ecosystems. This has revealed consistent seasonal development of eﬃcient subglacial drainage in the lower ablation area. While some hydrological models show qualitative agreement with ﬁeld data, conﬂicting evidence (both ﬁeld- and model-based) maintains uncertainty in the extent and rate of eﬃcient drainage development under thick ( ∼ 1 km) ice. Here, we present the ﬁrst simultaneous time series of directly-observed subglacial drainage evolution, supraglacial hydrology and ice dynamics over 11 weeks in a large GrIS catchment. We demonstrate development of a fast/eﬃcient subglacial drainage system extending from the margin to beneath ice > 900 m thick, which then persisted with little response to highly variable moulin inputs including extreme melt events and extended periods (2 weeks) of low melt input. This eﬃcient system evolved within ∼ 3 weeks at a moulin initiated when a fracture intersected a supraglacial river (rather than hydrofracture and lake drainage). Ice ﬂow response to surface melt inputs at this site follows a pattern commonly observed in the lower GrIS ablation area, and by assuming a strong relationship between ice dynamics and subglacial hydrology, we infer that eﬃcient subglacial drainage evolution is widespread under 900 m-thick ice in west Greenland. This time series of tracer transit characteristics through a developing and then persistent eﬃcient drainage system provides a unique data set with which to validate and constrain existing numerical drainage system models, extending their capability for simulating drainage system evolution under current and future conditions.


Introduction
The Greenland Ice Sheet's (GrIS) subglacial drainage system has been intensively studied over the last two decades, both to assess ice sheet response to increasing melt (e.g. Zwally et al., 2002;Sundal et al., 2011), and more recently to quantify nutrient * Corresponding author.
E-mail address: dcha@norceresearch.no (D.M. Chandler). and sediment export from the ice sheet to coastal/fjord environments (e.g. Sejr et al., 2014;Hawkings et al., 2015;Hopwood et al., 2020). Some aspects of the GrIS subglacial drainage system are now understood in considerable detail, with field and remotesensing data revealing a persistent, widespread seasonal pattern in the ablation area, characterised by inland-propagation of evolution from slow/inefficient to fast/efficient drainage (Nienow et al., 2017). Similar seasonal behaviour is observed on many other ice masses at smaller scales (e.g. Willis, 1995;Bingham et al., 2006). In Greenland, the inland propagation of efficient drainage is linked to the establishment of hydrological connections between the surface and subglacial drainage pathways under increasingly thick ice, enabled by fracturing and subsequent moulin development, sometimes associated with supraglacial lake drainage events (Nienow et al., 2017;Hoffman et al., 2018). Transient periods of extensional longitudinal stresses enable fracturing even in ice normally under compressive longitudinal stress (Christoffersen et al., 2018). The inland extent of efficient subglacial drainage is currently unclear Chandler et al., 2013;Meierbachtol et al., 2013;Hoffman et al., 2016;Christoffersen et al., 2018), but theory suggests this is limited by increasing ice thickness, decreasing hydraulic potential gradient, and decreasing ablation rate further inland (e.g. Röthlisberger, 1972;Schoof, 2010;Dow et al., 2014). Fast, efficient summer drainage reverts to slow, inefficient drainage during winter except very near the ice margin, as ice creep gradually closes subglacial conduits once melt inputs have ceased.
Owing to the inaccessibility of the GrIS subglacial drainage system, much of our current hydrological understanding is based on indirect observations. For example, borehole and moulin water pressure records, and corresponding diurnal to seasonal changes in ice surface motion, can reflect drainage system evolution via relationships between drainage system efficiency, subglacial water pressure, basal traction and ice motion (see references in Nienow et al., 2017). The observed response to melt inputs in Greenland is very similar to that in smaller ice masses (e.g. Hubbard et al., 1995). However, the impact of meltwater-driven ice acceleration has an uncertain impact on net annual ice displacement, particularly in the upper ablation area of the GrIS (Sundal et al., 2011;Sole et al., 2013;van de Wal et al., 2015;Doyle et al., 2014).
Predictions of how interactions between ice dynamics and hydrology will impact GrIS mass balance and nutrient/sediment export, especially under global warming scenarios, require numerical subglacial drainage models. Such models can qualitatively replicate the seasonal characteristics described above (e.g. Schoof, 2010;Hewitt, 2013;Werder et al., 2013;Banwell et al., 2016;Hoffman et al., 2016;de Fleurian et al., 2018;Koziol and Arnold, 2018), but quantitative evaluation of these models is severely limited by a lack of relevant field data with which to compare simulations. At present, a detailed time series of subglacial drainage evolution based on observations of water transit through the drainage system is not available from any part of the GrIS.
Here, we present the first simultaneous time series of drainage system tracing, supraglacial hydrology and ice dynamics during the melt season in a large GrIS catchment (Leverett Glacier, west Greenland). We undertook a series of traces at a moulin in ∼900 m thick ice , located 41 km from where its draining waters emerge at the GrIS margin (Fig. 1). We also collected simultaneous observations of ice surface motion, surface ablation, and moulin melt water input. Finally, discharge from the catchment was measured where the proglacial river emerges at the ice margin. The moulin chosen for this experiment is ideally situated in a region of the GrIS where there is currently disagreement regarding the seasonal extent of efficient drainage evolution (Hewitt, 2013;Werder et al., 2013;Banwell et al., 2016;Hoffman et al., 2016;de Fleurian et al., 2018;Koziol and Arnold, 2018). Our comprehensive time series of hydrology and ice dynamics provides valuable insights into how an efficient drainage system develops and persists under thick ice, and will provide an ideal dataset with which to validate subglacial drainage models.

Field site
Leverett Glacier is a land-terminating catchment of the GrIS at approximately 67 • N in west Greenland, with a theoretical hydro- water sampling site near Leverett Glacier terminus (yellow triangle), bed topography (shading) and surface topography (blue dashes, at 200 m intervals) from Lindbäck et al. (2014); and estimated hydrological catchment (white) and subglacial water routing (green) following Lindbäck et al. (2015) with k = 1 (pressure potential equal to ice overburden pressure) in Eq. 2 of their supplementary material. The estimated melt water routing between the moulin and terminus is marked in orange. The inset shows locations of the new Moulin L41A traced in this study, inactive Moulin L41 traced in the previous year by Chandler et al. (2013), and locations of GPS receivers. (b) Ice surface and bed topography along the estimated melt water route from the moulin to terminus as shown in (a). (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.) logical catchment estimated as ∼900 to 1200 km 2 Lindbäck et al., 2015). The hydrology and ice dynamics in this catchment have been studied for over a decade (e.g. van de Wal et al., 2008;Bartholomew et al., 2011;Palmer et al., 2011;Sundal et al., 2011;Chandler et al., 2013;van As et al., 2017) and follow the seasonal patterns common to other major catchments in Greenland, including the neighbouring Issunguata Sermia (Meierbachtol et al., 2013), the Paakitsoq region Andrews et al., 2014), south Greenland  and northeast Greenland (Neckel et al., 2020). Therefore, the hydrology and ice dynamics at Leverett Glacier can be considered to be widely representative of many of Greenland's outlet glaciers, as well as many smaller ice masses in the Arctic (e.g. Bingham et al., 2006). Data were collected during the 2012 melt season at a moulin (L41A) approximately 700 m south of Moulin L41 traced in 2011 (Chandler et al., 2013) (Fig. 1a). Moulin L41 reopened in 2012 but was inactive as the main supraglacial river draining that part of the ice sheet was captured by L41A. Ice thickness at L41A is estimated from radio-echo sounding  as 800 m, with an uncertainty of 18 m, and increases to 910 m within 2 km of the moulin along the estimated subglacial flow path (Fig. 1b). The ice thickness is between 800 and 920 m for the first 7 km of the estimated flow path.

Tracer experiments
We used sulphur hexafluoride (SF 6 ) as the tracer, which has been applied in a variety of environments (including in ocean-scale experiments) with high dilution (Watson et al., 1987;Clark et al., 1996;Chandler et al., 2013). This tracer can be detected at very low concentrations, does not adsorb onto sediments, and is nonreactive; therefore, it is advected at the same rate as the water in water-filled drainage systems as previously validated by compar-ison with dye tracing in the Leverett catchment (Chandler et al., 2013). We note the increasing environmental concerns regarding SF 6 , and while the amount required for each tracer experiment (3 kg) was very small in the context of ∼9 gG global emissions (Engel et al., 2018), alternatives should be sought if similar tracer experiments are conducted in the future.
A total of 11 traces were carried out between 16 June and 6 August, 2012. Tracers were injected at ∼19:00 hrs, shortly after the peak daily moulin melt water input (typically between 15:00 and 19:00 hrs). Tracer concentrations were monitored by collecting water samples from the proglacial river draining Leverett Glacier and were analysed by gas chromatography. The solubility of SF 6 in water is low (74 mg/l at 0 • C) (Bullister et al., 2002), and the choice of injection method is critical to minimise gas loss at the injection site. Following trials in previous years, the method adopted in 2012 involved direct injection through a 150 m, 5 mm internal diameter hose lowered into the moulin. The hose was connected to the gas cylinder via a precalibrated flow meter and regulator. By monitoring the flow rate, a known amount of gas was injected. The bottom of the hose was left open, and the presence of residual pressure in the hose when disconnecting the apparatus following the injection confirmed that the injection was into water rather than into air.
Proglacial river samples were collected in 120 ml gas-tight Wheaton vials (in duplicate or triplicate) sealed with rubber stoppers and a crimp top. Care was taken to ensure no gas bubbles remained in the bottle once it had been stoppered. N 2 headspace was added to the bottles in the field by piercing the stopper with a long needle pushed through to the bottom of the bottle; then, a 60 ml gas-tight syringe filled with N 2 at ambient temperature and pressure was used to add 20 ml of N 2 through a short needle just protruding through the stopper. Displaced water was ejected through the long needle. Bottles were shaken and left for >2 hours to allow headspace equilibration.
The SF 6 concentration in the headspace was analysed in the LOWTEX laboratories, Bristol University, UK, using a portable Cambridge Scientific 300-series gas chromatograph (GC) fitted with an electron capture detector (ECD). The headspace in the sample bottle was pressurised by injecting water into the sealed bottle through a long needle; a sharpened gas inlet tube connected to the GC was then used to pierce the stopper, thus injecting the pressurised headspace into the GC sample loop. Direct injection rather than transfer via a gas tight syringe achieved greater sensitivity than our previous method (Chandler et al., 2013). The GC used a carrier gas of CP-or ECD-grade N 2 at a flow rate of ∼18 ml/minute. Water and oxygen traps were installed in the carrier gas supply. The injector, column and detector temperatures were 150 • C, 100 • C and 200 • C, respectively. Two sample loops were used (1 ml or ∼0.05 ml), with the smaller loop being fitted when peak SF 6 concentration was expected to saturate the detector with the 1 ml loop. Samples were analysed in duplicate, and a 1 parts per billion by volume (ppbv) standard was used to correct for instrumental drift every 5 samples. All samples for an individual trace were analysed with the same sample loop. Specifically, traces 1, 2, 4, 5, 7, 8 were analysed with the small loop and traces 3, 6, 9, 10, 11 with the large loop. Bottles from the same trace were analysed together as a batch but in a random order. Peak area was converted to headspace SF 6 concentration by calibration with the 1 ppbv standard. Calibration with just one standard assumed a linear ECD response over the range of interest, as confirmed via dilution experiments (Fig. 2). Calculation of the corresponding SF 6 concentration in the river water sample prior to adding the headspace followed Chandler et al. (2013).
Errors in the tracer concentrations estimated empirically using duplicate bottles yielded mean percentage errors of 15% for the small sample loop and 11% for the large sample loop. The limit of detection (LOD) of the instrument (expressed as headspace concentration) was 0.8 parts per trillion by volume (pptv) for the large sample loop and 3 pptv for the small sample loop, equivalent to respective SF 6 concentrations of 1 × 10 -3 parts per trillion by mass (pptm) and 5 × 10 -3 pptm in the water samples.
A key characteristic of the tracer return curves is the tracer velocity, which quantifies the efficiency of the subglacial drainage system. Several measures of velocity are available including the maximum velocity and various measures of mean velocity (Chandler et al., 2013). Here we report both the maximum velocity and mean velocity for all traces, calculated using the straight-line distance between moulin and sampling site divided by the time when 5% or 50% of the recovered tracer has emerged (v 05 and v 50 , respectively). Due to the assumption of 'straight-line distance', the estimated flow velocities are minima as sinuosity in the flow path will cause the calculated velocity to be lower than the actual flow velocity. Errors in v 05 and v 50 were estimated in a Monte-Carlo simulation (with n = 1000) using the empirically-derived analytical errors.

Ice surface ablation
Surface ablation rates were measured daily using changes in surface height at five ablation stakes, arranged in a cross configuration at ∼2 m separation. The stakes were located ∼700 m north of Moulin L41A, and were in the supraglacial hydrological catchment feeding that moulin. The stakes were installed in holes deeper than the length of the stake (so each measurement was from the ice surface down to the top of the stake), to avoid the problem of enhanced surface melting caused by solar radiation absorbed by the stake . All ablation measurements were carried out by the same observer to ensure consistency, for example in the interpretation of the level of a rough ice surface.

Ice surface motion
Ice surface motion was recorded by five dual-frequency Leica SR520 GPS receivers deployed on poles drilled 2 m into the ice surface (Fig. 1). GPS data were post-processed kinematically (King, 2004) with Track v.1.27 software (Chen, 1998) against bedrockmounted reference stations using a precise ephemeris from the International GNSS Service (Dow et al., 2009). Reference stations were located 1 km from the terminus of Russell Glacier and at Kellyville, giving baseline lengths less than 41 km. Due to gaps in the time series caused by power outage, we averaged the horizon- tal velocities recorded at the five stations with the fewest gaps to give a single record. Positions were recorded at 30 s intervals; 1-hr means were then smoothed using a 5-point binomial filter. Since there was generally little difference in velocity between the stakes (Fig. 3), the mean velocity across the network gives a better indication of the seasonal pattern of ice motion with fewer gaps than in the individual records. Velocities are centred differences of hourly displacements, i.e. the velocity at time T i is (X i+1 -X i-1 )/(T i+1 -T i-1 ). GPS stakes required periodic re-drilling as they gradually melted out. Since re-drilling sometimes coincided with gaps in the data, a full season of vertical displacements is not available; the reliable uplift record is restricted to the early part of the season.

Supraglacial and proglacial river discharge
Supraglacial river discharge was measured by monitoring water depth with a custom-built pressure sensor, at 1-minute intervals, and converting depth to discharge with a rating curve established from salt dilution gauging. This method is well suited to supraglacial river gauging as the electrical conductivity of supraglacial melt water is very low. The pressure sensor was installed on the bottom of the channel and was able to move downwards with the ice surface as the channel gradually incised. For further details, see Wadham et al. (2016).
Proglacial river discharge was monitored using stage measurements (collected by a HOBO pressure sensor) and dye dilution gauging at a stable bedrock section near the terminus of Leverett Glacier, following Bartholomew et al. (2011) and Tedstone et al. (2013).

Results
The 11-week field campaign on the ice sheet commenced on 28 May 2012, prior to any visible evidence of melt water ponding. The surface at this time was characterised by seasonal snow and patches of bare ice . Snow melt over the subsequent week created extensive areas of slush and eventually melt water ponding in surface depressions and relict moulins, including Moulin L41 traced in 2011. A lake developed approximately 200 m up-glacier from L41, with an estimated maximum lateral extent of ∼500 m and an unknown depth. This lake was >200 m up-glacier from the new Moulin L41A, but the exact extent is unknown as it was not accessible on foot. On 31 May, the westwards ice surface velocity increased from 100-120 m yr -1 to 140-190 m yr -1 (Fig. 3). A period of frequent, audible cracking and the development of extensive surface fractures (generally oriented N-S and reaching ∼2 cm width) commenced on 3 June (Fig. 4). Subsequent inspection of this fracture zone showed that it extended from the vicinity of Moulin L41, southwards to Moulin L41A (Fig. 1a), and continued for a least 1 km further south of Moulin L41A.
We consider 3 June as marking the onset of the hydrologically active part of the season at L41A, which we divide into three distinct periods as described below.

Period 1: surface-bed connection (3-5 June)
During Period 1, surface melt water remained ponded and the audible cracking continued. The westwards ice motion accelerated further from 140-170 m yr -1 to 230 m yr -1 by 23:00 on 5 June. At this time, the water that had ponded in relic Moulin L41 and in other smaller relic moulins drained abruptly (Fig. 3). Drainage was accompanied by ∼0.3 m of vertical uplift. Both vertical uplift and westwards ice velocity peaked 12-24 hr after moulin L41 drainage.

Period 2: evolution from slow to fast drainage (5 June -10 July)
Moulin L41 received very little melt water input despite its proximity to extensive ponding, as the supraglacial channel feed- ing L41 was blocked with snow. Instead, the ponding melt water drained southwest along another supraglacial channel. We were unable to access the new moulin (L41A) fed by this channel until 15 June, as surface conditions prevented access on foot. Therefore, moulin L41A was first traced on 16 June, likely 11-14 days after ice-bed hydrological connection. L41, L41A and several other smaller moulins were all located in the linear fracture zone described above, and at least 200 m down-glacier from the lake (Figs. 1, 4, 5). There was no rapid lake drainage; instead, the lake gradually decreased in extent over the following month, presumably due to ongoing supraglacial channel incision. Trace #1 emerged very slowly over 3-4 days, with respective v 05 and v 50 velocities of 0.27 and 0.19 m s -1 ; subsequent traces (until Trace #6 on 10 July) emerged progressively quicker and with less dispersion, with v 05 and v 50 reaching 1.27 and 0.94 m s -1 , respectively (Figs. 6, 7).
Melt water input into moulin L41A showed strong diurnal cyclicity, with morning minima <3 m 3 s -1 and evening peaks reaching 5 -12 m 3 s -1 . An increase in melt input (most significantly, an increase in the diurnal minima) was associated with warm, cloudy conditions and corresponding rapid ablation exceeding 50 mm day -1 (Fig. 7).
Westwards ice velocity peaked at ∼260 m yr -1 on 6 June, shortly after L41 drained, then declined to ∼130 m yr -1 by 9th June during a period of low ablation and new snow. Diurnal cyclicity commenced abruptly on 16 June (Fig. 7). Surface motion generally decreased from 16 June -10 July, except for two relatively strong diurnal peaks on 20 June and 4 July, coincident with peaks in ablation and moulin melt input (Fig. 7). The relationship between moulin melt water input and ice velocity showed strong anticlockwise hysteresis with faster ice velocities on the falling limb of the diurnal melt input cycle (e.g. 16 -19 June, Fig. 8).  6. Tracer returns from all 11 traces of moulin L41A. Errors were empirically estimated from duplicate samples as 15% for traces 1, 2, 4, 5, 7 and 8, and 11% for traces 3, 6, 9, 10 and 11. This period ended with the extreme melt event of 9-10 July reported by Nghiem et al. (2012), when rapid ablation during warm, humid and cloudy conditions caused 2-3 days of sustained high moulin melt water inputs of 8-10 m 3 s -1 (with no clear night-time minimum) and temporary loss of the supraglacial gauging station. Discharge in the Watson River (of which Leverett is a tributary) Fig. 7. Times series of ice surface ablation, moulin melt water input, tracer velocity, east-west ice surface motion (hourly averages, black line; daily averages, red line) and Leverett river discharge. Ablation errors are 1 s.d. of 5 measurements; errors in melt water input (grey shading) are 95% confidence intervals for the rating curve; 95% confidence intervals for tracer velocities (calculated by Monte-Carlo simulations) are smaller than the symbols; errors in Leverett river discharge (grey shading) follow Tedstone et al. (2013). Fig. 8. Anticlockwise hysteresis between ice velocity and moulin melt water input for three periods: 16-19 June, shortly after the onset of diurnal cyclicity; 15-23 July, during the prolonged cold period after the first July melt event; and 29 July -1 August, during a period of relatively high melt after the second July melt event.
Points are marked at 1 hour intervals. reached a record high and damaged the bridge in Kangerlussuaq on 11 July (Mikkelsen et al., 2016).

Period 3: persistent fast drainage (10 July onwards)
The 9-10 July melt event was followed by ∼13 days of low ablation (<20 mm day -1 ), low but still diurnally cyclic moulin melt water input (average ∼2 m 3 s -1 , dropping to near-zero at night) and slower ice velocity (90-130 m yr -1 ). However, trace #7 on 16 July emerged at a similar velocity to trace #6, and all subsequent traces continued to emerge rapidly with low dispersion (Fig. 6). Another intense melt peak on 26-27 July (ablation reaching 91 mm day -1 ) was accompanied by corresponding peaks in moulin melt water input and surface motion. Interestingly, there was little response in Leverett river discharge.
Anticlockwise hysteresis between moulin melt water input and ice velocity persisted during Period 3, remaining particularly stable during the cold conditions of mid July (Fig. 8), and with greater lag between velocity and melt input than during 16-19 June. Despite similar diurnal cycles of melt water input during 16-19 June and 29 July -1 August, velocity in the latter period (Fig. 8) dropped to much lower levels during the rising limb of the melt water input cycle.

Discussion
The patterns of ice velocity and hydrological changes described above are consistent with evolution from slow/inefficient to fast/efficient drainage, qualitatively similar to many previous studies of ice dynamics in west Greenland and elsewhere (see references in Hubbard et al., 1995;Nienow et al., 2017). Importantly, the similarity between our results and those collected over several seasons on this and other land-and marine-terminating GrIS outlet glaciers (e.g. Sole et al., 2011;Banwell et al., 2013;Meierbachtol et al., 2013;Neckel et al., 2020) shows that our observations are widely representative rather than site-specific. Given this similarity, our discussion here focuses only on the additional insights gained from this detailed hydrological time series.

Surface-bed hydrological connection
Ice acceleration commenced on 31 May, prior to the onset of surface fracturing (3 June), surface uplift (4 June) and Moulin L41 drainage (5 June) (Fig. 3). Therefore, we associate this initial acceleration with longitudinal coupling to downstream ice, which would itself have accelerated as ice-bed hydrological connections became established down-glacier during the period of relatively high melt (∼30 mm day -1 ) in late May. As this zone of new hydrological connections propagated closer to the field site, increasing tensile stress caused the observed surface fracturing on 3 June (Fig. 4). We assume these new fractures initiated the hydrological connection and subsequent moulin development. This upglacier propagation of acceleration, fracturing, and new hydrological connections was simulated by Christoffersen et al. (2018). However, we note that surface uplift at L41A preceded the observed sudden drainage of the relic moulins by about 12 hours, suggesting this uplift was associated with non-local meltwater supply, perhaps from supraglacial lake drainages higher in the catchment. Similarly, Stevens et al. (2015) reported lake hydrofracture triggered by tensile fracturing when meltwater inputs from other moulins flooded the local subglacial drainage system.
Rapid supraglacial lake drainage events release large volumes of water over timescales of a few hours (e.g. Doyle et al., 2013;Tedesco et al., 2013;Stevens et al., 2015). In contrast, moulins such as L41A which open under a supraglacial channel supplied by a lake are hydrologically distinct in their 'slow but steady' release of stored supraglacial water; lake drainage is limited by downwards incision of the supraglacial channel between the lake and moulin. This occurs over timescales of days to weeks, leading to a longer but weaker ice dynamic response (Tedesco et al., 2013). These moulins can open by reactivation of a previous season's moulin (Tedesco et al., 2013); however, seasonal snow remnants in a former supraglacial channel 'downstream' of L41A (Fig. 4) suggest that L41A formed after 3 June when the observed surface fractures intersected a supraglacial channel. This is not certain as we did not access L41A until 15 June. Considering these several different patterns of hydrological connection, moulin development at L41A was likely initiated by the fracturing associated with upglacier propagation of tensile stresses (Christoffersen et al., 2018), while the observed dynamic response likely reflects a combination of local and non-local meltwater inputs. From a hydrological viewpoint, the main distinction between sites is in the fast versus slow drainage modes of moulins opening under lakes and channels, respectively. However, from a modelling viewpoint, the somewhat different ice dynamic triggers and responses between several sites (Doyle et al., 2013;Tedesco et al., 2013;Stevens et al., 2015, and here at L41A) suggest the conditions and processes driving this initial connection need further investigation, since accurately simulating the timing and distribution of moulin development is a vital aspect of subglacial hydrological models (Banwell et al., 2016).

Evolution and persistence of fast/efficient drainage
During the relatively cold conditions and new snow in the week immediately following establishment of moulins on 3-5 June, ablation decreased (average 9 mm day -1 ) and surface velocity decelerated. The surface-bed hydrological connection was presumably maintained primarily by gradual release of stored supraglacial melt water, as supraglacial channel incision released additional areas of ponding. Relatively weak hysteresis during 16-19 June (red cycles, Fig. 8) with sustained high velocities, despite only moderate melt water input, suggests that surface motion was responding rapidly to local melt inputs into a pressurised drainage system. This is consistent with initially inefficient drainage indicated by Trace #1.
The rapid increase in tracer velocity and decreasing dispersion during 16 June to 10 July indicated development of fast/efficient drainage consistent with the limited tracer experiments at L41 in 2011 (Chandler et al., 2013). However, the three-fold increase in v 05 from 0.27 to 0.61 m s -1 in 7 days (16 -23 June) at L41A can be compared with the much slower increase from 0.25 to 0.38 m s -1 over 8 days (26 June -4 July) in 2011 at L41. This difference is likely related to the much lower melt inputs in 2011: supraglacial stream discharge into L41A was ∼4 m 3 s -1 during that period in 2012 compared with ∼2 m 3 s -1 into L41 during the respective period in 2011 (Chandler et al., 2013); the corresponding Leverett catchment discharges averaged 310 and 180 m 3 s -1 . Decreasing tracer transit times from L41A are in good agreement with a decrease in calculated melt routing delay from approximately 4 days (June average) to 1 day (July average) for melt at 1000 m elevation in 2012 (Mikkelsen et al., 2016).
Each tracer return represents an integrated signal from all parts of the system through which it has transited. The consistently short transit times and low dispersion of all traces from #3 onward show that the entire flow path between Moulin L41A and the proglacial river remained efficient, consistent with gas and dye tracing of several other moulins in previous melt seasons (Chandler et al., 2013) and supported by GPS observations of surface velocity (e.g. Sole et al., 2013). Efficient drainage development within 2-3 weeks of surface-bed connection is currently underestimated in hydrological models, including those that do eventually develop efficient drainage under similar conditions (e.g. Werder et al., 2013). The efficient drainage persisted through extended periods of low surface melt and moulin melt water input (e.g. 12-24 July), perhaps aided by gradual release of stored melt water higher in the catchment (3-4 days delay in melt routing from 1500 m elevation in July 2012: van As et al., 2017).
In a channelised drainage system with rapid evacuation of melt water, we would expect local subglacial water pressure and ice velocity (near the moulin) to closely follow melt input, with little lag. In contrast, the anticlockwise hysteresis between melt water input and ice velocity during low melt in mid-July, and again during a period of higher melt (29 July -1 Aug) (Fig. 8), reveals higher ice velocity on the falling limb of the local moulin water input cycle. Peak ice velocity lags peak moulin input by ∼6 hr, compared with a time scale of only 12 to 15 hr for tracers to transit the entire >40 km flow path. The hysteresis suggests that local subglacial water pressure at L41A is driven partly by a non-local source, consistent with an efficient drainage system receiving melt water inputs from upstream regions (L41A is close to the predicted subglacial pathway draining a considerable portion of the upper catchment: Fig. 1). This non-local up-glacier melt water supply will have maintained high local subglacial water pressure beneath the moulin (and, therefore, higher ice velocity) for several hours after the peak of local water input, due to the travel times required to deliver melt water from more distant upstream sources.
The potentially complex relationship between local ice velocity and local melt input, due to non-local subglacial hydrology, supports Covington et al. (2012) who note caution when using proglacial hydrographs to infer drainage system evolution. This clearly highlights the benefit of our simultaneous time-series of moulin input, water transit and ice dynamics.

Implications for subglacial hydrological models
Subglacial hydrology models have invoked a variety of flowpath geometries for melt water, including sheets, linked cavities, channels and porous till (de Fleurian et al., 2018). Rapid tracer returns with low dispersion are not consistent with flow through a thin sheet, linked cavities or porous till, as these would be strongly dispersive; therefore, our tracer experiments provide strong evidence for the development of channels at the ice-bed interface. Such channels could be incised upwards into the ice (R-channels: Röthlisberger, 1972), downwards into the substrate ('Nye channels' in bedrock, or 'canals' in sediments), or of course both. Our data are supportive of those models in the SHMIP experiment that develop channels in the idealised ice sheet geometry (de Fleurian et al., 2018).
Neither the tracing experiments nor the other data sets can at present identify whether subglacial conduits are incised upwards or downwards; however, future comparison of this data set with drainage system models may enable such distinction if the hydraulic parameters are found to be sufficiently different (e.g., if dispersion in rough-walled 'canals' is much greater than in smooth-walled R-channels). Few in-situ measurements of conduit roughness characteristics are available, and are limited to accessible conduits under relatively thin ice (e.g., Gulley et al., 2012). Importantly, the rate of conduit development may be very different for incision into sediments, perhaps faster if sediment erosion requires less energy per unit volume than ice melt. When channelised drainage persists for several weeks in regions where the orientation of the conduit is parallel to that of the basal sliding direction, a conduit incised into the ice would likely erode a conduit in any underlying sediments (or vice-versa); this is less likely if the channel orientation and sliding direction are not aligned, as a conduit in the ice would be moving laterally relative to the bed. Downwards incision of channels into sediments may also be limited by the increasingly coarse substrate, as only the finer material can be removed and the channel floor becomes 'boulder armoured' (Gulley et al., 2013), though lateral development of canals could be unhindered. These boulders could then reduce the rate of channel closure by ice creep during periods of lower subglacial water pressure. Downwards incision of conduits aligned parallel to sliding direction could nevertheless be one means by which the drainage system can remain efficient through periods of low melt input, if these conduits take longer to close. In reality, we suspect that multiple drainage forms may exist simultaneously.

Conclusions
We find clear evidence for the development of efficient subglacial drainage extending from the margin to beneath ice at least ∼900 m thick, in a region of the ablation zone where the extent of seasonal drainage evolution has remained uncertain in previous field studies and numerical simulations. We have also demonstrated rapid development of efficient drainage (by late June, following estimated surface-bed connection on 3 June) from a moulin initiated by a fracture intersecting a supraglacial channel, rather than by lake hydrofracture. The efficient drainage system persisted and showed little response to very variable moulin inputs includ-ing extreme melt events and a 2-week period of subdued melt (Figs. 6, 7). Given the many similarities in seasonal evolution observed in this catchment and elsewhere in Greenland -including in tidewater glaciers upstream of the influence of the calving front Nienow et al., 2017) -we consider our data to be generally representative of similar sectors of the GrIS rather than site-specific.
Despite considerable advances in the numerical and physical complexity of subglacial drainage models, their development and validation remain limited by a lack of observations. This comprehensive data set will now provide a unique reference against which to test model configurations and parameters (for example, in future SHMIP phases: de Fleurian et al., 2018). The data are particularly relevant as they represent a region of the ablation zone where conflicting results in previous work show that calculated subglacial drainage characteristics are sensitive to the choice of model configuration.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.