Greenland subglacial drainage evolution regulated by weakly connected regions of the bed

Abstract

Penetration of surface meltwater to the bed of the Greenland Ice Sheet each summer causes an initial increase in ice speed due to elevated basal water pressure, followed by slowdown in late summer that continues into fall and winter. While this seasonal pattern is commonly explained by an evolution of the subglacial drainage system from an inefficient distributed to efficient channelized configuration, mounting evidence indicates that subglacial channels are unable to explain important aspects of hydrodynamic coupling in late summer and fall. Here we use numerical models of subglacial drainage and ice flow to show that limited, gradual leakage of water and lowering of water pressure in weakly connected regions of the bed can explain the dominant features in late and post melt season ice dynamics. These results suggest that a third weakly connected drainage component should be included in the conceptual model of subglacial hydrology.

ARTICLE

Received 10 Jun 2016 | Accepted 10 Nov 2016 | Published 19 Dec 2016

DOI: 10.1038/ncomms13903 OPEN

Greenland subglacial drainage evolution regulated by weakly connected regions of the bed

Matthew J. Hoffman1, Lauren C. Andrews2, Stephen A. Price1, Ginny A. Catania3,4, Thomas A. Neumann2, Martin P. Lthi5, Jason Gulley6, Claudia Ryser7, Robert L. Hawley8 & Blaine Morriss9

Penetration of surface meltwater to the bed of the Greenland Ice Sheet each summer causes an initial increase in ice speed due to elevated basal water pressure, followed by slowdown in late summer that continues into fall and winter. While this seasonal pattern is commonly explained by an evolution of the subglacial drainage system from an inefcient distributed to efcient channelized conguration, mounting evidence indicates that subglacial channels are unable to explain important aspects of hydrodynamic coupling in late summer and fall. Here we use numerical models of subglacial drainage and ice ow to show that limited, gradual leakage of water and lowering of water pressure in weakly connected regions of the bed can explain the dominant features in late and post melt season ice dynamics. These results suggest that a third weakly connected drainage component should be included in the conceptual model of subglacial hydrology.

1 Fluid Dynamics and Solid Mechanics Group, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA. 2 Cryospheric Sciences Laboratory, NASA Goddard Space Flight Center, Greenbelt, Maryland 20771, USA. 3 Institute for Geophysics, Jackson School of Geosciences, The University of Texas at Austin, Austin, Texas 78758, USA. 4 Department of Geological Sciences, Jackson School of Geosciences, The University of Texas at Austin, Austin,Texas 78758, USA. 5 Glaciology and Geomorphodynamics Group, Department of Geography, University of Zrich, 8057 Zrich, Switzerland. 6 School of Geosciences, University of South Florida, Tampa, Florida 33620, USA. 7 Laboratory of Hydraulics, Hydrology and Glaciology, Swiss Federal Institute of Technology (ETH) Zrich, 8093 Zrich, Switzerland. 8 Department of Earth Sciences, Dartmouth College, Hanover, New Hampshire 03755, USA.

9 Cold Regions Research and Engineering Laboratory, Hanover, New Hampshire 03755, USA. Correspondence and requests for materials should be addressed to M.J.H. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13903

a

In the ablation zone of the Greenland Ice Sheet (GrIS), the drainage of surface melt to the ice sheet bed via moulins and crevasses causes ice ow acceleration every summer14. The

inux of surface melt overwhelms the capacity of the subglacial drainage system, increasing subglacial water pressure and reducing basal traction of the ice sheet, inducing enhanced sliding5. However, ice speed subsequently lowers over the summer despite sustained meltwater input14,6, which is generally explained by increasing efciency of the subglacial drainage system2,3,6.

This seasonal evolution of subglacial drainage beneath the GrIS is currently interpreted in the context of traditional theory of a two-component subglacial drainage system consisting of distributed and channelized drainage14,6. Theory suggests that at low subglacial discharge, drainage occurs through inefcient, distributed pathwayssuch as linked cavities formed in the lee of bedrock bumps as the glacier slides over the bed or pathways eroded into basal sedimentsfor which increasing water ux leads to increased water pressure and sliding711. It is thought that when a critical discharge is reached in the distributed drainage system, dissipation of heat within the water ow causes a positive feedback between melting of the ice roof and cavity growth, leading to the formation of discrete, efcient channels incised into the ice above5,8,9,11,12. Such channels would then rapidly evacuate water from the distributed drainage system and lower the water pressures over a large region, terminating a sliding event despite sustained meltwater inputs to the drainage system5,1113. We use a model to illustrate the need for an additional type of drainage system, here termed weakly connected (Fig. 1), and that evolution within this system is responsible for previously unexplained seasonal adjustments to ice velocity.

Observations and modelling suggest that the current conceptual model of the subglacial hydrologic system overemphasizes the role of channelization in controlling GrIS subglacial drainage system capacity and water pressure. Depressed ice speeds persist through fall and much of winter1,6,14,15, which is inconsistent with a timescale of hours to days for channel collapse under the thick GrIS ice3,6,16. Additionally, modelling has indicated that gentle surface slopes on the ice sheet should suppress channel formation16,17, and low water pressures18 and depressed summer ice speeds19 characteristic of highly efcient channels are only observed near the ice sheet margin. Therefore, in interior regions there may be unrecognized drainage capacity elsewhere in the system. Similarly, Andrews et al.4 recently described direct observational evidence that even where channelization occurs, it is unable to explain lowering ice speed during late summer. At the extensively studied drill site FOXX in western Greenland4,14,2023 (Fig. 2a), water pressure in moulins feeding subglacial channels showed little change during the latter part of the melt season, yet velocities were observed to decrease over this same time period (Fig. 3a). Simultaneous observations of declining borehole water pressures sampling poorly connected regions of the bed suggested that weak drainage out of these isolated regions could potentially account for the unexpected increasing system efciency.

Here, we use a subglacial hydrology model (see Methods: Subglacial hydrology model overview) to demonstrate that the observations at FOXX can be explained by gradual evacuation of water from weakly connected, but spatially extensive, areas of the bed. We suggest that these areas exert the dominant control on the large-scale subglacial water pressure and basal resistance felt by the ice sheet, which we show by using the modelled ice effective pressure to force an ice dynamics model13,24 (see Methods: ice velocity calculations) to reproduce observed summer ice speed changes. Our subglacial hydrology model includes coupled components for distributed drainage, channelized drainage and drainage from weakly connected

b

c

Figure 1 | Conceptual model of three-component subglacial hydrologic system for the GrIS. (a) Onset of the melt season: a large fraction of the bed is composed of weakly connected cavities at a higher water pressure than the surrounding distributed system. (b) Middle of melt season: meltwater draining from the surface through moulins is largely accommodated by the formation of efcient channels (dashed grey outline). Concurrently, some of the weakly connected cavities have leaked water, lowering their water pressure, due to increasing connectivity with the rest of the system initiated during periods of pressurization. (c) End of melt season: channels collapse within days after melt inputs cease, but the partially drained weakly connected cavities take months to recharge by basal melting, leaving higher integrated basal traction than before summer began.

2 NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13903 ARTICLE

a

1,500

1,200

900

600

300

0

GULL

MODEL SITE

Elevation (m)

FOXX

300

55 50 45 40 35 30 25 20 15 10 5 0

x-position (km)

b

1,500

1,000

Elevation (m)

500

0

100 80 60 40 20

0 20

y

x-position (km)

Figure 2 | Study area and model domain. (a) Ice surface and bed longitudinal prole for GrIS owline at FOXX study site (thin lines) compared with those for idealized model domain (thick dashed lines). (b) Schematic of model domain showing ice surface (light blue) and bed topography (orange). The model study site is indicated by red vertical line. The channel added during the summer simulations is shown by the dark blue line. Purple line on ice sheet surface indicates equilibrium line above which no runoff occurs. Note that the gure is not to scale and the actual grid spacing used in the model is

Dx Dy 200 m. Inset: schematic of how weakly connected system is represented in the numerical model. Within each grid cell (black box), a fraction of

