Examining the effect of ice dynamic changes on subglacial hydrology through modelling of a synthetic Antarctic glacier

Abstract

Hydrologic pathways beneath ice sheets and glaciers play an important role in regulating ice flow. Antarctica has experienced, and will continue to experience, changes in ice dynamics and geometry, but the associated changes in subglacial hydrology have received less attention. Here, we use the GlaDS subglacial hydrology model to examine drainage evolution beneath an idealised Antarctic glacier in response to steepening ice surface slopes, accelerating ice velocities and subglacial lake drainages. Ice surface slope changes exerted a dominant influence, redirecting basal water to different outlet locations and substantially increasing channelised discharge crossing the grounding line. Faster ice velocities had comparatively negligible effects. Subglacial lake drainage results indicated that lake refilling times play a key role in drainage system evolution, with lake flux more readily accommodated following shorter refilling times. Our findings are significant for vulnerable Antarctic regions currently experiencing dynamic thinning since subglacial water re-routing could destabilise ice shelves through enhanced sub-shelf melting, potentially hastening irreversible retreat. These changes could also affect subglacial lake activity. We, therefore, emphasise that including a nuanced and complex representation of subglacial hydrology in ice-sheet models could provide critical information on the timing and magnitude of sea-level change contributions from Antarctica.

1. Introduction

Pressurised subglacial water plays a central role in the motion of Antarctic glaciers (Bell and others, 2007; Stearns and others, 2008; Smith and others, 2009; Wright and Siegert, 2012; Bougamont and others, 2019). Dynamic (i.e. velocity) and geometric (i.e. surface slope) changes impact the development and evolution of subglacial drainage systems and changes in subglacial drainage feed back into glacier dynamics and geometry. Understanding the interactions between ice dynamics and subglacial hydrology is especially important to Antarctic regions that are vulnerable to ocean warming. For instance, beneath the ice shelves of Pine Island, Thwaites and Totten glaciers, warm, modified Circumpolar Deep Water is accessing ice shelf cavities and grounding lines, where it is enhancing basal melt rates (Depoorter and others, 2013; Rignot and others, 2013; Roberts and others, 2018). Sub-shelf melting thins the ice shelf, which increases the surface slope between the grounded ice and the ice shelf, resulting in accelerated flow velocities and an increase in the amount of ice discharged into the ocean, a process known as dynamic thinning (Pritchard and others, 2008). Recent modelling efforts to reproduce Antarctic glaciers observed dynamic thinning behaviour (Joughin and others, 2019) and predict their future evolution (Brondex and others, 2017, 2019) have highlighted that the representation of basal conditions strongly influences the obtained results. However, in these simulations, the subglacial drainage system was coarsely represented. Features such as basal channels and cavities, which have contrasting impacts on the rate of basal sliding (Flowers, 2015), were not included.

Geophysical data have provided unparalleled detail of the Antarctic subglacial landscape and informed modelling studies. For example, subglacial channels extending hundreds of kilometres inland have been inferred beneath Thwaites (Schroeder and others, 2013; Hager and others, 2022) and Totten glaciers (Dow and others, 2020) and draining across the grounding lines of the Getz (Wei and others, 2020) and Filcher-Ronne ice shelves (Dow and others, 2022). Modelling studies suggest that Antarctic subglacial channels operate at relatively high-pressure, are stable over seasonal timescales and are coincident with regions of observed fast flow (e.g. Dow and others, 2018, 2020, 2022; Hager and others, 2022). All these characteristics contrast with Greenland's low-pressure, seasonally evolving channels that can considerably slow ice flow (Flowers, 2015; Davison and others, 2019). Furthermore, the geophysically derived data show elevated basal ice shelf melt rates that coincide with the modelled subglacial channel outlet locations (Wei and others, 2020; Dow and others, 2022). The high basal melt rate, combined with the fact that the channels can influence water pressure within a 100 km radius (Dow and others, 2022), suggests that changes in basal hydrology could have critical and widespread impacts on ice dynamics. This finding is especially salient since the regions in the abovementioned studies are candidates for the marine ice-sheet instability (Schoof, 2007; Ross and otherd, 2011; Joughin and others, 2014; Greenbaum and others, 2015; Seroussi and others, 2017), which describes the runaway retreat of glaciers resting on reverse sloping beds when buttressing provided by ice shelves is removed (Schoof, 2007).

Changes in either subglacial hydrology or ice dynamics can initiate changes in the other, with consequences for glacier stability. In Antarctica, increased flow velocities of these glaciers would lead to increased volumes of basal meltwater production, potentially leading to greater water fluxes across the grounding line that further enhance the already-elevated ice shelf melt rates. However, ice-sheet modelling experiments that produce time-varying ice geometries, dynamics and basal melt rates in response to climatic changes generally use simplified representations of subglacial drainage, not differentiating between efficient and inefficient systems (e.g. Paxman and others, 2020; DeConto and others, 2021).

In addition to increasing meltwater production, accelerated flow velocities would cause dynamic thinning at the terminus, which increases surface slopes and driving stress (Scott and others, 2009; Cuffey and Paterson, 2010; Flament and Rémy, 2012). Changes in hydropotential gradients because of steeper surface slopes could affect the amount of ice–bed coupling and the stability of ice flow. Although projections indicate that dynamic thinning will likely continue in the future (e.g. Golledge and others, 2015, 2021; Seroussi and others, 2014, 2017), it remains to be determined whether velocity changes or surface slope changes will exert the dominant control on subglacial hydrology. Here, we perform experiments investigating how increases in ice velocity and surface slope impact subglacial drainage evolution. We then consider how these processes alter frictional melt rates and the resulting consequences for the subglacial drainage network.

Subglacial lake drainage can also trigger variability in ice dynamics. In an Antarctic lake inventory updated by Livingstone and others (2022), an additional 10 subglacial lakes were classified as ‘active’, bringing the total number of known systems to 140. Active subglacial lakes are defined as lakes that drain and fill repeatedly and have been detected using satellite altimetry data (e.g. Fricker and others, 2007; Siegfried and Fricker, 2018). During lake filling, the ice surface elevation rises. As the lake drains, the ice surface elevation lowers and, in some case the downstream ice flow speed increases (e.g. Stearns and others, 2008; Siegfried and Fricker, 2018). Numerous studies have indicated that when an upglacier lake ‘activates’, or drains, it can trigger a cascade of lake drainages and catchment-wide flow acceleration (Bell and others, 2007; Stearns and others, 2008; Hoffman and others, 2020; Malczyk and others, 2020). Two cycles of cascading lake drainage and refilling were recently detected on Thwaites Glacier (Hoffman and others, 2020; Malczyk and others, 2020) and the drainage period of these water bodies could have important consequences for the future flow stability of this vulnerable glacier. If lake drainage cycles occur more frequently, it could promote ice flow self-regulation and flow stability as lake fluxes appropriate remnant channels formed during previous lake drainage events (Livingstone and others, 2022).

