Ice-water-sediment interactions beneath RG

The physical properties in the subglacial sediment model integrates recent geophysical observations, revealing where available that RG is underlain by a porous, mechanically weak sediment34,35, of similar character to tills produced by glaciers in Canada and Svalbard37,45 (see Methods). The subglacial sediment model was coupled to a hydrological model, which has previously been used to calculate the routing and fluxes of water associated with the episodic drainage of subglacial lakes in Antarctica46, and is well suited for analysis of SGL drainage events. The explicit inclusion of interacting models of subglacial sediment and water is new, yet consistent with the extremely high suspended sediment loads observed in proglacial streams draining RG catchment47,48. The latter equates to bulk catchment erosion rates of 4.8 mm a−1 (ref. 48), which is far greater than other rates estimated for regions where glaciers override a hard crystalline rock (0.004–0.1 mm a−1, refs 49, 50).

Current understanding of SGL dynamics suggest that a rapid, high-magnitude influx of water to the bed (considering peak fluxes as high as 5,000 m3 s−1, ref. 10) cannot be instantaneously accommodated by expansion of a channelized basal drainage system51. Rather, the accelerated flow observed in late summer, when an efficient drainage system had already developed12,28,51, points to subglacial evacuation of water in a high-pressure system5,30. Here, we assume that meltwater delivered to the bed is transported down the hydraulic potential surface in an efficient basal water system, which coexists and interacts with a hydraulically inefficient subglacial sediment layer. The physics controlling the rate of water infiltrating the sediment layer is likely complex, but with a paucity of data to constrain this interaction, we use a simple parameterization in which basal water is transferred into the underlying sediment in proportion to the magnitude of the horizontal water flux associated with the hydrological forcing (See Methods). This perturbation induces excess pore pressure and vertical hydraulic gradients, and thus flow of water within the sediment, which increases its porosity accordingly. This causes sediment shear strength to drop, along with basal traction and resistance to ice flow (equation (1)). The expansion of sediment pore space is represented in our model through compressibility, a well-established material property, which has been determined for a variety of subglacial tills37,52.

In this manner, we assessed the time-varying hydrological impact of seasonal meltwater delivery on subglacial sediment shear strength, thereby defining patterns of basal traction across the model domain at a horizontal resolution of 1 km (see Methods). The model was initialized by performing a data inversion on a composite image of the winter 2010 observed velocity (Fig. 1; Table 1, see also Methods). We then performed a series of forward experiments for a 6-month period starting on the 15 May 2010 (Methods). First, we investigated the impact of SGL drainage events by forcing the model with the 2010 record of drained SGL volumes10 (Fig. 1a; Supplementary Movie 1). Second, we examined the effect of the runoff produced daily in each SGL sub-catchment for that same year43. Third, we tested the possibility that ice flow responds to hydrological perturbations defined by the variability in runoff together with SGL drainage events (Supplementary Movie 2).

Figure 1: Supraglacial lakes and seasonal ice flow in the Russell Glacier catchment. (a) Colour-coded summer 2010 lake-drainage map, with elevation contours (m). (b) Composite winter 2010 observed velocity map (m per year). (c–e) Velocity maps (m per year) derived from TerraSAR-X satellite data, showing the average velocity over 11-day periods centred over the date indicated above each panel. (f) Modelled initial ice flow for winter 2010 (m per year). (g–i) Modelled ice flow (m per year) averaged over the same 11-day periods used to generate TerraSAR-X velocity maps as shown in c–e. The locations of Russell Glacier (RG), Isunngata Sermia Glacier (ISG), Orkendalen Glacier (OG) and SHR site are indicated in f. Full size image

Table 1 Comparison of remotely sensed data and model output. Full size table

Seasonal ice flow driven by observed SGL drainage

The 2010 record of SGL drainage volumes was first used to drive the model. Patterns of surface velocity derived from TerraSAR-X satellite image pairs acquired with 11-day separation and centred on 19 June, 22 July and 11 November53,54 are shown in Fig. 1 together with maps of modelled surface velocity averaged over the exact same periods. Flow acceleration, both observed and modelled, is most apparent on 19 June, ~30 km from the margin of RG, and for Orkendalen Glacier, south of RG (Fig. 1). On 22 July, modelled and observed surface velocities have declined substantially, and by 11 November, both have approached the previous winter mean. Modelled seasonal flow evolution is in overall good agreement with the TerraSAR-X velocity snapshots, with net errors of 10% and correlation coefficients of 0.79–0.94 (Table 1).