the area, fw, is assumed to be covered by patches of the weakly connected system (dark regions), with the remaining area being composed of the distributed system (light regions). For simplicity, the shape of the weakly connected system is assumed to be many small circular regions, though in reality we expect it to be highly irregular.

regions of the bed (see Methods: weakly connected drainage model). While the former two components have been routinely included in subglacial hydrology system models11, this is the rst application of the latter component. The simulations use an idealized ice geometry based on drill site FOXX and are forced by surface meltwater input into the channelized system based on observed melt rates and by observed ice sheet sliding14 (Fig. 4a; see Methods: model setup). We conceptualize the weakly connected regions as discrete patches of linked cavities (Fig. 2b), similar to the distributed drainage component (Supplementary Methods), but with a much lower hydraulic connectivity (see Methods: weakly connected drainage model; Supplementary Table 1). Based on observed diurnal and seasonal changes in water pressure in moulins and boreholes4,20, we assume these patches cover roughly two-thirds of the area of the bed (see Methods: model sensitivity to weakly connected area fraction). There is no through-ow between individual weakly

connected patches; instead, water movement occurs as a leaky exchange with the surrounding distributed system in each grid cell of the model and which is prescribed to become more transmissive over summer.

ResultsModelled channelized drainage. Hydraulic head (a measure of water pressure, see equation (8)) in the modelled channel demonstrates correspondence to hydraulic head measured in a moulin at FOXX (Fig. 4b), which was interpreted as representing the channelized drainage system4. Our model results show that an efcient channel remains in approximate equilibrium with melt inputs for the second half of summer (Figs 4 and 5). This is consistent with observations and channel modelling performed by Andrews et al.4 and in contrast to channel modelling performed for locations further inland on the GrIS where atter

NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13903

a b c

180

180

180

Ice surface speed (m yr1 )

160

Observed

160

Modelled: static

235

160

230

140

140

140

225

120

120

120

220

100

100

Day of 2012

100

215

210

80

80

80

205

60

60

60

200

40

40

40

195

200 150 100 50 0

200 150 100 50 0 200 150 100 50 0

Channel hydraulic head (m) (relative to floatation elevation)

Figure 3 | Seasonal evolution of the relationship between subglacial water pressure and ice speed. Each plot shows the minimum and maximum daily values of hydraulic head in the moulin-channel system and ice surface speed for the second half of the 2012 summer. (a) Observed relationship showing seasonal hysteresis of lowering ice speed for the same moulin head as summer progresses. Modied from Andrews et al.4 (b) Modelled relationship with increasing permeability of the weakly connected regions of the bed showing similar relationship. (c) Modelled relationship from control simulation with static weakly connected system showing lack of coherent seasonal evolution.

Ablation rate, m d1

0.40.20.0 50

300

200

100

0

Slidingspeed, m yr1

Hydraulic head (m),

relative to floatation elevation

0

50

100

150

200

250

Modelled BH4

BH6 BH7

Modelled Moulin3

150 160 170 180 190 200 210 220 230 240

Day of year 2012

Figure 4 | Model results for subglacial hydraulic head. (a) Observed melt rate (black) and ice sliding speed (grey) used to force the subglacial hydrology model. (b) Modelled and observed subglacial hydraulic head in the weakly connected and channelized systems. Modelled channel hydraulic head (blue) reproduces most features of the measured moulin hydraulic head (green). Modelled hydraulic head in the weakly connected system (red) reproduces most features of the hydraulic head measured in boreholes (pink, orange and maroon). BH borehole. The datum for both observations and

model results is the elevation corresponding to local oatation pressure.

slopes and thicker ice have been proposed to prevent channelization16,17.

A relatively large initial channel area is required in our model to accommodate surface melt draining to the bed quickly enough for pressures below oatation to develop within a few days, as found in previous models during large pulses of melt25,26. Previous modelling studies16,17,27 have highlighted the unrealistically long time scales required for such large channels to develop. However, our results and those of Andrews et al.4 demonstrate that if channels are able to grow large enough to accommodate surface melt inputs, they can explain the pressure record observed during the second half of summer. Therefore we consider the possibility that channel formation is preconditioned26,28. For example, extensive ooding of the bed

at the start of the melt season from supraglacial lake drainage3,25,29,30 or other moulin development opens cavity space that facilitates channel formation. This is supported by model results (Fig. 5) showing that channel growth is restricted until the cavity space (represented as water layer thickness) in the distributed system has grown to its maximum seasonal value. These results are consistent with studies suggesting an important role of distributed drainage in developing drainage efciency13,17. Alternatively, year to year persistent moulin locations31 may facilitate repeated occupation by channels of the same locations and rapid channel growth through cumulative erosion of basal sediments creating preferential ow pathways26,28,32. Finally, the prescription of a single channel in our model rather than an anastomosing network of channels likely hinders our models

4 NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13903 ARTICLE

a

Channel area (m2 )

10

9

8

7

6

5

4

3

2

1 140 160 180 200 220 240

b

Water thickness (m)

0.20

0.18

0.16

0.14

0.12

0.10

0.08

0.06

140 160 180 200 220 240

Time (day)

Figure 5 | Modelled subglacial channel area and water layer thicknesses at study site. (a) Modelled channel area. (b) Modelled water layer thickness in distributed (blue) and weakly connected (red) systems.

ability to realistically grow channels while they are small, as the most dominant channels in a network grow in large part by capturing drainage from smaller channels5,26,33.

Modelled weakly connected drainage. Similar to the results for channelized drainage, the modelled hydraulic head in the weakly connected system reproduces the seasonal changes in hydraulic head observed in the boreholes (Fig. 4b). In the model and in measured borehole 4 (and to a lesser extent, borehole 6, see Andrews et al.4 and Fig. 3d), hydraulic head at the beginning of summer is above that corresponding to ice overburden pressure and then gradually decreases over the course of the season. The trends in hydraulic head in the modelled and measured weakly connected system follow those in the modelled channel during the rst part of summer as the channel grows to equilibrium, supporting our hypothesis that dropping summer borehole pressures are caused by slow, down-gradient leakage towards well-connected portions of the drainage system. After the channel reaches its equilibrium size (day B200220; Fig. 5), hydraulic head in the weakly connected system continues to drop due to the enhanced connectivity between the weakly connected system and the rest of the drainage system.

This seasonal pattern is overlain by short-term variations in the weakly connected system that contrast with behaviour in the well-connected drainage system. On almost all days, modelled hydraulic head in the weakly connected system is out of phase

with channel pressure, with diurnal amplitudes of a few percent of overburden, matching borehole observations. This out of phase behaviour, observed for isolated or weakly connected boreholes on both the GrIS and mountain glaciers4,20,3437, has been explained as the transfer of normal stress from hydraulically well-connected regions of the bed20,34,36,37, which at our site occurs over kilometers20. Because our model does not include diurnal variations in normal stress transfer, our results indicate that diurnal variations in cavity opening rate associated with changes in ice sliding are a possible alternative or additional mechanism for inducing pressure variations that are out of phase with the well-connected drainage system4,13,37. We note that the modelled diurnal amplitude of these pressure variations grows unrealistically near the end of summer, which may indicate that a more sophisticated description for cavity opening than our simple linear parameterization (equation (2)) is required to explain all of the observations.

Modelled ice velocity. Having validated model water pressure results, we consider implications of the inclusion of the weakly connected drainage component on ice dynamics. The ice dynamics model (see Methods: Ice velocity calculations) reproduces the basic features of the ice velocity observations: while higher channel pressure results in higher ice speed on each day, there is a drop in ice speed for the same channel pressure as the summer progresses (Fig. 3a,b). To eliminate the possibility

NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13903

that the seasonal hysteresis in the relationship between channel water pressure and ice velocity could be caused by the channelized and distributed components evolving in capacity as the summer progresses, we perform an additional control simulation where those two components are free to evolve as in the rst simulation, but the weakly connected component is static with a prescribed, xed water pressure. In this control simulation, there is limited seasonal evolution (Fig. 3c) supporting our conclusion that the weakly connected system is controlling the late summer slowdown and that inclusion of the weakly connected system is necessary to reproduce the observed changes in ice dynamics.