This study investigates the roles of velocity changes, glacier surface slope changes, their combined effects on basal meltwater production and subglacial lake drainages on the development of the subglacial drainage network in an idealised Antarctic glacier using the Glacier Drainage System model (GlaDS; Werder and others, 2013), which includes co-evolving channelised and distributed drainage. Determining the level of control that ice dynamic changes have on the hydrology of a simulated glacier system is an essential first step towards understanding their influence on real subglacial drainage networks.

2. Methods

2.1. GlaDS model

The GlaDS subglacial hydrology model is a 2D finite-element model that incorporates both channelised and distributed drainage and captures the time-evolving interactions between these two forms of drainage systems. Fundamental model equations are summarised here; a complete description of the model is available in Werder and others (2013).

Distributed drainage approximates a linked cavity system and is represented as a continuous water sheet that occupies the domain elements. A Darcy–Weisbach turbulent flow relation represents discharge (q, in m²s^-1) through the distributed system:

(1)$$q = -kh^\alpha\vert \nabla{\phi}\vert ^{\beta}\nabla\phi$$

where k is the distributed system's conductivity (in m^7/4kg^-1/2), h is the thickness of the water sheet in metres, α and β are turbulent flow parameters (set to 5/4 and 3/2, respectively) and $\nabla \phi$ is the hydraulic potential gradient.

The evolution of the sheet over time, t, is given by:

(2)$${\partial h\over \partial t} = {u_{\rm b}( h_{\rm r} - h) \over l_{\rm r}} - Ah\vert N\vert ^{n-1}N + {\Xi_{\rm s}-\Gamma_{\rm s}\over \rho_{\rm i}L}$$

where u _b is the basal sliding speed in ms^-1. h _r is the typical bedrock obstacles size and l _r is the typical cavity spacing; both in metres. A is an ice flow constant (in Pa^-ns^-1), N is the effective pressure (Pa) and n is Glen's flow law exponent (unitless). $\Xi _{\rm s}$ and Γ_s represent the viscous energy dissipation and pressure melt term for the sheet, respectively, and are in units of Wm^-1. ρ _i is the density of ice (in kgm^-3) and L is the latent heat of fusion (units: Jkg^-1).

The first term in Eqn (2) describes the cavity opening rate, while the second term describes creep closure. The last term in Eqn (2) represents changes in viscous dissipation of heat and sensible heat changes for the sheet and were not included in original sheet evolution equation in Werder and others’ (2013) formulation of GlaDS due to the potentially unstable growth of the sheet concentrated over a narrow region. However, Dow and others (2018a) found that runaway growth does not occur and therefore incorporated these effects.

Our simulations were performed using an idealised glacier setup (Fig. 1), consisting of ~20 000 nodes with a domain length of 1500 km and width of 600 km, similar to the synthetic glacier dimensions adopted by Dow and others (2016). We use ‘Relative Northing’ and ‘Relative Easting’ to represent the domain coordinates of our idealised glacier. The glacier lies on a curved, inclined bed with a square root ice surface profile that thickens from a 1000 m ice cliff at the terminus to a maximum ice thickness of 3500 m at the upper boundary. The ice thicknesses at the upper and lower boundaries resemble that of Dow and others’ (2016) idealised Recovery Glacier, and our square-root surface profile is similar to theoretical ice surface profiles discussed in Cuffey and Paterson (2010). Channels are not permitted to form on the lateral boundaries of the domain. We impose a Dirichlet boundary condition at the terminus, which describes the outlet pressure condition, specifying that the ice is at overburden because it is floating. We apply Neumann flux conditions at the upper and lateral boundaries, with water supplied to the distributed system at the upper boundary, assumed to be generated from the wider upglacier catchment, and zero-flux conditions at the lateral boundaries.

Figure 1. Synthetic mesh (a) and topography (b) used in this study designed to emulate an Antarctic outlet glacier.

Channelised drainage develops on element edges, with water exchange between channels occurring at the nodes, while water exchange with the sheet occurs over the edges. The evolution of the channel cross-sectional area (S, in m²) is given by:

(3)$${\partial S\over \partial t} = {\Xi-\Gamma\over \rho_{\rm i}L}- AS\vert N\vert ^{n-1}N$$

Similar to above, but for the channelised system, $\Xi$ represents the viscous dissipation of potential energy and Γ is the pressure melt term. Channels are not prescribed at the beginning of the model run but instead develop and co-evolve in concert with changes in the water pressure and flux of the distributed system.

For each of the scenarios described below, we first run the model to steady state, close to overburden pressure (mean pressure for each scenario >84%). Consistent with previous GlaDS simulations (Dow and others, 2016, 2018b, 2020), we consider steady state to have been reached when the difference in distributed sheet thickness at the final two timesteps is less than 3.33 × 10⁻⁷ mday^-1. We run the model for another 100 years from this state for the slope and velocity change runs, and 150 years for the lake drainage runs, performing sensitivity tests of internal GlaDS parameters to assess the variables responsible for transmitting ice dynamic changes to the subglacial network. The range of values in our sensitivity tests (Table 1) straddles the typical values used in previous observation-informed GlaDS modelling (Dow and others, 2016, 2018b, 2020, 2022; Poinar and others, 2019; Indrigo and others, 2021). We focus on the sheet and channel conductivities, bedrock obstacle size and water input into the distributed system. Channel conductivity is related to the roughness of the channel walls, with smaller channel conductivities representing rougher channels that inhibit channel growth. Distributed system conductivity measures the ease with which water can flow through the linked cavity system and is constrained to a smaller range than the other parameters, outside of which the model fails to run to completion. Setting the conductivity value too low results in the system becoming over-pressurised. If conductivity is too high, the system is prevented from fully pressurising and remains in steady state (Dow and others, 2016, 2020; Poinar and others, 2019). The typical bedrock obstacle size, h _r, controls the evolution of the distributed system by influencing cavity opening rates (first term in Eqn (2)), with larger (smaller) bedrock bumps resulting in bigger (smaller) cavities. The bedrock bumps’ size also impacts the system's pressurisation rate, as more (less) water is required to pressurise larger (smaller) obstacle sizes. It follows that larger (smaller) bedrock bumps prolong (shorten) the time required to pressurise the system. We also test the water input rate into the distributed system, representing meltwater produced by geothermal and frictional heat. We include this test since large uncertainties in geothermal heat flux, which exceeds 50% in some areas (Fox-Maule and others, 2005; Burton-Johnson and others, 2020), result in poorly constrained basal water production rates. The results presented here in the main text are from the water input sensitivity test. Model outputs from the bedrock obstacle size, and the distributed and channelised conductivity tests are included in the supplementary material (Figs S1–S6).

