The land-terminating margin of the Greenland Ice Sheet has slowed down in recent decades, although the causes and implications for future ice flow are unclear. Explained originally by a self-regulating mechanism where basal slip reduces as drainage evolves from low to high efficiency, recent numerical modeling invokes a sedimentary control of ice sheet flow as an alternative hypothesis. Although both hypotheses can explain the recent slowdown, their respective forecasts of a long-term deceleration versus an acceleration of ice flow are contradictory. We present amplitude-versus-angle seismic data as the first observational test of the alternative hypothesis. We document transient modifications of basal sediment strengths by rapid subglacial drainages of supraglacial lakes, the primary current control on summer ice sheet flow according to our numerical model. Our observations agree with simulations of initial postdrainage sediment weakening and ice flow accelerations, and subsequent sediment restrengthening and ice flow decelerations, and thus confirm the alternative hypothesis. Although simulated melt season acceleration of ice flow due to weakening of subglacial sediments does not currently outweigh winter slowdown forced by self-regulation, they could dominate over the longer term. Subglacial sediments beneath the Greenland Ice Sheet must therefore be mapped and characterized, and a sedimentary control of ice flow must be evaluated against competing self-regulation mechanisms.

It is clear that the opposing implications of the two hypotheses for the future evolution of ice flow and mass loss from the GrIS must be reconciled. Direct observations that assert the simulated modifications of subglacial sediment strengths, and thus test the hypothesis of a sedimentary control of ice flow, are a requirement in achieving this reconciliation. Seismic reflection surveys are particularly well suited for hypothesis testing because they have a successful history of delineating and characterizing the substrate beneath the Antarctic Ice Sheet, and increasingly so beneath the GrIS ( 13 , 16 – 18 ). On the basis of an in situ study in the Kangerlussuaq sector in West Greenland, we present here the first seismic observations of subglacial sediments affected by supraglacial lake drainages, and explain them in a catchment-wide context using numerical modeling.

Either hypothesis of self-regulation or widespread sediment strengthening can therefore explain the observed decadal slowdown of GrIS’s land-terminating margin. However, a straight contradiction arises in their implications for the long-term evolution of ice flow, as forced not only by increased meltwater volume but also by increased variability in its delivery in a warming climate. Self-regulation and therefore slowdown of ice flow would endure or even magnify in the long term according to the first hypothesis ( 7 – 12 ). This conflicts directly with a modeled increase in the future significance of daily meltwater variability in weakening subglacial sediments and forcing enhancements of melt season ice flow ( 3 ). According to the second hypothesis, mean annual acceleration, and not deceleration, of ice sheet flow would result in the long term, as melt season accelerations progressively outweigh late-season and winter slowdown ( 3 , 27 ).

The second hypothesis builds on the growing recognition that large swathes of the GrIS are underlain by soft sediments that control basal traction and ice flow so that the assumption of negligible hydrodynamic effects may be flawed ( 13 – 21 ). Fluctuating subglacial water pressure gradients are proposed to drive water into or out of a basal sediment layer, respectively weakening or strengthening it by decreasing or increasing its pore water pressure and shear resistance. Sediment weakening or strengthening would correspondingly reduce or enhance basal traction, thereby causing the ice flow of the GrIS to accelerate or decelerate ( 19 , 22 ). This mechanism has been suggested by numerical modeling to be greatly enhanced when supraglacial lakes drain and large quantities of meltwater are discharged rapidly into the subglacial environment, resulting in an overall net strengthening of sediment ( 3 ). Moreover, the cumulative effect of hundreds of melt season supraglacial lake drainages ( 23 – 26 ) on mean annual sediment traction and ice flow is currently modeled to outweigh that of aggregated meltwater variability caused by daily changes in subglacial runoff ( 3 ). Supraglacial lake drainages thus emerge as the primary predicted control on present-day seasonal ice flow, and widespread sediment strengthening by cumulative supraglacial lake drainages is modeled to effect the observed slowdown of the GrIS in recent decades ( 3 ). To date, in situ observations that assert the modeled sedimentary control of ice flow have not been available.

