# Threshold response to melt drives large-scale bed weakening in Greenland

https://doi.org/10.1038/s41586-022-04927-3

Nathan Maier<fnref n="1" />,², Florent Gimbert<fnref n="1" /> & Fabien Gillet-Chaulet<fnref n="1" />

Received: 25 October 2021

Accepted: 1 June 2022

Published online: 27 July 2022

<figure><img src="image_1.png" /></figure>

Check for updates

Ice speeds in Greenland are largely set by basal motion<fnref n="1" />, which is modulated by meltwater delivery to the ice base²⁻⁴. Evidence suggests that increasing melt rates enhance the subglacial drainage network's capacity to evacuate basal water, increasing bed friction and causing the ice to slow⁵⁻¹⁰. This limits the potential of melt forcing to increase mass loss as temperatures increase¹¹. Here we show that melt forcing has a pronounced influence on dynamics, but factors besides melt rates primarily control its impact. Using a method to examine friction variability across the entirety of western Greenland, we show that the main impact of melt forcing is an abrupt north-to-south change in bed strength that cannot be explained by changes in melt production. The southern ablation zone is weakened by 20–40 per cent compared with regions with no melt, whereas in northern Greenland the ablation zone is strengthened. We show that the weakening is consistent with persistent basal water storage and that the threshold is linked to differences in sliding and hydropotential gradients, which exert primary control on the pressures within drainage pathways that dewater the bed. These characteristics are mainly set by whether a margin is land or marine terminating, suggesting that dynamic changes that increase mass loss are likely to occur in northern Greenland as temperatures increase. Our results point to physical representations of these findings that will improve simulated ice-sheet evolution at centennial scales.

Meltwater delivered to the bed of glaciers impacts water pressure and storage at the ice base²⁻⁴,¹², which modifies the basal friction by weakening or strengthening the coupling between the ice and its bed⁹,¹³,¹⁴. Current evidence from the Greenland Ice Sheet—the world’s second-largest ice sheet—suggests that increased meltwater delivery to the bed increases the prevalence of channelized subglacial drainage pathways⁵⁻⁹ that efficiently transport meltwater beneath the ice sheet¹²,¹⁵. Subglacial water flow through channels occurs at pressures lower than in the other two proposed drainage components¹²—the distributed drainage network composed of interconnecting cavities, and the semi-isolated network composed of either cavities or porous sediment that is only weakly connected to active drainage pathways⁹,¹². As summers get hotter and longer¹⁶, it is thought that the increased occurrence and extent of channels will draw additional water out of the weakly connected drainage system, which is spatially expansive and is estimated to cover between 60% and 90% of the bed⁹,¹²,¹⁷. This strengthens the ice–bed coupling and decreases ice speeds, resulting in a reduction of direct ice discharge into the ocean and dynamic thickening that decreases surface melting¹¹. This response leads to a negative feedback between increased melting and mass loss.

The framework for how warmer temperatures will influence dynamics into the future is mainly based on decadal-scale observations that show an inverse relationship between melt production and flow speeds⁵⁻⁷,¹⁰. However, the scope of these observations, which all concentrate on the same land-terminating sector, is limited, and the high variability of the velocity time series makes it difficult to directly infer the processes that

control ice speeds⁵⁻⁷,¹⁰,¹⁸. Thus, these observations alone are insufficient to assess the impact of melt across the entire ice sheet that spans a vast range of melt conditions and settings. This is particularly the case for marine-terminating regions¹⁹, which are significantly different in terms of dynamics, and inland areas where channelization is inferred to have a small or inconsistent role in seasonality²⁰. Correspondingly, other melt responses, such as threshold changes, could have a major role in marginal dynamics and yet remain unobserved across the limited temporal and spatial scales examined⁵. Improved approaches that can reveal the time-integrated impact of melt across expansive areas of the ice sheet are necessary to fully assess the dynamic response to melt in Greenland. At present, because of the assumed negative feedback and complexity of implementing subglacial hydrology into centennial-scale models²¹, the impact of melt on dynamics is often not considered in upper-limit sea-level projections²².

Here we assess the long-term impact of melt forcing at the ice-sheet scale by analysing basal-friction variations of seven ice catchments that span the entirety of the western margin of Greenland (Fig. 1). To isolate friction variations that reflect changes in hydrology, we normalize the friction field estimated using a state-of-the-art ice-sheet model (Methods) by the mean rate-strengthening friction relationship derived from high-elevation regions where no summer melting occurs (Fig. 1i and Methods). This normalization removes the interdependency between sliding speed and friction, and allows us to focus on friction variations resulting from other controls²³. As the friction relationship is derived from regions with no summer melt, changes in bed strength where melt

<footnote n="1">¹Université Grenoble Alpes, CNRS, IGE, Grenoble, France. ²e-mail: ntmaier@gmail.com</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1 | Summary of bed-strength variations across western Greenland.**
a–g, Along-flow summary of normalized friction, mean JJA SMB (2001–2015), and flotation fraction for the seven catchments shown in h. For each catchment, the left axis shows the mean normalized friction for 20 evenly spaced elevation bins (vertical bars indicate the 95% confidence interval). The colour of the square marker corresponds to the catchment. The colour of the outer marker denotes the ice flotation fraction in relation to sea level. The right axes (red axis and red-dashed line) show JJA SMB (in metres water equivalent (m.w.e.)) for each elevation bin. The vertical dashed lines show the mean snowline elevation⁴³. h, Normalized friction mapped across western Greenland. The coloured</figcaption></figure>

forcing has an active role in dynamics will be clearly identifiable. The
normalized friction is reported as the percentage difference from the
mean rate-strengthening relationship (Methods), where negative and

<figure><img src="image_3.png" /></figure>

boundaries demarcate the seven drainage catchments⁴⁴ examined. The black line
shows the multiyear averaged snowline (2001–2017)⁴³, and the maroon line shows
where the mean JJA SMB is zero (2001–2015). The yellow dots show the locations of
marine-terminating outlet glaciers. i, Rate-strengthening relationship (black line)
from which the normalized friction was derived (Methods). The magenta dots
show the binned friction through the velocity range. The interquartile range and
90th percentile for each bin are indicated with vertical black lines and black dots,
respectively. j, Multiyear velocity mosaic of western Greenland²⁴. Overlays are the
same as in h. All panels were produced using MATLAB; the background image on
the maps is from a land-ice classification⁴⁵,⁴⁶.

positive values indicate relatively weak and strong beds, respectively.
Unlike velocity time series, which show high interannual variability ⁵,⁶,¹⁰,¹⁸
and can have seasonal biases ⁶,¹⁰,¹⁸, our approach uses velocities averaged
over 20 years²⁴, which includes both winter and summer flow, and thus the normalized friction field reflects the time-averaged effect of melt forcing.

The current model of hydrology-dynamic coupling implies differ-
ences in bed strength should arise where surface melt is significantly
different⁸. With the greater than 2,000-km north-south span of the
study reach, which results in an approximately fourfold difference
in summer melt along the margins (Fig. 1), we can examine how bed
strength varies between regions with higher and lower melt rates
across western Greenland where friction variations are expected to
adhere to the following model. At the transition between the firn and
the bare ice, which demarcates the beginning of the ablation zone,
summer melt can begin to access the basal interface either through
moulins²⁵ or more sporadically through lake drainages⁴,²⁶,²⁷. These
high-elevation regions have lower melt rates²⁸,²⁹ and shorter sum-
mers²⁸,²⁹ compared with the lower-elevation regions near the margin.
The development of channelized drainage is unlikely outside of short
time periods¹⁵,²⁰,³⁰,³¹, and melt is likely to be predominantly routed at
high pressure through distributed drainage²⁰, leading to a relatively
lower basal friction (weakened beds). Moving towards the margin,
increasing melt rates²⁸,²⁹, melt season duration²⁸,²⁹ and upgradient
drainage area³² create conditions that are increasingly favourable to
channelized drainage development that can more efficiently dewater
the distributed and isolated regions of the bed⁹,¹². This would result in
a relative friction increase towards the margins (strengthened beds).
Correspondingly, regional differences in channelization and friction
should arise where melt rates are significantly different³³. The large
and continuous extent of the study region also permits us to look at
variations owing to other hydrology-dynamic factors such as hydro-
potential gradients and sliding speed. These aspects are theorized to
have a significant impact on drainage evolution³³,³⁴, but have yet to be
investigated at the ice-sheet scale.