Table 1. Model parameters and values used in GlaDS sensitivity tests

Constraints on all the tested model parameters require high-resolution observations of the basal environment. The lower and upper values for each parameter represent the minimum and maximum values that permit the simulation to run to completion. For instance, spatially uniform water production rates above 8 mma^-1 were not computationally feasible with our synthetic setup, and we acknowledge that actual basal water volumes in Antarctica may exceed this value in some regions, particularly towards grounding lines. To represent spatially variable melt rates, we also include simulations where meltwater production is frictionally derived (see ‘combination run’, discussed below). However, since we are modelling an idealised system, we cannot state which of the tested parameter values would be most appropriate under the imposed dynamic changes. Nevertheless, these sensitivity tests allow us to investigate the interactions between model parameters and ice dynamic changes and their effects on the subglacial drainage network, and we include a discussion of these interactions in section 3.4 of the results section, titled "Sensitivity tests".

2.2. Dynamic forcing experiments

Our model simulations represent a one-way dynamic forcing (velocity and slope changes) on the subglacial drainage system, and are not coupled to an ice-sheet model. We refer to the collection of these forcings as ‘dynamics runs’ due to changes in boundary conditions (described below), which differ from the baseline run in which the boundary conditions are time-invariant.

2.2.1. Baseline run

To facilitate the comparison between dynamics runs, we consider a baseline run. The outlet pressure is always at 100% of ice overburden pressure, indicative of fully floating ice, and the basal velocity field is spatially- and temporally constant and set to 100  ma^-1. With this setup, the ice surface profile described above remains unchanged throughout the simulation time.

2.2.2. Linear velocity

Basal sliding speed is a primary control on the opening rates of the inefficient cavity system in GlaDS and in other subglacial hydrology models (De Fleurian and others, 2018). Basal velocity impacts the size of the cavities and the pressurisation rate. With low velocities, cavity opening rates are low, resulting in smaller cavities that require less water to become pressurised. In the opposite scenario, high velocities result in larger cavities that require more water to become pressurised. To determine how sensitive subglacial drainage networks are to increases in basal velocity, we represent the velocity as a spatially variable field: the initial velocity is 100 ma^-1 and is highest at the centre of the grounding line, then decreases over the domain following a Gaussian distribution from this centre point. We increase the initial spatially variable velocity by 4%a^-1, such that, at the end of our 100-year simulation, the velocity at the grounding line is 500 ma^-1. We hereafter refer to these simulations as the ‘linear velocity’ runs. We do not consider any seasonal changes in basal velocity since surface meltwater does not presently contribute to seasonally modulated changes in ice velocity in Antarctica.

2.2.3. Surface slope changes

Determining the temporal and spatial development of efficient and inefficient drainage systems in response to changes in surface slope has not yet been explored in GlaDS. Here, we apply spatially variable thinning rates to our idealised glacier. We adopt a 1 ma^-1 thinning rate at the centre point of the grounding line, which decreases linearly with distance from the grounding line. This rate of thinning falls within the range observed on Totten Glacier (Flament and Rémy, 2012; Li and others, 2015) and glaciers on the Antarctic Peninsula (Wouters and others, 2015). The grounding line position is fixed in all the simulations mentioned above, hereafter referred to as the ‘slope change’ model runs.

2.2.4. Combination run

Our next set of simulations assesses how slope and velocity changes impact basal meltwater production. We base our melt rate calculations on the basal mass balance, E _b, equation from Cuffey and Paterson (2010):

(4)$$E_{\rm b} = G + u_{\rm b}\tau_{\rm b} + K_{\rm T}{\partial T\over \partial Z}$$

where G is the geothermal heat flux (set to 45  Wm^-2), u _b is the basal velocity and τ _b is the basal shear stress. The last term contains K _T, the thermal conductivity of ice and ∂T/∂Z, the temperature gradient in the ice. In our idealised setup, we omit this last term representing heat conduction and assume that the ice is everywhere at the pressure melting point.

Dividing E _b by the product of ice density, ρ _i and latent heat of fusion L _f for ice gives the basal melt rate, which, for these simulations, serves as the distributed water input into the system, replacing the meltwater production sensitivity test:

(5)$$m = {E_{\rm b}\over \rho_{\rm i}L_{\rm f}}$$

We consider two scenarios for temporally and spatially variable melt rates. The first incorporates changes in velocity only, with u _b increasing at 4 %a^-1 as in the linear velocity runs described above. The second simulation includes increases in velocity and surface slope changes, with the latter the same as those used for our slope change runs.

2.2.5. Lakes run

In our lakes simulations, we focus on lake drainages of subglacial lake Thw₁₄₂, located beneath Thwaites Glacier in West Antarctica. Thw₁₄₂ is an ideal test case as it occurs at an ice thickness similar to the thickness halfway up our domain along the centreline. Since Thw₁₄₂ experienced two distinct drainage cycles with different lengths and magnitudes, we can examine the short- and long-term impacts of the drainage interval lengths and fluxes on the drainage system. This GlaDS domain does not easily emulate lake drainage and filling on controlled timelines that we require here, so we apply the flux from Thw₁₄₂ into a channel placed to represent the downstream lip of the lake. We do not model the refilling cycle of the lake. Instead, we set the outlet channel point source input to zero when the lake refills.

For the first 114 years, our model setup for lakes has the same configuration as the baseline run with the standard distributed water input. We started drainage at 114 years since, before this, a high-pressure region occurred at our lake input site. This high-pressure region occurs because the flux from the upper boundary (at 1500 km Relative Easting) passes through the domain, but the water has not yet emptied into the lower-pressure channels. Allowing the high-pressure bulge associated with the upper boundary flux to first pass through the domain first prevents us from drawing erroneous conclusions about the lake's effects on pressure and channel discharge in its immediate vicinity and over the domain.

To investigate the impact of changing draining and refilling times, we perform sensitivity tests in which we halve or double the drainage and refilling times. We also include experiments in which we alter the drainage time but keep the refilling time the same as that measured by Malczyk and others (2020). Figure 9 shows our lake input rate and refilling time sensitivity tests. While we do not attribute any of these cycle length changes to any specific process, they could reasonably be caused by velocity and slope changes (Livingstone and others, 2022), such as those investigated in our other experiments.

3. Results

3.1. Linear velocity

In Figures 2a–c, we present the channel discharge differences associated with the linear velocity run at the final timestep. The mean channel discharge difference over the entire simulation time under standard water input between the baseline and the linear velocity run is −0.0020 m³s^-1. At the grounding line, the difference in the total channelised discharge ranges from −1.85 m³s^-1 in the low water input run to −3.02 m³s^-1 in the high water input run. Less water is therefore carried in the channels in the linear velocity run than in baseline run. Channel growth in the linear velocity run does not keep pace with that of the baseline run, and this becomes more evident as the water input rate increases. Channels in the linear velocity run reach a maximum distance from the grounding line of 435 km in the low water input run and 728  km in the high water input run. There are 9 and 14 outlet nodes with at least 1  m³s^-1 of channelised discharge for the low and high water input runs, respectively, in both the baseline run and the linear velocity run.

