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.
Introduction
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 flow acceleration every summer^{1,2,3,4}. The influx 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 sliding^{5}. However, ice speed subsequently lowers over the summer despite sustained meltwater input^{1,2,3,4,6}, which is generally explained by increasing efficiency of the subglacial drainage system^{2,3,6}.
This seasonal evolution of subglacial drainage beneath the GrIS is currently interpreted in the context of traditional theory of a twocomponent subglacial drainage system consisting of distributed and channelized drainage^{1,2,3,4,6}. Theory suggests that at low subglacial discharge, drainage occurs through inefficient, distributed pathways—such as linked cavities formed in the lee of bedrock bumps as the glacier slides over the bed or pathways eroded into basal sediments—for which increasing water flux leads to increased water pressure and sliding^{7,8,9,10,11}. It is thought that when a critical discharge is reached in the distributed drainage system, dissipation of heat within the water flow causes a positive feedback between melting of the ice roof and cavity growth, leading to the formation of discrete, efficient channels incised into the ice above^{5,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 system^{5,11,12,13}. 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 winter^{1,6,14,15}, which is inconsistent with a timescale of hours to days for channel collapse under the thick GrIS ice^{3,6,16}. Additionally, modelling has indicated that gentle surface slopes on the ice sheet should suppress channel formation^{16,17}, and low water pressures^{18} and depressed summer ice speeds^{19} characteristic of highly efficient 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 Greenland^{4,14,20,21,22,23} (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 efficiency.
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 largescale 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 model^{13,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 regions of the bed (see Methods: weakly connected drainage model). While the former two components have been routinely included in subglacial hydrology system models^{11}, this is the first 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 sliding^{14} (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 boreholes^{4,20}, we assume these patches cover roughly twothirds of the area of the bed (see Methods: model sensitivity to weakly connected area fraction). There is no throughflow 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.
Results
Modelled 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 system^{4}. Our model results show that an efficient 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 flatter slopes and thicker ice have been proposed to prevent channelization^{16,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 floatation to develop within a few days, as found in previous models during large pulses of melt^{25,26}. Previous modelling studies^{16,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 preconditioned^{26,28}. For example, extensive flooding of the bed at the start of the melt season from supraglacial lake drainage^{3,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 efficiency^{13,17}. Alternatively, year to year persistent moulin locations^{31} may facilitate repeated occupation by channels of the same locations and rapid channel growth through cumulative erosion of basal sediments creating preferential flow pathways^{26,28,32}. Finally, the prescription of a single channel in our model rather than an anastomosing network of channels likely hinders our model’s 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 channels^{5,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 first part of summer as the channel grows to equilibrium, supporting our hypothesis that dropping summer borehole pressures are caused by slow, downgradient leakage towards wellconnected portions of the drainage system. After the channel reaches its equilibrium size (day ∼200–220; 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 shortterm variations in the weakly connected system that contrast with behaviour in the wellconnected 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 glaciers^{4,20,34,35,36,37}, has been explained as the transfer of normal stress from hydraulically wellconnected regions of the bed^{20,34,36,37}, which at our site occurs over kilometers^{20}. 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 wellconnected drainage system^{4,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 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 first simulation, but the weakly connected component is static with a prescribed, fixed 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 defining 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 theory^{38,39,40} and observations^{4,41}, and we confirm 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 sliding—the same water pressure induces less sliding later in the season. To quantify this behaviour, we calculate the Pearson product–moment correlation coefficient 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 significant positive correlation with the observations (r=+0.36, P=0.01), while the model with the static weakly connected system is not significantly 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 system^{34,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 find 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} find 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 presummer 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 ∼2.0 years to return to its original water thickness. This estimate is a minimum value because water would continue to drain out as melt refills 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 significant variation, and not all cavities demonstrate significant drainage^{4,20}. Based on these estimates, varying degrees of drainage from the weakly connected system during each summer provide a plausible explanation for the multiyear selfregulation of annual velocity hypothesized by Tedstone et al.^{42}
Conceptual model
While subglacial drainage has generally been categorized into distributed and channelized components^{5,10,11,12}, 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 identified 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 speed^{4,20,34,35,36,37,43}. In many areas, these inactive regions cover large fractions of the bed; Hodge^{44} 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 site^{4,20}, and boreholes only 30 m apart showed no obvious connection and large hydraulic gradients^{20}, 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 confirm 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 wellconnected and poorly connected regions of the bed, and poorly connected regions therefore moderate active drainage regions^{40,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 thicknesses^{45,46,47} (>∼km 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 moulinconnected system increases or inactive patches switch to active behaviour^{4,34,35,36}.
Based on this new understanding, we propose a threecomponent conceptual model for subglacial hydrology that adds a ‘weakly connected system’ to the traditional distributed and channelized forms of subglacial drainage (Fig. 1). We find this third component necessary to model the ice sheet velocity response properly. Meltwater input of sufficient magnitude drives the formation of channels which are the primary control on overall system efficiency. Channels, as inherently linear features, 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 channeldriven 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 system^{43}, 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 orifices or the presence of low hydraulic conductivity till. Indeed, extensive subglacial till has been identified at FOXX^{21}, and stickslip ice motion observed seismically is associated with spatially differing connectivity of till^{23}. More generally, both modelling^{48,49} and observations^{21,50} have argued for the existence of significant amounts of till underlying marginal regions of the GrIS. We find 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 drainage^{10,11,12}.
Discussion
The 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 observations^{4} 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 reequilibration of low pressure conduits after the cessation of melt is hours to days^{3,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 years^{6,15,42}. This points to a slowreacting component of the basal drainage system. The gradual dewatering of the weakly connected system would take longer to recharge than the other systems^{42}, which is confirmed 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 melt^{42}, 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 flow resistance^{21}. 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 Price^{13} 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 (∼10^{2} 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 onedimensional 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 described^{13}, 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, f_{w}, is assumed to be covered by the weakly connected cavity system, with the remaining fraction, 1−f_{w}, covered by the throughflowing 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 filled 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 twothirds 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 FOXX^{4,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 (h_{w}) is described by the balance between locally generated melt (m_{w}) and exchange of water with the sheet (γ_{w}),
where ρ_{w} is the density of water, A_{w} is the area of the weakly connected system within each grid cell, defined as f_{w}ΔxΔy.
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:
where u_{b} is the sliding velocity, A is the temperaturedependent rate factor for ice deformation, N_{w} is the effective pressure in the weakly connected system, and h_{r} and l_{r} 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 flow modelling performed by Ryser et al.^{14} where the basal layer is temperate and an enhancement factor of 2.5–4 (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, m_{w}, also matches that for the distributed system:
where G is the geothermal heat flux, is the basal traction vector and L is the latent heat of fusion of water.
For all components, hydraulic potential is defined as
where the subscript _{{d,c,w}} indicates one of the distributed, channel, or weakly connected systems, respectively. g is the acceleration due to gravity (m^{2} s^{−1}), z_{b} is the bed elevation (m), p is water pressure (Pa) and p_{ice} is the ice overburden pressure (Pa). For the distributed and channel systems the ice overburden pressure is calculated using a hydrostatic assumption, p_{ice}=ρ_{i}gH, where ρ_{i} 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 floatation, as observed, we include a simple representation for normal stress transfer^{34,36} in the weakly connected system by increasing the ice overburden pressure according to
We set r_{w} to a constant value of 0.13, because this corresponds to the maximum value of borehole water pressure observed at our study site^{20}, and thus suggests a limiting value when N_{w}=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 flux, γ_{w}, between the surrounding distributed system and the weakly connected system within each grid cell based on a Darcy flow law:
where and are the hydraulic potential in the weakly connected and distributed systems, respectively, k_{0w} is the permeability between the weakly connected and distributed systems, P_{w} 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 P_{w} and Δs. Specifically, 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 elsewhere^{4,34,35,52}) and set Δs=10 m. For the chosen values of the fractional area of the weakly connected system (f_{w}=0.67) and grid spacing (Δx=Δy=200 m) we can calculate P_{w}≈5,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 k_{0w}, 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 (γ_{c}; 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, k_{0w}, is many times smaller than the permeability within the distributed system itself, k_{0} (by ∼9 orders of magnitude in our model configuration; Supplementary Table 1), defining 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 k_{0w} to increase over the course of the summer by
where is a base value during winter, and k_{rate} is a rate at which the permeability increases beginning at the start of summer, t_{s}. We choose the value of k_{rate} to reproduce observations of borehole water pressure (Fig. 4b). With our parameter choices, k_{0w} increases by a factor of ∼70 × over summer; the end of summer value of k_{0w} remains ∼8 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 FOXX^{4,20} and at mountain glaciers^{34,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 k_{0w} 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 fields, 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 fine grid, while the subgrid parameterization allows the weakly connected system to be represented at a coarser resolution, making this approach transferrable to largescale 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 summer^{4}. Similarly, some studies have observed weakly connected cavities becoming fully ‘connected’ during periods of high water pressure in the surrounding drainage system^{34,35}, which would correspond to temporal changes in our f_{w} 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 (h_{r}, l_{r}) 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 reflect unique measurements from specific boreholes.
Model setup
We generate a simplified 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 subglacial catchment in this region of GrIS (estimated moulin density in this region is 0.2–0.25 km^{−1} (refs 1, 53)). The ice sheet geometry is a flat bed with a ‘plastic’ glacier shape^{51,54} assuming a constant basal shear stress of 10^{5} 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 farfield constraints on the model solution at the study site. There is a zero subglacial water flux boundary condition for the distributed system at x=−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 simplified geometry over the more complex geometry of the field 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 distances^{14} and that ‘active’ and ‘passive’ subglacial regions might be affected by bed topography^{20}. However, the twodimensional basal topography of our study region is only known approximately. Only a handful of flight lines of airborne ice thickness measurements exist in close proximity to our study area, and for those that do, differences in ice thickness at flightline crossover points exceed 100 m in some locations. Additionally, the commonly used BedMachine product^{55}, which uses mass conservation constraints to improve estimates of ice thickness between radar flight lines, does not cover this area. Previous subglacial hydrologic modelling has shown that valleys in the bedrock topography can concentrate subglacial drainage^{33}, so the true bedrock topography is likely to play an important role in defining the largescale 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 stability^{26,28,33}. Because both our observations and simulations are focused on a site at a moulin, we expect the simplified 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 h_{r}, l_{r} and k_{0} (Supplementary Table 1) are chosen to yield a hydraulic head in the distributed system at our study site of about 200 m below local floatation elevation, and is chosen to yield a hydraulic head in the weakly connected system of about 40 m above local floatation elevation.
Model forcing data
The subglacial hydrology model uses time series of two forcing fields, both of which are based on observations in the vicinity of the FOXX drill site. Surface melt input to the subglacial drainage system (ω in Supplementary Equation (5)) represents the volume flux 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 (u_{b} 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 6h average ablation rates^{4} measured using a pressure transducer installed below the surface connected to a surface reservoir^{56}. Because the measured ablation rates do not have the precision necessary to directly generate a timeseries with the required time step of 2h, and the 6h average rates reduce the amplitude of the diurnal cycle, we generate an idealized diurnal cycle by fitting a sine curve to the 6h timeseries where we adjust the amplitude on each day to maintain the measured 6h average. The ablation sensor began to fail after the large melt event on days 229–230, so we estimate the ablation rate for days 231–240 using a linear fit between mean daily air temperature and ablation rate for the 10 days prior to the large melt event (Pearson correlation coefficient r=0.702, P<0.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 (ω). As discrete moulins can result in the generation of kinematic waves within our channel model, this linear distribution of melt input improves model 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 season^{57}, our melt forcing ignores supraglacial and englacial routing and storage processes that affect the timing and magnitude of melt delivery to the bed^{58,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 boreholederived measurements of ice internal deformation from Global Positioning Systemderived measurements of ice surface velocity to calculate basal slip. Four gaps of about a day are filled 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 constant^{14}. 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, is held constant at 10^{5} 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 spunup 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 timeseries). 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 nucleation^{28,33}. The channel initial condition is an area of 10.0 m^{2} at the margin decreasing linearly to zero at 55 km inland, which yields an area of 5.45 m^{2} at the study site. This value was chosen to allow water pressure to drop below floatation 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 moulins^{4} and is broadly consistent with previous modelling efforts.
Hydraulic head observations
We compare model results to hydraulic head measured at site FOXX^{4,20} in 2012, where hydraulic head, d, relates to hydropotential as
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 simplified 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 d_{f}, the elevation corresponding to floatation pressure, as the vertical datum:
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, threedimensional, firstorder Stokes approximation^{60,61,62,63} momentum balance solver in the Community Ice Sheet Model v2.0^{24} as described by Hoffman and Price^{13}. The magnitude of basal traction, τ_{b}, is defined by a physically based basal friction law for sliding over hard beds that allows for cavitation and bounded basal drag^{39,64},
where C is a Coulomb friction constant, λ_{max} and m_{max} are the wavelength (m) and maximum slope, respectively, of the dominant bedrock bumps, and n is the exponent in Glen’s flow law. is an addition to the original formulation added here to make basal traction calculated for our simplified model domain more realistic. Because our domain lacks the rough, heterogeneous topography of the real ice sheet where resistance to flow is likely to be concentrated in our study area^{20}, we use to represent these missing ‘sticky’ spots^{40}. We apply a constant value of 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 to the lefthand side to clarify the form)
Alternatively, at large effective pressures (low water pressure) the friction law takes a power law form
The basal traction appropriate for ice dynamics is the integrated basal traction of both connected and isolated regions of the bed^{43},
where and 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 areaweighted average in each grid cell, N_{int}, of the effective pressure in the distributed and weakly connected systems,
This simplification is justified 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 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 first simulation, equation (14) is calculated from N_{d} and N_{w} 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 N_{w} field is held steady while N_{d} comes from the summer subglacial hydrology simulations (Fig. 3c). This is a control simulation that confirms the lack of seasonal hysteresis in the effective pressure–velocity 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 temperaturedependent rate factor, A, is a function of ice temperature^{40}, and the vertical ice temperature profile in the model is taken equal to that measured at our study site^{14,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 fixed 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 ydirection. 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 field of 100 m a^{−1}. We confirm 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, >∼4–10 ice thicknesses^{65}).
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. Specifically, and C are varied in combination to approximately match the observed diurnal range of surface speeds. The ratio 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 defining 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 flatter 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 creep^{38,39,40}. This insensitivity of velocity to effective pressure at large effective pressures has been observed in mountain glaciers^{41} 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 (<125 m), we confirm that the observations and both model versions have a similar sensitivity to changes in effective pressure 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 seasonalscale 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 study’s main conclusions are affected by the choice of key parameter, f_{w}, 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 f_{w}=0.67. Because a total of six boreholes at two different drill sites exhibited outofphase behaviour for the majority of both summers measured^{4,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 f_{w} can be found based on the observations^{4,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 inphase with the moulin pressure variations^{4} that are driving the system, we can assume that the areaweighted diurnal variations in the distributed system connected to the moulin must be larger than the areaweighted variations in the weakly connected system:
This provides an upper bound of f_{w}<0.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 wave^{32}), then the upper bound for f_{w} should be smaller than the calculated value.
As anticipated from equation (15), if we prescribe f_{w}=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 areaweighted diurnal variations of the weakly connected system roughly cancel out the areaweighted 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 f_{w} is less constrained than its upper bound, but based on the weakly connected behaviour of all boreholes in the study area, we assume f_{w}=0.50 forms a reasonable lower constraint. Model results with that value (Supplementary Fig. 6) are somewhat similar to the baseline value of f_{w}=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 find that f_{w} values of ∼0.60–0.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 f_{w}, we also consider the possibility that f_{w} 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 f_{w} 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 f_{w} declines is the more observationally supported change—observations on mountain glaciers have shown that ‘isolated’ boreholes can become connected during periods of high water pressure in the active drainage system^{34,35}. However, in the run where f_{w} 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 f_{w} 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/). Model code, processing scripts, input datasets and model output are all available on request from the corresponding author.
Additional information
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).
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Change history
07 February 2017
The original version of this Article contained a typographical error in the spelling of the author Stephen F. Price, which was incorrectly given as Stephen A. Price. This has now been corrected in both the PDF and HTML versions of the Article.
References
 1
Zwally, H. J. et al. Surface meltinduced acceleration of Greenland icesheet flow. Science 297, 218–222 (2002).
 2
Bartholomew, I. et al. Seasonal evolution of subglacial drainage and acceleration in a Greenland outlet glacier. Nat. Geosci. 3, 408–411 (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, 80–83 (2014).
 5
Schoof, C. Icesheet acceleration driven by melt supply variability. Nature 468, 803–806 (2010).
 6
Sole, A. et al. Winter motion mediates dynamic response of the Greenland Ice Sheet to warmer summers. Geophys. Res. Lett. 40, 3940–3944 (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, 407–421 (1981).
 8
Walder, J. S. Hydraulics of subglacial cavities. J. Glaciol. 32, 439–445 (1986).
 9
Kamb, B. Glacier surge mechanism based on linked cavity configuration of the basal water conduit system. J. Geophys. Res. 92, 9083–9100 (1987).
 10
Fountain, A. G. & Walder, J. S. Water flow through temperate glaciers. Rev. Geophys. 36, 299–328 (1998).
 11
Flowers, G. E. Modelling water flow 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, 302–314 (2011).
 13
Hoffman, M. & Price, S. Feedbacks between coupled subglacial hydrology and glacier dynamics. J. Geophys. Res. Earth Surf. 119, 1–23 (2014).
 14
Ryser, C. et al. Sustained high basal motion of the Greenland ice sheet revealed by borehole deformation. J. Glaciol. 60, 647–660 (2014).
 15
van de Wal, R. S. W. et al. Selfregulation of ice flow varies across the ablation area in southwest Greenland. Cryosphere 9, 603–611 (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, 1044–1052 (2014).
 17
Meierbachtol, T., Harper, J. & Humphrey, N. Basal drainage system response to increasing surface melt on the Greenland Ice Sheet. Science 341, 777–779 (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, 1–14 (2016).
 19
Sundal, A. V. et al. Meltinduced speedup of Greenland ice sheet offset by efficient subglacial drainage. Nature 469, 521–524 (2011).
 20
Ryser, C. et al. Caterpillarlike ice motion in the ablation zone of the Greenland ice sheet. J. Geophys. Res. Earth Surf. 119, 1–14 (2014).
 21
Walter, F., Chaput, J. & Luthi, M. P. Thick sediments beneath Greenland’s ablation zone and their potential role in future ice sheet dynamics. Geology 42, 487–490 (2014).
 22
Lüthi, M. P. et al. Heat sources within the Greenland Ice Sheet: dissipation, temperate paleofirn and cryohydrologic warming. Cryosphere 9, 245–253 (2015).
 23
Roeoesli, C., Helmstetter, A., Walter, F. & Kissling, E. Meltwater influences on deep stickslip icequakes near the base of the Greenland Ice Sheet. J. Geophys. Res. Earth Surf. 121, 1–18 (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 flooding. Proc. R. Soc. A 467, 537–558 (2010).
 26
Hewitt, I. J. Seasonal changes in ice sheet motion due to melt water lubrication. Earth Planet. Sci. Lett. 371–372, 16–25 (2013).
 27
Dow, C. F. et al. Modeling of subglacial hydrological development following rapid supraglacial lake drainage. J. Geophys. Res. Earth Surf. 120, 1127–1147 (2015).
 28
Gulley, J. et al. The effect of discrete recharge by moulins and heterogeneity in flowpath efficiency at glacier beds on subglacial hydrology. J. Glaciol. 58, 926–940 (2012).
 29
Das, S. B. et al. Fracture propagation to the base of the Greenland Ice Sheet during supraglacial lake drainage. Science 320, 778–781 (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, 1–5 (2010).
 32
Hubbard, B., Sharp, M., Willis, I., Nielsen, M. & Smart, C. Borehole waterlevel variations and the structure of the subglacial hydrological system of Haut Glacier d’Arolla, Valais, Switzerland. J. Glaciol. 41, 572–583 (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, 2169–9011 (2013).
 34
Murray, T. & Clarke, G. K. C. Blackbox modeling of the subglacial water system. J. Geophys. Res. 100, 10231–10245 (1995).
 35
Gordon, S. et al. Seasonal reorganization of subglacial drainage inferred from measurements in boreholes. Hydrol. Process. 12, 105–133 (1998).
 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, 959–976 (2014).
 37
Meierbachtol, T., Harper, J., Humphrey, N. & Wright, P. Mechanical forcing on water pressure in a hydrologically isolated reach beneath Western Greenland’s ablation zone. Ann. Glaciol. 57, 62–70 (2016).
 38
Weertman, J. On the sliding of glaciers. J. Glaciol. 3, 33–38 (1957).
 39
Schoof, C. The effect of cavitation on glacier sliding. Proc. R. Soc. A 461, 609–627 (2005).
 40
Cuffey, K. & Paterson, W. S. B. The Physics of Glaciers 4th edn ButterworthHeinneman (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, 101–119 (1986).
 42
Tedstone, A. J. et al. Decadal slowdown of a landterminating sector of the Greenland Ice Sheet despite warming. Nature 526, 692–695 (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, 328–338 (1997).
 44
Hodge, S. M. Direct measurement of basal water pressures: progress and problems. J. Glaciol. 23, 309–319 (1979).
 45
Balise, M. J. & Raymond, C. F. Transfer of basal sliding variations to the surface of a linearly viscous glacier. J. Glaciol. 31, 308–318 (1985).
 46
Gudmundsson, G. H. Transmission of basal variability to a glacier surface. J. Geophys. Res. 108, 1–19 (2003).
 47
Armstrong, W., Anderson, R. S., Allen, J. & Rajaram, H. Modeling the WorldViewderived seasonal velocity evolution of Kennicott Glacier, Alaska. J. Glaciol. 62, 763–777 (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. & Gilletchaulet, F. Basal resistance for three of the largest Greenland outlet glaciers. J. Geophys. Res. Earth Surf. 121, 168–180 (2016).
 50
Dow, C. et al. Seismic evidence of mechanically weak sediments underlying Russell Glacier, West Greenland. Ann. Glaciol. 54, 135–141 (2013).
 51
Schoof, C., Hewitt, I. J. & Werder, M. A. Flotation and free surface flow in a model for subglacial drainage. Part 1. Distributed drainage. J. Fluid Mech. 702, 126–156 (2012).
 52
Fudge, T. J., Humphrey, N. F., Harper, J. T. & Pfeffer, W. T. Diurnal fluctuations in borehole water levels: configuration of the drainage system beneath Bench Glacier, Alaska, USA. J. Glaciol. 54, 297–306 (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, 1282–1295 (2013).
 54
Nye, J. F. The flow of glaciers and icesheets as a problem in plasticity. Proc. R. Soc. A 207, 554–572 (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, 418–422 (2014).
 56
Bøggild, C. E., Olesen, O. B., Ahlstrøm, A. P. & Jòrgensen, P. Automatic glacier ablation measurements using pressure transducers. J. Glaciol. 50, 303–304 (2004).
 57
Smith, L. C. et al. Efficient meltwater drainage through supraglacial streams and rivers on the southwest Greenland ice sheet. Proc. Natl Acad. Sci. 112, 1001–1006 (2015).
 58
Nienow, P., Sharp, M. & Willis, I. Seasonal changes in the morphology of the subglacial drainage system, Haut Glacier d’Arolla, Switzerland. Earth Surf. Process. Landforms 23, 825–843 (1998).
 59
Clason, C. C. et al. Modelling the transfer of supraglacial meltwater to the bed of Leverett Glacier, Southwest Greenland. Cryosphere 9, 123–138 (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.) 185–202 (Springer, 1987).
 61
Blatter, H. Velocity and stressfields in grounded glaciers—a simple algorithm for including deviatoric stress gradients. J. Glaciol. 41, 333–344 (1995).
 62
Pattyn, F. A new threedimensional higherorder thermomechanical ice sheet model: basic sensitivity, ice stream development, and ice flow across subglacial lakes. J. Geophys. Res. 108, 1–15 (2003).
 63
Dukowicz, J. K., Price, S. F. & Lipscomb, W. H. Consistent approximations and boundary conditions for icesheet dynamics from a principle of least action. J. Glaciol. 56, 480–496 (2010).
 64
Gagliardini, O., Cohen, D., Råback, P. & Zwinger, T. Finiteelement modeling of subglacial cavities and related friction law. J. Geophys. Res. 112, F02027 (2007).
 65
Kamb, B. & Echelmeyer, K. A. Stressgradient coupling in glacier flow: I. Longitudinal averaging of the influence of ice thickness and surface slope. J. Glaciol. 32, 267–284 (1986).
Acknowledgements
This work was supported by a grant to M.J.H. from the Laboratory Directed Research and Development Early Career Research Program (LDRDECR) at Los Alamos National Laboratory, Climate Modeling Programs within the U.S. Department of Energy Office of Science, and by the National Science Foundation, under grant ANT0424589 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 EwingWorzel 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 OPP0908156 and OPP0909454, Swiss National Science Foundation grant 200021_127197, National Geographic Society grant 906712 and NASA Cryospheric Sciences.
Author information
Affiliations
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.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing financial interests.
Supplementary information
Supplementary Information
Supplementary Figures, Supplementary Tables, Supplementary Methods and Supplementary References. (PDF 1026 kb)
Rights and permissions
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 article’s 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/
About this article
Cite this article
Hoffman, M., Andrews, L., Price, S. et al. Greenland subglacial drainage evolution regulated by weakly connected regions of the bed. Nat Commun 7, 13903 (2016). https://doi.org/10.1038/ncomms13903
Received:
Accepted:
Published:
Further reading

Hydraulic transmissivity inferred from icesheet relaxation following Greenland supraglacial lake drainages
Nature Communications (2021)

Rapid and sensitive response of Greenland’s groundwater system to ice sheet change
Nature Geoscience (2021)

Dynamic response of the Greenland ice sheet to recent cooling
Scientific Reports (2020)

Methane beneath Greenland’s ice sheet is being released
Nature (2019)

Cascading lake drainage on the Greenland Ice Sheet triggered by tensile shock and fracture
Nature Communications (2018)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.