Along-flow friction variations

We begin testing the model outlined above by examining along-flow
changes in bed strength (binned by elevation; Fig. 1a–g) moving from
the interior to the margin of each catchment. Transitioning from the
highest elevations in the firn down into the ablation zone, we find
abrupt bed weakening in catchments 1–3. This occurs near the onset of
negative June, July and August (JJA) surface mass balance (SMB; red lines
Fig. 1), which is near the average snowline elevation (vertical dashed
lines in Fig. 1a–g) in catchments 2 and 3 and at higher elevations within
catchment 1. This suggests that surface melt reaches the basal interface
and increases water pressure and/or storage compared with the up-ice
regions, correspondingly weakening the bed by approximately 20–40%
(Figs. 1 and 2). We examine other interpretations not related to melt
forcing (that is, changes in geology, roughness, deforming till cover,
ceice deformation and geothermal heating) and determine that they are
unlikely to produce the close correspondence between the onset of
weak beds and melting observed in multiple catchments (Methods and
Extended Data Figs. 1 and 2), but acknowledge that these aspects have
limited constraints. In stark contrast, catchments 4–7 do not show large
friction changes near the snowline, except for modest strengthening in
catchment 5, indicating that surface melt delivered to the base at high
elevations in these regions⁴,²⁷ has a limited influence on bed friction.

Continuing through the ablation zone towards the margin, except
for catchment 2, the catchments that show pronounced bed weakening
near the snowline do not have a prominent rebound to higher fric-
tion. This indicates that channelization, and its expected increasing
prevalence towards the margin as melt increases and accumulates,
has a limited strengthening effect and thus net impact on the basal
pressure field and/or quantity of water stored at the bed in these catch-
ments. Catchment 2 shows an abrupt strengthening in the along-flow
friction directly adjacent to the margin (Fig. 1f). This occurs due to

marginal strengthening in the northern part of the catchment, which
is marine terminating, whereas the southern land-terminating portion
of the catchment remains weak-bedded (Fig. 1h). In northern Greenland
(catchments 4–6) where a limited friction response was observed at
high elevations, a small increase in bed strength is observed, particu-
larly near the margin (small humps beginning at about 1,000 m above
sea level). This is consistent within an increase in efficient, channelized
drainage pathways that is expected moving towards the margin.

Finishing the marginal transect, the start of the tidewater region (in the catchments where present) corresponds to pronounced weakening, which occurs as the flotation fraction owing to sea level begins to increase above zero (outer markers, Fig. 1a–g). This weakening is consistent with our understanding of marine-terminating glaciers³⁸—when the flotation fraction is above zero sea level sets the base pressure and pressure gradients, and weakening can occur independent from conditions relating to evolving drainage. Notably, this serves as further evidence that friction variations observed using our methodology are interpretable in terms of changes in the basal hydrology.

Regional friction variations

Examining the along-flow friction variations for each catchment together shows that the bed strength is similar above the ablation zone but sharply deviates at the onset of surface melting (Fig. 2a). This correspondence implies that the differences in ablation-zone bed strength are mainly related to melt forcing. Interestingly, we observe a strong latitudinal gradient in ablation-zone bed strength. The weakest ablation zones occur in the south and the strongest in the north, a pattern that resembles a latitudinal wave in the along-flow friction presented in Fig. 1. Closer examination of latitudinal variations in ablation-zone bed strength shows the weak-bed regions abruptly transition to stronger beds moving from south to north near Jakobshavn (Figs. 1–3), which separates regions of western Greenland that are predominantly land terminating from those that are predominantly marine terminating. This threshold does not correspond to an abrupt change in melt rates but does closely correspond to an abrupt increase in hydropotential gradients and flow speeds, suggesting that hydrology-dynamic characteristics other than melt may have a primary role in the bed-strength transition (Fig. 3 and Extended Data Fig. 3).

Controls on bed weakening

We propose that the weakening observed in the southern catchments is primarily a result of time-persistent (multiyear) water storage at the ice-bed interface. We support this assertion by contextualizing the observations using hard-bed friction laws (Fig. 2b), an approach backed by previous work that found friction relationships in western Greenland mainly show rate strengthening consistent with sliding over non-deforming beds²³ (we also demonstrate that our observations could be consistent with deforming-bed friction laws (Methods and Extended Data Fig. 4)). We fit a rate-strengthening relationship that does not explicitly account for water storage to the elevation bins that occur above the ablation zone where the influence of melt forcing on dynamics is assumed to be minimal (solid line in Fig. 2b; Methods). In doing so, we find that the elevation bins for the ablation zones of the northern catchments continue to lie close or just above this reference scaling, indicating no increase in water storage and cavitation compared with the upgradient regions. By contrast, the elevation bins from the southern catchments that show melt weakening are shifted downward from this scaling. Plotting a constrained pressure-based friction relationship (dotted lines in Fig. 2b; Methods) shows that this shift is consistent with increased water storage and cavitation compared with upgradient regions. Notably, the points shifted to the right form a parallel trend, indicating that a similar cavitation state (that is, cavity size) occurs within the southern catchments such that their
<figure><img src="image_4.png" /><figcaption>**Fig. 2 |** Threshold change in ablation-zone bed strength. a, Along-flow summary of binned normalized friction for all catchments plotted together (colours and bins are the same as in Fig. 1; vertical bars indicate the 95% confidence interval). Elevation is normalized to maximum elevation with negative JJA SMB, which separates marginal regions subject to net summer melt (negative elevations) from inland regions that do not experience net summer melt (positive elevations). Markers with a lighter edge colour denote elevation bins near the margin where flotation fractions are greater than zero. b, Mean basal velocity and friction for elevation bins from each catchment are plotted alongside hard-bed friction laws<sup>47,48</sup>. The circle markers are interior elevation bins that experience no summer melt (from all catchments), the square markers are ablation-zone bins for northern Greenland (catchments 5-7) and the stars are ablation-zone bins for southern Greenland (catchments 1-4). The horizontal and vertical bars indicate the 95% confidence interval.</figcaption></figure>

<figure><img src="image_5.png" /></figure>

Reference scaling (solid line, $r^2 = 0.85$) is constrained using bins with no melt forcing (circles; Methods). The ablation-zone bins for northern Greenland (squares) adhere to this relationship whereas the ablation-zone bins for southern Greenland (stars) are shifted downwards, indicating cavity expansion compared with interior regions. This trend is captured by a cavitation scaling (dashed line, $r^2 = 0.76$) (Methods). A constrained, effective-pressure-based friction law<sup>47</sup> (dotted lines, effective pressure: $N = P_{\text{ice}} - P_{\text{water}}$, where $P$ is pressure and $N_{\text{ref}}$ is an arbitrary effective pressure dependent on $C$ (Methods)) is plotted for various $C N$ assuming that the cavitation scaling is consistent with Iken's bound (Methods). Deviation from the reference scaling (that is, the curvature in dotted lines) is due to decreasing effective pressure, which increases basal cavitation. Bins where flotation fractions arising from sea level are greater than zero are excluded from the plot but shown for each catchment in Extended Data Fig. 5.

dynamics can be captured by a single cavitation scaling (dashed line in Fig. 2b). This finding is similar to that obtained for a well constrained alpine glacier, where differences in winter and summer sliding rates also resulted in two distinct parallel friction scalings inferred to be related to a feedback between sliding, cavity size and cavity drainage<sup>36</sup>. By comparing the reference and cavity scaling parameters we provide an additional estimate of the weakening magnitude, which here is 33%.