Figure 2. Channel discharge, water sheet thickness and per cent overburden differences between the linear velocity and the baseline runs at the final timestep under different volumes of water input into the distributed system. The panels show (a)–(c) channel discharge difference; (d)–(f) water sheet thickness difference; (g)–(i) per cent overburden difference. Positive differences in blue indicate where the value is larger in the linear velocity run, while negative differences are shown in red and indicate where the value is larger in the baseline run.

The water sheet in the suite of linear velocity runs is generally thicker than that of the baseline run (Figs 2d–f). This overall pattern is represented by positive sheet thickness differences and is consistent across simulations, albeit the magnitude and spatial extent of change become more pronounced as water input increases. The resulting pattern occurs since the evolution of the distributed sheet in GlaDS depends on the cavity opening rate, which itself is a function of basal velocity, the bedrock obstacle size and the spacing of the bedrock obstacles. Increases in velocity result in a corresponding increase in the cavity opening rate, with larger cavities requiring more water to become pressurised. However, once the bedrock obstacles are flooded, the evolution of the sheet relies on changes in water pressure alone rather than the cavity opening rate. The cavity opening rate plays a central role in the evolution of the sheet in the linear velocity run. In this set of simulations, the highest velocities occur at the outlet and decrease upglacier. This spatial trend is mirrored in the cavity opening rates, and higher water volumes are needed to pressurise the cavities at the grounding line than further inland. As the linear velocity simulation progresses, the 100 ma^-1 velocity contour, which represents the spatially and temporally fixed velocity in the baseline run, migrates further upglacier, and so too do larger cavity opening rates. Since the velocity in the baseline run is, in general, slower than those in the linear velocity run, cavity opening rates are also smaller, which means that the flooded condition is met earlier in the baseline run. In all simulations, negative sheet thickness differences (and, therefore, thicker sheets in the baseline run) occur at the top of the domain where the upper boundary flux enters the domain.

The banded blue regions near 400, 500 and 600 km in Figures 2d–f highlight the impact of spatially variable cavity opening rates on sheet thickness differences. Here, the sheet thickness in the linear velocity run is greater than that of the baseline run by >1 cm. These regions coincide with the initiation point of channels in the baseline run, highlighting that the faster pressurisation of the baseline run's smaller cavities leads to a sufficiently strong hydraulic gradient to allow for channel formation. This stronger hydraulic gradient drains more water from the sheet into the channels, resulting in a thinner sheet in the baseline run compared to the linear velocity run. Therefore, by increasing the velocity, and hence, the cavity opening rate, channel locations are altered, which in turn affects distributed sheet growth and shrinkage.

In Figures 2g–i, we show the difference in per cent overburden between the baseline and linear velocity run. While the overburden differences are less than 1%, there are notable spatial trends present in all simulations. Within the channelised region, the positive differences occur where channels in the linear velocity run have smaller discharges than in the baseline run, and vice versa for negative differences. Behind the channelised region, one might expect the difference in overburden pressure to be negative where the sheet is thicker in the baseline run. However, this is not the case, and the differences are generally positive. This indicates that the linear velocity run is closer to overburden relative to the baseline run. Positive per cent overburden differences occur because the larger channels in the baseline run exert the dominant influence on effective pressure. The impact of channels on effective pressure and drainage efficiency is particularly evident in Figure 2i, where patterns in the overburden differences closely match the locations where channels are smaller in the linear velocity run.

3.2. Slope change

Channels in the slope change runs carry larger discharges and are substantially longer than those in the baseline runs, as shown in Figures 4a–c. The blue regions in the figure panels represent greater discharges in the slope change run than in the baseline run. These discharge differences are up to an order of magnitude larger than those associated with the linear velocity (Fig. 2) runs. Within the suite of slope change runs, there is considerable variability in the discharge differences. For example, under standard water input volumes, the mean channel discharge difference between the baseline and the slope change run over the entire domain at the final timestep is 0.0175 m³s^-1. At the grounding line, there are also appreciable discharge differences between simulations. Here, at the last timestep, the unaveraged discharge differences range from 0.980 m³s^-1 under low water input to +15.17 m³s^-1 under high water input. Another notable result is that the channels in the slope change run extend the furthest from the grounding line of all the forcing experiments tested. The maximum lateral distance these channels reach is 546 km for the low water input run and 1001 km for the high water input run. By increasing the water input from low to high volumes, the number of grounding line nodes with channel discharges of 1 m³s^-1 and above also increases, going from 9 nodes for the low water input run to 14 for the high water input run.

Figure 3. Comparison of velocity and ice thickness fields relative to the start of their respective runs. Panels (a–c) show the velocity difference between the linear velocity run at different time snapshots relative to the start of the transient run. Panels (d) and (e) respectively show the basal velocity and ice thickness evolution along the centreline of the domain (i.e. at 300 km Northing). Panels (f–h) show the ice thickness difference for the slope change run at different timesteps relative to the start of the transient run.

Figure 4. Channel discharge, water sheet thickness and per cent overburden differences between the slope change and baseline runs at the final timestep under different volumes of water input into the distributed system. The panels show (a)–(c) channel discharge difference; (d)–(f) water sheet thickness difference; (g)–(i) per cent overburden difference. Positive differences in blue indicate where the value is larger in the slope change run, while negative differences are shown in red and indicate where the value is larger in the baseline run. The central channel discussed in the main text is labelled in panel (c).

Steepening the ice surface slope steepens the hydraulic potential gradient and affects where the dominant channels form. Figure 4c illustrates a clear example of this: at the centre of the grounding line in the slope change run, a channel (hereafter the central channel; labelled in Fig. 4c) develops. The central channel carries discharges that are more than 10 m³s^-1 larger than those in the baseline run. Within the slope change run itself, the hydraulic gradient is steepest in the direction of the central channel and therefore redirects water towards the central channel (Fig S11) from the surrounding higher-pressure channels.

The per cent overburden difference plot (Fig. 4i) reflects the siphoning effect of water from the higher pressure channels to the central channel. Where channels carry smaller discharges in the slope change run, there are negative per cent overburden differences, and vice versa for positive per cent overburden differences (see also Fig. S10). In contrast to the slope change run's central channel, the largest positive discharge difference occurs near ~238 km Northing, where a channel in the baseline run carries discharges that are 13, and up to 30 m³s^-1 larger than the slope change run.

