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

(1)

Contents lists available atScienceDirect

Earth and Planetary Science Letters

www.elsevier.com/locate/epsl

Rapid development and persistence of eﬃcient subglacial drainage under 900 m-thick ice in Greenland

David M. Chandler

^a^,^b^,^∗

, Jemma L. Wadham

^b^,^c

, Peter W. Nienow

^d

, Samuel H. Doyle

^e

, Andrew J. Tedstone

^f

, Jon Telling

^g

, Jonathan Hawkings

^h^,ⁱ

, Jonathan D. Alcock

^b

, Benjamin Linhoff

^j

, Alun Hubbard

^c^,^k

aNORCENorwegianResearchCentreandBjerknesCentreforClimateResearch,Bergen,Norway bBristolGlaciologyCentre,UniversityofBristol,Bristol,UK

cCAGE- CentreforArcticGasHydrate,EnvironmentandClimate,DepartmentofGeosciences,UiTtheArcticUniversityofNorway,9010,Tromsø,Norway dDepartmentofGeosciences,UniversityofEdinburgh,Edinburgh,UK

eCentreforGlaciology,DepartmentofGeographyandEarthSciences,AberystwythUniversity,Aberystwyth,UK fDepartmentofGeosciences,UniversityofFribourg,Fribourg,Switzerland

gSchoolofNaturalandEnvironmentalSciences,NewcastleUniversity,Newcastle,UK

hNationalHighMagneticFieldLaboratoryGeochemistryGroupandtheDepartmentofEarth,OceanandAtmosphericSciences,FloridaStateUniversity,USA iGermanResearchCentreforGeosciencesGFZ,Potsdam,Germany

jUSGeologicalSurvey,NewMexicoWaterScienceCenter,Albuquerque,NM,USA kKvantumInstitute,UniversityofOulu,90014Oulu,Finland

a rt i c l e i n f o a b s t r a c t

Articlehistory:

Received21July2020

Receivedinrevisedform2April2021 Accepted28April2021

Availableonlinexxxx Editor:J.-P.Avouac

Keywords:

GreenlandIceSheet hydrology subglacialdrainage icefracture moulin

Intensive study of the Greenland Ice Sheet’s (GrIS) subglacial drainage has been motivated by its importanceforicedynamicsandfornutrient/sedimentexport tocoastalecosystems.Thishasrevealed consistentseasonaldevelopmentofefficientsubglacialdrainageinthelowerablationarea.Whilesome hydrological modelsshow qualitative agreement with fielddata, conflicting evidence (bothfield- and model-based) maintains uncertainty in the extent and rate ofefficient drainage development under thick(∼¹ ^km)îce.^Here,^we^present ^the^first simultaneoustimeseriesofdirectly-observedsubglacial drainageevolution,supraglacial hydrologyand icedynamicsover11 weeksinalargeGrIScatchment.

Wedemonstratedevelopmentofafast/efficientsubglacial drainagesystemextending fromthemargin tobeneathice>900mthick,whichthenpersistedwithlittleresponsetohighlyvariablemoulininputs includingextrememelteventsandextendedperiods(2weeks)oflowmeltinput.Thisefficientsystem evolvedwithin∼³^weeksâtâ^moulinînitiated^whenâ^fractureintersectedasupraglacialriver(rather thanhydrofractureand lake drainage). Ice flowresponse tosurfacemeltinputs atthissitefollows a patterncommonly observed inthe lower GrIS ablation area,and by assuming a strong relationship betweenicedynamicsandsubglacial hydrology,weinferthatefficient subglacialdrainageevolutionis widespreadunder900m-thick iceinwest Greenland.Thistimeseriesoftracer transitcharacteristics througha developing and then persistent efficient drainage system providesa unique data set with whichtovalidateand constrainexistingnumericaldrainagesystemmodels,extending theircapability forsimulatingdrainagesystemevolutionundercurrentandfutureconditions.

©²⁰²¹^TheÂuthors.^Published^byÊlsevier^B.V.^Thisîsânôpenâccessârticleûnder^the^CC^BY^license (http://creativecommons.org/licenses/by/4.0/).

1. 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;Sundaletal.,2011),andmorerecentlytoquantifynutrient

*

Correspondingauthor.

E-mailaddress:[email protected](D.M. Chandler).

andsediment export fromthe ice sheet to coastal/fjordenviron- ments (e.g. Sejret al., 2014; Hawkings etal., 2015; Hopwood et al., 2020). Some aspects of the GrIS subglacial drainage system arenowunderstoodinconsiderabledetail,withfield andremote- sensing data revealing a persistent, widespread seasonal pattern intheablationarea,characterisedbyinland-propagationofevolu- tionfromslow/inefficienttofast/efficientdrainage(Nienowetal., 2017). Similar seasonal behaviour is observedon manyother ice massesatsmallerscales(e.g.Willis,1995;Binghametal.,2006).

https://doi.org/10.1016/j.epsl.2021.116982