It is notable that the observations and both model versions exhibit variations in the slope of the lines dening the minimum and maximum channel hydraulic head and ice surface speed on each day (Fig. 3). These varying slopes represent differences in the sensitivity of sliding to changes in effective pressure as effective pressure changes. This changing relation is expected from theory3840 and observations4,41, and we conrm that both models share a similar sensitivity with the observations (see Methods: ice velocity calculations; Supplementary Fig. 1). Thus, the varying slope of the lines in Fig. 3a is expected from having a different range of effective pressure on each day and is not associated with evolution within the weakly connected system.

In contrast, evolution within the weakly connected system is required to explain the lowering of the lines in Fig. 3a as the summer progresses (the downward propagating rainbow pattern in the plot). This downward trend represents seasonal changes in the relationship between moulin water pressure and slidingthe same water pressure induces less sliding later in the season. To quantify this behaviour, we calculate the Pearson product moment correlation coefcient between the observations and model results for the intercept of each line with hydraulic head of

75 m (this value chosen as approximately the centre of the range of hydraulic head values). The seasonal evolution in the model with evolving weakly connected system has a signicant positive correlation with the observations (r 0.36, P 0.01),

while the model with the static weakly connected system is not signicantly correlated with the observations (r 0.12,

P 0.43).

Because some observations on mountain glaciers have shown that isolated boreholes can become connected during periods of high water pressure in the active drainage system34,35, we consider an alternative hypothesis that it is changes to the areal extent of the weakly connected system, and not changes in its connectivity, that cause declining ice speed. However, we nd that such a parameterization primarily affects the diurnal range of ice speed and results in minimal changes to ice speed at the seasonal scale (see Methods: model sensitivity to weakly connected area fraction; Supplementary Figs 2 and 3). While we cannot rule out the possibility that the area fraction of the weakly connected system changes modestly during summer, our model results indicate it is not the primary mechanism causing ice speed to drop.

Recharge time scale. Tedstone et al.42 nd annual ice motion is more strongly correlated to summer melt volume from the previous 1 to 4 years than summer melt from the corresponding year, and they suggest that this multiyear response is due to gradual net drainage of water stored in unchannelized regions. Our model predicts that water layer thickness in the distributed system remains elevated over the entire summer relative to its pre-summer value, while the water layer thickness in the weakly connected system gradually lowers during summer (Fig. 5). Thus, based on our results, net summer drainage only occurs from the weakly connected system, suggesting that it is these regions that

are controlling the multiyear changes in ice motion. A simple calculation of the time scale of recharge of the weakly connected system made using the basal melt rate and assuming no water drains out suggests that it would take the weakly connected system B2.0 years to return to its original water thickness. This estimate is a minimum value because water would continue to drain out as melt rells the system unless exchange completely stopped, which is unlikely. If exchange between the weakly connected and distributed systems is included, assuming the hydraulic gradient at the end of summer and that the permeability immediately returns to its winter value at the end of summer, the recharge time scale increases to 5.1 years. Of course, this is an illustrative time scale because the extent to which weakly connected cavities drain during summer is expected to have signicant variation, and not all cavities demonstrate signicant drainage4,20. Based on these estimates, varying degrees of drainage from the weakly connected system during each summer provide a plausible explanation for the multiyear self-regulation of annual velocity hypothesized by Tedstone et al.42

Conceptual model. While subglacial drainage has generally been categorized into distributed and channelized components5,1012, there is ample evidence that large fractions of the beds of glaciers are inactive and largely disconnected from the pathways directly conveying surface melt that has drained to the bed. These regions are identied by boreholes that drain slowly when reaching the bed and exhibit diurnal water pressure variations that are out of phase with meltwater input and ice speed4,20,3437,43. In many areas, these inactive regions cover large fractions of the bed; Hodge44 estimated 90% of the bed of South Cascade Glacier in Washington, USA, was hydraulically isolated. This behaviour was observed at all three of the FOXX boreholes, as well as three boreholes at an additional study site4,20, and boreholes only 30 m apart showed no obvious connection and large hydraulic gradients20, suggesting that weakly connected drainage likely covers most of the GrIS bed in this region. The implications of large fractions of the bed being relatively isolated at the ice sheet scale has not been previously considered.

Our modelling results conrm that these isolated or weakly connected regions of the bed can play an important role where present by virtue of their large area fraction; modest changes in effective pressure in these extensive regions will have a large impact on basal traction. Importantly, ice dynamics respond to the integrated basal traction over both well-connected and poorly connected regions of the bed, and poorly connected regions therefore moderate active drainage regions40,43. The heterogeneous nature of basal drainage (active regions interwoven with weakly connected regions on length scales of metres to tens of metres) is smoothed out because ice dynamics responds to basal traction at the length scale of a few ice thicknesses4547 (4Bkm for GrIS). Both observations and our modelling efforts indicate that, despite the apparent isolation of these areas, they are not entirely static. Often, water pressure changes occur in these areas when water pressure in the active moulin-connected system increases or inactive patches switch to active behaviour4,3436.

Based on this new understanding, we propose a three-component conceptual model for subglacial hydrology that adds a weakly connected system to the traditional distributed and channelized forms of subglacial drainage (Fig. 1). We nd this third component necessary to model the ice sheet velocity response properly. Meltwater input of sufcient magnitude drives the formation of channels which are the primary control on overall system efciency. Channels, as inherently linear features,

6 NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13903 ARTICLE

have little direct effect on ice sheet basal traction, but indirectly control basal traction by acting to lower subglacial water pressure within the spatially expansive distributed system. The strong connectivity between these two components acts to couple the channel-driven drainage state to the ice dynamics. The new qualitative behaviour we are proposing occurs in regions of the distributed system that are so weakly connected to the rest of the bed that they act as a regulator on the active system43, meaning changes there affect integrated basal traction.

We have envisaged the weakly connected system being composed of cavities formed in the lee of bedrock bumps or clasts in lodged till. The low hydraulic connectivity that differentiates these regions from the rest of the distributed system may be due to local bed geometry causing smaller cavity orices or the presence of low hydraulic conductivity till. Indeed, extensive subglacial till has been identied at FOXX21, and stick-slip ice motion observed seismically is associated with spatially differing connectivity of till23. More generally, both modelling48,49 and observations21,50 have argued for the existence of signicant amounts of till underlying marginal regions of the GrIS. We nd an alternative formulation of the weakly connected system being composed entirely of till yields similar results (Supplementary Methods, Supplementary Table 2 and Supplementary Fig. 4), and we hypothesize that the weakly connected system is composed of a mixture of cavities and till. We emphasize that the weakly connected system, as conceptualized here, is a subset of distributed drainage with extremely low hydraulic conductivity, and an alternative modelling approach would be to model a distributed system with spatially variable conductivity at high grid resolution. However, we propose it as a separate drainage component due to the qualitatively different behaviour observed there; a broad parameterization of weakly connected drainage is consistent with the parameterized descriptions typically employed for distributed and channelized drainage1012.

DiscussionThe impact of weakly connected drainage on the seasonal evolution of hydrodynamic coupling of the GrIS explains apparent contradictions in a less complex conceptual model. First, our model can explain recent observations4 of water pressure in moulins that show subglacial channels quickly equilibrate with meltwater input while ice velocity continues to decline over the summer. Second, the time scale for re-equilibration of low pressure conduits after the cessation of melt is hours to days3,6,16, yet the impact of summer melt on GrIS ice dynamics is clearly sustained through most of the winter and may extend for multiple years6,15,42. This points to a slow-reacting component of the basal drainage system. The gradual dewatering of the weakly connected system would take longer to recharge than the other systems42, which is conrmed by a recharge time scale of years in our model.