We posit that water storage within a spatially expansive weakly connected drainage system is consistent with the widespread bed weakening observed in southern Greenland. During early summer, water can be stored in weakly connected regions when drainage efficiency and effective pressures (ice pressure minus the water pressure) are low, and drain out during late summer as drainage efficiency and effective pressures in active drainage pathways increase<sup>2–4</sup>. We therefore propose that the net bed strength and basal storage is dependent on either the occurrence or the effective pressure within efficient drainage pathways that can drain weakly connected regions during late summer, although we cannot rule out groundwater transport as another possible mechanism to remove water from the base. Modelling using a framework that represents the physics of both cavities and channels indicates little difference in the ability to form channelized drainage pathways along the margin (Methods and Extended Data Fig. 5), yet finds effective pressures within active drainage pathways are approximately twice as high in northern Greenland compared with southern Greenland during summer (Fig. 4a). This arises not because of the change in melt rates as typically envisioned to control the long-term hydrologic response

on the ice sheet<sup>5,6</sup>, but because drainage pressures also have a strong dependency on hydropotential gradients<sup>33</sup> and sliding speed<sup>34</sup>, which can enlarge drainage pathways independent of discharge. Hydropotential gradients are approximately 1.5-times higher and sliding speeds nearly an order-of-magnitude higher moving north of Jakobshavn (Fig. 3). With this context, our observations can be explained using current theory. Higher effective pressures in active drainage pathways in northern Greenland increase the pressure gradients between the active and the weakly connected drainage systems. This increases the rate at which water can be removed from weakly connected regions of the bed. Moreover, because of the sensitivity of effective pressures on sliding, concentrated areas of fast flow in north Greenland act as glacier-wide ‘funnels’ where pressure gradients within active drainage pathways regionally drive water to a main hydrologic axis (Fig. 4a). Although these effects allow the bed to remain strong in northern Greenland, we find the opposite conditions in southern Greenland. Low effective pressures within active drainage systems limit pressure gradients between, and drawdown from, isolated regions of the bed. The lack of concentrated fast-flow areas in the south also inhibits concentrating drainage, potentially explaining why high effective pressures at the ice edge in southern Greenland generally do not correspond to strengthening, contrary to that observed in the north.

Our observation and analysis indicate a direct link between margin type and the net impact of melt forcing. The transition between weak and strong beds corresponds to the switch between predominantly land-terminating and predominantly marine-terminating regions, as
<figure><img src="image_6.png" /><figcaption>**Fig. 3** | Latitudinal variations in bed strength and drainage-system variables. **a–d**, Latitudinal relationship between normalized friction (**a**), JJA SMB (**b**), basal velocity (**c**) and hydropotential gradient (**d**) is shown using 100</figcaption></figure>

evenly spaced latitudinal bins for the upper and lower ablation zone and upgradient regions with no JJA melt. The markers represent the mean of each bin.

we conceptualize in Fig. 4b. The high speeds and steep surface slopes in the north, which keep effective pressures high within active drainage pathways, are necessary to accommodate ice discharge rates at the terminus of marine-outlet glaciers³⁷, which contrasts with land-terminating margins in the south that have no stress imbalance at the terminus. Moreover, the link between hydrology and sliding yields a physically consistent mechanism to explain why fast-flowing regions maintain a similar friction scaling to their upgradient regions even though cavitation is generally expected to occur at high flow speeds³⁸. Cavitation can be suppressed if the effective pressures increase with sliding speed, as we show in Fig. 4a. This implies that the general character of land- and marine-terminating regions is strongly influenced by hydrology, and vice versa. Steep surfaces and narrow ablation zones in northern Greenland are in part maintained by hydrologic strengthening, whereas in southern Greenland, the flat and broad ablation zone reflects a hydrologically weakened bed.

## Implications for future mass changes

The abrupt north–south bed-strength change (Figs. 1–3) and side-by-side friction scalings (Fig. 2b) indicate the presence of a previously unrecognized hydrologic threshold that we link to sliding and hydropotential gradients. This implies a significant change in dynamics is possible, particularly in northern Greenland where the ablation zone has the potential to shift from strong to weak beds. With the insight from our analysis, we can hypothesize how the future evolution and mass loss of the ice sheet will be modified. As air and ocean temperatures increase and outlet glaciers in northern Greenland recede and become land terminating as many are expected to do over the next century¹¹,³⁷,³⁹, decreasing ice speeds and mass flux at the termini will cause the upgradient regions to slow and flatten. This will create conditions where effective pressures within active drainage pathways decrease, decreasing pressure gradients between drainage-system components and drawdown of

water out of weakly connected regions of the bed. This will most likely begin at higher-elevation regions where velocities and hydropotential gradients are at their lowest. The associated acceleration would impact the surface geometry resulting in upgradient thinning, amplifying the SMB–elevation feedback already thought to be critical to mass-loss rates over the coming century¹¹. It also implies an increase in mass flux towards the margin, possibly forcing outlet glaciers to re-advance to the marine margin resulting in another period of ice discharge.

Our results show weakening in the ablation zone has the potential to reach about 30%, far exceeding known hydrology feedbacks linked to changes in melt rates. Decadal-scale velocity variability attributed to melt-induced drainage efficiency changes shows maximum variability on the order of about 20%, which we estimate would correspond to a friction change of only 4–6% (Methods). This suggests that previously hypothesized changes in drainage efficiency linked to melt quantity are second order to drainage changes set by sliding and hydropotential gradients. Importantly, this indicates that subglacial hydrology, which is commonly assumed unimportant to long-term projections¹¹, has a fundamental role in ice-sheet evolution.

Our findings indicate that a three-component drainage system and sliding feedbacks are necessary to model the hydrologic threshold, neither of which are currently incorporated into the most commonly used drainage model⁴⁰. Our observations also closely match behaviours expected with current friction laws and we demonstrate that a physically based parameterization can be used to model the weakening observed by way of an effective pressure threshold that determines the rate-strengthening scaling to apply (Extended Data Fig. 6 and Methods). Owing to the sensitivity of sea-level projections to the sliding-law choice⁴¹,⁴² and the uncertainty with regard to how to represent hydrology-dynamic coupling in long-term projections, these constraints represent a significant step towards understanding and accounting for complex hydrologic processes in ice-sheet models.
<figure><img src="image_7.png" /></figure>

**Fig. 4** | *Controls on bed-strength threshold.* *a*, Modelled effective pressure (N) in active drainage pathways using accumulated discharge derived from JJA SMB (Methods). Patterns are consistent with the north-to-south threshold observed in normalized friction. The grey and maroon contours show the time-average snowline⁴³ and a JJS SMB of zero, respectively. The background image is from a land-ice classification map⁴⁵,⁴⁶ plotted using MATLAB. *b*, Conceptual model detailing the link between margin type and friction variability from the ice edge to the equilibrium line altitude (ELA). Actively draining distributed and channelized regions are lumped together (centre flow path) and are surrounded by a weakly connected drainage system (isolated patches). Late-summer conditions are depicted where efficient drainage

<figure><img src="image_8.png" /></figure>

pathways have developed. In marine-terminating regions, steep
hydropotential gradients and high sliding velocities lead to high effective
pressures within active drainage pathways, resulting in strong pressure
gradients that draw water out of the weakly connected drainage system,
producing strong beds. Where marine glaciers terminate, flow speeds are high,
but hydropotential gradients are low owing to the influence of sea level, which
restricts the ability to draw out water stored in isolated regions. In
land-terminating regions, low sliding and hydropotential gradients result in
low effective pressures within the actively draining pathways. This reduces
pressures gradients between the active and isolated drainage systems, which
limits the drawdown of stored water and produces weak beds.

Online content

Any methods, additional references, Nature Research reporting
summaries, source data, extended data, supplementary information,
acknowledgements, peer review information; details of author con-
tributions and competing interests; and statements of data and code
availability are available at https://doi.org/10.1038/s41586-022-04927-3.