0012-821X/©²⁰²¹^TheÂuthors.^Published^byÊlsevier^B.V.^Thisîsânôpenâccessârticleûnder^the^CC^BY^license(http://creativecommons.org/licenses/by/4.0/).

(2)

In Greenland, the inland propagation of eﬃcient drainage is linked to the establishment of hydrological connections between the surface and subglacial drainage pathways under increasingly thick ice, enabled by fracturing andsubsequent moulin development,sometimesassociatedwithsupraglaciallakedrainageevents (Nienow et al., 2017; Hoffman etal., 2018). Transient periods of extensionallongitudinalstressesenablefracturingeveninicenor- mally undercompressive longitudinalstress (Christoffersenetal., 2018). The inland extent of eﬃcient subglacial drainage is cur- rently unclear (Bartholomew et al., 2011; Chandler et al., 2013;

Meierbachtol etal., 2013; Hoffman etal., 2016; Christoffersen et al., 2018), but theory suggests this is limited by increasing ice thickness, decreasing hydraulic potential gradient,and decreasing ablationratefurtherinland(e.g.Röthlisberger,1972;Schoof,2010;

Dow etal., 2014). Fast,eﬃcientsummerdrainagereverts toslow, ineﬃcientdrainageduringwinterexceptveryneartheicemargin, asicecreepgradually closessubglacialconduitsoncemelt inputs haveceased.

OwingtotheinaccessibilityoftheGrISsubglacialdrainagesys- tem, muchofour currenthydrological understandingisbased on indirect observations. For example, borehole and moulin water pressure records,and corresponding diurnal to seasonal changes in ice surface motion, can reﬂect drainage system evolution via relationships betweendrainagesystemeﬃciency,subglacialwater pressure,basal tractionandicemotion(see referencesinNienow et al., 2017). The observed response to melt inputsin Greenland is very similar tothat in smallerice masses(e.g.Hubbard etal., 1995). However, the impact of meltwater-driven ice acceleration hasan uncertainimpacton netannualicedisplacement,particularly inthe upper ablationarea ofthe GrIS (Sundal et al., 2011;

Soleetal.,2013;vandeWaletal.,2015;Doyleetal.,2014).

Predictionsofhowinteractions betweenicedynamicsandhy- drology willimpactGrIS massbalance andnutrient/sedimentexport,especiallyunderglobalwarmingscenarios,requirenumerical subglacialdrainagemodels.Suchmodelscanqualitativelyreplicate theseasonalcharacteristicsdescribedabove(e.g.Schoof,2010;He- witt, 2013; Werder et al., 2013; Banwell et al., 2016; Hoffman et al., 2016; de Fleurian et al., 2018; Koziol and Arnold, 2018), butquantitativeevaluationofthesemodelsisseverelylimitedby a lack of relevant ﬁeld data withwhich to compare simulations.

At present,adetailedtime seriesofsubglacialdrainageevolution based onobservations ofwater transitthrough the drainagesys- temisnotavailablefromanypartoftheGrIS.

Here,wepresentthefirstsimultaneoustimeseriesofdrainage system tracing, supraglacial hydrology and ice dynamics during the meltseason ina largeGrIS catchment (LeverettGlacier, west Greenland). Weundertooka seriesoftracesata moulinin∼⁹⁰⁰ m thick ice (Lindbäck et al., 2014), located 41 km from where its draining waters emerge at the GrIS margin (Fig. 1). We also collected simultaneousobservationsoficesurface motion,surface ablation,andmoulinmelt waterinput.Finally,dischargefromthe catchmentwasmeasuredwheretheproglacialriveremergesatthe icemargin.Themoulin chosenforthisexperimentisideallysitu- atedinaregionoftheGrISwherethereiscurrentlydisagreement regarding theseasonalextent ofefficientdrainage evolution(He- witt,2013;Werder etal., 2013; Banwelletal., 2016;Hoffman et al., 2016; de Fleurian etal., 2018; Kozioland Arnold, 2018). Our comprehensivetimeseriesofhydrologyandicedynamicsprovides valuable insights into how an efficient drainagesystem develops andpersistsunderthickice,andwillprovideanidealdatasetwith whichtovalidatesubglacialdrainagemodels.

2. Fieldsite

LeverettGlacierisa land-terminatingcatchmentoftheGrISat approximately67^◦N inwestGreenland, witha theoretical hydro-

Fig. 1.(a)MapoftheﬁeldsiteshowinglocationofMoulinL41A(yellowcircle);

watersamplingsitenearLeverettGlacierterminus(yellowtriangle),bedtopography (shading)andsurfacetopography(bluedashes,at200mintervals)fromLindbäck etal.(2014);andestimatedhydrologicalcatchment(white) andsubglacialwater routing(green)followingLindbäcketal.(2015) withk=1(pressurepotentialequal toiceoverburdenpressure)inEq.2oftheirsupplementarymaterial.Theestimated meltwaterrouting betweenthemoulin andterminusis markedinorange.The insetshowslocationsofthenewMoulinL41Atracedinthisstudy,inactiveMoulin L41tracedinthe previousyearbyChandleret al.(2013),and locationsofGPS receivers.(b)Icesurfaceandbedtopographyalongtheestimatedmeltwaterroute fromthemoulintoterminusasshownin(a).(Forinterpretationofthecoloursin theﬁgure(s),thereaderisreferredtothewebversionofthisarticle.)

logicalcatchmentestimatedas∼⁹⁰⁰^to¹²⁰⁰^km² ^(Palmer^et^al., 2011; Lindbäck etal., 2015). The hydrology andice dynamicsin thiscatchment have beenstudied forover a decade (e.g.van de Walet al., 2008; Bartholomew etal., 2011; Palmer et al., 2011;

Sundal et al., 2011; Chandler et al., 2013; van As et al., 2017) andfollow the seasonal patterns common to other major catch- mentsinGreenland,includingtheneighbouringIssunguataSermia (Meierbachtol et al., 2013), the Paakitsoq region (Banwell et al., 2013; Andrews et al., 2014), south Greenland (Sole et al., 2011) andnortheast Greenland(Neckel et al., 2020). Therefore,the hy- drologyandicedynamicsatLeverettGlaciercanbeconsideredto bewidelyrepresentativeofmanyofGreenland’soutletglaciers,as well as many smaller ice masses in the Arctic (e.g. Bingham et al., 2006). Data were collected during the2012 melt season ata moulin(L41A)approximately700msouthofMoulinL41tracedin 2011(Chandleretal.,2013)(Fig.1a).MoulinL41reopenedin2012 butwasinactive asthemainsupraglacialriverdraining thatpart oftheicesheetwascapturedbyL41A.IcethicknessatL41Aises- timatedfromradio-echosounding(Lindbäcketal.,2014)as800m, withanuncertaintyof18m,andincreasesto910mwithin2km ofthe moulin along the estimatedsubglacialﬂow path (Fig. 1b).

Theicethicknessisbetween800and920mfortheﬁrst7kmof theestimatedﬂowpath.

3. Methods

3.1. Tracerexperiments

We used sulphur hexaﬂuoride (SF6) as the tracer, which has beenappliedinavarietyofenvironments(includinginocean-scale experiments)withhighdilution(Watsonetal., 1987; Clarketal., 1996; Chandler etal., 2013). This tracer can be detected at very low concentrations,doesnot adsorb ontosediments,andis non- reactive;therefore,itisadvectedatthesamerateasthewaterin water-ﬁlleddrainage systemsas previously validatedby compar-

(3)

ison withdye tracing in theLeverett catchment (Chandler etal., 2013). We note the increasing environmental concerns regarding SF6,andwhiletheamountrequiredforeachtracerexperiment(3 kg)wasverysmallinthecontextof∼⁹^gG^global^emissions^(Engel etal.,2018),alternativesshouldbesoughtifsimilartracerexperi- mentsareconductedinthefuture.

A total of11 traces were carried out between16 June and 6 August, 2012. Tracers were injected at ∼^19:00 ^hrs, ^shortly ^after the peakdaily moulinmelt waterinput (typicallybetween15:00 and19:00hrs).Tracerconcentrationsweremonitoredbycollecting water samplesfrom theproglacial river draining LeverettGlacier andwereanalysedbygaschromatography.

The solubilityofSF₆ inwateris low(74 mg/l at0^◦C)(Bullis- teretal., 2002), andthechoice ofinjectionmethod iscriticalto minimise gas loss at the injectionsite. Following trials in previous years, the method adopted in 2012involved direct injection through a150m, 5mm internaldiameterhoseloweredinto the moulin. The hose was connected to the gas cylinder via a pre- calibratedflow meter andregulator. Bymonitoringthe flowrate, aknownamountofgaswasinjected.Thebottomofthehosewas leftopen, andthepresenceofresidualpressureinthehosewhen disconnectingtheapparatusfollowingtheinjectionconfirmedthat theinjectionwasintowaterratherthanintoair.

Proglacial river samples were collected in 120 ml gas-tight Wheatonvials(induplicateortriplicate)sealedwithrubberstop- pers and a crimp top. Care was taken to ensure no gasbubbles remainedinthebottleonceit hadbeenstoppered.N2 headspace wasaddedtothebottlesintheﬁeldbypiercingthestopperwith a long needlepushed through to the bottomof thebottle; then, a 60 ml gas-tightsyringe ﬁlled withN2 atambient temperature andpressurewasusedtoadd20mlofN2 throughashortneedle just protruding through thestopper. Displacedwaterwas ejected throughthelongneedle.Bottleswereshakenandleftfor>2hours toallowheadspaceequilibration.

The SF6 concentration in the headspace was analysed in the LOWTEXlaboratories,BristolUniversity,UK,usingaportableCam- bridge Scientific 300-series gas chromatograph (GC) fitted with an electron capturedetector(ECD). The headspace inthe 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 thantransfer via a gastight syringe achievedgreater sen- sitivitythan ourprevious method(Chandler etal.,2013). TheGC used a carriergasof CP- orECD-grade N2 ata flow rateof∼¹⁸ ml/minute. Waterand oxygen traps were installed in the carrier gas supply. Theinjector, columnanddetector temperatureswere 150^◦C, 100^◦C and 200^◦C, respectively. Two sample loops were used(1mlor∼^0.05 ^ml),^with^the^smaller^loop^being^fitted^when peakSF6 concentrationwasexpectedtosaturatethedetectorwith the1mlloop.Sampleswere analysedinduplicate,anda1 parts per billionbyvolume(ppbv) standardwasusedtocorrectforin- strumental drift every 5 samples. All samples for an individual tracewereanalysedwiththesamesampleloop.Specifically,traces 1,2,4,5,7,8 wereanalysedwiththesmallloop andtraces3,6, 9, 10, 11 withthe large loop. Bottles from the same trace were analysedtogetherasabatchbutinarandomorder.Peakareawas convertedtoheadspaceSF6concentrationbycalibrationwiththe1 ppbvstandard.Calibrationwithjustonestandardassumedalinear ECDresponseovertherangeofinterest,asconfirmedviadilution experiments(Fig.2).Calculationofthecorresponding SF6 concentration in the river water sample prior to adding the headspace followedChandleretal.(2013).

Errors inthetracerconcentrationsestimatedempiricallyusing duplicate bottles yielded mean percentage errors of 15% for the small sample loop and 11% for the large sample loop. The limit

Fig. 2.Electron-capturedetector(ECD)linearityandlevelofdetectionforGCanaly- siswiththetwosampleloops.EachSF6concentrationwasobtainedbysuccessive dilutionofthe1ppbstandard,andanalysedintriplicate.

ofdetection(LOD)oftheinstrument(expressedasheadspacecon- centration)was0.8partspertrillionbyvolume(pptv)forthelarge sampleloopand3pptvforthesmallsampleloop,equivalenttore- spectiveSF₆ concentrationsof1×¹⁰^–3^parts^per^trillion^by^mass (pptm)and5×¹⁰^–3^pptmⁱⁿ^the^water^samples.

Akeycharacteristicofthetracerreturncurvesisthetracerve- locity, which quantifies the efficiency of the subglacial drainage system. Several measures of velocity are available including the maximumvelocityandvariousmeasures ofmeanvelocity(Chan- dleretal.,2013).Herewereportboththemaximumvelocityand meanvelocityforall traces,calculatedusingthestraight-linedis- tancebetweenmoulinandsamplingsitedividedbythetimewhen 5% or50% ofthe recovered tracer hasemerged (v05 andv50,respectively). Due to the assumption of ‘straight-line distance’, the estimatedflowvelocitiesareminimaassinuosityintheflowpath willcausethecalculatedvelocitytobelowerthantheactualflow velocity. Errors in v05 and v50 were estimated in a Monte-Carlo simulation(withn =¹⁰⁰⁰⁾ûsing^theempirically-derived analyti- calerrors.

3.2.Icesurfaceablation

Surface ablation rates were measured daily using changes in surfaceheightatfiveablationstakes,arrangedinacrossconfigura- tionat∼²^mseparation.Thestakeswerelocated∼⁷⁰⁰^m^northôf MoulinL41A,andwere inthesupraglacialhydrologicalcatchment feedingthatmoulin.Thestakeswereinstalledinholesdeeperthan the length of the stake(so each measurement was from the ice surfacedowntothetopofthestake),toavoidtheproblemofen- hancedsurfacemeltingcausedby solarradiationabsorbedbythe stake(Chandleretal.,2015).Allablationmeasurementswerecar- riedoutby thesameobservertoensureconsistency, forexample intheinterpretationofthelevelofaroughicesurface.

3.3.Icesurfacemotion

Icesurface motion was recorded by ﬁve dual-frequencyLeica SR520 GPS receivers deployed on poles drilled 2 m into the ice surface(Fig.1).GPSdatawerepost-processedkinematically(King, 2004) with Track v.1.27 software (Chen, 1998) against bedrock- mounted reference stations using a precise ephemeris from the InternationalGNSS Service (Dowet al., 2009). Reference stations were located 1 km from the terminus of Russell Glacier and at Kellyville,givingbaselinelengthslessthan41km.Duetogapsin thetimeseriescausedbypoweroutage,weaveragedthehorizon-

(4)

Fig. 3.Ice surface motion at the establishment of a hydrological connection between the supraglacial and subglacial drainage systems (‘spring event’).

talvelocitiesrecordedattheﬁvestationswiththefewestgapsto giveasinglerecord.Positionswererecordedat30sintervals;1-hr means were thensmoothed usinga 5-pointbinomialﬁlter. Since therewasgenerallylittledifferenceinvelocitybetweenthestakes (Fig.3),themeanvelocityacrossthenetworkgivesabetterindica- tionoftheseasonalpatternoficemotionwithfewergapsthanin theindividualrecords.Velocitiesare centreddifferencesofhourly displacements,i.e.thevelocityattimeTiis(Xi+1- Xi-1)/(Ti+1- Ti-1).

GPS stakesrequired periodic re-drilling asthey gradually melted out. Since re-drillingsometimes coincided withgaps inthe data, afullseasonofverticaldisplacementsisnotavailable;thereliable upliftrecordisrestrictedtotheearlypartoftheseason.

3.4. Supraglacialandproglacialriverdischarge

Supraglacial river discharge was measured by monitoring water depth with a custom-built pressure sensor, at 1-minute in- tervals, and converting depth to discharge with a rating curve established fromsaltdilutiongauging.Thismethodiswell suited to supraglacial river gauging as the electrical conductivity of supraglacial melt wateris very low.The pressure sensor was in- stalledonthebottomofthechannelandwasabletomovedown- wards withthe ice surface asthe channel gradually incised. For furtherdetails,seeWadhametal.(2016).

Proglacialriverdischarge was monitoredusingstage measurements (collected by a HOBO pressure sensor) and dye dilution gauging atastablebedrocksection neartheterminusofLeverett Glacier, following Bartholomewet al. (2011) and Tedstone et al.

(2013).

4. Results

The 11-week ﬁeld campaignon the ice sheet commenced on 28May2012,priortoanyvisibleevidenceofmeltwaterponding.

The surface atthistime was characterisedby seasonalsnow and patchesofbareice(Chandleretal.,2015).Snowmeltoverthesub- sequentweekcreatedextensiveareasofslushandeventuallymelt waterpondinginsurfacedepressionsandrelictmoulins,including MoulinL41tracedin2011.Alakedevelopedapproximately200m up-glacierfromL41,withanestimatedmaximumlateralextentof

∼⁵⁰⁰^mândânûnknown^depth.^This^lake^was>200mup-glacier fromthenewMoulinL41A,buttheexactextentisunknownasit wasnotaccessibleonfoot.

Fig. 4.Left:linearfracturesopenedon3Juneextendingseveralhundredmetres N-S(transversetoﬂow),andintersectedexistingchannelsanddepressions.This photowastakenbetweentheactivemoulinL41AandinactivemoulinL41.Drillfor scale.Right:MoulinL41Aon19June,showingdeepsnowﬁllingtheformerchannel downstreamofthemoulin.

On31May,the westwardsicesurface velocity increasedfrom 100-120m yr^-1to140-190m yr^-1(Fig.3).Aperiodoffrequent,au- diblecracking andthedevelopment ofextensivesurface fractures (generallyorientedN-Sandreaching∼²^cm^width)^commenced^on 3June(Fig.4).Subsequentinspectionofthisfracturezoneshowed that it extended from the vicinity of Moulin L41, southwards to MoulinL41A(Fig.1a),andcontinuedforaleast1kmfurthersouth ofMoulinL41A.

Weconsider3Juneasmarkingtheonsetofthehydrologically activepartoftheseasonatL41A,whichwedivideintothreedis- tinctperiodsasdescribedbelow.

4.1. Period1:surface-bedconnection(3-5June)

DuringPeriod1,surface melt waterremainedpondedandthe audiblecrackingcontinued.Thewestwards icemotionaccelerated furtherfrom140-170m yr^-1 to230m yr^-1by 23:00on5June.At thistime, the waterthat hadpondedin relic MoulinL41andin othersmallerrelicmoulinsdrainedabruptly(Fig.3).Drainagewas accompaniedby∼^0.3^môf^verticalûplift.^Both ^verticalûpliftând westwardsicevelocitypeaked12-24hraftermoulinL41drainage.

4.2.Period2:evolutionfromslowtofastdrainage(5June- 10July) Moulin L41 received very little melt water input despite its proximityto extensiveponding, asthesupraglacial channel feed-

(5)

Fig. 5.Top:Photolookingeast(up-glacier)on19Junefromnearthefracturesin Fig.4.Thesurfacedepressionwithasmalllakegraduallydrainedviashortrivers intomoulinL41in2011andmoulinL41Ain2012.Fromthisviewpoint,inactive moulinL41istothe left(north)and activemoulinL41Aistothesouth(right).

Bottom:thenextlargemoulin ∼¹^km^south^of^L41A^situatedⁱⁿ^the^same^linear fracturezoneasL41/L41Aandmoreclearlydisplayingthefracturesintersectingthe previousseason’smeltchannel.

ing L41was blockedwithsnow.Instead, thepondingmelt water drained southwest along another supraglacial channel. We were unable toaccessthenewmoulin(L41A)fedby thischannel until 15June,assurfaceconditionspreventedaccessonfoot.Therefore, moulin L41A was ﬁrst traced on 16 June, likely 11-14 days after ice-bed hydrological connection. L41, L41Aand several other smaller moulins were all located in the linearfracture zone de- scribed above, and at least 200 m down-glacier from the lake (Figs. 1, 4, 5). There was no rapid lake drainage; instead, the lakegraduallydecreasedinextentoverthefollowingmonth,presumably due to ongoing supraglacial channel incision. Trace #1 emerged very slowly over 3-4days, with respective v05 and v50 velocitiesof0.27and0.19m s^-1;subsequenttraces(untilTrace#6 on 10 July) emerged progressively quicker and with lessdisper- sion, withv05 andv50 reaching 1.27and 0.94m s^-1, respectively (Figs.6,7).

Melt water input into moulin L41A showed strong diurnal cyclicity, with morning minima <3 m³s^-1 and evening peaks reaching 5 - 12 m³s^-1. An increase in melt input (most signif- icantly, an increase in the diurnal minima) was associated with warm,cloudyconditionsandcorrespondingrapidablationexceed- ing50mm day^-1(Fig.7).

Westwards ice velocity peaked at ∼²⁶⁰ ^{m yr}^-1 ^on ⁶ ^June, shortlyafterL41drained,thendeclinedto∼¹³⁰^{m yr}^-1^by^9th^June during a period oflow ablation andnew snow.Diurnal cyclicity commencedabruptlyon16June(Fig.7).Surfacemotiongenerally decreasedfrom16June - 10July, exceptfortwo relativelystrong diurnalpeakson20Juneand4July,coincidentwithpeaksinab- lation and moulin melt input (Fig. 7). The relationship between moulinmeltwaterinputandicevelocityshowedstronganticlock- wisehysteresis withfastericevelocitiesonthefallinglimbofthe diurnalmeltinputcycle(e.g.16- 19June,Fig.8).

Fig. 6.Tracerreturnsfromall11tracesofmoulinL41A. Errorswereempirically estimatedfromduplicatesamplesas15%fortraces1,2,4,5,7and8,and11%for traces3,6,9,10and11.

Thisperiodendedwiththeextrememelteventof9-10Julyre- portedbyNghiemetal.(2012),whenrapidablationduringwarm, humid and cloudy conditions caused 2-3 days of sustainedhigh moulinmeltwaterinputsof8-10m³s^-1(withnoclearnight-time minimum)andtemporarylossofthesupraglacialgaugingstation.

Discharge in the Watson River (of which Leverett is a tributary)

(6)

Fig. 7.Timesseriesoficesurfaceablation,moulin meltwaterinput,tracervelocity, east-westicesurfacemotion(hourlyaverages,blackline;dailyaverages,redline) andLeverettriverdischarge.Ablationerrorsare1s.d. of5measurements;errorsinmeltwaterinput(greyshading)are95%conﬁdenceintervalsfortheratingcurve;

95%conﬁdenceintervalsfortracervelocities(calculatedbyMonte-Carlosimulations)aresmallerthanthesymbols;errorsinLeverettriverdischarge(greyshading)follow Tedstoneetal.(2013).

Fig. 8.Anticlockwisehysteresisbetweenicevelocityandmoulinmeltwaterinput forthreeperiods:16-19June,shortlyaftertheonsetofdiurnalcyclicity;15-23July, duringtheprolongedcoldperiodafterthe ﬁrstJulymeltevent;and 29July- 1 August,duringaperiodofrelativelyhighmelt afterthesecondJulymeltevent.

Pointsaremarkedat1hourintervals.

reacheda record highand damaged thebridge in Kangerlussuaq on11July(Mikkelsenetal.,2016).

4.3.Period3:persistentfastdrainage(10Julyonwards)

The9-10Julymelteventwasfollowedby∼¹³^daysôf^lowâb- lation(<20 mm day^-1), low butstill diurnally cyclicmoulin melt water input (average ∼² ^m³^s^-1^, ^dropping ^to ^near-zero ât ^night) andslowericevelocity(90-130m yr^-1).However, trace #7on16 Julyemergedatasimilarvelocity totrace#6,andallsubsequent 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 meltwaterinputandsurfacemotion.Interestingly,therewaslittle responseinLeverettriverdischarge.

Anticlockwisehysteresisbetweenmoulinmeltwaterinput and icevelocitypersistedduringPeriod3,remainingparticularlystable duringthecoldconditionsofmidJuly(Fig.8),andwithgreaterlag betweenvelocity andmeltinput thanduring 16-19June. Despite similardiurnal cyclesofmelt waterinput during 16-19June and 29July- 1August,velocityinthelatterperiod(Fig.8)droppedto muchlower levelsduring therisinglimbofthe meltwaterinput cycle.

(7)

5. Discussion

Thepatternsoficevelocityandhydrologicalchangesdescribed aboveareconsistentwithevolutionfromslow/inefficienttofast/efficient drainage, qualitatively similar to manyprevious studies of icedynamicsinwest Greenlandandelsewhere(see referencesin Hubbardetal., 1995;Nienowetal.,2017). Importantly,thesimi- laritybetweenourresultsandthosecollectedoverseveralseasons onthisandotherland- andmarine-terminatingGrISoutletglaciers (e.g. Sole et al., 2011; Banwell et al., 2013; Meierbachtol et al., 2013;Neckeletal., 2020) showsthat ourobservationsare widely representative rather than site-specific. Given this similarity, our discussionherefocusesonlyontheadditionalinsightsgainedfrom thisdetailedhydrologicaltimeseries.

5.1. Surface-bedhydrologicalconnection

Ice accelerationcommenced on 31 May,prior tothe onset of surface fracturing(3June),surface uplift(4June)andMoulinL41 drainage (5 June) (Fig. 3). Therefore, we associate this initial acceleration with longitudinal coupling to downstream ice, which would itselfhave acceleratedasice-bed hydrologicalconnections became established down-glacier during the period of relatively highmelt(∼³⁰^{mm day}^-1⁾ⁱⁿ^late^May.Âs^this^zoneôf^new^hydrological connections propagated closerto the field site, increasing tensile stress caused the observed surface fracturing on 3 June (Fig. 4). We assume these new fracturesinitiated the hydrologi- calconnectionandsubsequentmoulindevelopment.Thisupglacier propagationofacceleration,fracturing, andnewhydrologicalcon- nections was simulated by Christoffersen et al. (2018). However, wenotethatsurfaceupliftatL41Aprecededtheobservedsudden drainage ofthe relic moulins by about 12hours, suggestingthis uplift was associated with non-local meltwater supply, perhaps fromsupraglaciallakedrainageshigherinthecatchment.Similarly, Stevensetal.(2015) reportedlakehydrofracturetriggeredbyten- silefracturingwhenmeltwaterinputsfromother moulinsflooded thelocalsubglacialdrainagesystem.

Rapid supraglacial lakedrainage events release large volumes of water over timescalesof a few hours (e.g. Doyle et al., 2013;

Tedescoetal.,2013;Stevensetal.,2015).Incontrast,moulinssuch as L41A which open under a supraglacial channel supplied by a lake are hydrologically distinct in their ‘slow butsteady’ release of stored supraglacial water; lake drainage is limited by downwards incision ofthe supraglacial channel between the lakeand moulin. This occursover timescales ofdays to weeks,leading to a longer butweakerice dynamic response(Tedesco etal., 2013).

These moulins can open by reactivation of a previous season’s moulin(Tedescoetal.,2013);however,seasonalsnowremnantsin aformersupraglacialchannel‘downstream’ofL41A(Fig.4)suggest thatL41Aformedafter3Junewhentheobservedsurfacefractures intersectedasupraglacialchannel.Thisisnotcertainaswedidnot accessL41Auntil15June.Consideringtheseseveraldifferentpat- ternsofhydrologicalconnection,moulindevelopmentatL41Awas likely initiated by the fracturing associated with upglacier propagation of tensilestresses (Christoffersen et al., 2018), while the observed dynamic response likelyreﬂects a combinationof local andnon-localmeltwaterinputs.Fromahydrologicalviewpoint,the maindistinctionbetweensitesisinthefastversusslowdrainage modesofmoulinsopeningunderlakesandchannels, respectively.

However, fromamodellingviewpoint,thesomewhatdifferentice dynamictriggersandresponsesbetweenseveralsites(Doyleetal., 2013;Tedescoetal.,2013;Stevensetal.,2015,andhereatL41A) suggesttheconditionsandprocessesdrivingthisinitialconnection need furtherinvestigation, since accuratelysimulating thetiming and distributionof moulin development isa vital aspect ofsub- glacialhydrologicalmodels(Banwelletal.,2016).

5.2.Evolutionandpersistenceoffast/eﬃcientdrainage

Duringtherelativelycoldconditionsandnewsnowintheweek immediatelyfollowingestablishmentof moulinson3-5June, ab- lationdecreased(average9 mm day^-1) andsurface velocitydecel- erated. The surface-bed hydrological connection was presumably maintainedprimarilybygradualreleaseofstoredsupraglacialmelt water, as supraglacial channel incision released additional areas ofponding.Relativelyweakhysteresis during 16-19June (redcy- cles,Fig.8) withsustainedhighvelocities,despiteonlymoderate melt water input, suggests that surface motion was responding rapidly to local melt inputs into a pressurised drainage system.

This is consistent with initially ineﬃcient drainage indicated by Trace#1.

Therapidincreaseintracer velocityanddecreasingdispersion during 16 June to 10 Julyindicated development offast/eﬃcient drainageconsistentwiththelimitedtracerexperimentsatL41in 2011(Chandler et al., 2013). However, the three-fold increase in v05 from0.27to 0.61m s^-1 in7 days(16- 23 June)atL41Acan 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 dif- ference is likely relatedto the much lower melt inputs in2011:

supraglacialstreamdischargeintoL41Awas∼⁴^m³^s^-1^during^that periodin2012comparedwith∼² ^m³^s^-1 ^into^L41^during ^the^re- spectiveperiodin2011(Chandleretal.,2013);thecorresponding Leverettcatchment dischargesaveraged 310 and 180m³s^-1. De- creasingtracertransittimesfromL41Aareingoodagreementwith adecreaseincalculatedmelt routing delayfromapproximately4 days(Juneaverage)to1day(Julyaverage)formeltat1000mel- evationin2012(Mikkelsenetal.,2016).

Eachtracerreturnrepresentsanintegratedsignalfromallparts of the system through which it has transited. The consistently shorttransit times andlow dispersion of all traces from#3 on- ward show that the entire flow path between Moulin L41A and theproglacialriverremainedefficient,consistentwithgasanddye tracingofseveralothermoulinsinpreviousmeltseasons(Chandler etal.,2013)andsupportedbyGPSobservationsofsurfacevelocity (e.g. Sole et al., 2013). Efficient drainage development within 2- 3weeksofsurface-bedconnectioniscurrentlyunderestimatedin hydrologicalmodels,includingthosethatdoeventuallydevelopef- ficientdrainageundersimilarconditions(e.g.Werderetal.,2013).

The eﬃcientdrainage persisted through extendedperiods of low surface melt and moulin melt water input (e.g. 12-24July), per- hapsaided by gradual releaseofstoredmelt waterhigherinthe catchment(3-4daysdelay inmelt routingfrom1500melevation inJuly2012:vanAsetal.,2017).

Inachanneliseddrainagesystemwithrapidevacuationofmelt water,wewouldexpectlocalsubglacialwaterpressureandiceve- locity(nearthemoulin)tocloselyfollowmeltinput,withlittlelag.

Incontrast,theanticlockwise hysteresisbetweenmeltwaterinput andicevelocity during low melt inmid-July, andagainduring 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 ∼⁶ ^hr, ^com- paredwithatime scaleofonly 12to15hr fortracerstotransit the entire >40 km ﬂow path. The hysteresis suggests that local subglacialwaterpressure atL41A isdriven partlyby a non-local source,consistentwithaneﬃcientdrainagesystemreceivingmelt waterinputsfromupstreamregions(L41Aisclosetothepredicted subglacial pathwaydraining a considerable portion of the upper catchment:Fig.1).Thisnon-localup-glaciermeltwatersupplywill havemaintainedhighlocalsubglacialwaterpressurebeneaththe moulin(and,therefore,highericevelocity)forseveralhours after thepeakoflocalwaterinput, dueto thetraveltimesrequiredto delivermeltwaterfrommoredistantupstreamsources.

(8)

The potentially complex relationship betweenlocal 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 beneﬁt of our simultaneous time-series of moulininput,watertransitandicedynamics.

5.3. Implicationsforsubglacialhydrologicalmodels

Subglacial hydrology models have invoked a variety of ﬂow- path geometries for melt water, including sheets,linked cavities, channelsandporoustill(de Fleurianetal.,2018).Rapidtracerre- turns withlow dispersionare notconsistent withﬂowthrougha thinsheet,linkedcavitiesorporoustill,asthesewouldbestrongly dispersive; therefore, our tracer experiments provide strong evidence for the development of channels at the ice-bed interface.

Such channelscould beincisedupwardsintotheice(R-channels:

Röthlisberger,1972),downwardsintothesubstrate(‘Nyechannels’

in bedrock,or ‘canals’insediments), orofcourse both. Ourdata aresupportiveofthosemodelsintheSHMIPexperimentthatde- velopchannelsintheidealisedicesheetgeometry(de Fleurianet al.,2018).

Neither the tracing experiments nor the other data sets can at present identify whether subglacial conduits are incised up- wardsordownwards;however,futurecomparisonofthisdataset with drainagesystem models may enable such distinction if the hydraulic parameters are found to be suﬃciently different (e.g., if dispersion in rough-walled ‘canals’ is much greater than in smooth-walled R-channels).Few in-situ measurementsof conduit roughness characteristics are available, and are limited to acces- sible conduits under relatively thinice (e.g., Gulley et al., 2012).

Importantly, therateofconduitdevelopment maybevery differ- entforincisionintosediments,perhapsfasterifsedimenterosion requires lessenergy per unit volume than icemelt. When chan- nelised drainagepersists forseveral weeks inregions where the orientation of the conduit is parallel to that ofthe basal sliding direction, aconduitincisedintotheicewouldlikelyerodeacon- duitinanyunderlyingsediments(orvice-versa);thisislesslikely ifthechannelorientationandslidingdirectionarenotaligned,as aconduitintheicewouldbemovinglaterallyrelativetothebed.

Downwardsincisionofchannelsinto sedimentsmayalsobelim- itedbytheincreasinglycoarsesubstrate,asonlytheﬁnermaterial canberemovedandthechannelﬂoorbecomes‘boulderarmoured’

(Gulleyetal.,2013),thoughlateraldevelopmentofcanalscouldbe unhindered.Theseboulderscouldthenreduce therateofchannel closurebyicecreepduringperiodsoflowersubglacialwaterpres- sure. Downwards incision of conduits aligned parallel to sliding directioncouldneverthelessbe onemeansbywhichthedrainage systemcan remaineﬃcientthrough periodsoflow melt input,if theseconduitstakelongertoclose.Inreality,wesuspectthatmul- tipledrainageformsmayexistsimultaneously.

6. Conclusions

We ﬁnd clear evidence for the development of eﬃcient sub- glacialdrainageextendingfromthemargintobeneathiceatleast

∼⁹⁰⁰^m^thick, ⁱⁿâ ^regionôf^the âblation^zone^where ^theêxtent of seasonal drainage evolution has remained uncertain in previ- ousfieldstudiesandnumericalsimulations.Wehavealsodemon- strated rapiddevelopment ofefficient drainage(bylate June,fol- lowingestimatedsurface-bedconnectionon3June)fromamoulin initiated by a fracture intersecting a supraglacial channel, rather thanbylakehydrofracture.Theefficientdrainagesystempersisted andshowedlittleresponsetoveryvariablemoulin inputsinclud-

ing extreme melt events and a 2-week period of subdued melt (Figs.6,7).Giventhe manysimilarities inseasonalevolutionob- served inthiscatchment andelsewherein Greenland- including intidewaterglaciersupstreamoftheinﬂuenceofthecalvingfront (Soleetal.,2011;Nienowetal.,2017)- weconsiderourdatatobe generallyrepresentativeofsimilar sectorsoftheGrIS ratherthan site-speciﬁc.

Despite considerable advances in the numerical and physical complexity ofsubglacial drainagemodels,their development and validationremainlimitedby a lackofobservations.Thiscompre- hensivedatasetwillnowprovideauniquereferenceagainstwhich totestmodelconfigurationsandparameters(forexample,infuture SHMIPphases:de Fleurianetal., 2018).The dataare particularly relevantastheyrepresentaregionoftheablationzonewherecon- flicting results in previous work show that calculated subglacial drainagecharacteristicsare sensitive to thechoice ofmodelcon- figuration.

CRediTauthorshipcontributionstatement

Conceptualization:DC, JW, PN, AH. Methodology: DC, JW, PN, SD,AT, JT,AH.Formalanalysis:DC.Investigation,writing,editing:

Allauthors.

Declarationofcompetinginterest

Theauthorsdeclarethattheyhavenoknowncompetingﬁnan- cialinterestsorpersonalrelationshipsthatcouldhaveappearedto inﬂuencetheworkreportedinthispaper.

Acknowledgements

Theprojectwas fundedby UKNERCgrant nos.NE/H023879/1 (JW), NE/F021380/1 (PN) and NE/G005796/1 (AH). JW acknowledges support via a Royal Society Wolfson Merit Award and a LeverhulmeTrust ResearchGrant (RPG-2016-439)andfellowship.

AH gratefully acknowledges support by the Research Council of NorwaythroughitsCentresofExcellencefundingscheme,project number223259 andan Academy ofFinland ArcI visiting fellow- shiptotheUniversityofOulu.Thankstothemanyﬁeldassistants whocontributedtotheﬁeld logisticsandwatersampling. Thanks alsoforthecommentsfromtwoanonymousreviewers.

References

Andrews,L.C.,Catania,G.A.,Hoffman,M.J.,Gulley,J.D.,Lüthi,M.P.,Ryser,C.,Hawley, R.L.,Neumann,T.A.,2014.Directobservationsofevolvingsubglacialdrainage beneaththe GreenlandIceSheet.Nature 514,80–83.https://doi.org/10.1038/ nature13796.

vanAs,D.,BechMikkelsen,A.,HoltegaardNielsen,M.,Box,J.E.,ClaessonLiljedahl, L.,Lindbäck,K., Pitcher,L.,Hasholt, B., 2017.Hypsometric ampliﬁcation and routingmoderationofGreenlandicesheetmeltwaterrelease. Cryosphere 11, 1371–1386.https://doi.org/10.5194/tc-11-1371-2017.

Banwell,A.F.,Hewitt,I.,Willis,I.,Arnold,N.,2016.Moulindensitycontrolsdrainage developmentbeneaththeGreenlandicesheet.J.Geophys.Res.,EarthSurf. 121, 2248–2269.https://doi.org/10.1002/2015JF003801.

Banwell,A.F.,Willis,I.C., Arnold,N.S.,2013.Modelingsubglacialwaterroutingat Paakitsoq,WGreenland.J.Geophys.Res.,EarthSurf. 118,1282–1295.https://

doi.org/10.1002/jgrf.20093.

Bartholomew, I., Nienow, P., Sole, A., Mair, D., Cowton, T., Palmer, S., Wad- ham,J.,2011.Supraglacial forcingofsubglacialdrainageintheablationzone of the Greenland ice sheet. Geophys. Res. Lett. 38. https://doi.org/10.1029/ 2011GL047063.

Bingham,R.G.,Nienow,P.W.,Sharp,M.J.,Copland,L.,2006.Hydrologyanddynamics ofapolythermal(mostlycold)HighArcticglacier.EarthSurf.Process.Landf. 31, 1463–1479.https://doi.org/10.1002/esp.1374.

Bullister,J.L.,Wisegarver,D.P.,Menzia,F.A.,2002.Thesolubilityofsulfurhexaﬂu- orideinwater and seawater. Deep-Sea Res., PartI: Oceanogr. Res. Pap. 49, 175–187.https://doi.org/10.1016/S0967-0637(01)00051-6.