While the annual cycle of dewatering and recharge of weakly connected regions of the bed suggests a resilience of GrIS ice dynamics to increasing melt42, the likely importance of erodible till for governing the connections of the weakly connected system to more active regions of the drainage system could allow rapid changes in ow resistance21. Furthermore, the spacing of channels, and the moulins that feed them, may control the extent to which weakly connected areas are tapped during the melt season. Future work should attempt to improve direct observations of these weakly connected regions of the bed and improve their description in models to better constrain future changes of hydrodynamic coupling of the GrIS and its effects on ice sheet mass balance.

Methods

Subglacial hydrology model overview. The subglacial hydrology model is based on the model described by Hoffman and Price13 for coupled distributed and channelized drainage. The model takes a continuum approach to subglacial drainage that captures the bulk behaviour of the subglacial drainage system at spatial scales relevant to ice dynamics (B102 m) rather than attempting to model every basal conduit explicitly. The three modes of drainage, distributed, channelized and weakly connected, are modelled as separate components, each with appropriate physics, that are coupled together by exchanges of water driven by gradients in hydropotential between the systems. Distributed drainage is modelled as a continuous macroporous sheet in two dimensions on the primary model grid, while channelized drainage is represented by a single one-dimensional channel located along a line of grid cell edges of the primary model grid taking up no area within the primary model grid (Fig. 2b). Summaries of the previously described13, distributed and channelized drainage components are included in the Supplementary Methods for completeness. The weakly connected cavity system is represented as a subgrid element within the distributed system (Fig. 2b, inset). Within each grid cell of the model, a fraction of the bed, fw, is assumed to be covered by the weakly connected cavity system, with the remaining fraction, 1 fw,

covered by the through-owing distributed drainage system. Surface melt draining from the surface is delivered to the channelized system, which can exchange water with the distributed system. The distributed system can exchange water with both the channelized system and the weakly connected system, while the weakly connected system and the channelized system have no direct exchange. In all three systems, conduit space is assumed to be entirely lled with water at all times, a common approach in subglacial hydrology models (c.f. Schoof et al.51). All model parameters are listed in Supplementary Table 1 and are spatially uniform unless otherwise mentioned.

Weakly connected drainage model. The weakly connected component is implemented as a discontinuous reservoir that can exchange water with the surrounding distributed drainage system within each grid cell (Fig. 2b, inset). The weakly connected areas are prescribed to cover two-thirds of each grid cell, with the remainder covered by the distributed system. This fraction is chosen to capture the fact that all boreholes observed at drill site FOXX4,20 exhibited characteristics of hydraulic isolation from the active drainage system over most of the observing period while avoiding being overly prescriptive with the new component.

Using a macroporous sheet continuum formulation for the weakly connected system analogous to that used for the distributed system (Supplementary Methods), mass conservation of the water thickness in the weakly connected component (hw) is described by the balance between locally generated melt (mw)

and exchange of water with the sheet (gw),

@hw

@t

mw rw

gwAw ; 1

where rw is the density of water, Aw is the area of the weakly connected system within each grid cell, dened as fwDxDy.

Evolution of cavity space within the weakly connected component is the same as for the distributed system, a balance of cavity opening by sliding of the ice over bedrock bumps and close by creep of their ice roof:

@hw

@t ub

j j

hr hwlr

2A

27 hwN3w; 2

where ub is the sliding velocity, A is the temperature-dependent rate factor for ice deformation, Nw is the effective pressure in the weakly connected system, and hr and lr are parameters describing the height and wavelength, respectively, of bumps on the bed. The value for A for the basal layer (Supplementary Table 1) is approximated from results of borehole temperature and deformation observations and ice ow modelling performed by Ryser et al.14 where the basal layer is temperate and an enhancement factor of 2.54 (taken here as 2.7) is appropriate for ice from the late Wisconsin found at this depth.

Energy for local melting within the weakly connected component, mw, also matches that for the distributed system:

mwLG ub tb; 3

where G is the geothermal heat ux, tb is the basal traction vector and L is the latent heat of fusion of water.

For all components, hydraulic potential is dened as

f d;c;w

f grwgzb p d;c;w

f grwgzb pice N d;c;w

f g; 4 where the subscript {d,c,w} indicates one of the distributed, channel, or weakly connected systems, respectively. g is the acceleration due to gravity (m2 s 1), zb is the bed elevation (m), p is water pressure (Pa) and pice is the ice overburden

pressure (Pa). For the distributed and channel systems the ice overburden pressure is calculated using a hydrostatic assumption, pice rigH, where ri is the ice density

(kg m 3) and H is the ice thickness (m). To allow the weakly connected system to attain water pressure greater than oatation, as observed, we include a simple representation for normal stress transfer34,36 in the weakly connected system by

NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13903

increasing the ice overburden pressure according to

pice rigH 1 rw

; 5 We set rw to a constant value of 0.13, because this corresponds to the maximum value of borehole water pressure observed at our study site20, and thus suggests a limiting pice value when Nw 0. Note that this simple representation does not

include temporal variations in normal stress transfer.

Similar to the coupling between the channel and distributed system (Supplementary Methods), the weakly connected component is coupled to the distributed system by calculating a ux, gw, between the surrounding distributed system and the weakly connected system within each grid cell based on a Darcy ow law:

gw

k0wh3w

Zw

fw fd

Ds Pw; 6

where fw and fd are the hydraulic potential in the weakly connected and distributed systems, respectively, k0w is the permeability between the weakly connected and distributed systems, Pw represents the perimeter between the two systems within each grid cell and s is a characteristic spacing between the two systems.

While in reality the boundary between the two systems is likely to be complex, for simplicity in our parameterization we assume a simple geometry for the weakly connected patches within each grid cell for the purposes of generating reasonable values for Pw and Ds. Specically, we assume each weakly connected patch is a circle of radius 10 m (based on observations of differing connectivity on that scale in GrIS and elsewhere4,34,35,52) and set Ds 10 m. For the chosen values of the

fractional area of the weakly connected system (fw 0.67) and grid spacing

(Dx Dy 200 m) we can calculate PwE5,355 m from basic geometry. While

other choices could be made here, these details are not important from a practical standpoint without detailed knowledge of the bed, as these parameters, along with k0w, are free variables in equation (6). Essentially, equation (6) parameterizes the exchange of water between the two systems primarily as a function of the difference in hydraulic potential and a permeability constant. As with the exchange between the distributed and channelized components (gc; Supplementary Methods), the exchange between the weakly connected and distributed systems can occur in either direction, with water always moving from the component with higher hydraulic potential to that with lower.

The permeability between the weakly connected and distributed systems, k0w, is

many times smaller than the permeability within the distributed system itself, k0 (by B9 orders of magnitude in our model conguration; Supplementary Table 1), dening the weakly connected nature of this new component. To parameterize increases in connectivity between isolated regions of the bed and the active drainage system that are inferred to occur during summer, we allow k0w to increase over the course of the summer by

k0wk0w krate t ts

; 7 where k0w is a base value during winter, and krate is a rate at which the permeability increases beginning at the start of summer, ts. We choose the value of krate to reproduce observations of borehole water pressure (Fig. 4b). With our parameter choices, k0w increases by a factor of B70 over summer; the end

of summer value of k0w remains B8 orders of magnitude smaller than the permeability constant for the distributed system. While it changes in time, the permeability is spatially uniform. This is a simple parameterization for the increased connectivity in some boreholes observed during summer at site FOXX4,20 and at mountain glaciers34,35. It should be highlighted that there currently is little observational or physical basis by which to construct a governing relation for how permeability may evolve during periods of meltwater input, and improving on the simple linear relationship used here is a critical area for further research. Having explored a number of possible, more complicated ad hoc relations, we found that any relation that caused k0w to increase substantially during summer generated similar qualitative behaviour.