1. Maier, N., Humphrey, N., Harper, J. & Meierbachtol, T. Sliding dominates slow-flowing margin regions, Greenland Ice Sheet. *Sci. Adv.* 5, eaaw5406 (2019).

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., Catania, G. A., Neumann, T., Andrews, L. & Rumrill, J. Links between acceleration, melting, and supraglacial lake drainage of the western Greenland Ice Sheet. *J. Geophys. Res. Earth Surf.* **116**, F04035 (2011).

4. Andrews, L. C. et al. Seasonal evolution of the subglacial hydrologic system modified by supraglacial lake drainage in western Greenland. *J. Geophys. Res. Earth Surf.* **123**, 1479–1496 (2018).

5. Williams, J. J., Gourmelen, N. & Nienow, P. Dynamic response of the Greenland Ice Sheet to recent cooling. *Sci. Rep.* **10**, 1647 (2020).

6. Tedstone, A. J. et al. Decadal slowdown of a land-terminating sector of the Greenland Ice Sheet despite warming. *Nature* **526**, 692–695 (2015).

7. Van de Wal, R. et al. Self-regulation of ice flow varies across the ablation area in south-west Greenland. *Cryosphere* **9**, 603–611 (2015).

8. Davison, B. J., Sole, A. J., Livingstone, S. J., Cowton, T. R. & Nienow, P. W. The influence of hydrology on the dynamics of land-terminating sectors of the Greenland Ice Sheet. *Front. Earth Sci.* **7**, 10 (2019).

9. Hoffman, M. J. et al. Greenland subglacial drainage evolution regulated by weakly connected regions of the bed. *Nat. Commun.* **7**, 13903 (2016).

10. Stevens, L. A. et al. Greenland Ice Sheet flow response to runoff variability. *Geophys. Res. Lett.* **43**, 11295–11303 (2016).

11. Pattyn, F. et al. The Greenland and Antarctic ice sheets under 1.5 °C global warming. *Nat. Clim. Change* **8**, 1053–1061 (2018).

12. Andrews, L. C. et al. Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet. *Nature* **514**, 80–83 (2014).

13. Cowton, T., Nienow, P., Sole, A., Bartholomew, I. & Mair, D. Variability in ice motion at a land-terminating Greenlandic outlet glacier: the role of channelized and distributed drainage systems. *J. Glaciol.* **62**, 451–466 (2016).

14. Bougamont, M. et al. Sensitive response of the Greenland Ice Sheet to surface melt drainage over a soft bed. *Nat. Commun.* **5**, 5052 (2014).

15. Chandler, D. et al. Evolution of the subglacial drainage system beneath the Greenland Ice Sheet revealed by tracers. *Nat. Geosci.* **6**, 195–198 (2013).
16. Fettweis, X. et al. Estimating the Greenland Ice Sheet surface mass balance contribution to future sea level rise using the regional atmospheric climate model MAR. *TCryosphere* **7**, 469–489 (2013).

17. Mejia, J. et al. Isolated cavities dominate Greenland Ice Sheet dynamic response to lake drainage. *Geophys. Res. Lett.* **48**, e2021GL094762 (2021).

18. Joughin, I., Smith, B. E. & Howat, I. Greenland ice mapping project: ice flow velocity variation at sub-monthly to decadal time scales. *Cryosphere* **12**, 2211 (2018).

19. Moon, T. et al. Distinct patterns of seasonal Greenland glacier velocity. *Geophys. Res. Lett.* **41**, 7209–7216 (2014).

20. Doyle, S. H. et al. Persistent flow acceleration within the interior of the Greenland Ice Sheet. *Geophys. Res. Lett.* **41**, 899–905 (2014).

21. Gagliardini, O. & Werder, M. A. Influence of increasing surface melt over decadal timescales on land-terminating Greenland-type outlet glaciers. *J. Glaciol.* **64**, 700–710 (2018).

22. Goelzer, H. et al. The future sea-level contribution of the Greenland Ice Sheet: a multi-model ensemble study of ISMIP6. *Cryosphere* **14**, 3071–3096 (2020).

23. Maier, N., Gimbert, F., Gillet-Chaulet, F. & Gilbert, A. Basal traction mainly dictated by hard-bed physics over grounded regions of Greenland. *Cryosphere* **15**, 1435–1451 (2021).

24. Joughin, I., Smith, B., Howat, I. & Scambos, T. MEaSURES Multi-year Greenland Ice Sheet Velocity Mosaic, Version 1. NSIDC https://doi.org/10.5067/QUA5Q9SVMSJG (2016).

25. Smith, L. C. et al. Efficient meltwater drainage through supraglacial streams and rivers on the southwest Greenland Ice Sheet. *Proc. Natl Acad. Sci. USA* **112**, 1001–1006 (2015).

26. Stevens, L. A. et al. Greenland supraglacial lake drainages triggered by hydrologically induced basal slip. *Nature* **522**, 73–76 (2015).

27. Selmes, N., Murray, T. & James, T. Fast draining lakes on the Greenland Ice Sheet. *Geophys. Res. Lett.* **38**, L15501 (2011).

28. Tedstone, A. J. et al. Greenland Ice Sheet motion insensitive to exceptional meltwater forcing. *Proc. Natl Acad. Sci. USA* **110**, 19719–19724 (2013).

29. Fettweis, X. et al. Reconstructions of the 1900–2015 Greenland Ice Sheet surface mass balance using the regional climate MAR model. *Cryosphere* **11**, 1015–1033 (2017).

30. Sole, A. et al. Winter motion mediates dynamic response of the Greenland Ice Sheet to warmer summers. *Geophys. Res. Lett.* **40**, 3940–3944 (2013).

31. Meierbachtol, T., Harper, J. & Humphrey, N. Basal drainage system response to increasing surface melt on the Greenland Ice Sheet. *Science* **341**, 777–779 (2013).

32. Banwell, A., Hewitt, I., Willis, I. & Arnold, N. Moulin density controls drainage development beneath the Greenland Ice Sheet. *J. Geophys. Res. Earth Surf.* **121**, 2248–2269 (2016).

33. Röthlisberger, H. Water pressure in intra-and subglacial channels. *J. Glaciol.* **11**, 177–203 (1972).

34. Schoof, C. Ice-sheet acceleration driven by melt supply variability. *Nature* **468**, 803–806 (2010).

35. Stearns, L. A. & van der Veen, C. J. Friction at the bed does not control fast glacier flow. *Science* **361**, 273–277 (2018).

36. Gimbert, F., Gilbert, A., Gagliardini, O., Vincent, C. & Moreau, L. Do existing theories explain seasonal to multi-decadal changes in glacier basal sliding speed? *Geophys. Res. Lett.* **43**, e2021GL092858 (2021).

37. Catania, G., Stearns, L., Moon, T., Enderlin, E. & Jackson, R. Future evolution of Greenland's marine-terminating outlet glaciers. *J. Geophys. Res. Earth Surf.* **125**, e2018JF004873 (2020).

38. Helanow, C., Iverson, N. R., Woodard, J. B. & Zoet, L. K. A slip law for hard-bedded glaciers derived from observed bed topography. *Sci. Adv.* **7**, eabe7798 (2021).

39. Moon, T. A., Gardner, A. S., Csatho, B., Parmuzin, I. & Fahnestock, M. A. Rapid reconfiguration of the Greenland Ice Sheet coastal margin. *J. Geophys. Res. Earth Surf.* **125**, e2020JF005585 (2020).

40. 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**, 2140–2158 (2013).

41. Brondex, J., Gagliardini, O., Gillet-Chaulet, F. & Durand, G. Sensitivity of grounding line dynamics to the choice of the friction law. *J. Glaciol.* **63**, 854–866 (2017).

42. Ritz, C. et al. Potential sea-level rise from Antarctic ice-sheet instability constrained by observations. *Nature* **528**, 115–118 (2015).