The driving forces are currently subject to debate, and two alternative hypotheses exist. The first of these considers self-regulation of mean annual ice flow ( 7 – 11 ). Early in the melt season, subglacial channel networks are poorly developed and easily overwhelmed by incipient meltwaters, reducing basal traction and accelerating ice flow. Because highly efficient networks of low-pressure subglacial channels develop through the boreal summer, basal traction increases progressively by either stored waters being drawn into the channels ( 7 – 11 ) or complex lateral stress transfers ( 12 ). The resulting slowdown of ice flow through the melt season and the subsequent winter is then proposed to overcompensate for early-season accelerations, sustaining mean annual ice flow reductions over recent decades as increasingly larger channel networks persisted through the winter ( 8 , 10 ). Any hydrodynamic effects of the ice sheet substrate are assumed to be negligible in this self-regulation hypothesis.

It was proposed initially that increased meltwater volumes in a warming climate would scale inversely with basal traction and therefore directly with future ice flow speed ( 4 ). Since then, ice flow accelerations in the summer have been observed to correlate with meltwater variability ( 5 – 7 ). This shorter-term variability arises principally from the drainage of thousands of surface (supraglacial) lakes, diurnal melt cycles, and rainfall, affecting subglacial water pressure gradients, hydrological network development, and thus basal traction ( 6 ). On multiannual time scales, an entirely different picture has been emerging, where enhanced meltwater availability may lead to annually averaged ice flow deceleration, rather than longer-term acceleration ( 7 – 11 ). Perhaps most strikingly, the land-terminating southwestern margin of the GrIS appears to have slowed down by ~12% between 1985–1994 and 2007–2014, despite a 50% increase in surface meltwater production ( 10 ). It is therefore important that the processes driving this observed decadal deceleration, and the future evolution of ice flow, are fully understood and embedded within ice sheet models.

The contribution of the Greenland Ice Sheet (GrIS) to eustatic sea level rise has increased more than sixfold since 1992, due in almost equal measure to increased surface melt runoff and dynamic ice discharge ( 1 ). Current understanding of dynamic discharge processes is, however, disproportionally poor compared to their significance in conveying interior ice into the ocean. Pressurized subglacial water reduces traction and thus enables well-lubricated basal slip, and as such is modeled to exert a strong control on the form and flow of the GrIS ( 2 ). The impacts of basal water on the ice sheet’s mass balance are less certain, however, because observations of the substrate and spatial and temporal variations in its traction are sparse ( 2 , 3 ). Our study directly addresses this problem.

If a thin layer of subglacial water was perched above either lodged (4.26 ± 0.59 × 10 6 kg m −2 s −1 ) or dilatant (3.0 × 10 6 to 3.4 × 10 6 kg m −2 s −1 ) sediment, our existing models of subglacial AVA responses ( 16 ) would predict strongly negative zero-incidence reflectivities with respective values of −0.22 and −0.26 to −0.28. These models are incompatible with the weakly positive zero-incidence reflectivities measured in both profiles S1 ( Fig. 3C ) and S2 ( Fig. 3D ), thus excluding the possibility of a thin water layer beneath the ice-bed interface.