Acknowledging that the weakly connected system may be a subset of distributed drainage but with very low permeability, an alternative implementation would be to directly model a distributed system with spatially varying permeability, avoiding the need for a new component. From a practical perspective, the subgrid parameterization of a third component used here is advantageous for two reasons. First, because GPS and satellite measurements of ice surface velocity indicate smooth velocity elds, it can be inferred that spatial variability in effective pressure and associated drainage conditions at the bed is at the length scale of an ice thickness or less. Representing such spatial heterogeneity at the grid scale would require a very ne grid, while the subgrid parameterization allows the weakly connected system to be represented at a coarser resolution, making this approach transferrable to large-scale ice sheet models. Second, explicitly modelling variable permeability of the distributed system would require additional assumptions and parameters about the spatial distribution of these variations.

Our parsimonious parameterization of the weakly connected system leaves room for additional complexity as empirical knowledge of the system increases. We have assumed in our model that all parts of the weakly connected system have changing permeability during summer (equation (7)), while borehole observations suggest that only some weakly connected cavities become leakier during summer4. Similarly, some studies have observed weakly connected cavities becoming fully

connected during periods of high water pressure in the surrounding drainage system34,35, which would correspond to temporal changes in our fw parameter that we have kept steady in time (see Methods: model sensitivity to weakly connected area fraction for additional discussion). Additionally, the bedrock geometry for the weakly connected system may have different characteristics (hr, lr) than for the distributed system. Though we acknowledge the inherent complexity of the subglacial system, we have chosen the simplest formulation that includes the dominant processes inferred to occur. Our parameterization of the weakly connected system is meant to represent the mean conditions of the bed and thus will not necessarily directly reect unique measurements from specic boreholes.

Model setup. We generate a simplied model domain that is consistent with the basic geometry of our study site (Fig. 2). The domain represents a 100 km long sector of the GrIS from margin to lower accumulation zone, with our study site located 25 km inland from the ice margin. The domain is 5 km wide with periodic lateral boundaries to approximate the typical width of a supra- and sub-glacial catchment in this region of GrIS (estimated moulin density in this region is0.20.25 km 1 (refs 1,53)). The ice sheet geometry is a at bed with a plastic glacier shape51,54 assuming a constant basal shear stress of 105 Pa. This geometry provides an idealized but consistent setting with our study site (Fig. 2 and Supplementary Table 3). The subglacial hydrology model is applied for the entire domain, but model results are primarily analysed at the study site location (Fig. 2 and Supplementary Table 3); the larger domain is used only to generate realistic far-eld constraints on the model solution at the study site. There is a zero subglacial water ux boundary condition for the distributed system atx 100 km, and water pressure in the distributed system is assumed to

be at atmospheric pressure at the ice sheet margin, x 0 km.

We choose this simplied geometry over the more complex geometry of the eld site to simplify our interpretation of model results, given uncertainty in the true ice sheet geometry. Previous work in our study area has shown that complex topography can affect sliding and deformation over short distances14 and that active and passive subglacial regions might be affected by bed topography20. However, the two-dimensional basal topography of our study region is only known approximately. Only a handful of ight lines of airborne ice thickness measurements exist in close proximity to our study area, and for those that do, differences in ice thickness at ightline crossover points exceed 100 m in some locations. Additionally, the commonly used BedMachine product55, which uses mass conservation constraints to improve estimates of ice thickness between radar ight lines, does not cover this area. Previous subglacial hydrologic modelling has shown that valleys in the bedrock topography can concentrate subglacial drainage33, so the true bedrock topography is likely to play an important role in dening the large-scale drainage network and the initial rate of subglacial channelization. However, the location of moulins that input water to the basal drainage system is the key control on channel initiation and stability26,28,33. Because both our observations and simulations are focused on a site at a moulin, we expect the simplied model geometry to have little effect on the local subglacial drainage conditions at the study site.

To generate a model initial condition for the summer simulations, we spin up the coupled distributed and weakly connected systems to steady state, assuming no channelized drainage, to represent late winter conditions. We use the winter ice sliding speed to force the model and apply no meltwater forcing. Model parameters hr, lr and k0 (Supplementary Table 1) are chosen to yield a hydraulic head in the distributed system at our study site of about 200 m below local oatation elevation, and k0w is chosen to yield a hydraulic head in the weakly connected system of about 40 m above local oatation elevation.

Model forcing data. The subglacial hydrology model uses time series of two forcing elds, both of which are based on observations in the vicinity of the FOXX drill site. Surface melt input to the subglacial drainage system (o in Supplementary

Equation (5)) represents the volume ux of water drained to the bed through moulins and crevasses, and is the primary source of new water to the subglacial drainage system during summer (basal melting being of much smaller magnitude). Ice sliding speed (ub in equation (2) and Supplementary Equation (2)) controls the rate at which new cavity space opens in the distributed and weakly connected systems, and, less importantly, the frictional melt rate (equation (3) and Supplementary Equation (4)).

We base the melt forcing (Fig. 4a) off of 6-h average ablation rates4 measured using a pressure transducer installed below the surface connected to a surface reservoir56. Because the measured ablation rates do not have the precision necessary to directly generate a time-series with the required time step of 2-h, and the 6-h average rates reduce the amplitude of the diurnal cycle, we generate an idealized diurnal cycle by tting a sine curve to the 6-h time-series where we adjust the amplitude on each day to maintain the measured 6-h average. The ablation sensor began to fail after the large melt event on days 229230, so we estimate the ablation rate for days 231240 using a linear t between mean daily air temperature and ablation rate for the 10 days prior to the large melt event (Pearson correlation coefcient r 0.702, Po0.01). We scale the point ablation rate by the width of our

model domain and apply it as a linear source term along the length of the channel model (o). As discrete moulins can result in the generation of kinematic waves within our channel model, this linear distribution of melt input improves model

8 NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13903 ARTICLE

stability, while reproducing synchronous changes in moulin water level as observed by Andrews et al.4 We apply a linear lapse rate such that runoff is zero at the equilibrium line altitude of 1,100 m. This provides plausible runoff rates for the entire domain, recognizing that model results will be most sensitive to runoff at the location of our study site which are well constrained. Though supraglacial storage on the GrIS is known to change over the melt season57, our melt forcing ignores supraglacial and englacial routing and storage processes that affect the timing and magnitude of melt delivery to the bed58,59, and is meant to be an idealized representation of the diurnal and seasonal variations in melt forcing at our study site.

The ice sliding speed forcing (Fig. 4a) comes from the results of Ryser et al.14 for site FOXX, which subtracts borehole-derived measurements of ice internal deformation from Global Positioning System-derived measurements of ice surface velocity to calculate basal slip. Four gaps of about a day are lled by averaging the diurnal cycles on either side of the gap and smoothing the edges to avoid any sharp transitions. Basal slip is 73% of motion during winter and, though speeds vary dramatically over summer, the contribution of deformation is roughly constant14. We apply this sliding forcing uniformly over our entire model domain. Though ice speed is not spatially uniform in reality, this simplifying assumption will not directly affect our conclusions because we only analyse the model results at our study site. The magnitude of basal traction, tb

j j is held constant at 105 Pa for the

entire simulation; though it should vary in reality, we do not have the information to provide a more accurate value, and, in any case, the results are not sensitive to this choice as it only affects frictional melting, which is orders of magnitude smaller than summer surface melt input.

Summer subglacial hydrology simulations. Using this spun-up state as an initial condition, we model the 2012 summer for days 150 (onset of summer speedup in the GPS record) to 250 (end of GPS and melt forcing observation time-series). For the summer simulation we add a single active channel along the centreline of the domain. While this precludes the formation of a network of channels(c.f. refs 26,33), we expect this approach to resolve the dominant channel effects near our study site, as moulin location strongly controls channel nucleation28,33. The channel initial condition is an area of 10.0 m2 at the margin decreasing linearly to zero at 55 km inland, which yields an area of 5.45 m2 at the study site. This value was chosen to allow water pressure to drop below oatation and diurnal variations to develop within days of melt onset, as seen in the ice velocity record. The channel roughness parameter F is chosen to yield diurnal hydraulic head minima in late summer comparable to that measured in moulins4 and is broadly consistent with previous modelling efforts.