43. Vandecrux, B. et al. Firn data compilation reveals widespread decrease of firn air content in western Greenland. *Cryosphere* **13**, 845–859 (2019).

44. Zwally, H. J., Giovinetto, M. B., Beckley, M. A. & Saba, J. L. *Antarctic and Greenland Drainage Systems* (GSFC Cryospheric Sciences Laboratory, 2012); https://earth.gsfc.nasa.gov/cryo/data/polar-altimetry/antarctic-and-greenland-drainage-systems

45. Howat, I., Negrete, A. & Smith, B. The Greenland Ice Mapping Project (GIMP) land classification and surface elevation data sets. *Cryosphere* **8**, 1509–1518 (2014).

46. Howat, I., Negrete, A. & Smith, B. MEaSURES Greenland Ice Mapping Project (GIMP) digital elevation model from GeoEye and WorldView Imagery, Version 1. NSIDC https://doi.org/10.5067/HOKUYVF53Q8M (2017).

47. Gagliardini, O., Cohen, D., Råback, P. & Zwinger, T. Finite-element modeling of subglacial cavities and related friction law. *J. Geophys. Res. Earth Surf.* **112**, F02027 (2007).

48. Weertman, J. The theory of glacier sliding. *J. Glaciol.* **5**, 287–303 (1964).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2022
# Methods

## Quantifying the basal-friction field

We use publicly available elevation and surface velocity datasets to invert for the friction field ($\tau_b$) and basal velocity ($u_b$) across the entirety of Greenland. Surface elevations are from the Greenland Ice Mapping Project surface-elevation model$^{45,46}$, bed elevations are from BedMachine v4$^{49,50}$ and surface velocities are from a multiyear velocity mosaic (MEaSURES program) with low absolute error$^{24,51}$. The surface velocities (Fig. 1) represent decadal-scale averages generated from observations collected over 1995–2015, but are mostly derived from observations over 2005–2015. The time averaging of velocities incorporates both winter and summer flow, and thus they are not expected to reflect a strong seasonality although the exact time period of the measurements is not explicitly quantified$^{51}$.

The inverted friction and basal velocity field (Extended Data Fig. 8) was generated using a state-of-the-art, high-resolution ice-sheet model implemented in Elmer/Ice, a three-dimensional fine-element ice-flow modelling software$^{52,53}$. We used a full-Stokes representation of the momentum balance, which optimizes the friction field in full consideration of three-dimensional flow and Glen's ice rheology (using $n = 3$) with a temperature-dependent viscosity derived from a palaeo-spin-up of the SICOPOLISv5.1 (Simulation Code for Polythermal Ice Sheets) model$^{53}$. The inversion domain consists of a triangular-element mesh with resolution ranging from 800 m in regions of complex flow near the margins to 3 km through to the interior of the ice sheet with a vertical resolution of 16 layers (Extended Data Fig. 8). The basal friction was inverted using the control-inverse method where the 'effective' linear friction coefficient is optimized through 100 iterations to minimize the misfit between the observed and modelled velocity$^{52,54}$ using an initial friction field estimated using the shallow stream approximation (SSA). Even though a linear friction coefficient is used, the inverted friction field has been shown to be insensitive to the type of friction law used$^{55}$. Sensitivity analyses indicate that different regularization, lower mesh resolution, bed topography (BedMachine v3$^{50}$), ice rheology (enhancement factor = 2, temperatures from GRISLlv2 (GRenoble Ice Sheet and Land Ice model)$^{56}$), averaged speed map$^{57}$ or representation of the momentum balance (SSA, Shallow Ice Approximation) yield very similar results to those presented in the main text. Overall, the consistency of the results stems from the analysis being focused on the large-scale features of the friction field that are mainly set by the ice geometry.

## Defining a mean friction relationship

Following the approach of ref. 23, we average and regrid the modelled basal velocity and friction over 3 km by 3 km grid cells to yield spatially independent grid cells of constant dimension across western Greenland. Grid-cell averaging over length scales greater than several ice thicknesses results in a more local force balance, which minimizes the sensitivity of the friction field to inversion parameters$^{23}$ (that is, ice rheology and regularization). We then mask the inverted basal velocity and friction field so that only the likely thawed-bed regions of the ice sheet$^{58}$ where basal motion is expected to occur are included, and remove any grid cells with ice-free areas. We also mask regions with high thickness uncertainty (greater than 1H (ice thickness)) that could impact the inverted friction field. Although, in general, the bed topography in western Greenland is well constrained using radar and resolved using mass conservation$^{49}$, radar flight lines are less dense in southern Greenland and the uncertainty mask mainly impacts spatial coverage in catchment 1 (Fig. 1). We then fit a hard-bed-style friction relationship (rate-strengthening, power-law relationship between basal velocity and friction) to the data from all catchments in western Greenland that experience no summer surface melting using the MATLAB 2019b curve-fitting toolbox (robust nonlinear least squares) with the form:

$$ \tau_{b,RS} = \left( \frac{u_b}{A_s} \right)^{\frac{1}{m}} \qquad (1) $$

where $A_s$ is a parameter related to sliding physics, $u_b$ is the basal velocity, $\tau_b$ is the basal friction and $m$ is the apparent sliding exponent. The model ($A_s = 2.29 [2.09-2.49] \times 10^6 \text{ m yr}^{-1} \text{ MPa}^m$, $m = 4.56 [4.52-4.59]$) shows good adherence to the data ($r^2 = 0.71$), and the binned velocity-friction relationship (median friction for 10 logarithmically spaced velocity bins) and the interquartile range and 90th percentiles indicate that most of the data fall close to the rate-strengthening relationship (Fig. 1).

## Normalizing the friction field to the mean friction relationship

We normalize the friction field to the mean rate-strengthening relationship derived from regions with no surface melt to identify regions that are weaker or stronger than expected. The normalized friction is the percentage difference of basal-friction field from the mean rate-strengthening relationship:

$$ \tau_{b,\text{norm}} = 100 \times \frac{\tau_b - \tau_{b,RS}}{\tau_{b,RS}} \qquad (2) $$

where $\tau_{b,RS}$ is the friction predicted using equation (1) and $\tau_b$ is the friction estimated using Elmer/Ice and the grid-averaging scheme described above. This normalization removes the rate-strengthening friction variability tied to sliding speeds$^{23}$ to infer beds that are weaker or stronger than expected given the basal velocity in each grid cell.

## Comparison with surface mass balance and snowline data

We compare the spatial and along-flow variations in the normalized friction to the time-averaged JJA SMB (2001–2015) and time-averaged snowline$^{43}$ (Fig. 1). The time-averaged JJA SMB is derived from the RACMO 2.3 p2 monthly time series$^{59}$ (Extended Data Fig. 8). RACMO is a regional climate model that is coupled to a surface-energy-balance model that computes the SMB$^{59}$. The output was regridded from 1-km resolution to 3-km resolution for our analysis. The multiyear averaged snowline (2001–2017) is used to delineate the bare ice from the firn. The snowline closely approximates the time-averaged equilibrium line altitude$^{43}$ and is used to roughly separate the accumulation and ablation zones.

## Consideration of non-hydrologic factors on friction variability