The northwestern subglacial valley’s substrate (profile N, Fig. 2A ) is characterized by a positive zero-incidence reflectivity and a negative AVA gradient ( Fig. 3B ) and has an acoustic impedance of 4.17 ± 0.11 × 10 6 kg m −2 s −1 and a Poisson’s ratio of 0.06 ± 0.05. These AVA attributes are compatible with a very stiff substrate of low water content ( 16 ) and are interpreted as lodged subglacial sediment ( 34 ). The positive AVA gradients observed for profiles S1 ( Fig. 3C ) and S2 ( Fig. 3D ) are conceptually indicative of either dilatant sediment or water beneath the ice base ( Fig. 3A ), although observed zero-incidence reflectivities are positive ( Fig. 3 , C and D) instead of negative, as would be expected from theory in these cases ( Fig. 3A ). We showed previously that this apparent contradiction is diagnostic of a composite seismic reflection, where lodged sediment is overlain by a seismically thin cap of dilatant sediment ( 16 ). On the basis of our established models of “thin-layer” seismic responses, we thus infer the presence of a layer of dilatant sediment within the subglacial basins in profiles S1 and S2 ( Figs. 2 , B and C, and 3 , C and D) that is less than 2 m thick with an acoustic impedance range of 3.0 × 10 6 to 3.4 × 10 6 kg m −2 s −1 and a Poisson’s ratio approaching 0.5 ( 16 ). This thin dilatant sediment layer is interpreted to be underlain by stiffer, lodged sediment with acoustic impedance of 4.26 ± 0.59 × 10 6 kg m −2 s −1 ( 16 ). The inferred acoustic impedances of the lodged sediments identified by profiles N, S1, and S2 are therefore indistinguishable ( Fig. 3 , B to D).

( A ) Conceptual AVA curves for possible ice substrates ( 48 ), with range of panels (B to D) highlighted. ( B to D ) AVA data and uncertainty ranges derived respectively for profiles N, S1, and S2, where solid red lines are best fits from Bayesian modeling, as explained in Materials and Methods. The presence of lodged subglacial sediment is revealed in profile N, and an additional layer of thin dilatant sediment is indicated by composite reflections in profiles S1 and S2 ( 16 ).

All three seismic profiles imaged undulating glacier beds characterized by topographic variations of up to 150 m ( Fig. 2 ). Profile N shows a continuous reflector beneath, and subparallel to, the basal reflector between positions 1.0 and 2.3 km ( Fig. 2A ). Profiles S1 and S2 show planar sections of reflectivity, dipping northeast at ~1° and northwest at ~3.5°, respectively, beneath which concave reflections suggest subglacial basins, up to 700 m wide ( Fig. 2 , B and C). These interlake basins are distinct from the ice-bed parallel reflector dominant in profile N ( Fig. 2A ). Amplitude-versus-angle (AVA) seismic analysis ( Fig. 3A ) revealed that the AVA gradient for profile N ( Fig. 3B ) is negative, whereas the AVA gradients for profiles S1 ( Fig. 3C ) and S2 ( Fig. 3D ) are both positive, and all three profiles are characterized by positive reflection coefficients at incidence angles of 0° (“zero incidence”).

( A to C ) Respective seismic structure of the ice sheet base along profiles N, S1, and S2 within a 300-m depth window. Automatic gain control, with a 300-ns window, was applied for display, and the yellow dashed lines show the intersections between profiles S1 and S2 in (B) and (C). Major subglacial sediment basins and the range of ice substrate reflectors supplied to our AVA analyses are indicated.

We acquired three two-dimensional (2D) seismic reflection profiles, including profile N in the northwestern subglacial valley between 3 and 5 July 2010, and profiles S1 and S2 beyond the southern margin of Lake F on 10 July and 16 to 17 July 2010, respectively ( Figs. 1C and 2 ). Profile N was located immediately adjacent to the input point of Lake F waters into the northwestern subglacial valley ( Fig. 1C ), and therefore surveyed a basal environment that experienced extreme discharge rates during its rapid drainage some 4 to 6 days before. Calculated hydrological gradients reveal that a hydraulic barrier separated the interlake ice sheet bed beneath profiles S1 and S2 from the subglacial input point of Lake F waters ( Fig. 1C ). Seismic data acquisition and processing are explained in Materials and Methods.