The slope change runs demonstrate the largest differences in sheet thickness of all the simulations (Figs 4d–f). Sheet thickness differences are generally on the order of ±0.20 m, alternating between positive and negative values close to the grounding line, co-located with major channel locations. For instance, the red regions in panels (d–f) indicate negative sheet thickness differences, occurring where the channels are larger in the slope change run. As in the linear velocity minus baseline comparison, band features occur. However, the bands have the opposite sign in the slope change minus baseline runs. Here, negative sheet thickness differences of up to −0.40 m occur and extend further back across the domain as water input increases. The steeper ice surface slopes near the grounding line (Fig. 3e) in the slope change run cause steeper hydraulic gradients. As a result, the rate of water drainage from the sheet in the slope change run is faster compared to the baseline run, which allows channels to form more easily. Thus, the red bands in Figures 4d–f indicate the channel initiation points in the slope change simulations.

Approaching the top of the domain, ice surface slopes are shallower than downglacier. Since ice surface slope exerts the dominant control on hydraulic gradients, a shallower ice surface means that the hydraulic potential gradient is weaker. As a result, water depths and pressures build up. However, since drainage in the slope change run is more efficient, water is removed more quickly than in the baseline run, and this affects sheet thickness and per cent overburden over the entire domain. For instance, behind the band feature, sheet thickness differences are small (~10^-2 m), but negative, indicating that the sheet is thinner in the slope change run. At the top of the domain, positive overburden differences show the accumulation of water pressure in the slope change run. These positive overburden differences represent the interplay between drainage efficiency, the distributed water input and the upper boundary flux.

As with the other dynamics runs, the greatest per cent overburden differences correspond with the locations of greatest channel discharge differences (Figs 4g–i). Positive overburden differences occur where channels in the slope change run carry less discharge than in the baseline run, and vice versa for negative overburden differences. This spatial pattern intensifies as the water input volume is increased. The differences are the largest of all individual dynamics runs, generally falling within ±1%.

3.3. Combination run

We next describe the results from our combination experiments, in which velocity and surface slope both change over time and determine basal meltwater production rates. Figure 5 compares the results of our combination experiments and the baseline run with standard volumes of water input. This difference plot shows similar spatial patterns to those in the individual dynamics runs, albeit the differences are greater in magnitude. This is unsurprising since the velocity-dependent melt produces melt rates ranging from 0.012 to 0.771  ma^-1, up to two orders of magnitude larger than the constant melt input of 0.005 ma^-1 in the baseline run. Melt rates of 0.012–0.721 ma^-1 occur in the run where both velocity and slope increases determine melt rates.

Figure 5. Channel discharge, water sheet thickness, and per cent overburden differences between the combination runs and the baseline run at the final timestep. Panels (a)–(c) show the channel discharge difference, sheet thickness difference, and per cent overburden difference for the simulation where velocity increases alone determine melt rates. Panels (d)–(f) show the results from the experiment where slope and velocity changes dictate melt rates. Positive differences in blue indicate where the value is larger in the combination experiment run, while negative differences are shown in red and indicate where the value is larger in the baseline run.

Model outputs from the run in which melt rates are based solely on velocity increases closely resemble the simulation in which slope and velocity changes determine melt rates. When differenced (i.e. the velocity and slope change combination run minus the run where melt depends solely on velocity), channel discharges are generally larger in the run that includes slope changes, up to 73 m³s^-1 larger than the velocity-only run. In contrast, the distributed sheet is generally thicker in the velocity-only run, up to 0.5 m, compared to the full combination experiment. Here, we focus primarily on the results from the simulation incorporating both velocity and slope changes since this simulation deviated the most from the baseline run.

When ice dynamic and geometric changes dictate the basal melt rate, channels carry up to 365 m³s^-1 more discharge than when the melt rate is spatially and temporally invariant, as in the baseline run. The mean channel discharge difference across the domain and the entire simulation time is 1.32 m³s^-1. These channels extend slightly further across the domain than in the slope change run, reaching 1058 km upglacier. At the grounding line, the difference in total channelised discharge is 1041 m³s^-1. Overall, the sheet is thicker in the combination run than in the baseline run, with a mean difference of 0.075 m at the final timestep. However, like in the individual dynamics runs, well-defined bands in the difference plot show the channel initiation points. Per cent overburden differences are also predominantly negative over the domain where channels in the combination run carry larger discharges, similar to the slope change run results. Together, these results point to the dominant role of surface slope in steering water flow (Shreve, 1972) and that as Antarctic outlet glaciers respond to changes in climate, surface slope changes will be more influential than velocity accelerations in controlling drainage evolution, with the caveat that these slope changes could be caused by velocity increases via dynamic thinning.

3.4. Sensitivity tests

We next compare the results of our sensitivity tests for the individual dynamics runs (see also Figs S1–S6). We focus on the general differences between the parameter tests and their impacts on the drainage systems of the dynamics runs. Sensitivity tests for the combination run were not computationally feasible.

The horizontal extent of channels varied depending on the sensitivity test. For all dynamics runs, the maximum horizontal distance of channels under the high k _c value is more than twice that of the low k _c value and is the largest factor increase in channel length of all parameter tests. This wide range in channel length occurs since lower channel conductivities represent rougher channel walls that suppress channel growth, whereas high conductivities promote it. Raising the water input rate resulted in maximum channel lengths ~85% longer than those under low water input rates. Increasing distributed system conductivity decreased channel lengths, but the channels under high k were longer than those under low k _c. In contrast, changing the bedrock obstacle size had negligible impacts on the maximum channel length.

Lowering k _c, or increasing the value of k or h _r, tends to inhibit channel growth. Slope and velocity changes either amplify or modulate this effect. For instance, differences in total channelised discharge across the grounding line lessen as k _c decreases. For the slope change minus baseline runs, the steeper ice surface gradients in the slope change run that promote channel growth partially offset the effects of lowering k _c. In the linear velocity run, lowering k _c modulates the impact of the run's faster cavity opening rates that prolong the time for the cavities to pressurise, which hinders channel growth. Raising the distributed system conductivity, k, inhibits channel growth as water flux moves more easily through the distributed system. In the slope change runs, this means that the steeper hydraulic gradients play a less dominant role. In the linear velocity runs, faster velocities cause the system to be further from overburden and raising k reinforces this effect. When the bedrock obstacle size (h _r) increases, the time required for the cavities to pressurise increases – for the slope change suite of runs, raising the obstacle size results in more minor differences in total channelised grounding line discharge. In contrast, for the linear velocity runs, the differences become larger.

3.5. Lakes

To quantify the effect of lakes on subglacial drainage, we compare the drainage system during the peak lake input to the timestep before lake drainage starts (hereafter referred to as pre-event). We adopt the convention that positive differences indicate larger quantities during peak lake input, while negative differences indicate that the pre-event value was larger.