To assess model efficacy and its ability to capture the dynamic events observed, continuous model output was compared with daily mean surface velocity measured by a Global Positioning System (GPS) receiver at the SHR site (Fig. 1), ~15 km from the margin of RG (Fig. 2a). All the main velocity features at SHR are captured in the model; throughout June, July and August, modelled mean daily velocity at SHR was generally well estimated (within 16%), with a good correlation to observed values (r2=0.83, Supplementary Table 1). On 10 June, modelled velocity at SHR attained a maximum of 326 m per year, within 6% of that measured by GPS. Likewise, on 25 June modelled velocity peaked at 248 m per year, within 15% of that observed. Noteworthy however, is that the model did not initially reproduce the first ‘spring-event’ acceleration in late May, as fluxes calculated from SGL drainage events were insufficient to perturb subglacial conditions at SHR (Supplementary Fig. 1). To reproduce this first spring event, the lake volume loss on 24–27 May needed to be three times those shown in Fig. 2b. The model may be unable to sufficiently respond to spring-event stimuli possibly because the ice flow may be more sensitive to water input in the lower ablation zone early in the melt season, when the basal system is not yet fully developed12,51. Another plausible explanation is that basal meltwater produced in situ by frictional and geothermal heating accumulates at the bed over the course of winter and is released together with the first SGL drainage events.

Figure 2: Modelled versus observed daily mean ice flow speeds at site SHR and daily water forcing estimated for RG catchment. (a) Timeseries of observed mean daily velocity acquired with Global Positioning System at the SHR site (black), and comparison with model output at the same location, when forced with SGL-only volumes (red), with SGL volumes plus runoff rates (blue) and with absolute runoff volumes (grey, RHS scale). The speed-up event on 24–27 May is reproduced assuming a threefold increase in SGL volume loss (red dashed line). This can be explained, as in reality, the ice flow may be more sensitive to water input early in the season, when the basal water system is not yet fully developed. Alternately, the observed ice flow may have responded to SGL water loss combined with the basal meltwater produced over winter and released early in the ablation season. The shaded red zone corresponds to uncertainties in observing the timing of SGL drainage on 10 June and 25 June, due to cloud cover. (b) SGL volume loss data used to drive the model. For periods of time where no satellite data were available (horizontal bars), the timing of drainage is centred over that period10. (c) Daily runoff rates (blue bars), calculated from the total runoff volume estimated for 2010 (ref. 43) (grey line). The daily mean water input rates (RHS scale, cyan dots in b,c) are highest for SGL volume loss. Full size image

Seasonal ice flow driven by observed SGL drainage and runoff

We isolated and tested the extent to which SGL drainage controls seasonal ice flow by performing experiments in which the model was driven by local runoff production delivered to the bed in each SGL sub-catchment (see Methods). Using total daily runoff volumes43, we found that modelled ice flow was fastest in late July (Fig. 2a; Supplementary Figs 2 and 3), an outcome inconsistent with the relatively slow flow observed at this time (Figs 1 and 2a)12,54. This discrepancy is not surprising since both theory and observations demonstrate that it is the variability rather than the absolute volume of meltwater delivery to the bed that drives ice flow dynamics19,27,28,51,55. Indeed, a reduction in the daily variability of surface melt at this peak period of runoff feasibly explains why observed ice flow remained slower from early July and onwards28, a reasoning consistent with the existence of an efficient system capable of evacuating large volumes of water. Our total runoff experiment supports this interpretation. The model was unable to reproduce seasonal ice motion when assuming that the sediment layer takes in all of the runoff, which suggests that some of this water is instead integrated in an existing efficient hydrological system. Accordingly, we performed a final suite of experiments in which the model was forced with daily meltwater perturbations, including the estimated daily difference in runoff in addition to the SGL drainage volumes. The underlying assumption for these experiments is that it is the sum of positive daily increases in runoff at each sub-catchment and the SGL volume that collectively represents the net hydrological forcing to which the ice flow responds (equation (7)), and that the remaining quantity of water is routed to the margin without effect (as previously inferred18,19,28,51,55). Although the forcing volume of water was over three times greater than in the SGL-only drainage experiments (Fig. 2c, see Methods), the model was still able to successfully yield the observed flow structure (Fig. 2a; Supplementary Figs 2 and 3), showing no appreciable difference to the SGL-only runs (Table 1; Supplementary Table 1). This limited impact stems from the fact that daily differences in runoff calculated at each lake site (henceforth referred to as runoff rates) were smaller than the typical volume of water contained in lakes observed to drain (Fig. 2b,c). These results suggest that the hydrological forcing associated with the runoff variability alone was not capable of inducing a substantial ice flow response in 2010 (Supplementary Fig. 3).

