Research ArticleGEOLOGY

Seismic evidence for complex sedimentary control of Greenland Ice Sheet flow

See allHide authors and affiliations

Science Advances  16 Aug 2017:
Vol. 3, no. 8, e1603071
DOI: 10.1126/sciadv.1603071


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.


Short-term versus long-term changes of ice sheet flow in a warming climate

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.

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 (57). 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 (711). 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.

Hypotheses of subglacial hydrological control of ice sheet flow

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 (711). 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 (711) 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.

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 (1321). 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 (2326) 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.

Contradictory forecasts of long-term ice sheet flow

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 (712). 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).

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, 1618). 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.


Study site and data acquisition

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 (1517, 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).

Fig. 1 Kangerlussuaq sector study site.

(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.

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

The rapid drainage of Lake F was characterized by peak discharge and uplift rates—due to elastic hydraulic jacking—of ~3300 m3 s−1 and 0.8 m hour−1 attained ~1 hour after initiation (25). The lake had a volume of ~7.4 × 106 m3 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).

Data acquisition

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.

Fig. 2 Seismic images of the ice sheet base.

(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.

Nature and properties of the ice sheet substrate

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

Fig. 3 AVA seismic analysis.

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

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 × 106 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 × 106 to 3.4 × 106 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 × 106 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).