Hydraulic head observations. We compare model results to hydraulic head measured at site FOXX4,20 in 2012, where hydraulic head, d, relates to hydropotential as

df= rwg

; 8 Modelled hydraulic head in the channelized system is compared with hydraulic head measurements in moulin 3 and modelled hydraulic head in the weakly connected system to hydraulic head measurements in boreholes 4, 6 and 7. Because of the simplied model geometry and our desire to make the comparison of our model data general to both FOXX and neighbouring GULL study sites (which both exhibited similar borehole behaviour)4,20, there is no direct correspondence in ice thickness and bed elevation between the model and the observational study sites (Supplementary Table 3). Therefore, we plot both the model and measured hydraulic head using df, the elevation corresponding to oatation pressure, as the vertical datum:

df rwgzb rigH

= rwg

; 9 which is slightly different at each measurement site (see ref. 4, Extended Data

Table 1) and in the model (Supplementary Table 3) due to modest differences in bed elevation and ice thickness.

Ice velocity calculations. Ice velocity for Fig. 3b,c is modelled using the thermomechanical, three-dimensional, rst-order Stokes approximation6063 momentum balance solver in the Community Ice Sheet Model v2.024 as described by Hoffman and Price13. The magnitude of basal traction, tb, is dened by a physically based basal friction law for sliding over hard beds that allows for cavitation and bounded basal drag39,64,

tb

1=nN; L

lmaxA

mmax ; 10

where C is a Coulomb friction constant, lmax and mmax are the wavelength (m) and maximum slope, respectively, of the dominant bedrock bumps, and n is the exponent in Glens ow law. t b is an addition to the original formulation added here to make basal traction calculated for our simplied model domain more realistic. Because our domain lacks the rough, heterogeneous topography of the real ice sheet where resistance to ow is likely to be concentrated in our study area20, we use tob to represent these missing sticky spots40. We apply a constant value of t b

across our study domain (Supplementary Table 1), chosen to reproduce a realistic range of model speeds across the range of effective pressure forcing applied.

With widespread cavitation (i.e., when effective pressure approaches 0 at high water pressure), the friction law becomes a Coulomb friction law of the form (moving t b to the left-hand side to clarify the form)

tb

t bCN: 11 Alternatively, at large effective pressures (low water pressure) the friction law takes a power law form

tb

t b / u1=nb: 12 The basal traction appropriate for ice dynamics is the integrated basal traction of both connected and isolated regions of the bed43,

tb

t btbd 1 fw

tbwfw; 13 where tbd and tbw are the basal traction in the distributed and weakly connected systems, respectively.

We approximate equation (13) by assuming the effective pressure used in equation (10) is an area-weighted average in each grid cell, Nint, of the effective

pressure in the distributed and weakly connected systems,

NintNd 1 fw

Nwfw; 14 This simplication is justied by the fact that for the effective pressure and parameter values in the simulations, equation (10) remains primarily in the Coulomb friction law regime (equation (11)) where the use of equation (14) yields the appropriate description of basal traction, equation (13). (This approximation becomes increasingly inaccurate as the basal friction law transitions into the form of equation (12) where tb is not directly proportional to N).

Two standalone simulations of ice velocity are performed, both forced by effective pressure generated by the summer subglacial hydrology simulations. In the rst simulation, equation (14) is calculated from Nd and Nw calculated in the summer subglacial hydrology simulation. This simulation demonstrates the impact of weakly connected drainage on ice dynamics (Fig. 3b). In the second simulation, the Nw eld is held steady while Nd comes from the summer subglacial hydrology simulations (Fig. 3c). This is a control simulation that conrms the lack of seasonal hysteresis in the effective pressurevelocity relationship when the evolution of weakly connected drainage is not included in the model. Comparison of the velocity versus pressure relationship from these two simulations eliminates the possibility that the observed seasonal hysteresis is due to changing conditions in the distributed and/or channelized systems.

In both ice velocity simulations, the temperature-dependent rate factor, A, is a function of ice temperature40, and the vertical ice temperature prole in the model is taken equal to that measured at our study site14,22 using 11 uniform vertical levels. Based on ice deformation calculations from Ryser et al.14 we apply an enhancement factor to A of 2.7 within the deepest modelled layer. Ice geometry

is held steady for the duration of these summer simulations; ice velocity is simply calculated diagnostically at each time step based on the xed geometry and changing basal boundary condition, which is a function of our modelled subglacial hydrological evolution. As for the subglacial hydrology model, the domain is periodic in the y-direction. Because we are only concerned with the ice speed at the study site, the model domain for the ice dynamics calculations is subset to span the area 10 km upstream and downstream of the study site to make the calculations less expensive. At the upstream and downstream boundaries we apply a vertically uniform Dirichlet boundary condition on the velocity eld of 100 m a 1. We conrm that the velocity solution at the study site is independent of the choice of boundary condition value (as expected when the boundaries are far enough from the study site, 4B410 ice thicknesses65).

Parameters in equation (10) (Supplementary Table 1) are tuned to to yield model velocity of the observed magnitude, and the same values are used for both model versions. Specically, t b and C are varied in combination to approximately match the observed diurnal range of surface speeds. The ratio lm is tuned to achieve the observed variation in the sensitivity of ice speed to effective pressure. This is assessed by an analysis of the slopes of the lines dening the minimum and maximum channel hydraulic head and ice surface speed on each day in Fig. 3. This slope represents the relationship between water pressure and sliding within a single day. A key feature in the observations (Fig. 3a) is a tendency for lines restricted to lower hydraulic head values to have atter slopes than those at higher hydraulic head values (lines to the left tend to be more horizontal). This behaviour is the result of a reduction in sliding sensitivity to water pressure at low water pressures (high effective pressures). It is predicted by theory and is a key feature of the basal friction law used (equation (10)). At large effective pressure when subglacial cavities are small, sliding is controlled by regelation and enhanced creep3840. This insensitivity of velocity to effective pressure at large effective pressures has been observed in mountain glaciers41 and can be clearly seen in high temporal resolution data at our study site (see ref. 4; Extended Data Fig. 4b). To ensure the models are calibrated to correctly represent this changing sensitivity of sliding to effective pressure, we calculate a linear regression between the minimum channel hydraulic head and the slope of the channel hydraulic head/ice speed relationship on each day. Restricting the regression to the range of hydraulic head minimum values in common between all three data sets (o125 m), we conrm that the observations and both model versions have a similar sensitivity to changes in effective pressure

t b C

ub ub N

nL

NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13903

for the parameter values used (Supplementary Fig. 1). The differing slopes of the lines in Fig. 3c represent this variable sensitivity of sliding to effective pressure and are what cause the modest seasonal-scale changes exhibited by the model with a static weakly connected system.

Model sensitivity to weakly connected area fraction. While a complete analysis of the sensitivity of all parameters used in the models is beyond the scope of this study, we assess how the studys main conclusions are affected by the choice of key parameter, fw, the area fraction of the weakly connected system, and the possibility that it could change in time.

In the main text, we present results for fw 0.67. Because a total of six boreholes

at two different drill sites exhibited out-of-phase behaviour for the majority of both summers measured4,20, we assume the fraction must be large, but we want to avoid being overly prescriptive with the new component. In fact, an upper bound for fw can be found based on the observations4,20 that hydraulic head in the moulin system has typical diurnal variations of about 20% of overburden pressure, while diurnal variations in the boreholes are about 2%. Because the observed surface velocity is precisely in-phase with the moulin pressure variations4 that are driving the system, we can assume that the area-weighted diurnal variations in the distributed system connected to the moulin must be larger than the area-weighted variations in the weakly connected system:

0:20 1 fw

40:02fw; 15

This provides an upper bound of fwo0.91. Note that if diurnal variations in the moulin are in actuality larger than those across representative regions of the distributed system connected to it (as might be expected for a diffusive pressure wave32), then the upper bound for fw should be smaller than the calculated value.

As anticipated from equation (15), if we prescribe fw 0.90 within the model