Impact of SGL discharge on basal sediment properties

In our model, water at the ice–sediment interface travels 20–70 km down glacier in 1 day, an efficient routing that is consistent with SF6 gas tracing experiments that estimate subglacial water velocities from 22 to >86 km per day (ref. 22). As water travels down the hydraulic potential surface, it interacts with the hydraulically inefficient subglacial sediment layer below it. On 10 June, modelled water fluxes locally attain 570 m3 s−1 (Fig. 3a) and average 57 m3 s−1 over a distance of ~65 km. The volume of water entering the sediment was 3.8 × 106 m3, which is ~10% of the lake water available on that date. The remaining ~90% is routed away by the basal hydrological system. This is consistent with observations, which demonstrate that the majority of water delivered to the bed travels to the ice margin within a few days10,22,28,43. Nevertheless, local sediment shear strength and basal traction fell to about a third (~50 kPa) of the pre-discharge value (Fig. 3c), as the porosity of the sediment increased with the water intake (equation (1) and Supplementary Fig. 4). This led to a substantial (up to 200%) surface flow acceleration in a region considerably larger than that of the local meltwater production and input, and extending as far as the ice sheet’s margin (Fig. 3b,d). Furthermore, an equivalent sediment strengthening took place over the following days (12–13 June, Fig. 3e), as pore-water pressure within the temporarily expanded sediment layer returned to its post-pulse equilibrium driving a net flow deceleration of up to 75% (Fig. 3f).

Figure 3: Modelled dynamic impact of SGL drainage on 10 June 2010. (a) Bed elevation (m a.s.l., grey scale) overlain with modelled water flux (m3 s−1, colour scale) from lakes draining on 10 June (red and white dots, v loss ). (b) Modelled (pre-discharge) mean daily velocity on 10 June (m per year). (c) Modelled basal shear stress (kPa, grey scale) overlain with basal stress reduction calculated on 11 June (kPa, colour scale) and (d) associated speedup (%). (e) Modelled basal shear stress (kPa, grey scale) overlain with increase in basal stress calculated on 13 June (kPa, colour scale), relative to that calculated on 11 June, and (f) associated slowdown (%). The location of transects (T lower , T upper ) used to calculate velocities in Fig. 4 is shown in b. Full size image

Sensitivity to surface water inputs

In all our experiments (Supplementary Fig. 3), modelled ice flow at the end of the season was slower than that at the start (winter 2009/10). Although this seasonal velocity reduction may be subtle (Fig. 2a; Supplementary Fig. 3), it is a consistent feature of the model’s sensitivity to surface meltwater delivery to the bed (Fig. 4; Supplementary Fig. 5). This modelled self-regulation of ice flow is consistent with recent findings from GPS measurements across the ablation area of RG, revealing that flow enhancement from increased surface melting during warm summers is negated by reduced winter velocities17,18. To date, the latter has been exclusively attributed to increased effective pressure provided by summer expansion of subglacial channels over a hard bed19,20,21,28,51. Here, model results provide an alternative explanation for the soft-bed condition. Summer ice flow is greater with higher SGL discharge because higher hydraulic gradients drive a proportionally larger volume of water into the subglacial sediment layer (Methods, equations (2) and (3)). Upon evacuation of surface meltwater, reversed but equally high hydraulic gradients develop and drive water out of the sediment. With sufficiently high gradients, the new state of water pressure equilibrium in the sediment is attained with an overall sediment strengthening compared with its pre-summer state, and therefore lower velocities in autumn.

Figure 4: Modelled ice flow sensitivity to variation in surface water inputs. Modelled velocity along the transects T lower (red) and T upper (black), using SGL-only volumes (solid lines) and runoff rates in addition to SGL volumes (dashed lines). See Fig. 3 for transects location. (a) Change (%) in mean summer velocity (open squares) and winter velocity (coloured squares). The latter is assumed to be represented by the end-of-run (15 November) value. (b) Mean annual velocity change (%). Full size image