Geology, sediment, bed roughness, unaccounted-for ice warming and geothermal heating are factors other than melt forcing that could result in friction variations in Greenland. We assess each of these factors in terms of their likelihood to cause the abrupt transition to weak beds at the snowline in southern Greenland and the north-south threshold in the ablation-zone bed strength. Transitions between bed geologies with different mechanical properties under the ice sheet could influence the friction field via differences in roughness$^{60}$ or a greater propensity for sediment generation$^{61}$. The underlying geology beneath the ice sheet is poorly constrained, but initial approximations of major geologic boundaries$^{62}$ do not show geologic transitions coinciding where major changes in friction are observed (Extended Data Fig. 1). Radar-derived roughness proxies allow us to examine roughness variability at the topographic (>200 m, ref. 63) and sub-topographic (<100 m) scale in relation to the normalized friction variations (Extended Data Fig. 2). Along-flow topographic scale roughness, which influences sliding speeds, shows only limited along-flow variability (shown in ref. 63), does not show a sharp transition between north and south Greenland, or abrupt along-flow changes in roughness in southern Greenland, and is thus unlikely to produce the friction variations we interpret$^{63}$. At sub-topographic length scales, scattering-derived roughness shows a pattern that resembles the normalized friction patterns. This roughness proxy is sensitive to both bed roughness and scattering due to
stored water, which drowns roughness and creates a more specular bed reflector⁶³. Given the close correspondence to the snowline, geo-physical surveys that have found extensive areas of basal water during winter in this region⁶⁴ and the fact that pressure transients common in the ablation zone should generate rougher beds compared with their upgradient regions⁶⁵,⁶⁶, this transition is most likely related to water rather than an abrupt change to smooth beds. Thick, temperate ice layers owing to cryohydrologic warming are not accounted for in the SICOPOLIS temperature field and could also result in an apparent weakening through an unaccounted-for increase in the deformational velocity. Cryohydrologic warming due to water-filled surface crevasses is shown to have a negligible impact on deformation motion⁶⁷, and thus rapid change in rheology that significantly effects the deformational velocity can probably be explained only through widespread basal crevassing, which produces thick temperate basal ice⁶⁸. This has not been found in inland regions⁶⁸ and it is unclear why such a warming process would operate only in southern Greenland. Finally, geothermal heating at the ice base could produce spatial variations in basal melt that influence the bed strength. This aspect is poorly constrained, but modelled geothermal heating⁶⁹ does not show an abrupt marginal change, or a north-south threshold change consistent with the nor- malized friction variations.

## Estimating friction variability associated with decadal-scale velocity variations

We use equation (1) to estimate the magnitude of the friction change (that is, the friction-change equivalent) required to account for the magnitude of decadal velocity variations observed in Greenland⁵,⁶. We rearrange equation (1), a rate-strengthening friction law in terms of velocity:

$$u_b = A_s \tau_b^m \quad (3)$$

Now we define $\Delta u_b/u_b$ in terms of friction changes using equation (3):

$$\frac{\Delta u_b}{u_b} = \frac{A_s \tau_b^m - A_s (\tau_b + \Delta \tau_b)^m}{A_s \tau_b^m} \quad (4)$$

This simplifies to:

$$\frac{\Delta u_b}{u_b} = \left(1 + \frac{\Delta \tau_b}{\tau_b}\right)^m - 1 \quad (5)$$

which is an expression that can be used to approximate the friction-change equivalent of velocity changes. Assuming $m = 3$, which would be expected for a rheology with $n = 3$ and sliding where regulation is negligible⁴⁸, we estimate that the 20% variability in velocities ($\Delta u_b/u_b = 0.2$) observed over decadal timescales⁵,⁶ would reflect friction changes of 6%, whereas assuming $m \approx 4.6$ as found in our study, the friction variability would decrease to 4%.

## Comparison with pressure-based friction relationship

We fit a rate-strengthening relationship from equation (1) to interior elevation bins where no summer melting occurs (that is, non-negative JJA SMB; circles on Fig. 2b) to constrain $A_s$ (8.23 [-1.9-1.83] x 10⁵ m yr⁻¹ MPa⁻¹) and $m$ (4.10 [3.59-4.61]) using a nonlinear least-squares regression with the MATLAB curve-fitting toolbox (solid line, Fig. 2b). Interior bins from catchment 1 are excluded in this fit as Fig. 2a (darkest blue markers) shows that these bins are not consistent with the data from the other catchments. This reference relationship is not pressure based and does not explicitly include cavitation. However, the parameter $A_s$ can implicitly include the effect of cavities under the consideration that cavitation results in a typical 'apparent bed roughness' set by both the geometry of the bed and cavity system. The reference relationship shows a close fit

to the interior bins ($r^2 = 0.85$). The ablation-zone bins from catchments in northern Greenland (squares on Fig. 2b) also follow this scaling, indicating no large change in $A_s$ compared with upgradient regions. By contrast, the bins from southern Greenland (stars on Fig. 2b) are shifted down from the reference scaling and form an offset parallel trend.

To further understand conditions resulting in the offset parallel trend (stars), we plot a frequently used pressure-based hard-bed friction relationship⁴⁷ (dotted lines, Fig. 2b) of the form:

$$\tau_b = CN \left( \frac{X}{1 + \alpha X^q} \right)^{\frac{1}{m}} \quad (6)$$

where $\alpha$ is a function of $q$ and $X = \frac{u_b}{C m^m N^m A_s}$. Here $q$ is a bed-shape parameter, $C$ is a parameter related to the steepness of the bed roughness and $N$ is the effective pressure at the bed. $A_s$ and $m$ are constrained using values from the reference fit to the interior bins. For realistic beds, $q$ is expected to be between 1 and 3 (refs. ³⁶,⁴⁷). In Fig. 2b, $q$ is taken as 1.44, and is derived by assuming that the rightward shift is near Iken's limit (discussed below). The relation is plotted for CN values varying from 0.02 MPa to 0.2 MPa (range of dotted lines, Fig. 2b). The curvature away from the reference relationship is due to a reduction of apparent bed roughness due to the expansion of basal cavities. It follows that the rightward shift of the southern Greenland bins (stars) suggests that the bed is weakened compared with upgradient regions owing to an expansion of basal cavities. Moreover, the linear trend suggests that the apparent bed roughness owing to this cavity expansion converges on a common value that can be captured by replacing the coefficient $A_s$ in equation (1) with $A_{s,eff}$, which incorporates the effective bed roughness owing to both bed bumps and expanded cavitation. Fitting the offset data (stars) using equation (1) with $A_{s,eff}$ and the $m$ from the reference scaling yields a good fit between the model and the offset data ($r^2 = 0.76$; dashed line in Fig. 2b). $A_{s,eff}$ is found to be 2.6 [2.23-3.13] x 10⁶ m yr⁻¹ MPa⁻¹, which is 3.3-times higher than $A_s$. In terms of friction, this suggests that the bed ablation zones of southern Greenland are weakened by 33% compared with the interior bins.

If the assumption is made that the cavitation state lies near Iken's bound, which describes the linear trend at the apex of the pressure-based curves when $\tau_b = CN$ in equation (6) then⁴⁷

$$\tau_b = \left( \left( \frac{q-1}{q} \right) \frac{u_b}{A_s} \right)^{\frac{1}{m}} \quad (7)$$

can be used to constrain $q$ (1.44 [1.33-1.55]). We note that although we assume the cavitation scaling is at Iken's limit to retrieve $q$, this might represent an upper bound. Sliding beyond Iken's limit enters the rate-weakening regime and can lead to unstable ice flow⁷⁰. As our interpretation focuses on dynamics associated with cavitation, rather than dynamics associated with Iken's bound, the value of $q$ does not impact the analysis or implications presented.

## Comparison with other friction relationships

We also contextualized the data using a till-bed friction law from ref. ⁷¹, which is modified to accommodate empirically derived constraints (Extended Data Fig. 4). In this context, scaling for the weakened regions lies in-between the rate strengthening and Coulomb portion of the friction curve and would suggest that the bed is in a mixed state where some regions are deforming and some regions are not⁷¹.

We describe the derivation of the modified till-bed friction law used for Extended Data Fig. 4. The till-bed model in ref. ⁷¹ incorporates a rate-strengthening range, where the ice deforms around clasts before shear failure and plowing occur, and a rate-independent range when plowing and till deformation initiate:
$$ \tau_b = N \tan(\varphi) \left( \frac{u_b}{u_b + u_t} \right)^{1/p} ; \qquad (8) $$