(and retuning the parameters in the basal friction law), diurnal variations in ice speed are almost entirely absent (Supplementary Fig. 5). This is because the area-weighted diurnal variations of the weakly connected system roughly cancel out the area-weighted diurnal variations in the distributed system because the two systems are out of phase from each other. The few days with substantial diurnal range in ice speed seen in these model results occur when conditions within the model temporarily shift the phasing of diurnal variations in the weakly connected system. Note that ice speed still drops over the summer when the weakly connected system is allowed to drain (Supplementary Fig. 5b), but overall the results differ markedly from the observations and are largely unphysical.

The lower bound on fw is less constrained than its upper bound, but based on the weakly connected behaviour of all boreholes in the study area, we assumefw 0.50 forms a reasonable lower constraint. Model results with that value

(Supplementary Fig. 6) are somewhat similar to the baseline value of fw 0.67 used

in the main text (Fig. 3), but the ice speed does not drop as substantially; with the weakly connected system covering a smaller fraction of the bed, changes there have less impact on the integrated basal traction. Noting that the modelled changes in the weakly connected system have been constrained by the borehole observations, we nd that fw values of B0.600.80 can give results that provide a reasonable match to the ice speed measurements, allowing for modest adjustments to the other parameters in the model.

In addition to assessing the baseline value of fw, we also consider the possibility that fw could change during the summer. Certainly, rather than existing areas of weakly connected drainage becoming more strongly connected to the active drainage system, a reasonable hypothesis would be a change in the area fraction of weakly connected regions. We test this by an additional model run where evolution within the weakly connected system occurs by fw declining (Supplementary Fig. 2)

or increasing (Supplementary Fig. 3) rather than changes to the permeability. These parameterizations do a worse job at reproducing the observed behaviour. Of these two runs, the situation where fw declines is the more observationally supported changeobservations on mountain glaciers have shown that isolated boreholes can become connected during periods of high water pressure in the active drainage system34,35. However, in the run where fw declines, little seasonal evolution in the ice speed occurs (Supplementary Fig. 2b). Instead, the diurnal range in ice speed increases as the summer progresses, an effect that is not seen in the observations. Similarly, the primary effect of increasing fw during summer (Supplementary

Fig. 3b) is a decrease in the diurnal range of ice speed. While we cannot rule out the possibility that the area fraction of the weakly connected system changes during summer, our model results indicate it is not the primary mechanism causing ice speed to drop.