We first focus on the original drainage pattern from Malczyk and others (2020) and the evolution of channel discharge and per cent overburden. We then discuss how using the same volumes of water input from Malczyk and others (2020) but changing the drainage interval length (i.e. changing the flux into the system) and refilling interval length affects the drainage system.

3.5.1. Original drainage pattern

Figures 6a,c show the channel discharge and differences between the peak of the first lake drainage event and the pre-event drainage system. For all the drainage patterns considered, the first lake drainage event had relatively minor effects on hydrology, and all simulations closely resembled each other. Thus, we only discuss the results of the first drainage event for the original drainage cycle simulation as it represents the trends observed in all runs. Over several channels, discharges at the peak of the first event increase by up to 1 m³s^-1 compared to the pre-event (Fig. 6).

Figure 6. Peak minus pre-event channel discharge (a, b) and per cent overburden differences (c, d) showing the first (a, c) and second (b, d) lake drainage events with the original drainage and refilling cycle. Positive differences (in blue) indicate that the channel discharge/per cent overburden is larger during the peak of lake drainage, while negative differences (in red) show where discharge/per cent overburden is larger before lake drainage. The yellow star indicates the lake site.

Prior to lake drainage, the average water pressure in channels with discharges at or above 1 m³s^-1 is 97.86%. The pressure in these channels increases by 0.01% during peak drainage. The greatest increase in pressure occurs at the lake site during the peak of this event. Here, water pressure reaches 105.78% of overburden, which represents a 6.49% increase relative to pre-event water pressure. This is demonstrated by the saturated blue area in Figure 6c. The region above overburden extends 36.49 km downglacier of the lake site. Mean catchment-wide water pressure is above 98%.

The second drainage event (Figs 6b,d) has a more substantial impact on the drainage. The flux through the channels increases by up to 13  m³s^-1 at the peak of drainage, relative to pre-event channel discharge. The largest increases occur immediately downglacier of the lake, with more moderate increases moving further downglacier towards the grounding line. As with the first event, there is a slight increase in pressure in the channels relative to the pre-event. The peak in water pressure occurs at the lake during drainage but is slightly lower than that associated with the first event, with water pressure reaching 104.11% of overburden. Water pressures at or above ice overburden pressure extend 100 km downglacier of the lake site (Fig. 6d), more than double the extent of the high-pressure region in the first event.

3.5.2. Short drainage cycle

Figure 7 shows the results from the short-drainage suite of experiments. In all simulations, the mean water pressure in channels with at least 1 m³s^-1 of discharge is above 97%, and the catchment-wide mean water pressure is above 98%. Shortening the drainage and refilling time leads to maximum increase in channel discharge of 8.49 m³s^-1 relative to pre-event channels for this simulation. This maximum difference is comparable to the maximum difference for the first event in our simulations that use the original drainage cycle from Malczyk and others (2020). Water pressure at the lake site during peak drainage is 104.71%. The above-overburden region extends 43.13 km downglacier of the lake site. With the original refilling time, channel discharge differences reach 13  m³s^-1, similar to those associated with the original drainage cycle. During peak drainage, water pressure is 105.11% of overburden at the lake site. The high pressure associated with lake drainage reaches 85 km downglacier of the lake, slightly shorter than the original drainage cycle. Doubling the lake refilling interval leads to the largest changes in peak minus pre-event differences in channel discharge and water pressure. Channel discharges increase by 25 m³s^-1 at the lake site, but this increase is localised. Smaller increases on the order of ~5  m³s^-1 are found moving towards the grounding line. Lake water pressure is 105.50% of overburden at peak drainage, and the high pressure extends 147.67 km downglacier of the lake, the furthest of any simulation. By the end of the second event, per cent overburden differences across the domain relative to the peak of lake drainage are predominantly negative, indicating that the spike in water pressure from the peak lake drainage dissipates quickly, within 150 days.

Figure 7. Peak minus pre-event channel discharge (a–c) and per cent overburden differences (d–f) showing the second lake drainage event with a shorter-than-original drainage cycle and varying refilling cycle lengths. Panels (a)–(c) show the channel discharge differences for experiments with (a) short drainage and a short refilling interval length, (b) short drainage and the original refilling interval length, and (c) short drainage and a longer refilling interval. Panels (d)–(f) show the per cent overburden differences for the short drainage runs with the following refilling interval lengths: (d) a short refilling time, (e) the original refilling time and (f) a longer refilling time. Positive differences (in blue) indicate that the channel discharge/per cent overburden is larger during the peak of lake drainage, while negative differences (in red) show where discharge/per cent overburden is larger prior to lake drainage. The yellow star indicates the lake site.

3.5.3. Long drainage cycle

When the refilling cycle length deviates from the original duration, there are considerable changes in the flux through the channelised system for the second event (Fig. 8). For the short refilling run, channel discharges at the peak of lake drainage increase by up to 8  m³s^-1 immediately downglacier of the lake site, then decrease slightly to 5  m³s^-1 moving towards the grounding line (Fig. 8a). When the refilling time is the original duration, changes in channel discharge are minor and do not exceed 3 m³s^-1 (Fig. 8b). In contrast, channel discharges for the long refilling run during peak lake drainage are up to 25  m³s^-1 larger than before lake drainage. Notably, this increase in channel discharge nearly five times larger than the increase associated with the peak of the first event.

Figure 8. Peak minus pre-event channel discharge (a–c) and per cent overburden differences (d–f) showing the second lake drainage event with a longer-than-original drainage cycle and varying refilling cycle lengths. Panels (a)–(c) show the channel discharge differences for experiments with (a) long drainage and a short refilling interval length, (b) long drainage and the original refilling interval length, and (c) long drainage and a longer refilling interval. Panels (d)–(f) show the per cent overburden differences for the long drainage runs with the following refilling interval lengths: (d) a short refilling time, (e) the original refilling time and (f) a longer refilling time. Positive differences (in blue) indicate that the channel discharge/per cent overburden is larger during the peak of lake drainage, while negative differences (in red) show where discharge/per cent overburden is larger prior to lake drainage. The yellow star indicates the lake site.

There are also more considerable differences in the water pressure for the second event in the pre- versus peak lake drainage compared to the first event. Water pressure at the lake site is 102.9% of overburden for the long drainage, short refilling simulation and 103.3% of overburden for the long drainage, long refilling simulation. At the lake site, lake drainage causes a more than 4% increase in water pressure relative to pre-event overburden pressure for all refilling interval lengths. This localised spike in water pressure at the lake site has catchment-wide impacts on the water pressure in the drainage system. However, the effects of this localised ≥4% increase in water pressure vary between simulations and depend on the configuration of the pre-existing drainage network. That is, the extent to which the high-pressure region associated with lake drainage spreads into downglacier areas depends on where the drainage system was below overburden before lake drainage.