where N is the effective pressure, φ is the internal angle of friction of the till, u_t is the threshold velocity at which the transition from rate-strengthening to rate-independent behaviour begins, and p is an empirically derived sliding exponent. u_t is the velocity at which the stresses imparted by the plowing clasts exceed the strength of the till and is derived from two formulations. The first describes the averaged stress on clasts protruding from the basal ice derived from a hard-bed friction relationship for linear ice flowing over a wavy bed⁷²):

$$ u_b = \left( \frac{\tau_c}{\eta(Ra)^2 k_o^3} + \frac{\tau_c 4 C_1}{(Ra)^2 k_o} \right) \qquad (9) $$

Here, η is the effective viscosity of the ice, R is the clast radius, a is the fraction of the clast protruding from ice, C_1 is a constant related to regulation, k_o = 2π/4R, and τ_c is the averaged stress on each clast. The second formulation describes the stresses that resist plowing:

$$ \tau_p = \frac{N_f N}{2 + N_f k}, \qquad (10) $$

where k accounts for a till-strength reduction owing to clast shadows and N_f is the bearing capacity factor of the till. Plowing occurs when τ_c = τ_p. We can substitute our data-derived parameters derived from the reference scaling (Fig. 2b and equation (1)) for equation (9), which yields u_b as a function of our constrained parameters and a characteristic clast roughness R_c, which is necessary to estimate the clast-averaged stress:

$$ u_b = \frac{A_s}{R_c} (\tau_b)^m \qquad (11) $$

Here A_s and m are derived from equation (1). By setting τ_b = τ_p, we find:

$$ u_t = \frac{A_s}{R_c} \left( \frac{N_f N}{2 + N_f k} \right)^m \qquad (12) $$

which is roughly equivalent to that in ref. ⁷¹. Parameter values come from the ranges presented in ref. ⁷¹, where N_f = 40, k = 0.45 and R_c = 0.535. The latter is tuned to fit the data in reference scaling shown in Extended Data Fig. 4.

## Drainage modelling

We estimate the likelihood for channelization to occur in summer by conservatively assuming channels form if the mean JJA discharge exceeds the critical discharge following ref. ³⁴:

$$ Q_c = \frac{u_b h}{c_1 (\alpha - 1) \psi} \qquad (13) $$

where h is the characteristic height of bed roughness, ψ is the hydro-potential gradient calculated assuming ice overburden pressure, and c_1 and α are constants related to viscous closure and the Darcy-Weisbach friction factor, respectively. We estimate the mean discharge (Q_m) using a D8 water-routing algorithm⁷³ that calculates upgradient flow accumulation weighted by the JJA SMB and basal melt⁶⁹. Using a step height of 0.1, we show that Q_c is exceeded near the margin in both northern and southern Greenland (Extended Data Fig. 5).

To approximate the effective pressure within actively draining pathways and examine how it changes in space, we again follow ref. ³⁴, where a steady-state pressure formulation is given for active drainage pathways consisting of either cavities or channels. This formulation includes

feedbacks related to discharge as well as hydropotential gradients and sliding opening:

$$ N^n = \frac{c_1 Q \psi + u_b h}{c_2 c_3^{-1/\alpha} Q^{1/\alpha} \psi^{-1/2\alpha}} \qquad (14) $$

where c_2 and c_3 are related to the ice viscosity and the Darcy-Weisbach friction factor, respectively. We apply this formulation to the accumulated flow estimate (Extended Data Fig. 5). It is important to note that although these equations apply to only discrete drainage pathways, melt inputs are discrete in nature²⁵ and will concentrate along low-potential flow pathways; thus, this formulation is a reasonable approximation to estimate the effective pressure within active drainage pathways and examine how it varies over large scales. Although steady-state conditions are not reached owing to rapidly fluctuating discharge during summer, a steady range of moulin values is consistently observed (that is, background plus diurnal variations) during late summer¹²,⁷⁴ and thus quasi-steady-state behaviour is suggested by observations. The parameters used in equations (13) and (14) follow those in ref. ³⁴.

## Parameterized implementation of threshold

We use the reference and cavity friction scalings identified in Fig. 2b to model the north-south friction threshold observed in the normalized friction field in Fig. 1h. We use the effective pressure within the active drainage pathways as modelled using the mean JJA SMB, ice geometry and sliding speeds to determine which friction scaling to apply. We apply the cavity scaling to N < 1.3 MPa and the reference scaling to N > 1.7 MPa. For effective pressures between the two scalings, we linearly interpolate between the two friction solutions to create a smooth transition in the friction field. The results show that the main north-south threshold can be reproduced (Extended Data Fig. 6), yielding a parameterization that can be coupled with future changes in mass balance and surface geometry to investigate the feedback observed here.

This approach has some limitations. We do not explicitly account for the limited strengthening observed in northern Greenland. In addition, although the model captures the north-south threshold, it overestimates the friction in the southern-most catchment, and underestimates it in the northern-most catchment. This suggests that although sliding and hydropotential gradients are the primary factors controlling the observed threshold, the large difference in melt conditions in northern and southern Greenland might also have a role that is not resolved in our analysis.

## Data availability

All data used in this study are currently archived and publicly available for download. Multiyear surface velocities: https://doi.org/10.5067/QUA5Q9SVMSJG. Ice surface elevation: https://doi.org/10.5067/H0KUYVF53Q8M. Land/Ice classification: https://doi.org/10.5067/B8X58MQBFUPA. Bed topography: https://doi.org/10.5067/2CIX82HUV88Y. Basal thermal state map: https://doi.org/10.5067/R4MWDWWUWQF9. Snowline elevation: https://doi.org/10.18739/A2V40JZ6C. Basal melt: https://doi.org/10.22008/FK2/PLNUEO. RACMO surface mass balance data can be requested from Brice Noël (b.p.y.noel@uu.nl).

## Code availability

Elmer/Ice files and code used to produce this paper are archived at https://doi.org/10.5281/zenodo.6535396.

49. Morlighem, M. et al. BedMachine v3: complete bed topography and ocean bathymetry mapping of Greenland from multibeam echo sounding combined with mass conservation. Geophys. Res. Lett. 44, 11051-11061 (2017).

50. Morlighem, M. IceBridge BedMachine Greenland, Version 3. NSIDC https://doi.org/10.5067/2CIX82HUV88Y (2018).
51. Joughin, I., Smith, B. E. & Howat, I. M. A complete map of Greenland ice velocity derived from satellite data collected over 20 years. *J. Glaciol.* **64**, 1–11 (2018).

52. Gagliardini, O. et al. Capabilities and performance of Elmer/Ice, a new-generation ice sheet model. *Geosci. Model Dev.* **6**, 1299–1318 (2013).

53. Goelzer, H., Robinson, A., Seroussi, H. & van de Wal, R. S. W. Recent progress in Greenland Ice Sheet modelling. *Curr. Clim. Change Rep.* **3**, 291–302 (2017).

54. Morlighem, M. et al. Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica. *Geophys. Res. Lett.* **37**, L14502 (2010).

55. Joughin, I., MacAyeal, D. R. & Tulaczyk, S. Basal shear stress of the Ross ice streams from control method inversions. *J. Geophys. Res. Solid Earth* **109**, B09405 (2004).

56. Quiquet, A. & Dumas, C. The GRISLI-LSCE contribution to the Ice Sheet Model Intercomparison Project for phase 6 of the Coupled Model Intercomparison Project (ISMIP6)—Part 1: projections of the Greenland Ice Sheet evolution by the end of the 21st century. *Cryosphere* **15**, 1015–1030 (2021).

57. Mouginot, J., Rignot, E., Scheuchl, B. & Millan, R. Comprehensive annual ice sheet velocity mapping using Landsat-8, Sentinel-1, and RADARSAT-2 data. *Remote Sens.* **9**, 364 (2017).

58. MacGregor, J. et al. A synthesis of the basal thermal state of the Greenland Ice Sheet. *J. Geophys. Res. Earth Surf.* **121**, 1328–1350 (2016).

59. Noël, B. et al. Evaluation of the updated regional climate model RACMO2.3: summer snowfall impact on the Greenland Ice Sheet. *Cryosphere* **9**, 1831–1844 (2015).