Data and code availability. Model code is development code based off of the Community Ice Sheet Model (CISM) version 2.0.4 (http://oceans11.lanl.gov/cism/

Web End =http://oceans11.lanl.gov/cism/). Model code, processing scripts, input datasets and model output are all available on request from the corresponding author.

References

1. Zwally, H. J. et al. Surface melt-induced acceleration of Greenland ice-sheet ow. Science 297, 218222 (2002).

2. Bartholomew, I. et al. Seasonal evolution of subglacial drainage and acceleration in a Greenland outlet glacier. Nat. Geosci. 3, 408411 (2010).

3. Hoffman, M. J., Catania, G. A., Neumann, T. A., Andrews, L. C. & Rumrill, J. A. Links between acceleration, melting, and supraglacial lake drainage of the western Greenland Ice Sheet. J. Geophys. Res. 116, F04035 (2011).

4. Andrews, L. C. et al. Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet. Nature 514, 8083 (2014).

5. Schoof, C. Ice-sheet acceleration driven by melt supply variability. Nature 468, 803806 (2010).

6. Sole, A. et al. Winter motion mediates dynamic response of the Greenland Ice Sheet to warmer summers. Geophys. Res. Lett. 40, 39403944 (2013).

7. Iken, A. The effect of the subglacial water pressure on the sliding velocity of a glacier in an idealized numerical model. J. Glaciol. 27, 407421 (1981).

8. Walder, J. S. Hydraulics of subglacial cavities. J. Glaciol. 32, 439445 (1986).9. Kamb, B. Glacier surge mechanism based on linked cavity conguration of the basal water conduit system. J. Geophys. Res. 92, 90839100 (1987).

10. Fountain, A. G. & Walder, J. S. Water ow through temperate glaciers. Rev. Geophys. 36, 299328 (1998).

11. Flowers, G. E. Modelling water ow under glaciers and ice sheets. Proc. R. Soc. A 471, 20140907 (2015).

12. Hewitt, I. J. Modelling distributed and channelized subglacial drainage: the spacing of channels. J. Glaciol. 57, 302314 (2011).

13. Hoffman, M. & Price, S. Feedbacks between coupled subglacial hydrology and glacier dynamics. J. Geophys. Res. Earth Surf. 119, 123 (2014).

14. Ryser, C. et al. Sustained high basal motion of the Greenland ice sheet revealed by borehole deformation. J. Glaciol. 60, 647660 (2014).

15. van de Wal, R. S. W. et al. Self-regulation of ice ow varies across the ablation area in south-west Greenland. Cryosphere 9, 603611 (2015).

16. Dow, C. F., Kulessa, B., Rutt, I., Doyle, S. & Hubbard, A. Upper bounds on subglacial channel development for interior regions of the Greenland ice sheet.J. Glaciol. 60, 10441052 (2014).17. Meierbachtol, T., Harper, J. & Humphrey, N. Basal drainage system response to increasing surface melt on the Greenland Ice Sheet. Science 341, 777779 (2013).

18. Wright, P. J., Harper, J. T., Humphrey, N. F. & Meierbachtol, T. W. Measured basal water pressure variability of the western Greenland Ice Sheet: implications for hydraulic potential. J. Geophys. Res. Earth Surf. 121, 114 (2016).

19. Sundal, A. V. et al. Melt-induced speed-up of Greenland ice sheet offset by efcient subglacial drainage. Nature 469, 521524 (2011).

20. Ryser, C. et al. Caterpillar-like ice motion in the ablation zone of the Greenland ice sheet. J. Geophys. Res. Earth Surf. 119, 114 (2014).

21. Walter, F., Chaput, J. & Luthi, M. P. Thick sediments beneath Greenlands ablation zone and their potential role in future ice sheet dynamics. Geology 42, 487490 (2014).

22. Lthi, M. P. et al. Heat sources within the Greenland Ice Sheet: dissipation, temperate paleo-rn and cryo-hydrologic warming. Cryosphere 9, 245253 (2015).

23. Roeoesli, C., Helmstetter, A., Walter, F. & Kissling, E. Meltwater inuences on deep stick-slip icequakes near the base of the Greenland Ice Sheet. J. Geophys. Res. Earth Surf. 121, 118 (2016).

24. Price, S. et al. Community Ice Sheet Model (CISM) v2.0.5 Documentation. Tech. Rep., Los Alamos National Laboratory (2015).

25. Pimentel, S. & Flowers, G. E. A numerical study of hydrologically driven glacier dynamics and subglacial ooding. Proc. R. Soc. A 467, 537558 (2010).

26. Hewitt, I. J. Seasonal changes in ice sheet motion due to melt water lubrication. Earth Planet. Sci. Lett. 371372, 1625 (2013).

27. Dow, C. F. et al. Modeling of subglacial hydrological development following rapid supraglacial lake drainage. J. Geophys. Res. Earth Surf. 120, 11271147 (2015).

28. Gulley, J. et al. The effect of discrete recharge by moulins and heterogeneity in ow-path efciency at glacier beds on subglacial hydrology. J. Glaciol. 58, 926940 (2012).

29. Das, S. B. et al. Fracture propagation to the base of the Greenland Ice Sheet during supraglacial lake drainage. Science 320, 778781 (2008).

30. Tedesco, M. et al. Ice dynamic response to two modes of surface lake drainage on the Greenland Ice Sheet. Environ, Res. Lett. 8, 034007 (9pp) (2013).

31. Catania, G. A. & Neumann, T. A. Persistent englacial drainage features in the Greenland Ice Sheet. Geophys. Res. Lett. 37, 15 (2010).

32. Hubbard, B., Sharp, M., Willis, I., Nielsen, M. & Smart, C. Borehole water-level variations and the structure of the subglacial hydrological system of Haut Glacier dArolla, Valais, Switzerland. J. Glaciol. 41, 572583 (1995).

33. Werder, M. A., Hewitt, I. J., Schoof, C. G. & Flowers, G. E. Modeling channelized and distributed subglacial drainage in two dimensions. J. Geophys. Res. Earth Surf. 118 21699011 (2013).

34. Murray, T. & Clarke, G. K. C. Black-box modeling of the subglacial water system. J. Geophys. Res. 100, 1023110245 (1995).

35. Gordon, S. et al. Seasonal reorganization of subglacial drainage inferred from measurements in boreholes. Hydrol. Process. 12, 105133 (1998).

10 NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13903 ARTICLE

36. Schoof, C., Rada, C. A., Wilson, N. J., Flowers, G. E. & Haseloff, M. Oscillatory subglacial drainage in the absence of surface melt. Cryosphere 8, 959976 (2014).

37. Meierbachtol, T., Harper, J., Humphrey, N. & Wright, P. Mechanical forcing on water pressure in a hydrologically isolated reach beneath Western Greenlands ablation zone. Ann. Glaciol. 57, 6270 (2016).

38. Weertman, J. On the sliding of glaciers. J. Glaciol. 3, 3338 (1957).39. Schoof, C. The effect of cavitation on glacier sliding. Proc. R. Soc. A 461, 609627 (2005).

40. Cuffey, K. & Paterson, W. S. B. The Physics of Glaciers 4th edn (Butterworth-Heinneman, 2010).

41. Iken, A. & Bindschadler, R. A. Combined measurements of subglacial water pressure and surface velocity of the Findelengletscher, Switzerland: conclusions about drainage system and sliding mechanism. J. Glaciol. 32, 101119 (1986).

42. Tedstone, A. J. et al. Decadal slowdown of a land-terminating sector of the Greenland Ice Sheet despite warming. Nature 526, 692695 (2015).

43. Iken, A. & Truffer, M. The relationship between subglacial water pressure and velocity of Findelengletscher, Switzerland, during its advance and retreat.J. Glaciol. 43, 328338 (1997).44. Hodge, S. M. Direct measurement of basal water pressures: progress and problems. J. Glaciol. 23, 309319 (1979).

45. Balise, M. J. & Raymond, C. F. Transfer of basal sliding variations to the surface of a linearly viscous glacier. J. Glaciol. 31, 308318 (1985).

46. Gudmundsson, G. H. Transmission of basal variability to a glacier surface.J. Geophys. Res. 108, 119 (2003).47. Armstrong, W., Anderson, R. S., Allen, J. & Rajaram, H. Modeling the WorldView-derived seasonal velocity evolution of Kennicott Glacier, Alaska.J. Glaciol. 62, 763777 (2016).48. Bougamont, M. et al. Sensitive response of the Greenland Ice Sheet to surface melt drainage over a soft bed. Nat. Commun. 5, 5052 (2014).

49. Shapero, D. R., Joughin, I. R., Poinar, K., Morlighem, M. & Gillet-chaulet, F. Basal resistance for three of the largest Greenland outlet glaciers. J. Geophys. Res. Earth Surf. 121, 168180 (2016).

50. Dow, C. et al. Seismic evidence of mechanically weak sediments underlying Russell Glacier, West Greenland. Ann. Glaciol. 54, 135141 (2013).

51. Schoof, C., Hewitt, I. J. & Werder, M. A. Flotation and free surface ow in a model for subglacial drainage. Part 1. Distributed drainage. J. Fluid Mech. 702, 126156 (2012).

52. Fudge, T. J., Humphrey, N. F., Harper, J. T. & Pfeffer, W. T. Diurnal uctuations in borehole water levels: conguration of the drainage system beneath Bench Glacier, Alaska, USA. J. Glaciol. 54, 297306 (2008).

53. Banwell, A. F., Willis, I. C. & Arnold, N. S. Modelling subglacial water routing at Paakitsoq, W Greenland. J. Geophys. Res. Earth Surf. 118, 12821295 (2013).

54. Nye, J. F. The ow of glaciers and ice-sheets as a problem in plasticity. Proc. R. Soc. A 207, 554572 (1951).

55. Morlighem, M., Rignot, E., Mouginot, J., Seroussi, H. & Larour, E. Deeply incised submarine glacial valleys beneath the Greenland ice sheet. Nat. Geosci. 7, 418422 (2014).

56. Bggild, C. E., Olesen, O. B., Ahlstrm, A. P. & Jrgensen, P. Automatic glacier ablation measurements using pressure transducers. J. Glaciol. 50, 303304 (2004).

57. Smith, L. C. et al. Efcient meltwater drainage through supraglacial streams and rivers on the southwest Greenland ice sheet. Proc. Natl Acad. Sci. 112, 10011006 (2015).

58. Nienow, P., Sharp, M. & Willis, I. Seasonal changes in the morphology of the subglacial drainage system, Haut Glacier dArolla, Switzerland. Earth Surf. Process. Landforms 23, 825843 (1998).

59. Clason, C. C. et al. Modelling the transfer of supraglacial meltwater to the bed of Leverett Glacier, Southwest Greenland. Cryosphere 9, 123138 (2015).

60. Herterich, K. in Dynamics of the West Antarctic Ice Sheet, Vol. 4 of Glaciology and Quaternary Geology (eds van der Veen, C. & Oerlemans, J.) 185202 (Springer, 1987).

61. Blatter, H. Velocity and stress-elds in grounded glaciersa simple algorithm for including deviatoric stress gradients. J. Glaciol. 41, 333344 (1995).

62. Pattyn, F. A new three-dimensional higher-order thermomechanical ice sheet model: basic sensitivity, ice stream development, and ice ow across subglacial lakes. J. Geophys. Res. 108, 115 (2003).

63. Dukowicz, J. K., Price, S. F. & Lipscomb, W. H. Consistent approximations and boundary conditions for ice-sheet dynamics from a principle of least action.J. Glaciol. 56, 480496 (2010).64. Gagliardini, O., Cohen, D., Rback, P. & Zwinger, T. Finite-element modeling of subglacial cavities and related friction law. J. Geophys. Res. 112, F02027 (2007).

65. Kamb, B. & Echelmeyer, K. A. Stress-gradient coupling in glacier ow:I. Longitudinal averaging of the inuence of ice thickness and surface slope.J. Glaciol. 32, 267284 (1986).

Acknowledgements

This work was supported by a grant to M.J.H. from the Laboratory Directed Research and Development Early Career Research Program (LDRD-ECR) at Los Alamos National Laboratory, Climate Modeling Programs within the U.S. Department of Energy Ofce of Science, and by the National Science Foundation, under grant ANT-0424589 to the Center for Remote Sensing of Ice Sheets (CReSIS). L.C.A. was supported by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by Universities Space Research Association under contract with NASA, and UTIG Ewing-Worzel and Gale White Graduate Student Fellowships. J.G. was supported by National Science Foundation Division of Earth Sciences (EAR) Postdoctoral Fellowship (No. 0946767). Fieldwork resulting in the presented observations was supported by United States National Science Foundation grants OPP-0908156 and OPP-0909454, Swiss National Science Foundation grant 200021_127197, National Geographic Society grant 9067-12 and NASA Cryospheric Sciences.

Author contributions

M.J.H. developed the model, designed and ran the experiments, and wrote the manuscript. L.C.A. analysed observational data and contributed to development of the conceptual model. S.A.P. consulted on model implementation, application and interpretation. All authors contributed to analysis and interpretation of observational data and participated in preparation of the manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Hoffman, M. J. et al. Greenland subglacial drainage evolution regulated by weakly connected regions of the bed. Nat. Commun. 7, 13903 doi: 10.1038/ncomms13903 (2016).

Publishers note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional afliations.

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2016

NATURE COMMUNICATIONS | 7:13903 | DOI: 10.1038/ncomms13903 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 11