High pressures at the lake site permeate into downstream channels. This effect is evident in the short refill and long refill simulations. In the short refill run, the drainage system has not efficiently evacuated the lake input from the first event, so water pressure builds up and only drops below overburden 43 km downglacier of the lake site. In the long refilling scenario, pressures stay at or above overburden pressure for 147.59 km downstream of the lake site. After lake drainage ends, water pressures remain elevated across most of the domain for the long drainage, original refilling run and the long drainage, long refilling run. In contrast, water pressures at the end of the second drainage event for the long drainage, short refilling simulation are lower than those before the lake drained.

3.5.4. Drainage cycles and grounding line discharge

In Figure 9, we show how changing draining and refilling interval impacts the flux of channelised discharge crossing the grounding line. After the first event, grounding line discharge rises to a peak that is similar for all simulations. Grounding line discharge then decreases until the onset of the second event. All runs show a rapid increase in grounding line discharge 100–150 days after the peak of the second lake drainage input, followed by a steady decline as the simulation progresses. For the short-drainage runs this timing also corresponds with the approximate timing for pressure to decrease over the domain from the peak of lake drainage to the end of the event.

Figure 9. Comparison of drain-refilling cycles with total channelised discharge crossing the grounding line. The y-axis on the left-hand side (in red) shows the total grounding line discharge (in m³s^-1) and the y-axis on the right-hand side (in black) shows the lake flux (in m³s^-1) into the drainage system. Panel (a) shows the original Thw₁₄₂ lake drainage pattern from Malczyk and others (2020), plotted against total channelised grounding line discharge from our GlaDS model outputs. Panels (b)–(g) show the results of our experiments in which drainage and refilling interval lengths are altered. Panels (b)–(d) show the channelised discharge and lake flux for shorter drainage cycles and (b) a shorter refilling interval, (c) the original refilling interval length and (d) a longer refilling interval. Panels (e)–(f) show the channelised discharge and lake flux for longer drainage cycles and (b) a shorter refilling interval, (c) the original refilling interval length and (d) a longer refilling interval. Note the different y-axis limits for the lake discharges between simulations.

4. Discussion and significance

Antarctic glaciers are speeding up, losing mass and retreating (Meredith and others, 2019). However, studies of recent and future changes in Antarctic ice loss and flow speed have seldom focused on how these changes impact the various components of the ice sheet's diverse basal hydrological network, such as efficient channels, inefficient linked cavities and subglacial lakes. Here, we modelled 2D co-evolving distributed and channelised subglacial drainage and included the effects of increases in ice surface slopes, basal velocities, their combined impact on basal meltwater production and subglacial lake activity to show how changes in any of these could affect subglacial hydrology. Below, we discuss the implications of our results, primarily focusing on the consequences of subglacial hydrological changes for Antarctic ice dynamics. We focus on the implications of our results for Antarctic catchments that have experienced and will likely continue to experience geometric, dynamic and subglacial drainage changes. However, since the current picture of Antarctic subglacial drainage systems is incomplete, as is the understanding of which and how many catchments will undergo any of the abovementioned changes in the future, we cannot discount the potential impact that subglacial hydrologic variability may have on the continental-scale.

4.1. Velocity increases

Our results indicated that channel formation is inhibited when basal velocities continually increase. Faster velocities increase the cavity opening rate, prolonging the time required for the cavities to pressurise. Channel growth is hindered since there is insufficient hydraulic potential gradient to force the water out of the sheet and into the channels. Therefore, the channels grow larger in the baseline run with temporally constant velocity since there is more hydraulic potential drive than that associated with the linear velocity run. These findings are germane to vulnerable regions of the Antarctic, where glaciers have accelerated in recent decades. For instance, the velocity of Crane Glacier was 548 ma^-1 before the 2002 breakup of the Larsen B Ice Shelf on the Antarctic Peninsula but then increased to 2484 ma^-1 following collapse (Wuite and others, 2015). Velocities of the Crane Glacier and other tributary glaciers feeding the former ice shelf have remained above pre-collapse values ever since (Berthier and others, 2012; Wuite and others, 2015). While numerical limitations of the GlaDS model prevent us from incorporating such large velocity accelerations, our findings nonetheless indicate that persistent velocity increases result in a continually adjusting subglacial drainage system. The sustained, elevated velocities years after ice shelf collapse suggest that there could also be long-term effects on subglacial hydrology, which would in turn, feed back onto ice flow speeds.

4.2. Ice surface slope changes

We also considered changes to ice surface slopes, which are particularly relevant to the Amundsen Sea Sector of West Antarctica, where surface slopes have steepened in recent years (Scott and others, 2009; Wingham and others, 2009; Flament and Rémy, 2012). As the ice surface slope increases, our model results indicate that channels grow further upstream and are larger than those found in the baseline run where surface slopes are time-invariant. As these channels extend further across the domain, there are considerable impacts on overburden pressure and water sheet thickness. Although we do not include two-way coupling between ice dynamics and subglacial hydrology, we can infer the impact of slope changes on ice dynamics: where our results indicate pressure drops, one would expect the ice to slow down. Conversely, where pressure increases, the ice would speed up. These changes in velocity would, in turn, alter the ice surface slope, therefore feeding back on the subglacial drainage system. In our model simulations, the largest pressure drop occurs near the grounding line, where the increase in surface slope is greatest, which in turn drives the development of larger channels. Our implementation of thinning is qualitatively in agreement with observations of Pine Island Glacier, which show that the largest magnitude of thinning is concentrated at the grounding line and weakens moving inland (Scott and others, 2009; Wingham and others, 2009; Flament and Rémy, 2012). These studies indicated that the slope steepening of Pine Island Glacier in recent decades caused glacier acceleration through dynamic thinning while also redirecting subglacial drainage (Scott and others, 2009; Wingham and others, 2009; Flament and Rémy, 2012; Bougamont and others, 2019). We have also shown that changes in surface slope impact where dominant channels form. Numerical limitations of the model preclude the use of high thinning rates of 5 ma^-1 recently observed near the grounding line of Pine Island Glacier (Konrad and others, 2017), but our adopted 1 ma^-1 thinning rate is representative of thinning on other vulnerable glaciers, such as Totten Glacier in East Antarctica (Flament and Rémy, 2012; Li and others, 2015) and glaciers along the Southern Antarctic Peninsula (Wouters and others, 2015). Furthermore, our results corroborate the findings of Wright and others (2008), who, using a simplified water routing scheme, showed that ice surface elevation changes of 15 m and less are sufficient to redirect water flow paths in East Antarctica. This is well within the range of ice thickness changes implemented in our more complex modelling approach.

4.3. Channelisation and ice-sheet stability