If a thin layer of subglacial water was perched above either lodged (4.26 ± 0.59 × 106 kg m−2 s−1) or dilatant (3.0 × 106 to 3.4 × 106 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.


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.

Synthesis and implications for long-term evolution of GrIS ice flow

Our numerical simulations suggest that cumulative supraglacial lake drainages are currently the dominant control on subglacial sediment strength, and therefore of the ice flow of the land-terminating margin of the GrIS (3). The simulations also predict that in the long term, daily runoff variability will replace these lake drainages as the dominant control. This control comprises cumulative runoff events, causing widespread sediment weakening both in marginal areas already affected by melt season drainage and in interior regions of the GrIS that will be affected in the future (3, 27, 35, 36). Melt season accelerations of ice flow would then become progressively enhanced relative to the present day, to a level where they can no longer be compensated for by winter slowdown and thus leading to mean annual ice sheet acceleration (3, 27, 37). Sustained softening of subglacial sediments by continued drainage would allow greater water storage throughout the melt season, thus contributing further to sustained and, owing to potentially lagged release of water from storage, ice acceleration on time scales longer than those recognized so far.

Our seismic observations lend credibility to these predictions and therefore support the notion that the future response of the GrIS to climate warming will be more complex than the relatively simplistic view of long-term ice sheet deceleration due to self-regulation of ice flow (711). Complexity arises particularly because the self-regulation and sediment hypotheses are not mutually exclusive and will interplay to varying degrees both spatially and temporally. Our seismic data show, for example, that sediment tends to be focused in basal depressions, whereas its presence on elevated slopes is less clear (Fig. 2). A patchy distribution of subglacial sediment would lessen the long-term effects of sediment weakening and ice flow acceleration, in favor of self-regulation of ice flow and long-term deceleration. Complexity can be further compounded by the transient evacuation of supraglacial lake–derived waters through efficient turbulent sheets (31) or pressurized water layers at the ice-bed interface, by the temporary modulation of subglacial drainage by extreme melt and rainfall events (37), through the regulation of subglacial drainage by weakly connected regions of the bed composed of cavities and patchy subglacial sediments (38), through the systematic erosion of subglacial sediments by hydraulic interactions between channelized and distributed subglacial water systems (14, 15, 20, 3942), or through enhanced storage of water in and lagged release from subglacial sediments, thus modifying its ability to deform and further increasing complexity in ice flow response. To enable realistic simulations of future ice sheet flow, numerical models must therefore be able to reproduce the complex trade-offs between self-regulation and sedimentary control.

It follows that future work should focus on mapping the spatial extents, thicknesses, and physical properties of sediments beneath the GrIS, on appropriate model developments and implementations of subglacial sediments, and on testing rigorously the ensuing trade-offs between self-regulation and sedimentary control of ice flow. Our study has focused on the effects of supraglacial lake drainages in a prominent land-terminating area of the GrIS. Future work will therefore also need to ascertain the degree to which daily variability in subglacial runoff affects sediment strengths and ice flow and on how this variability and lake drainages act together to control present-day and future ice flows of both land- and marine-terminating glaciers.


Seismic data acquisition and basal imaging

We acquired three 2D seismic reflection profiles using a Geometrics Geode–based acquisition system. Seismic sources were 250-g Pentolite explosive charges, frozen in ~3-m-deep holes augered at 80-m intervals. Receivers of seismic energy were 48 × 100 Hz vertically orientated geophones mounted on predrilled 15-kg concrete slabs embedded at 10-m intervals in the ice surface. Although our data acquisitions were optimized for seismic analysis of the nature and physical properties of the substrate (below), they also imaged the structures of the basal environment along a central portion of each of the three profiles (Fig. 2). All data were prestack Kirchhoff-migrated using a 2D velocity model derived from integrated semblance-based (43) and migration velocity analysis (44) and were depth-converted using an average ice seismic velocity of 3800 ± 40 m s−1. Elevations, varying by less than 2% of local ice thickness, were recorded using a Leica SR250 differential GPS system (with vertical precision of ~0.1 m) and were used to define static corrections in seismic data processing (17).

Seismic AVA processing

AVA methods are increasingly applied to ice-substrate characterization (1618, 4548). Variations in the seismic reflectivity of a horizon with incident angle, described by the Knott-Zoeppritz equations, are diagnostic of the material properties on either side of that horizon. The AVA response (Fig. 3A) of a horizon is controlled by contrasts in two diagnostic quantities: acoustic impedance, the product of a material’s density and its compressional wave velocity, and Poisson’s ratio, a function of compressional and shear-wave velocities. A negative reflection coefficient observed at zero incidence implies that the acoustic impedance of the overburden exceeds that of the substrate. If the AVA response exhibits a positive gradient, the Poisson’s ratio is elevated in the substrate, suggesting greater water saturation. An increasingly negative zero-incidence reflection coefficient coupled with an increasingly positive AVA gradient thus suggests a more slippery ice substrate (34), associated with sediment becoming increasingly soft, porous, and eventually dilatant (16, 48). Acoustic impedances less than ~3.4 × 106 kg m−2 s−1 are consistent with high-porosity (>40%) dilatant sediment, nondeforming lodged sediment for values of ~3.6 × 106 to 5 × 106 kg m−2 s−1, and at least partially lithified rock for values greater than ~5 × 106 kg m−2 s−1 (34). Dilatant water-saturated subglacial sediment has a Poisson’s ratio close to the theoretical maximum of 0.5 for water; dry uncompacted sediment has a Poisson’s ratio of ~0.15 or potentially even lower, and that of ice and consolidated rock is between these limits (16). Seismic amplitudes were treated carefully during data processing to compensate for geometrical spreading and anelastic losses (49) and also to correct for differences in shot strength and receiver coupling (16, 17). Observed amplitudes were normalized by the reflectivity observed for seismic incidences less than 10° (49) and were collated into angle bins of 2° incidence to improve signal-to-noise ratios. The resulting AVA data are shown in Fig. 3 (B to D); error bars show the uncertainty in both angle and reflectivity within each bin. The best-fit AVA curve to each data set (red; Fig. 3, B to D) was evaluated using Bayesian modeling, assuming fixed ice density (920 kg m−3) and compressional (3800 ± 40 m s−1) and shear (1898 ± 60 m s−1) wave velocities (48). With a set of observations and their associated uncertainties, here the distribution of AVA responses in each 2° angle bin (16), Bayesian analysis assesses the probability that a given model is explained by those observations. The best-fit AVA curve was therefore weighted by the uncertainty distribution in the observed data and corresponded to the most probable model of acoustic impedance and Poisson's ratio based on the observed AVA responses. Substrate properties were assumed to be constant over the range of the glacier bed (200 to 400 m) that we investigated, but this is not a major limitation in practice because the horizontal resolution (first Fresnel zone) is ~165 m.

Description of the ice sheet model

The model experiment was conducted using the higher-order thermodynamic and 3D Community Ice Sheet Model (CISM), coupled to subglacial sediment and subglacial hydrology models (3). CISM and the sediment model were coupled through the porosity-controlled basal shear strength (50). The latter was used to calculate the basal stress in the force balance equation solved by CISM, assuming a plastic yield stress basal boundary condition (51). Perturbations to the ice flow velocity can therefore only occur in response to changes in sediment water content and strength. We hypothesized that the sediment shear strength evolved in space and time after supraglacial lake waters reached the ice sheet bed, whence they were routed in the subglacial hydrological model according to hydraulic potential gradients. Large water volumes transiently delivered to the basal environment caused perturbations and thus short-lived vertical hydraulic gradients at the ice-sediment interface, thus initiating subglacial water flow into the sediment layer. After the water drained away, vertical hydraulic gradients reversed and allowed water to flow out of the sediment layer. Only about 10% of the total available water was assumed to infiltrate into the sediment, whereas the remaining 90% were assumed to drain away in an efficient subglacial hydrological system, thus allowing for flow through subglacial channels or turbulent sheets together with any associated uplift of the ice. The evolving sediment shear strength was calculated according to changing water content and assuming, in the absence of direct observations beneath the GrIS, sediment properties similar to tills found beneath Trapridge Glacier (52) and glaciers from other Arctic regions (53).

Application of the ice sheet model

The model was run at a 1-km resolution. The surface geometry was prescribed using the 2008 SPOT surface digital elevation (54) and a 500-m bed digital elevation model (DEM) (28). Initial conditions were obtained by performing a model inversion (55) of a composite image of the 2009–2010 winter velocity, achieving a very good match between model and observations (r2 = 0.99) (3). The ice flow model was driven by the record of supraglacial lake drainages for the Russell Glacier catchment, derived from MODIS imagery and comprising 663 discrete discharge events occurring during summer 2010 (30). The model was validated by comparing outputs with available observations (3). Modeled seasonal ice flow evolution showed good agreement with TerraSAR-X velocity snapshots [acquired with 11-day separation and centered on 19 June, 22 July, and 11 November 2010 (5658)], with an overall net error of 10% and correlation coefficients of 0.79 to 0.94. In addition, the modeled daily variation in velocity at SHR site (located ~15 km from the margin; Fig. 1A) was within 16% of the velocity recorded from 1 June to 31 August 2010, with a correlation coefficient of 0.83.


Supplementary material for this article is available at

fig. S1. Location and dates of drainage of supraglacial lakes in the Kangerlussuaq sector of the GrIS in the 2010 melt season.

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: Landmark Strategic University Alliance grants ProMAX access to the University of Leeds, under agreement 2004-COM-024982. Shot boxes were loaned by S. Anandakrishnan (Pennsylvania State University) and SEIS-UK at the University of Leicester. The facilities of SEIS-UK are supported by Natural Environment Research Council (NERC) under agreement no. R8/H10/64. Knott-Zoeppritz equations were calculated using source code from CREWES ( Funding: We acknowledge major support by UK NERC grants NE/H0126889/1, NE/G007195/1, and NE/G00692X/1; Greenland Analogue Project: Sub-Project A, funded by SKB Posiva and Nuclear Waste Management Organization (Sweden); Leverhulme Trust Research Leadership project F/00391/J and the Climate Change Consortium of Wales (C3W) that funded A.D.B.; an Aberystwyth University doctoral scholarship that funded S.H.D.; and a NERC doctoral training grant awarded to Swansea University that funded C.F.D. A.L.H. acknowledges a professorial fellowship from the Research Council of Norway through its Centres of Excellence funding scheme (grant 223259). Author contributions: B.K. wrote the manuscript and led NERC project NE/H012869/1, and A.L.H. led NERC project NE/G007195/1 and Greenland Analogue Project: Sub-Project A. A.D.B., B.K., S.H.D., A.L.H., and G.A.J. acquired the seismic and associated GPS data, processed respectively by A.D.B. and S.H.D. M.B., P.C., and C.F.D. conducted the numerical and hydrological modeling. K.L., R.P., A.L.H., and S.H.D. acquired and processed the ground-penetrating radar data that allowed the generation of the basal DEM. A.A.W.F. calculated daily supraglacial lake volumes and drainage discharges and critically guided the field site selection. All authors contributed to the glaciological data interpretation and writing. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The seismic data are available on request from the corresponding author, B.K. (b.kulessa{at}
View Abstract

Stay Connected to Science Advances

Navigate This Article