The rapid drainage of Lake F was characterized by peak discharge and uplift rates—due to elastic hydraulic jacking—of ~3300 m 3 s −1 and 0.8 m hour −1 attained ~1 hour after initiation ( 25 ). The lake had a volume of ~7.4 × 10 6 m 3 immediately before drainage and emptied in ~2 hours. Numerical simulations constrained by Global Positioning System (GPS) observations of ice uplift and acceleration at the ice surface before, during, and after lake drainage were consistent with the subglacial dispersal of the lake waters via a transient turbulent sheet and the surrounding distributed subglacial water system ( 25 , 31 , 32 ). Basal ice motion was measured to be enhanced for several hours ( 25 ). Passive seismic monitoring of lake drainage ( 25 ) and icequake analysis using novel techniques ( 33 ) revealed that input of lake waters into the subglacial environment occurred below its western margin ( Fig. 1 ). The waters then spread out into the northwestern subglacial valley, as affirmed initially by reconstruction of ice tectonics using fracture patterns, by uplift and separation of GPS stations and rates of seismicity ( 25 ), and subsequently by numerical modeling ( 31 ). The assertions are consistent with the morphological dominance of the northwestern subglacial valley within the Kangerlussuaq sector ( Fig. 1A ). Several moulins then fed the subglacial hydrological system for the remainder of the melt season ( 25 ). Conduit closure modeling indicates that the subglacial system beneath Lake F remained distributed and at high pressure throughout the drainage event, potentially with a number of smaller subglacial conduits or cavities operating within it ( 25 , 31 , 32 ).

We focused our seismic surveys on supraglacial Lake F in the sector’s Russell Glacier catchment ( Fig. 1 and fig. S1), because it was large and had a high recurrence rate of annual rapid drainage since records began in 2002 ( 30 ). The lake is located on ~1200-m-thick ice and at the head of a broad subglacial valley through which ice and subglacial water flows are channeled toward the northwest ( Fig. 1 ) ( 28 ). The mean date of rapid drainage of Lake F between 2002 and 2009 was 14 July, but 2010 was an abnormally warm melt season and rapid drainage occurred some 2 weeks earlier during the night of 29 to 30 June (fig. S1) ( 25 ). Multiple rapid lake drainages tended to occur in distinct clusters within the same elevation band, migrating up-glacier during the course of the melt season (fig. S1). Accordingly, Lake F was part of a cluster of lake drainages that occurred between 30 June and 2 July 2010 ( 25 , 30 ). Before this, only little surface melt accessed the ice sheet bed via small moulins ( 25 ) so that the subglacial hydrological system below Lake F was likely distributed and inefficient ( 31 ).

( A ) Map of ice sheet bed elevation revealing the location of Lake F at the head of a major subglacial valley that ultimately extends to the western ice sheet margin. The outline of Lake F is shown by the solid black line, and the estimated input point of Lake F waters into the subglacial environment is shown by the filled white circle. The white star shows the location of site SHR referred to in the text. ( B ) Location of study region in Greenland. ( C ) Map of the hydraulic potential gradients and subglacial topography, where the length of arrows is proportional to the magnitude of flux. Seismic profiles N, S1, and S2 and their acquisition directions are indicated by thick black arrows. The dotted box indicates the 6 km × 6 km area over which sediment strengths and ice flow velocities in Fig. 5 were averaged.

The Kangerlussuaq sector is well suited for an investigation of subglacial conditions because it is land-terminating, and hydrological forcing of ice flow is therefore isolated from oceanic influences. Furthermore, a rich scientific baseline includes high-density radar mapping of the bed ( Fig. 1A ) ( 28 ) and compelling evidence for an abundance of subglacial sediments ( 15 – 17 , 29 ). The sector has also a high areal extent of supraglacial lakes ( 23 ), of which some 28% drain rapidly and in clusters to the bed every melt season ( 30 ).

DISCUSSION