60. Woodard, J., Zoet, L., Iverson, N. R. & Helanow, C. Linking bedrock discontinuities to glacial quarrying. *Ann. Glaciol.* **60**, 66–72 (2019).

61. Crompton, J. W. & Flowers, G. E. Correlations of suspended sediment size with bedrock lithology and glacier dynamics. *Ann. Glaciol.* **57**, 142–150 (2016).

62. Dawes, P. R. The bedrock geology under the Inland Ice: the next major challenge for Greenland mapping. *Geol. Surv. Den. Greenl. Bull.* **17**, 57–60 (2009).

63. Cooper, M. A. et al. Subglacial roughness of the Greenland Ice Sheet: relationship with contemporary ice velocity and geology. *Cryosphere* **13**, 3093–3115 (2019).

64. Chu, W. et al. Extensive winter subglacial water storage beneath the Greenland Ice Sheet. *Geophys. Res. Lett.* **43**, 12484–12492 (2016).

65. Cohen, D., Hooyer, T. S., Iverson, N. R., Thomason, J. & Jackson, M. Role of transient water pressure in quarrying: a subglacial experiment using acoustic emissions. *J. Geophys. Res. Earth Surf.* **111**, F03006 (2006).

66. Hallet, B. Glacial quarrying: a simple theoretical model. *Ann. Glaciol.* **22**, 1–8 (1996).

67. Poinar, K., Joughin, I., LENAERTS, J. T. & Van Den Broeke, M. R. Englacial latent-heat transfer has limited influence on seaward ice flux in western Greenland. *J. Glaciol.* **63**, 1–16 (2017).

68. Harrington, J. A., Humphrey, N. F. & Harper, J. T. Temperature distribution and thermal anomalies along a flowline of the Greenland Ice Sheet. *Ann. Glaciol.* **56**, 98–104 (2015).

69. Karlsson, N. B. et al. A first constraint on basal melt-water production of the Greenland Ice Sheet. *Nat. Commun.* **12**, 3461 (2021).

70. 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).

71. Zoet, L. K. & Iverson, N. R. A slip law for glaciers on deformable beds. *Science* **368**, 76–78 (2020).

72. Nye, J. F. A calculation on the sliding of ice over a wavy surface using a Newtonian viscous approximation. *Proc. R. Soc. Lond. A* **311**, 445–467 (1969).

73. Schwanghart, W. & Scherler, D. TopoToolbox 2—MATLAB-based software for topographic analysis and modeling in Earth surface sciences. *Earth Surf. Dyn.* **2**, 1–7 (2014).

74. Covington, M., Gulley, J., Trunz, C., Mejia, J. & Gadd, W. Moulin volumes regulate subglacial water pressure on the Greenland Ice Sheet. *Geophys. Res. Lett.* **47**, e2020GL088901 (2020).

**Acknowledgements** This work was funded by the French National Research Agency grant ANR-17-CE01-0008. We thank B. Noël for providing the RACMO surface mass balance data. The Elmer/Ice computations presented in this paper were performed using HPC resources of both CINES under the allocation 2021- A0060106066 made by GENCI and of the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities. We thank A. Gilbert and P. Christoffersen for discussions that improved this work.

**Author contributions** N.M. designed the study, analysed the data, produced the figures and is the primary author of the manuscript. F.G. contributed significantly to the data analysis and writing of the manuscript. F.G.-C. performed the friction inversions implemented in Elmer/Ice and contributed to the writing of the manuscript.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-022-04927-3.

**Correspondence and requests for materials** should be addressed to Nathan Maier.

**Peer review information** Nature thanks Stephen Livingstone and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_9.png" /></figure>

**Extended Data Fig. 1 | Relationship between normalized friction and**

**major geologic boundaries.** Geologic boundaries (dashed lines) were
adapted from Dawes et al. (2009). Unit descriptions are as follows:

1: Palaeoproterozoic Precambrian Shield, 2: Meso-Neoarchean Precambrian Shield, 3: Palaeoproterozoic with reworked Archaean Precambrian Shield, 4: Basalts and intrusions, Palaeogene Volcanic Province, 5: Mesoproterozoic Red Beds, 6: Mainly Archaean Precambrian Shield, 7: Palaeoproterozoic Precambrian Shield, 8: Mesoproterozoic Porphyries and Red Beds, 9: Phanerozoic Basin Formation from the Cambrian through Silurian.
<figure><img src="image_10.png" /><figcaption>**Extended Data Fig. 2** | Relationship between normalized friction and radar-derived roughness proxies. (a) Normalized friction plotted adjacent to (b) topographic roughness and (c) scattering-derived roughness proxies</figcaption></figure>

derived from radar from ref. ⁶³. Topographic roughness reflects roughness at the greater than 200 m length scale, while the scattering-derived roughness proxy is estimated to reflect roughness at length scales less than 100 m.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 3** | Regression between glaciological variables and normalized friction. Relationship between hydropotential gradients, sliding speed, JJA SMB, and ice thickness and the normalized friction for the upper and

<figure><img src="image_14.png" /></figure>

lower ablation zone. Individual markers are the mean of 100 evenly spaced latitudinal bins along the margin for the upper and lower ablation zone.
<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 4 | Observations contextualized with till-bed friction law.** Mean basal velocity and friction for elevation bins from each catchment are plotted against a till-bed friction law modified from ref. ⁷¹ (Methods). Circle markers are interior elevation bins that experience no summer melt (from all catchments), square markers are ablation-zone bins for north Greenland (catchments 5–7), stars are ablation-zone bins for south Greenland (catchments 1–4). Horizontal and vertical bars indicate 95% C.I. Reference scaling (solid line, r² = 0.85) is constrained using bins with no melt forcing

(circles, Methods). The ablation-zone bins for northern Greenland (squares) adhere to this relationship while the ablation-zone bins for southern Greenland (stars) are shifted downward indicating increased bed deformation compared to upgradient regions. This trend is captured by a deformation scaling (dashed line, r² = 0.76) (Methods). A constrained, effective-pressure-based till-bed friction law<sup>71</sup> (Methods, dotted lines) is plotted for various N. Deviation from the reference scaling (i.e. the curvature in dotted lines) results from increasing bed deformation as effective pressures decrease.
<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 5** | Drainage analysis. Panel (a) shows the summer mean accumulated discharge ($Q_{sum}$), while panel (b) shows the critical discharge ($Q_c$) required for channel formation following ref. ³⁴ (Methods). Panel (c) shows the regions where the accumulated discharge exceeds the critical discharge

indicating channel formation is likely. Panel (d) shows the modelled effective pressure within active drainage pathways using a formulation from ref. ³⁴ which incorporates cavities and conduits (Methods).
<figure><img src="image_17.png" /><figcaption>**Extended Data Fig. 6** | Parameterized implementation of hydrologic threshold. Hydrologic threshold is reproduced using a parameterization which implements the reference and cavity friction scaling based on the modelled effective pressure within active drainage pathways (Methods).</figcaption></figure>
<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 7** | Friction summary. Relationship between basal velocity and friction for 20 evenly spaced elevation bins is plotted for each catchment. Markers with red outline are elevation bins located within ablation zone, and

those with black outline are elevation bins where the mean flotation fraction is greater than zero. Relationships from all catchments are plotted together in bottom right panel. Horizontal and vertical bars indicate 95% C.I.
<figure><img src="image_19.png" /></figure>

**Extended Data Fig. 8** | **Data summary.** (a) Multiyear averaged surface velocities (1995–2015)²⁴,⁵¹, (b) element nodes for the full-Stokes friction inversions, (c) grid cells used in analysis, (d) mean JJA surface mass balance

(2001–2015) from RACMO 2.3p2, (e) modelled basal friction and (f) modelled basal velocity from Elmer/Ice inversion (Methods). Coloured lines show catchment boundaries. All maps are plotted using MATLAB.