We have shown that even when we combine both dynamic forcings to allow frictional melt rates to vary spatially and temporally, changes in surface slope dominate changes in the drainage system. As a result of steeper slopes and the corresponding increase in the efficiency of the subglacial drainage system, the flow regime of Antarctic glaciers could critically depend on the presence and evolution of channelised drainage systems. However, whether channels will stabilise ice flow or reinforce flow instabilities is unknown. On one hand, the more efficient channels increase the coupling between the ice and the bed, drawing water into them from the surrounding region. Basal friction would increase, stabilising flow and potentially slowing the grounding line retreat, attenuating ice loss. On the other hand, more efficient channels could weaken ice shelves by increasing sub-ice shelf melt rates. Outflows from subglacial channels are associated with high and localised sub-ice shelf melt rates in Antarctica (Le Brocq and others, 2013; Alley and others, 2016; Adusumilli and others, 2020; Wei and others, 2020; Dow and others, 2022). Due to slope steepening, new discharge locations from subglacial channels could introduce additional weaknesses in the ice shelf, and greater volumes of subglacial discharge crossing the grounding line could increase the already-high sub-ice shelf melt rate. As a result, the ice shelf could be destabilised, decreasing the amount of buttressing provided by the ice shelf to the grounded ice and potentially triggering irreversible retreat of the grounding line and rapid rates of ice loss.

Since the rate and quantity of ice discharged into the ocean, and therefore sea-level rise contributions, are intricately tied to the characteristics of the subglacial drainage network, it is of critical importance to determine which of the competing effects of channelisation will dominate ice-sheet evolution. While we do not model two-way coupling of an ice sheet and the subglacial drainage system, our findings nevertheless highlight the importance of resolving time-evolving subglacial hydrology in ice-sheet models, which often do not include hydrology at all. We, therefore, motivate the research community to build on the work presented here by studying the impacts of geometric changes on subglacial hydrology, and the corresponding impacts on ice shelf and ice-sheet stability, by coupling ice-sheet and subglacial hydrology models with ice shelf cavity circulation models. Capturing these interactions will increase our understanding of Antarctica's ice dynamic response to changes at the ice–bed and ice–ocean interfaces and will provide more representative projections of Antarctic ice dynamics and contributions to sea-level rise.

4.4. Lakes

4.4.1. Effect of altering refilling intervals

Recent work focusing on lakes beneath Thwaites Glacier hypothesised that a lake drainage event can pre-condition the drainage system for subsequent lake drainages through the formation of efficient channels (Malczyk and others, 2020). Our results suggest that the refilling time between lake drainage events exerts a dominant control on the response of the drainage system to multiple lake drainages as it determines the relative importance of creep closure. When the lake is refilling, there is no water supplied to the drainage system from the lake, and therefore the only water maintaining the channels is that produced in situ. If this volume of water is insufficient to maintain the channels at the same size as during the first event, the channels contract through creep closure. When the refilling time is shorter, there is still sufficient flux through the channels carved by the first lake drainage event to maintain them at a similar size. Consequently, when the lake drains a second time, its flux is more readily accommodated by the pre-existing channels, and the increases in pressure and channel discharge are more moderate relative to those associated with longer refilling times. In contrast, when the refilling time is longer, most of the water pumped into the system during the first event has been evacuated, so there is not enough water in the channels to keep them as large. As a result, the channels contract substantially. When the lake drains again, the flux travels through the channels established during the first event. However, since flux through these pre-existing channels was comparatively low before the lake drained, lake drainage causes a spike in channel flux. The drainage system is overwhelmed by the lake input, and water pressures exceeding ice overburden pressure nearly 150 km downglacier of the lake site.

4.4.2. Lake activity and ice dynamics

Our refilling time sensitivity tests have implications for ice dynamics. With infrequent lake drainages due to long refilling intervals, ice velocity increases will be large since relic channels from previous lake drainage events have mostly shut down. Thus, the drainage system cannot accommodate the extra water. The length of time that the lake takes to drain is of secondary importance, but if the drainage time is shortened at the same time as the refilling time is lengthened, the pressure increase (and hence, velocity response) will be amplified. However, it is possible that in a warming climate, lake activity becomes more frequent, as suggested by Livingstone and others (2022). The authors indicated that there could be fewer and smaller lakes as ice surface slopes steepen. These lakes would evacuate large volumes of water in short bursts through channels. Our findings support Livingstone and others’ (2022) hypothesis that short-duration drainage events can quickly shed lake flux through channelised discharge. For our idealised glacier, whose dimensions are similar to Recovery Glacier (Dow and others, 2016), we found total channelised grounding line discharge peaks 100–150 days after the lake drainage peak. Our results also illustrate that the drainage system recovers more quickly when the drainage system is subject to short-duration lake drainage events. There are large pressure drops at the end of the drainage event compared to the beginning of lake drainage. When the drainage is prolonged, water pressure decreases over the domain at the end of the event, but the spatial extent is limited compared to the short drainage event.

5. Conclusion

Using an idealised Antarctic glacier system, we modelled the effects of steepening ice surface slopes, accelerating velocities and subglacial lake drainages on the evolution of the subglacial drainage system. Our simplified approach allowed us to determine which processes exert the dominant control on basal hydrological networks. In agreement with previous work (Wright and others, 2008; Bougamont and others, 2019), our findings show that surface slope changes can substantially alter the discharges and locations of primary drainage pathways. Velocity increases had more attenuated effects on drainage, with the associated changes in subglacial drainage an order of magnitude smaller than those caused by surface slope changes.

Our results do not represent ice-ocean-subglacial hydrology feedback processes since we do not couple the subglacial hydrology model to an ice or ocean model. While this can be viewed as a limitation, our results suggest feedbacks that can be explored in future work using coupled models. For instance, we have shown that due to ice surface slope changes, subglacial channels can be redirected to new locations. We suggested that such changes could impact the buoyant circulation and melt rates beneath ice shelves, potentially weakening the ice shelf and driving further dynamic and geometric changes. Projections indicate committed dynamic thinning in vulnerable regions such as the Amundsen Sea Embayment (e.g. Golledge and others, 2019, 2021), but the treatment of subglacial hydrology in many ice-sheet models has been simpler than the 2D co-evolving inefficient and efficient representation adopted in this work. Therefore, we encourage and provide motivation to the glaciological community to use subglacial hydrology models incorporating both channelised and distributed drainage in coupled simulations of ice sheet–ice shelf–ocean interactions to test whether ice shelves are more or less stable due to dynamically driven changes in subglacial hydrology. Results from this line of work would provide valuable information about the volume of subglacial discharge above which ice shelves could collapse due to enhanced ice shelf basal melting. Determining such thresholds along with the temporally and spatially changing characteristics of subglacial hydrological systems would better constrain the likelihood of unstable retreat mechanisms and projections of sea-level rise.