Although our seismic observations are consistent with common sediment textures and depositional histories in all profile locations, the upper horizon of the subglacial sediments in the basin to the south of Lake F is interpreted to be soft and weak, whereas subglacial sediments beyond its northwestern margin are interpreted to be stiff and strong. These observations are counterintuitive because the latter sediments had been affected by rapid subglacial drainage of Lake F some 4 to 6 days before, whereas those in the southern basin had not. The hypothesis of a sedimentary control of ice flow that we aim to test here was suggested by numerical modeling (3). We therefore use herein the same model to understand our tantalizing observations within the broader context of the seasonal evolution of subglacial hydrology and ice flow on the catchment scale. The numerical model is described in Materials and Methods.

Numerical modeling and hypothesis testing Our model simulations highlight cumulative subglacial water fluxes due to multiple supraglacial lake drainage because their aggregate impact on subglacial sediment strength currently outweighs that of daily changes in runoff (above). Although these fluxes were negligible in our study area before the night of 29 to 30 June 2010, the drainage of Lake F during that night caused a distinct spike in cumulative simulated fluxes that are particularly concentrated in the northwestern subglacial valley (Fig. 4A). These simulations agree with previous observations (25), numerical modeling (31), and its morphological dominance in the Kangerlussuaq sector (Fig. 1A). The subsequent drainage of several supraglacial lakes then caused widespread subglacial water fluxes in the Russell Glacier catchment during the 3 to 5 July acquisition period of profile N (Fig. 4B). Many more supraglacial lakes then drained and further enhanced subglacial water fluxes simulated between 5 and 10 July (Fig. 4C), when profile S1 was acquired, and between 10 July and 16 to 17 July (Fig. 4D), when profile S2 was acquired. Fig. 4 Modeled cumulative fluxes due to supraglacial lake drainages in the 2010 melt season. (A) On 30 June, showing fluxes due to Lake F drainage (white circle). (B) On 5 July, showing drainages on 30 June (Lake F, white circle) and 3 to 5 July (pink circles). (C) On 10 July, showing drainages on 6 July (orange circles with white border) and 10 July (orange circles with black border). (D) On 17 July, showing drainages on 11 to 16 July (red circles with black border) and 17 July (red circles with white border). The gray box outlines the 6 km × 6 km area over which sediment strength and ice velocity were averaged in Fig. 5. To facilitate an integrated comparison of our AVA results and model simulations, we have averaged the latter over an area of 6 km × 6 km that is centered on Lake F and includes all three seismic profiles (Fig. 1C). The drainage of Lake F on 29 to 30 June 2010 (days 180 to 181) is simulated to weaken the subglacial sediments in the northwestern valley, but by the time we acquired our seismic profile N on 3 to 5 July (days 184 to 186), the sediments are modeled to have restrengthened again (Fig. 5). In contrast, simulated sediment strengths were low both on 10 July 2010 (day 191) when profile S1 was acquired and on 16 to 17 July 2010 (days 197 to 198) when profile S2 was acquired (Fig. 5). Fig. 5 Modeled subglacial sediment strength and ice flow velocity, averaged over a 6 km × 6 km area centered on Lake F (dotted box in Fig. 1C ). The gray shaded areas indicate the acquisition periods of our seismic profiles N, S1, and S2. Sediments are modeled to be strong and weak when profiles N and S1/S2, respectively, were acquired. Corresponding modeled subglacial water fluxes are shown in Fig. 4. Our seismic observations are therefore consistent with changes in subglacial sediment strength forced by supraglacial lake drainages on the GrIS, in full agreement with our numerical simulations (3). Furthermore, the model simulates pronounced accelerations of ice flow when sediments are measured to be weak (for example, on days 191 and 197 to 198) and decelerations when sediments are measured to be strong (for example, on days 184 to 186) (Fig. 5). Because simulated changes in ice flow velocities and sediment strengths generally scale inversely with each other (Fig. 5) (3), our seismic observations confirm the hypothesis of a sedimentary control of ice sheet flow by implication.