Subglacial hydrology regulates oscillations in marine ice streams

Marine ice stream dynamics are sensitive to conditions at the grounding line and basal shear stress. Variations in subglacial hydrology have been implicated in ice stream speed-up and shutdown. To investigate the interplay between marine ice stream flow and subglacial hydrology, we couple models of a marine ice stream and a subglacial drainage system. The coupled system evolves dynamically due to a positive feedback between ice flow, heat dissipation at the ice stream bed, and basal lubrication. Our results show that, depending on the hydraulic conductivity of the bed, distinct dynamic regimes can be identified. These regimes include steady streaming, hydraulically controlled oscillations, and thermally controlled oscillations. Periodic fast flow can be initiated by activation waves travelling upstream or downstream or quasi-simultaneously everywhere. Different dynamical regimes are characterized by large differences in grounded ice volume, even under modest changes of grounding line positions. These results imply a strong dependence of marine ice sheet dynamics on evolving hydrological conditions at the bed and highlight the importance of a better understanding of subglacial hydrology.

1 Introduction

Much of the West Antarctic Ice Sheet rests on a bed that is below sea level, making it a marine ice sheet. Most of its ice discharge occurs in regions of fast flow with velocities in excess of 100 m yr⁻¹ (called ice streams) embedded in regions of slow flow with velocities of only a few metres per year (Rignot et al., 2011). While most of West Antarctica has experienced dramatic mass loss in the last four decades, the Siple Coast region has a net positive mass balance (Rignot et al., 2019, 2022). This is due to the current inactivity of Kamb ice stream, which slowed down approximately 130 years ago (Retzlaff and Bentley, 1993) and the ongoing slow down of Whillans, Mercer, and MacAyeal ice streams (Stephenson and Bindschadler, 1988; Bindschadler and Vornberger, 1998; Joughin et al., 2002; Rignot et al., 2022). On longer timescales, observations suggest that marine ice stream flow speeds vary on timescales of hundreds to thousands of years and can undergo shut-down-and-reactivation cycles (Hulbe and Fahnestock, 2007; Catania et al., 2012). Our aim in this paper is to better understand what role subglacial hydrology plays in these processes.

Some studies have suggested that subglacial hydrology is responsible for the speed up and slow down of Antarctic ice streams (Anandakrishnan and Alley, 1997; Elsworth and Suckale, 2016). Observations support the existence of a dynamic subglacial drainage system beneath West Antarctica's ice streams (Engelhardt and Kamb, 1997; Fricker et al., 2007). However, existing observational work does not address the interplay between subglacial hydrology and marine ice sheet dynamics because of a lack of direct observations.

Conversely, theoretical studies of marine ice sheet dynamics usually focus on static basal boundary conditions (e.g. Weertman, 1974; Schoof, 2007a, b, 2012; Tsai et al., 2015; Joughin et al., 2019; Sergienko and Wingham, 2019; Sergienko and Haseloff, 2023; Gregov et al., 2023) or prescribe their evolution (Brondex et al., 2017; Sergienko and Wingham, 2024). These studies show that the dynamics of the grounding line, which delineates the transition between grounded and floating ice, are sensitive to the basal boundary conditions. However, without an evolving subglacial drainage system, these models do not investigate which dynamic feedbacks exist between marine ice stream dynamics and subglacial drainage.

Investigation of the interaction between ice stream flow and subglacial drainage raises the question of how regions of fast and slow flow can co-exist in close proximity. One answer to this question invokes a positive feedback between heat dissipation and ice flow. Different mechanisms for this feedback have been hypothesized, relating it either to viscous weakening or to basal heat dissipation (e.g. Clarke et al., 1977; Fowler and Johnson, 1996; Payne and Dongelmans, 1997; Bougamont et al., 2011; Robel et al., 2014; Kyrke-Smith et al., 2014, 2015; Brinkerhoff and Johnson, 2015; Schoof and Mantelli, 2021). In this study, we focus on the latter feedback.

If the basal temperature is at the melting point, rapid sliding is possible. The dominant positive feedback mechanism then involves melting at the base of the ice; faster sliding leads to more heat dissipation by basal friction, which in turn produces additional melt water that reduces basal drag and permits even faster sliding until other processes suppress further weakening of the bed. For this positive feedback to operate, melt water production must raise the water pressure at the bed, and hence reduce the effective pressure, defined as the difference between ice overburden and basal water pressure. This behaviour can be caused by two mechanisms. In the first, the effective pressure decreases through an increase in water storage at the bed. If water storage is in subglacial till and no drainage is possible due to the material's low permeability, this is called the undrained bed model (e.g. Tulaczyk et al., 2000b; Sayag and Tziperman, 2008). In the second mechanism, an increase in subglacial water discharge reduces the effective pressure (as is common for distributed hydrologic systems, e.g. Fowler and Johnson, 1996; Kyrke-Smith et al., 2014).

Several studies using the undrained bed model have shown that the positive feedback between fast flow and lubrication can explain oscillatory behaviour of ice streams when modulated by a thermodynamic feedback (MacAyeal, 1993; Fowler and Johnson, 1996; Payne and Dongelmans, 1997; Fowler and Schiavi, 1998; Bougamont et al., 2003; Christoffersen and Tulaczyk, 2003; Robel et al., 2013, 2014; Feldmann and Levermann, 2017). Faster ice flow results in dynamic thinning. When the ice thickness becomes too small, the amount of conductive cooling at the bed exceeds the heat dissipation due to flow, water freezes, and the ice stream shuts down. However, these studies neglect the transport of subglacial water and its effect on water content of the till.

Fewer theoretical studies investigate the role of the second process (i.e. of subglacial hydrology), and these studies do not account for the subglacial conditions at the grounding line (Fowler and Johnson, 1996; Kyrke-Smith et al., 2014). Nevertheless, large-scale models capable of including these dynamics are beginning to emerge (e.g. Kazmierczak et al., 2022, 2024), motivating a better understanding of these processes. An extreme example of a subglacial drainage model assumes that subglacial hydrology is so efficient that hydraulic gradients cannot develop. Steady grounding lines are possible in this case (Tsai et al., 2015).

The goal of this study is to investigate the simplest possible feedback between fast flow, heat dissipation, and basal lubrication, and its role in marine ice sheet dynamics. We use a model that includes the positive feedback between fast flow and basal lubrication both through storage of water in the subglacial sediment and evacuation of water through an active subglacial drainage system. Subglacial water flow is modelled with a Darcy-style flux law. The hydraulic conductivity is the quotient of a conductivity factor and the effective pressure, i.e. it is assumed to increase as the effective pressure decreases, as is common for distributed systems. The overall balance of storage and subglacial discharge is determined by the hydraulic conductivity factor. At very low hydraulic conductivity ($K_{\mathrm{d}}\ll {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹), the model reproduces the undrained bed model, and the water content of the bed is determined by the local energy balance. At very high hydraulic conductivities ($K_{\mathrm{d}}\gg {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹), hydraulic gradients cannot develop at the bed; in this limit the effective pressure at the bed is set by the ice and bed geometry alone, independent of the local melt rate. This is the limiting case considered in Tsai et al. (2015).

The main novelty of our study is the investigation of marine ice sheet dynamics over the entire spectrum of different hydraulic conductivities of the subglacial system. This puts the results of existing studies with vanishing (Robel et al., 2014) and infinite (Tsai et al., 2015) hydraulic conductivities into context. We find that at intermediate hydraulic conductivities of the bed, the interplay between local water storage and subglacial drainage can lead to oscillations with distinct character and frequency. These hydraulically induced dynamics differ from those previously reported (MacAyeal, 1993; Robel et al., 2014); they are controlled by the efficiency of the subglacial drainage system, rather than the local melt rate at the bed.

Figure 1 Model geometry and details of freezing. A marine ice sheet of thickness H and width W flows on a downwards sloping bed, becoming afloat at the grounding line. The underlying sediment is of thickness $(e+\mathrm{1})h_{\mathrm{s}}+(e_{\mathrm{c}}+\mathrm{1})(h_{\mathrm{0}}-h_{\mathrm{s}})$ with e the void ratio of the unfrozen sediment, e_c the void ratio at which the sediment freezes, h₀ the thickness of the sediment layer if all void space was removed from it (i.e. at e=0), and (e+1)h_s the thickness of the expanded, hydraulically active, unfrozen portion of the sediment. At effective pressures N below the critical effective pressure N_c, freezing is focused at the ice–bed contact (a). At N≥N_c, freezing can intrude into the void space between sediment grains, and a frozen fringe forms (b).

2 The model

We formulate the simplest possible model with an evolving hydraulic system and grounding line dynamics. Naturally, such a model must parameterize or neglect some other important aspects of ice sheet dynamics and we discuss the limitations of our model in Sect. 4.4.

The computational domain is of length L_x=1000 km with a downward sloping bed at elevation $b=z_{\mathrm{0}}-z_{\mathrm{1}}x$. All parameters are listed in Table 1. We solve for velocity u, ice thickness H, grounding line position x_g, water content of the bed h_w (given as water column thickness), and thickness of unexpanded, hydraulically active till layer h_s, which can evolve due to freezing and melting in the sediment, see Fig. 1.

2.1 Ice mechanics

We use a standard, depth- and width-integrated flowline model for ice sheet dynamics (Dupont and Alley, 2005; Nick et al., 2009; Robel et al., 2014; Haseloff and Sergienko, 2018). The ice thickness H evolves according to the equation for mass balance,

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial H}{\partial t}}+{\displaystyle \frac{\partial \left(Hu\right)}{\partial x}}=a,\end{array}$$

with u the ice velocity in the downstream x direction and a the surface mass balance term (positive for net accumulation). The velocity u is determined from stress balance (e.g. MacAyeal, 1989; Hindmarsh, 2012; Pegler, 2016),

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }\mathrm{2}{\displaystyle \frac{\partial }{\partial x}}\left[A^{-\mathrm{1}/n}H{\left[{\left|{\displaystyle \frac{\partial u}{\partial x}}\right|}^{\mathrm{2}}+{\mathit{\epsilon }}^{\mathrm{2}}\right]}^{\frac{\mathrm{1}-n}{\mathrm{2}n}}{\displaystyle \frac{\partial u}{\partial x}}\right]-{\mathit{\tau }}_{\mathrm{b}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{C_{\mathrm{w}}{A}^{-\mathrm{1}/n}}{W^{\mathrm{1}/n+\mathrm{1}}}}H|u{|}^{\mathrm{1}/n}{\displaystyle \frac{u}{\left|u\right|}}={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\partial (H+b)}{\partial x}},\end{array}\end{array}$$

with A the rate factor, n a rheological parameter, ε a parameter that regularizes regions with zero strain rates (e.g. Schoof, 2006; Bueler and Brown, 2009), and τ_b the basal shear stress. The third term parameterizes the lateral shear stress with the width of the outlet glacier W; C_w is a constant (see Hindmarsh, 2012, for details).

Table 1 Model parameters and their values. Geometric (W, z₀, z₁), environmental (a, T_s, A, q_geo), and basal parameters (h₀, C) were chosen to be representative of Siple Coast ice streams and for comparability with previous studies (e.g. Robel et al., 2013, 2014; Tsai et al., 2015). Subglacial drainage parameters (e_r, e_c, N_r) are based on Tulaczyk et al. (2000a).

In accordance with recent studies, we assume that the relationship between basal shear stress τ_b and velocity u is of a regularized-Coulomb form (Zoet and Iverson, 2020; Helanow et al., 2021). For large velocities, the basal shear stress becomes independent of the velocity and only depends on the effective pressure $N={\mathit{\rho }}_{\mathrm{i}}gH-p_{\mathrm{w}}$, the difference between ice overburden and water pressure. We use a regularized version of the slip law used in Tsai et al. (2015),

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{\mathrm{b}}={\displaystyle \frac{C|u{|}^m\mathit{\mu }N}{C|u{|}^m+\mathit{\mu }N}}{\displaystyle \frac{u}{\left|u\right|}},\end{array}$$

with $m=\mathrm{1}/n$ (e.g. Weertman, 1974; Pattyn et al., 2012; Brondex et al., 2017). For $\mathit{\mu }N\gg C|u{|}^m$ this reproduces a Weertman-style sliding law ${\mathit{\tau }}_{\mathrm{b}}\sim C|u{|}^m$ (Weertman, 1957), while for $\mathit{\mu }N\ll C|u{|}^m$ this mimics the Coulomb-plastic rheology τ_b∼μN expected to be applicable to West Antarctic subglacial tills (Zoet and Iverson, 2020).

We are interested in interactions of the subglacial drainage system with ice stream flow, which are coupled through the effective pressure N. Modelling the evolution of N requires a theory for the evolution of the subglacial drainage system. We turn to this next.

2.2 Subglacial drainage

Many models describing subglacial drainage exist (Flowers, 2015). All models share the assumption that water is conserved at the bed and flows in the direction of hydraulic gradients (Shreve, 1972; Flowers, 2015). However, only a few models are applicable to soft-bedded glaciers (Walder and Fowler, 1994; Ng, 2000a, b; van der Wel et al., 2013; Bougamont et al., 2014; Bueler and van Pelt, 2015). We write conservation of subglacial water mass as (Clarke, 2005; Flowers, 2015)

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial h_{\mathrm{w}}}{\partial t}}+{\displaystyle \frac{\partial q_{\mathrm{w}}}{\partial x}}=m_{\mathrm{b}},\end{array}$$

where h_w is the effective water depth at the bed, q_w is the water flux, and m_b is the basal melt rate.

We assume that the hydraulic conductivity only depends on the water content of the bed, as is common for distributed systems (Hewitt, 2011). This can also be modified to capture the formation of subglacial conduits (Sommers et al., 2018). However, aiming to keep our model as simple as possible, we write for the water flux

$$\begin{array}{}\text{(5)}& q_{\mathrm{w}}=-{\displaystyle \frac{K_{\mathrm{d},\text{eff}}{h}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{w}}g}}{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial x}}\text{ with }{K}_{\mathrm{d},\text{eff}}=K_{\mathrm{d}}{\displaystyle \frac{N_{\mathrm{c}}}{N}}.\end{array}$$

The effective hydraulic conductivity K_d,eff (in m s⁻¹) depends on the effective pressure (e.g. Haseloff et al., 2019); ρ_w is the density of water and N_c is a constant (for convenience we use the critical effective pressure for freezing, introduced below). The hydraulic potential Φ is given by

$$\begin{array}{}\text{(6)}& \mathrm{\Phi }={\mathit{\rho }}_{\mathrm{w}}gb-N+{\mathit{\rho }}_{\mathrm{i}}gH.\end{array}$$

Its gradient drives the water flux. Equation (5) is similar to a model of distributed water flow in a system of subglacial conduits (“canals”) eroded into soft beds (Walder and Fowler, 1994; Ng, 2000b).

The basal melt rate m_b satisfies an energy conservation equation,

$$\begin{array}{}\text{(7)}& m_{\mathrm{b}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}{L}_{\mathrm{h}}}}\left(q_{\text{geo}}+{\mathit{\tau }}_{\mathrm{b}}u+k{\left.{\displaystyle \frac{\partial T}{\partial z}}\right|}^{+}\right),\end{array}$$

with L_h the latent heat of fusion, q_geo the geothermal heat flux, k the thermal conductivity of ice, and ${\left.\partial T/\partial z\right|}^{+}\approx (T_{\mathrm{s}}-T_{\mathrm{m}})/H$ the temperature gradient into the ice, which we approximate as linear between surface and melting point temperatures. This assumes that the bed is temperate; this condition is satisfied for all results shown here.

Motivated by observations of thick sediments at the base of the West Antarctic Ice Sheet (Alley et al., 1986; Anandakrishnan et al., 1998), we assume that water is stored locally in the till, and h_w is proportional to the void ratio e of the till and the thickness of the drained hydraulically active sediment column h_s,

$$\begin{array}{}\text{(8)}& h_{\mathrm{w}}=e{h}_{\mathrm{s}}.\end{array}$$

The void ratio increases with decreasing effective pressure (Tulaczyk et al., 2000a)

$$\begin{array}{}\text{(9)}& e=e_{\mathrm{r}}-C_e\log \left({\displaystyle \frac{N+N_{\mathrm{0}}}{N_{\mathrm{r}}}}\right)\end{array}$$

where e_r, C_e, N_r, and N₀ are constants (see Table 1). This relationship is based on the empirical relationship $e=e_{\mathrm{r}}-C_e\log (N/N_{\mathrm{r}})$ between the void ratio of subglacial sediment and effective pressure (Tulaczyk et al., 2000a). The constant N₀≪N_r has been introduced to satisfy the boundary condition N=0 at the grounding line (Eq. 13), which now corresponds to a void ratio of e₀=1. Void ratios in excess of e₀ would correspond to an overpressured system (N<0), but with our choice of parameters this does not happen. We expect that additional model physics would be necessary to capture such a case.

Note that h_s is the thickness of the water-saturated sediment column if all pore space is removed from it (e=0). As water-filled pores expand the hydraulically active sediment column, its total thickness is $(\mathrm{1}+e)h_{\mathrm{s}}=h_{\mathrm{s}}+h_{\mathrm{w}}$ (see Fig. 1).

2.3 Freezing and melting in the till

Both e and h_s can evolve dynamically, and we require a condition in addition to Eq. (8) to determine them. This condition comes from the dynamics of freezing and melting in the sediment, which determines the thickness of the hydraulically active sediment layer h_s. We adopt the basal model of Robel et al. (2013, 2014) and Haseloff (2015),

$$\begin{array}{}\text{(10)}& {\displaystyle }e\left(N\right){\displaystyle \frac{\partial h_{\mathrm{s}}}{\partial t}}{\displaystyle }=\left\{\begin{array}{ll}{m}_{\mathrm{b}}& \begin{array}{rl}{\displaystyle }\text{ if }(& {\displaystyle }{m}_{\mathrm{b}}<\mathrm{0}\text{ and }{h}_{\mathrm{s}}>\mathrm{0}\\ {\displaystyle }& {\displaystyle }\text{ and }(N\ge N_{\mathrm{c}}\text{ or }{h}_{\mathrm{s}}\le h_{\mathrm{0}}))\end{array}\\ \frac{e\left(N\right)}{e\left(N_{\mathrm{c}}\right)}{m}_{\mathrm{b}}& \text{ if }(m_{\mathrm{b}}>\mathrm{0}\text{ and }\mathrm{0}<h_{\mathrm{s}}<h_{\mathrm{0}})\\ \mathrm{0}& \text{ otherwise.}\end{array}\right.\end{array}$$

This condition states that to overcome surface tension and extend the ice surface into the pore space between the grains (corresponding to h_s<h₀), the effective pressure must reach a critical value N_c (see, e.g., Rempel, 2008, 2009). Therefore, as long as N<N_c and h_s=h₀, the pore space expands or consolidates in response to water addition or removal. If freezing occurs while N<N_c, then water is removed from the sediment and frozen onto the base of the ice stream (Fig. 1a). If N≥N_c, freezing results in a frozen fringe, i.e. ice intrudes into the void spaces of the till (Fig. 1b).

In adopting Eq. (10), we make two assumptions that simplify the dynamics of basal freezing. First, once a frozen fringe has formed, any further freezing occurs inside the sediment, even if the effective pressure drops below N_c, as described by the first case in Eq. (10). This assumption excludes the formation of distinct ice lenses. Second, the void ratio of the frozen sediment is always e_c=e(N_c). It is conceivable that freezing could occur at smaller void ratios if subglacial drainage has removed water before the onset of freezing so that N>N_c. However, melting can readily occur at effective pressures smaller than N_c, accounting for the second case in Eq. (10).

Note that the thickness of the unexpanded, hydraulically active sediment layer h_s is bounded between 0 (completely frozen) and h₀ (completely unfrozen). Our parameter choices prevent the basal sediment from freezing entirely ($h_{\mathrm{w}}=h_{\mathrm{s}}=\mathrm{0}$). In this case the basal temperature could be below the melting point. This scenario might become important in some areas of Antarctica (Pattyn, 2010; Mantelli et al., 2019). We simultaneously solve for h_w and h_s, with the void ratio e determined from Eq. (8) and the effective pressure determined from Eq. (9).

2.4 Boundary and initial conditions

We assume an ice divide at x=0 with zero surface slope, velocity and water flux

$$\begin{array}{}\text{(11)}& {\displaystyle }u={\displaystyle \frac{\partial (H+b)}{\partial x}}={\displaystyle \frac{\partial \mathrm{\Phi }}{\partial x}}{\displaystyle }=\mathrm{0}\text{ at }x=\mathrm{0}.\end{array}$$

The position of the grounding line is given through the floatation condition and stress balance condition

$$\begin{array}{}\text{(12)}& \begin{array}{r}{\displaystyle }\left.\begin{array}{rl}{\displaystyle }H& {\displaystyle }=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}}b\\ {\displaystyle }\mathrm{2}{A}^{-\mathrm{1}/n}H{\left|{\displaystyle \frac{\partial u}{\partial x}}\right|}^{\mathrm{1}/n-\mathrm{1}}{\displaystyle \frac{\partial u}{\partial x}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{i}}g\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}\right)H^{\mathrm{2}}\end{array}\right\}\\ {\displaystyle }\text{ at }x=x_g.\end{array}\end{array}$$

In the numerical model, we assume that an unbuttressed ice shelf fills the domain from the grounding line to the domain boundary; that is, we apply Eq. (12) at x=L_x instead. This is mathematically equivalent to Eq. (12), which assumes that the ice sheet terminates at the grounding line (e.g. Schoof, 2007a).

Assuming ice overburden pressure equal to water pressure at the grounding line, we set the effective pressure to zero at the grounding line and assume that the till layer is completely unfrozen where it is in contact with the ocean,

$$\begin{array}{}\text{(13)}& N=\mathrm{0}\text{ and }{h}_{\mathrm{s}}=h_{\mathrm{0}}\text{ at }x=x_g.\end{array}$$

For most cases, the model is initialized with a linearly decreasing ice thickness of $H_{\text{initial}}=max(\mathrm{2000}\mathrm{m}-\mathrm{6}\times {\mathrm{10}}^{-\mathrm{3}}x,\mathrm{100}\mathrm{m})$, a linearly decreasing effective pressure $N_{\text{initial}}=(\mathrm{20.85}-{\mathrm{10}}^{-\mathrm{6}}\frac{x}{\mathrm{1}\mathrm{m}})$ MPa, and h_s=h₀. To promote convergence, some runs for $K_{\mathrm{d}}\ge \mathrm{3}\times {\mathrm{10}}^{-\mathrm{2}}$ m s⁻¹ are initialized from solutions previously obtained for other values of K_d.

Figure 2 Parameter space of modelled behaviour and sample profiles. Ice is coloured according to speed, and the underlying slab of hydraulically active sediment is coloured according to effective pressure. Note that the thickness of the hydraulically active sediment is vertically exaggerated by 750 % and does not account for expansion due to void-ratio changes. The red line indicates the extent of the frozen fringe into the sediment (if located at the ice–bed interface, then there is no frozen fringe). Note the upstream travelling activation wave of fast flow in column (a), the downstream travelling activation wave in column (b), and the quasi-simultaneous speed-up of the entire ice stream in columns (c) and (d). Panel (e) shows the steady-state solution. The displayed profiles are selected to illustrate the qualitative differences in ice stream behaviour and match the information in Fig. 4. Animations of these solutions are provided as a video supplement (Haseloff et al., 2025).

3 Results

We are primarily interested in the interaction of subglacial drainage and ice sheet flow, and therefore consider a suite of models over which hydraulic conductivity decreases from infinity to zero. We expect steady-state behaviour at quasi-infinite hydraulic conductivity (e.g. Tsai et al., 2015) and oscillations between slow and fast flow at vanishing hydraulic conductivity (e.g. Robel et al., 2014). Our interest lies in the transition between these limiting cases. In Appendix B, we non-dimensionalize the model to show that the non-dimensional group

$$\begin{array}{}\text{(14a)}& {\displaystyle }\mathit{\kappa }& {\displaystyle }={\displaystyle \frac{\left[q_{\mathrm{w}}\right]/h_{\mathrm{0}}}{\left[u\right]}}=K_{\mathrm{d}}\times {\displaystyle \frac{\mathit{\mu }{N}_{\mathrm{c}}}{{\mathit{\rho }}_{\mathrm{w}}ga{L}_x}}\text{(14b)}& {\displaystyle }& {\displaystyle }\approx K_{\mathrm{d}}\times {\mathrm{10}}^{\mathrm{4}}{\text{ s m}}^{-\mathrm{1}}\end{array}$$

determines the ice stream behaviour (see Fig. 2). It describes the ratio between the scale of water velocity $\left[q_{\mathrm{w}}\right]/h_{\mathrm{0}}$ and the ice velocity [u]. For κ≫1, water moves much faster through the subglacial system than the ice sheet evolves, and vice versa. Note that N_c enters Eq. (14a) through its appearance in Eq. (5), not because of freezing processes.

Ice stream behaviour spans steady streaming (κ→∞), hydraulically controlled oscillations between fast and slow flow (κ≫1), and thermally controlled oscillations (κ≪1). These regimes are distinguished by ice stream morphology, ice discharge, and, where applicable, oscillation frequencies. Notably, hydraulically controlled oscillations are characterized by a quasi-simultaneous speed-up of the entire ice stream in less than 2 years (Fig. 2d), while thermally controlled oscillations exhibit an activation wave, which travels upstream at about 4 km yr⁻¹ during the speed-up phase (Fig. 2a). In the regime between these two limiting cases (κ∼1), ice stream activation can occur by downstream propagation of an activation wave at about 2 km yr⁻¹ (Fig. 2b). We analyse the processes leading to these different regimes next, starting with the fully drained limit (κ→∞).

Figure 3 Profiles of hydraulically controlled marine ice sheets with varying values of the hydraulic conductivity of the bed K_d in Eq. (5). Top row: ice sheet profile and hydraulic potential Φ. As in Fig. 2, a vertically exaggerated sediment layer is coloured according to effective pressure. Second row: basal shear stress τ_b and effective pressure contribution μN to it (note that the vertical axis is linear below μN_c=1 MPa and logarithmic above). Third row: speed. Bottom row: basal melt rate. Note that the hydraulic potential is constant for κ=10⁵ (red line in a₁) and gradually increases for decreasing κ (panels a₂ and a₃). Note that κ=10² is the smallest value of κ for which steady states exist.

3.1 Steady streaming

In the limit of quasi-infinite conductivity of the hydraulic system, water fluxes are high enough to evacuate basal melt water efficiently, so that hydraulic gradients ($\partial \mathrm{\Phi }/\partial x$) cannot develop, as can be seen from rearranging Eq. (5) with Eq. (14a):

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial \mathrm{\Phi }}{\partial x}}=-{\displaystyle \frac{\mathit{\mu }N}{\mathit{\kappa }a{L}_x{h}_{\mathrm{s}}}}{q}_{\mathrm{w}}\to \mathrm{0}\text{ for }\mathit{\kappa }\to \mathrm{\infty }.\end{array}$$

Consequently, $\mathrm{\Phi }={\mathit{\rho }}_{\mathrm{w}}gb-N+{\mathit{\rho }}_{\mathrm{i}}gH=$ constant. At the grounding line, boundary conditions (12) require $\mathrm{\Phi }\left(x_g\right)=\mathrm{0}.$ This requires the effective pressure to adjust such that $N={\mathit{\rho }}_{\mathrm{i}}gH+{\mathit{\rho }}_{\mathrm{w}}gb$; that is, the subglacial water content is set by the ice and bed geometry alone, and is independent of the basal melt rate (which is positive throughout, see Fig. 3d₁). This limit is also referred to as height-above-buoyancy or height-above-floatation model (Van der Veen, 1987; Asay-Davis et al., 2016; Kazmierczak et al., 2022), and it leads to high effective pressure and low water content throughout most of the domain (Fig. 3b₁) apart from a small region near the grounding line. Stable steady-state solutions exist, characterized by high ice thickness and surface slopes (Fig. 3a₁). In steady state, the water flux is given by Eq. (4), i.e. $q_{\mathrm{w}}={\int }_{\mathrm{0}}^x{m}_{\mathrm{b}}\text{ d}x$. Note that effective pressures in excess of the critical pressure for freezing N_c are possible, as water is removed from the sediment by subglacial drainage, rather than freezing. Here, κ→∞ corresponds to K_d≈10 m s⁻¹.

As conductivities decrease from infinity, water is less efficiently evacuated, leading to overall higher water content and lower effective pressures, as shown in Fig. 3b₁₋₃. However, basal shear remains high and the average water content increase from ${\overline{h}}_{\mathrm{w}}\approx \mathrm{0.43}$ to ${\overline{h}}_{\mathrm{w}}\approx \mathrm{0.53}$ m leads to a negligible increase in ice discharge at the grounding line (Fig. 5c). At κ=10², a small region of larger effective stress forms just upstream of the grounding line (Fig. 3b₃), corresponding to a lower water content and marginally higher basal resistance. While this region of elevated effective pressure and basal shear is a subtle feature for the displayed steady state, it becomes critically important for the behaviour at smaller hydraulic conductivities, as steady states cease to exist below a critical value κ⪅100. In this range of κ, basal melt evacuation cannot keep up with melt rates. This leads to a local build up in the basal water content and the periodic disappearance of the region of elevated basal resistance, and hence leads to oscillations between fast and slow flow. We analyse the processes leading to these different dynamics next.

3.2 Hydraulically controlled oscillations

At moderately high hydraulic conductivity, $\mathrm{1}\ll \mathit{\kappa }⪅\mathrm{100}$, basal melt rates are positive in all of the domain, and most of the water is continuously evacuated (Fig. 4a_3,5). However, the hydraulic conductivity is not high enough to evacuate all melt water that is produced at the ice–bed interface, and the water content of the bed slowly increases. This shrinks the region of higher basal resistance close to the grounding line (Fig. 4a₄). Eventually the bed is sufficiently weakened for well-lubricated conditions to exist in all of the ice stream, leading to a brief episode of fast flow, enhanced heat dissipation, and sudden grounding line advance (Fig. 4a at ≈6600 years and Fig. 2d₃). However, a lower effective pressure increases the effective hydraulic conductivity of the bed as $K_{\mathrm{d},\text{eff}}\propto \mathrm{1}/N$, and during this phase more water is discharged than produced. This leads to a strengthening of the bed, eventually restoring the conditions at the beginning of the oscillation cycle. Overall, these oscillations are characterized by short peaks in ice discharge and average water content of the bed (right panels in Fig. 5h–i).

Figure 4 Hovmöller diagrams for four different hydraulic conductivities, decreasing from top to bottom with dashed magenta line marking the grounding-line position. Column 1 shows ice thickness H, column 2 ice velocity u, column 3 basal melt rate m_b, column 4 effective pressure N, and column 5 the scaled water flux $q_{\mathrm{w}}/K_{\mathrm{d}}$. Note that the colour bar of the water flux is the inverted colour bar of the effective pressure to facilitate comparison between the two. Arrows indicate times with profiles in Fig. 2, and animated videos of these solutions are provided as a video supplement (Haseloff et al., 2025).

Reduction of the hydraulic conductivity decreases the period and amplitude of these oscillations (Fig. 5). This is due to the region of elevated effective pressure and basal shear upstream of the grounding line becoming less pronounced (compare the regions of elevated effective pressure in Fig. 4a₄ and b₄). At smaller hydraulic conductivity, the subglacial drainage configuration becomes increasingly uncoupled from the hydraulic gradient while the importance of local melting and freezing processes increases (in fact, as the conductivity approaches zero, the subglacial water content is determined by local processes only, as discussed in Sect. 3.4). The period of hydraulically controlled oscillations decreases from ≈2000 years at κ=50 to a minimum of ≈300 years at κ=1; with further reduction of κ the oscillation period increases again up to a period of ≈800 years at $\mathit{\kappa }={\mathrm{10}}^{-\mathrm{2}}$ (Fig. 5e).

3.3 Intermediate oscillations

For κ∼1, freezing becomes possible, albeit only in a small area close to the grounding line at first (Fig. 4b₃). However, this has no noticeable effect on the effective pressure and hence basal shear stress as the water flux remains high throughout the domain (Fig. 4b₅). Therefore, for κ=1, the ice flow characteristics are similar to those of the hydraulically controlled oscillations, with the entire ice stream speeding up almost simultaneously.

For smaller $\mathit{\kappa }={\mathrm{10}}^{-\mathrm{1}}$, freezing extends up to the ice divide and leads to strengthening of the bed, which alters the oscillation characteristics. In this case, the onset of fast flow is triggered by the thickening of ice upstream and the ensuing reduction in conductive cooling. This behaviour is shown in Figs. 4c₁–c₅ and 2b. At the end of the fast-flowing phase, the ice is at its thinnest, basal freezing is prevalent, and a frozen fringe forms in the downstream half of the ice stream (Fig. 2b₄ to b₁). Gradual thickening of ice upstream leads to a transition from basal freezing to melting, decreasing the effective pressure. This leads to faster ice flow in these upstream areas, and thus the formation of a steep surface gradient at the transition between low and high effective stress. The transition between fast and slow flow travels downstream (Figs. 4c₂ and 2b₁–b₃), similar to surge fronts in mountain glaciers (Kamb et al., 1985). While subglacial water fluxes generally mirror the inverse of the effective pressure for larger values of κ, here, they are elevated ahead of the downstream travelling surge front (Fig. 4c₅). This indicates that meltwater is transported into the region just downstream of the barrier between low and high effective stress, lowering the effective pressure ahead of the surge front, and thereby allowing its propagation.

As in the hydraulically controlled oscillations, the region of elevated effective pressure and basal shear stress downstream of the surge front acts as a buttress to the overall flow. Once this buttress vanishes, the entire ice stream speeds up, leading to grounding line advance and thinning of the ice stream. The thinner ice stream experiences larger conductive cooling rates, leading to freezing, bed strengthening, and slower flow, closing the oscillation cycle.

Figure 5 Bifurcation diagrams for finite κ (a–e) and selected time series after an initial spin up of 6000 years (f–i). Shaded area in plots (a)–(d) correspond to the range covered during an oscillation period; points are minimum and maximum values of the oscillation. Note that hydraulically controlled oscillations are characterized by a sharp peak in ice discharge and subglacial water content (e.g. red curves in plots h and i), while thermally controlled oscillations are changing more gradually (e.g. green curve in plots h and i).

3.4 Thermally controlled oscillations

At vanishing hydraulic conductivity ($\mathit{\kappa }\lesssim {\mathrm{10}}^{-\mathrm{2}}$), the conditions at the bed are fully determined by the local energy balance (7), $\partial h_{\mathrm{w}}/\partial t=m_{\mathrm{b}}$. Steady states would require a locally zero basal melt rate (m_b=0), but this can only occur if conductive cooling exactly balances heat dissipation and geothermal heat flux everywhere. This condition is generally not satisfied. Therefore, the local effective pressure, and hence basal resistance to flow, change on timescales dependent on the factors controlling the basal melt rate – the heat dissipation (τ_bu term in Eq. 7) and the conductive cooling ($k\partial T/\partial z{|}^{+}$ term in Eq. 7), which in turn is controlled by the ice thickness. This leads to thermally controlled oscillations with a period of ≈800 years as previously reported in Robel et al. (2014), see Fig. 4d₁–d₅.

During quasi-stagnant flow phases, the ice is sufficiently thin and slow that melt rates at the bed are negative and part of the sediment freezes. The effective pressure is at the critical effective pressure for freezing N_c in all of the domain (apart from a small boundary layer close to the grounding line which vanishes for κ→0). Slow flow also means that less mass is evacuated than accumulated, and hence the ice sheet thickens over time. This leads to a gradual reduction of the conductive cooling; the energy balance at the bed becomes more positive and the sediment starts to thaw. Once the sediment is fully thawed, further melt produced at the bed is stored locally, leading to a reduction of the effective pressure. As the ice is thicker towards the ice divide, this process happens first there. However, the flow remains stabilized by the partially frozen sediment in the downstream half of the ice stream.

Eventually, the transition to the fast streaming phase is triggered at the grounding line, where an increase in ice thickness and driving stress leads to speed-up (Fig. 4d at 6700 years and Fig. 2a₁–a₂). This increases the heat dissipation, lowers the basal resistance, and thus triggers a positive feedback between fast flow, heat dissipation, and basal lubrication. Force balance dictates that loss of basal resistance to flow affects the flow in a small boundary layer immediately upstream through longitudinal stresses ($\mathit{\eta }\partial u/\partial x$, see, e.g., Fowler and Schiavi, 1998; Robel et al., 2014). At high effective stress, basal melt rates in this boundary layer are high (≈10 mm yr⁻¹), quickly lowering the effective pressure, which then speeds up the flow immediately upstream. This creates an activation wave that travels upstream, which is clearly visible in the basal melt rate (Fig. 4d₃); its properties are discussed in detail in Robel et al. (2014). In its wake the ice stream bed is well-lubricated, leading to discharge of ice in excess of accumulation, thinning of the ice, and consequently increasing conductive cooling at the bed. This eventually leads to strengthening of the bed, slow down of flow, and transition to the stagnant phase. The frequency of these oscillations is primarily determined by how fast the ice stream can replenish its mass (i.e. the accumulation rate) and the terms in the basal melt rate (Eq. 7) (i.e. k, T_s, and q_geo, MacAyeal, 1993; Robel et al., 2013).

4 Discussion

In this study, we investigate the role of an evolving subglacial drainage system on marine ice stream dynamics using the simplest model of the relevant physics. A positive feedback between fast flow and basal lubrication is incorporated through both increasing local storage of water with decreasing effective pressure and increasing drainage efficiency with decreasing effective pressure. Subglacial conditions affect marine ice sheet dynamics through the basal shear stress, which depends on the effective pressure.

We find that the efficiency of a distributed subglacial drainage system alters the mode of grounding line dynamics. Lower hydraulic conductivities result in more water storage, lower effective pressure, and consequently lower basal resistance and faster flow, i.e. generally higher slip rates, thinner ice, and less grounded ice volume (Fig. 5b). With lower hydraulic conductivities, the melt rate at the bed gradually drops, from positive at all points in space and time (due to thicker ice mitigating conductive cooling) to zero on average at the impermeable limit, see Fig. 4a₄-d₄. This behaviour is accompanied by a transition from a quasi-simultaneous speed-up of the entire ice stream to a gradual speed-up and slow down (Fig. 5h). This speed-up can be due to downstream and upstream-travelling activation waves.

4.1 Ice stream oscillations and glacier surges

Observations suggest that the Siple Coast ice streams show temporal variability on timescales of 300 to 500 years (e.g. Retzlaff and Bentley, 1993; Hulbe and Fahnestock, 2007; Catania et al., 2012). While previous models have been able to reproduce stagnation and activation cycles, the period of oscillations in these models is generally on the order of 1000 years or longer (e.g. MacAyeal, 1989; Robel et al., 2013, 2014). Here, we illustrate that more frequent oscillations might be modulated by subglacial hydrology.

The similarities of the marine ice streams modelled here and the observed behaviour of surging mountain glaciers corroborate the hypothesis that the same dynamics govern both of these systems, albeit at different spatial and temporal scales (e.g. Clarke et al., 1984). In the regime intermediate between hydraulic and thermal control, the oscillations share similarities with the observed classical surges of Variegated Glacier, Trapridge Glacier and others: a downstream-travelling surge front, acceleration by up to three orders of magnitude, and a basal transition between upstream temperate conditions and downstream near-frozen bed conditions (Clarke et al., 1984; Kamb et al., 1985). Recent observations have also shown that some surges are instead characterized by activation waves travelling upglacier (Murray et al., 2003; Sevestre et al., 2018), leading to the hypothesis that there are at least two distinct surge mechanisms (Murray et al., 2003). Our results suggest that these different surge dynamics could be due to small variations in effective basal conductivity of the bed.

Previous conceptual, non-spatial (“lumped”) models for surging glaciers (Fowler et al., 2001; Benn et al., 2019, 2022; Terleth et al., 2024) include the same ingredients as our model. These models show that the existence of oscillations also depends on climatic (accumulation, surface temperature) and geometric factors (bed slope). While our focus is on the dynamic behaviour resulting from variable hydraulic conductivity, it is reasonable to assume that variation of climatic and geometric factors will have an equally important effect on marine ice stream dynamics in our model. This is supported by the appearance of the accumulation rate a and the characteristic ice sheet extent L_x in the ratio of water velocity to ice velocity κ (Eq. 14a). Conversely, we supplement existing, zero-dimensional models for glacier surges with a spatially extended model, enabling identification and analysis of previously inaccessible complexities.

4.2 Marine ice sheet dynamics

Marine ice sheet dynamics are typically framed to be either stable or unstable, with instability equated to marine ice sheet collapse. Our results indicate that such classifications do not always apply when basal boundary conditions dynamically evolve. Instead, marine ice streams might undergo periodic cycles of retreat and re-advance. Previous studies investigating grounding line dynamics with evolving boundary conditions have also shown that thermal oscillations can temporarily stabilize grounding lines on retrograde slopes (Robel et al., 2016).

Most models simulating future ice loss under warming scenarios do not account for evolving basal boundary conditions. Our results show that basal conditions can change on decadal to centennial timescales. Faster oscillations are likely possible for different parameter choices. However, incorporation of subglacial drainage into predictive ice sheet models is still hampered by uncertainty in how subglacial drainage should be upscaled to large spatial and temporal scales, as well as the computational cost of incorporating these dynamics.

4.3 Subglacial hydrology of ice streams

Given the control that effective conductivity of the bed exerts on ice stream dynamics, determining appropriate parameterizations of the effective hydraulic conductivity is crucial. The existence of subglacial till beneath parts of West Antarctica is well established (Blankenship et al., 1986, 1987; Alley et al., 1986, 1987; Peters et al., 2006; Schroeder et al., 2014). Based on these observations, Walder and Fowler (1994) and Ng (1998, 2000b) propose the existence of subglacial canals partially eroded into the bed and ice. In settings typical for ice streams, the dynamics of these conduits mimic those typical for distributed subglacial systems where the effective pressure decreases with increasing water discharge. Assuming the same qualitative behaviour as Walder and Fowler (1994), we model subglacial water flow following a Darcian-style description (e.g. Hewitt, 2011; Flowers, 2015) with K_d,eff the effective hydraulic conductivity. In our model, κ=1 corresponds to $K_{\mathrm{d}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹ and at N=100 kPa we get $K_{\mathrm{d},\text{eff}}=K_{\mathrm{d}}{N}_{\mathrm{c}}/N=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹. This is consistent with in-situ measurements of the hydraulic conductivity at the base of Whillans Ice stream, where a minimum of K=(0.5 to $\mathrm{1.4})\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹ has been estimated (Engelhardt and Kamb, 1997).

However, such high hydraulic conductivities are inconsistent with subglacial water flow within the void space of the till alone. For example, for subglacial tills derived from Whillans ice stream, Leeman et al. (2016) find conductivities around 10⁻¹² m s⁻¹. Measurements in the field report values around $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{9}}$ m s⁻¹ (Engelhardt et al., 1990). Moreover, the Kozeny–Carman relationship $K\propto e^{\mathrm{3}}/(\mathrm{1}+e)$ is typically used for flow of water in the void space e of sediment (e.g. Lambe and Whitman, 1991). When combined with the measured dependence of the void ratio on effective pressure (Eq. 9), a different dependence of K on N arises than what we assumed here. This might also alter the ice stream dynamics. In our hydraulically-controlled oscillations, the increase of drainage efficiency with lower N causes the termination of the fast flow phase, as more water is evacuated as the bed becomes more lubricated. In the Kozeny–Carman relationship, the dependence on N is weaker than in our model, which might not permit the dynamics seen here.

Assuming that intermediate κ are representative for West Antarctic ice stream dynamics, alternative pathways must facilitate a more efficient water transport. Engelhardt and Kamb (1997) suggest that a distributed system of sediment-incised canals, as proposed by Walder and Fowler (1994) is most compatible with their observations. Between these canals, water might flow in a thin film at the ice bed-interface, rather than in the subglacial sediment, further enhancing the hydraulic conductivity (Creyts and Schoof, 2009).

For high water discharge conditions, it is plausible that a transition from distributed to channelized drainage might occur (Walder and Fowler, 1994; Röthlisberger, 1972; Hager et al., 2022; Dow et al., 2022). Under these conditions, we expect the net discharge to increase with effective pressure, rather than decrease. In Greenland, where abundant surface meltwater can drain to the bed, increasing channelization leads to a strengthening of the bed at the end of the melt season (Schoof, 2010; Sundal et al., 2011; Bartholomew et al., 2012). We also expect the conductivity to vary locally, so that we interpret K_d,eff as an effective hydraulic conductivity, averaged over spatial scales resolvable in ice sheet models.

We have assumed that the basal shear stress adjusts instantaneously to changes in subglacial water content, as we are interested in the annual to centennial evolution of marine ice sheets. However, at the grounding line, marine ice sheets are subject to tidal forcing on much shorter timescales. At such short timescales, transient dilatant strengthening of subglacial till can lead to time-dependent sliding laws and affect grounding line motion (Warburton et al., 2020, 2023) It has been hypothesized that time-dependent sliding laws might also modulate surging of glaciers (Minchew and Meyer, 2020); this behaviour might be important for the quasi-simultaneous speed-up we observed in the hydraulically controlled limit at moderately high hydraulic conductivities.

4.4 Model limitations

In the formulation of our model, we have made simplifying assumptions that limit its applicability. Most notably, we use a width-integrated model with a parameterized lateral drag term. Assuming that the width of the ice stream is constant allows us to ignore the complex dynamics required to capture its evolution. However, observations indicate that ice stream widths can change over time (e.g. Harrison et al., 1998; Echelmeyer and Harrison, 1999). Modelling also shows that substantial englacial melt is produced in ice stream margins, which could affect ice stream hydrology and dynamics (Raymond, 2000; Suckale et al., 2014; Haseloff et al., 2019). Future modelling efforts should build on our work to account for the hydraulic and thermal processes affecting ice stream width.

Another limitation of the width-averaged approach is that it precludes the formation of ice ridges, whose surface slopes affect the ice stream stability by altering the hydraulic gradient (Kyrke-Smith et al., 2014) and the interaction of neighbouring ice streams. For example, it has been hypothesized that changes in subglacial water pathways might have contributed to the shutdown of Kamb ice stream approximately 140 years ago (e.g. Anandakrishnan and Alley, 1997). The advantage is that we are left with a minimal model of a marine ice stream in which dissipation, drainage, and freezing processes compete to shape the ice stream thickness and velocity field, and we can study the effect of different physical processes in the model.

We also ignore the role of ice shelves here. Buttressing ice shelves can alter marine ice stream dynamics (Gudmundsson et al., 2012; Gudmundsson, 2013; Haseloff and Sergienko, 2018, 2022). In addition, at the grounding line, subglacial hydrology can interact with the ocean in multiple ways not taken into account here, for example through seawater infiltration and the initialization of subshelf meltwater plumes which might alter ice shelf buttressing (e.g. Jenkins, 2011; Robel et al., 2022; Ehrenfeucht et al., 2024). More generally, where marine ice streams are in contact with stochastically-varying ocean and atmospheric conditions, their dynamics can significantly diverge from purely deterministic results (e.g. Mantelli et al., 2016; Robel et al., 2018; Christian et al., 2022; Sergienko and Haseloff, 2023). Such effects should be addressed in future studies.

Another model limitation is that we do not allow for full freezing of the bed; in all simulations reported here the bed remains temperate. Observations suggest temperate ice stream beds (e.g. Engelhardt and Kamb, 1997) but the bed might be frozen under the adjacent ice ridges (e.g. Bentley et al., 1998). Frozen–temperate transitions might also play a role in ice stream onset (Mantelli et al., 2019), although in our model upstream regions usually have higher water content, as the thicker ice there reduces conductive cooling of the bed. The formation of temperate to frozen boundaries is not possible in our model, and additional model physics would be required to capture these (Haseloff et al., 2015, 2018; Mantelli and Schoof, 2019; Schoof and Mantelli, 2021).

5 Conclusions

In this paper we showed that subglacial hydrology affects the stability, extent, and grounded ice volume of marine ice streams. Different dynamical regimes, including steady streaming, hydraulically controlled, quasi-simultaneous advances, and periodic oscillations with advancing or retreating activation fronts are possible for different hydraulic conductivities of the bed. Our results illustrate that marine ice streams can undergo internal variability at periods from a few centuries to millennia. This variability may need to be taken into account for accurate prediction of future sea level contributions on these timescales.

Appendix A: Model implementation and verification

The model is implemented in PETSc (Balay et al., 2023). PETSc's SNES library for solution of nonlinear systems of algebraic equations is used to solve individual equations. The equations for momentum balance (Eq. 2) and subglacial drainage (Eqs. 4–9) are discretized using conservative finite differences with implicit time steps (Katz et al., 2007). The equation for mass balance (Eq. 1) is discretized with a third-order upwind scheme. We also use a linear subgrid interpolation at the grounding line (Pattyn, 2003).

At each time step, we use a segregated loop to solve our model. That is, the equations for mass balance (H), momentum balance (u), and subglacial water content (h_s and h_w) are solved in three independent steps. To ensure that the combined solution at the current time step is converged, the steps are iterated using a Picard scheme until convergence is achieved. That is, at each time step t_i we calculate $u\left(t_i^{\mathrm{1}}\right)$, $H(t_i^{\mathrm{1}},x)$, $h_{\mathrm{s}}(t_i^{\mathrm{1}},x)$, and $h_{\mathrm{w}}(t_i^{\mathrm{1}},x)$ in the first iteration and use these to continue to calculate $u(t_i^{\mathrm{2}},x)$, $H(t_i^{\mathrm{2}},x)$, etc. until further iteration does not alter the solution at the current time step. The size of the time step is adjusted to achieve convergence of the individual SNES solvers as well as convergence of the segregated solution in no more than 10 Picard iterations.

Figure A1 Numerical convergence study. Panels (a)–(b) show grounding line position x_g and average water content ${\overline{h}}_{\mathrm{w}}$ under grid refinement for steady-state solutions (κ=∞); panels (c)–(d) show time series of the same for oscillating solutions (κ=10).

The implemented scheme is of O(Δx) accuracy, and comparison of steady-state results and time series under grid refinement are shown in Fig. A1. While the qualitative behaviour is consistent for values of Δx≲1000 m, convergence requires resolution at 10s of metres. Solutions shown here were calculated at resolutions of Δx=100 m or finer. Simulation times depend on the hydraulic conductivity parameter and, of course, the model resolution. For κ≲1 and a grid spacing of Δx=100 m, 10⁴ model years are computed in about 24 h on a single core. Simulations with κ≫1 require significantly more time (up to two weeks) due to limitations on the time step.

Appendix B: Non-dimensional model

Here, we non-dimensionalize Eqs. (1)–(4) to highlight scaling behaviour in each balance region. There are two limiting cases. For hydraulically controlled ice streams, ice is relatively thick, velocities are slow, and basal resistance is high, with both the Weertman-like terms and the Coulomb-like terms in Eq. (3) contributing to the basal shear stress. Subglacial water flow is efficient and the sediment does not experience freezing. We consider this limit first in Sect. B1. For thermally controlled ice streams, ice is thin and basal resistance is primarily of the Coulomb-plastic type. The basal energy balance oscillates between melting and freezing; this limit is considered in Sect. B2.

B1 Hydraulically controlled ice streams

In this limit, we choose the following scales with characteristic values as indicated,

$$\begin{array}{}\text{(B1a)}& {\displaystyle }\left[x\right]& {\displaystyle }=L_x=\mathrm{1000}\text{ km},\text{(B1b)}& {\displaystyle }\left[u\right]& {\displaystyle }={\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}g{a}^{\mathrm{2}}\left[x\right]}{C}}\right)}^{n/(\mathrm{2}n+\mathrm{1})}\approx \mathrm{100}\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}},\text{(B1c)}& {\displaystyle }\left[H\right]& {\displaystyle }={\displaystyle \frac{a\left[x\right]}{\left[u\right]}}\approx \mathrm{5000}\text{m},\text{(B1d)}& {\displaystyle }\left[t\right]& {\displaystyle }={\displaystyle \frac{\left[H\right]}{a}}\approx {\mathrm{10}}^{\mathrm{4}}\text{years},\text{(B1e)}& {\displaystyle }\left[N\right]& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\mu }}}C[u{]}^{\mathrm{1}/n}\approx \mathrm{500}\text{kPa}.\end{array}$$

We set $x=\left[x\right]x^{*}$, $H=\left[H\right]H^{*}$, $u=\left[u\right]u^{*}$, $t=\left[t\right]t^{*}$, $N=\left[N\right]N^{*}$, $h_{\mathrm{s}}=h_{\mathrm{0}}{h}_{\mathrm{s}}^{*}$, and $\mathrm{\Phi }={\mathit{\rho }}_{\mathrm{i}}g\left[H\right]{\mathrm{\Phi }}^{*}$. Dropping asterisks immediately, we obtain the non-dimensional equations for momentum balance (Eq. 2), mass balance (Eq. 1), and water balance (Eq. 4),

$$\begin{array}{}\text{(B2a)}& {\displaystyle }& {\displaystyle }\begin{array}{r}{\displaystyle }\mathrm{2}{\mathit{\alpha }}_{\mathrm{1}}{\displaystyle \frac{\partial }{\partial x}}\left(H{\left|{\displaystyle \frac{\partial u}{\partial x}}\right|}^{\mathrm{1}/n-\mathrm{1}}{\displaystyle \frac{\partial u}{\partial x}}\right)-{\mathit{\alpha }}_{\mathrm{2}}H{u}^{\mathrm{1}/n}\\ {\displaystyle }-{\displaystyle \frac{N{u}^{\mathrm{1}/n}}{u^{\mathrm{1}/n}+N}}{\displaystyle \frac{u}{\left|u\right|}}=H{\displaystyle \frac{\partial H}{\partial x}},\end{array}\text{(B2b)}& {\displaystyle }& {\displaystyle \frac{\partial H}{\partial t}}+{\displaystyle \frac{\partial \left(Hu\right)}{\partial x}}=\mathrm{1},\text{(B2c)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \left(e{h}_{\mathrm{s}}\right)}{\partial t}}& {\displaystyle }+\mathit{\kappa }{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial x}}\right)\\ {\displaystyle }& {\displaystyle }=\left({\mathit{\alpha }}_{\mathrm{3}}+{\mathit{\alpha }}_{\mathrm{4}}{\displaystyle \frac{N{u}^{\mathrm{1}/n+\mathrm{1}}}{u^{\mathrm{1}/n}+N}}-{\mathit{\alpha }}_{\mathrm{5}}{\displaystyle \frac{\mathrm{1}}{H}}\right),\end{array}\end{array}$$

with the non-dimensional groups

$$\begin{array}{}\text{(B3a)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\left[H\right]}{A^{\mathrm{1}/n}C{L}_x^{\mathrm{1}/n+\mathrm{1}}}}\approx \mathrm{6}\times {\mathrm{10}}^{-\mathrm{4}},\text{(B3b)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{C_{\mathrm{w}}\left[H\right]}{A^{\mathrm{1}/n}C{W}^{\mathrm{1}/n+\mathrm{1}}}}\approx {\mathrm{10}}^{-\mathrm{1}},\text{(B3c)}& {\displaystyle }\mathit{\kappa }& {\displaystyle }=K_{\mathrm{d}}{\displaystyle \frac{\mathit{\mu }{N}_{\mathrm{c}}}{{\mathit{\rho }}_{\mathrm{w}}g{L}_{x}a}}\approx K_{\mathrm{d}}\times \mathrm{6}\times {\mathrm{10}}^{\mathrm{3}}{\text{ s m}}^{-\mathrm{1}},\text{(B3d)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{3}}& {\displaystyle }={\displaystyle \frac{q_{\text{geo}}\left[H\right]}{{\mathit{\rho }}_{\mathrm{w}}{L}_h{h}_{\mathrm{0}}a}}\approx \mathrm{60},\text{(B3e)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{4}}& {\displaystyle }={\displaystyle \frac{C\left[u{]}^{\mathrm{1}/n}\right[x]}{{\mathit{\rho }}_{\mathrm{w}}{L}_h{h}_{\mathrm{0}}}}\approx \mathrm{700},\text{(B3f)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{5}}& {\displaystyle }={\displaystyle \frac{k\mathrm{\Delta }T}{{\mathit{\rho }}_{\mathrm{w}}{L}_h{h}_{\mathrm{0}}a}}\approx \mathrm{8}.\end{array}$$

With ${\mathit{\alpha }}_{\mathrm{1}}\ll {\mathit{\alpha }}_{\mathrm{2}}\ll \mathrm{1}$, the leading-order momentum balance is between basal resistance and driving stress,

$$\begin{array}{}\text{(B4)}& {\displaystyle }-{\displaystyle \frac{N{u}^{\mathrm{1}/n}}{u^{\mathrm{1}/n}+N}}{\displaystyle \frac{u}{\left|u\right|}}{\displaystyle }\sim H{\displaystyle \frac{\partial H}{\partial x}}.\end{array}$$

In the basal water balance, the non-dimensional group of interest is κ, which controls the importance of subglacial hydrology. It can be understood at the ratio of subglacial water velocity to ice velocity, as

$$\begin{array}{}\text{(B5)}& \mathit{\kappa }={\displaystyle \frac{\left[q_{\mathrm{w}}\right]/h_{\mathrm{0}}}{\left[u\right]}}\text{ with }\left[q_{\mathrm{w}}\right]=K_{\mathrm{d}}{\displaystyle \frac{N_{\mathrm{c}}{h}_{\mathrm{0}}\left[\mathrm{\Phi }\right]}{{\mathit{\rho }}_{\mathrm{w}}g\left[N\right]\left[x\right]}}.\end{array}$$

If $K_{\mathrm{d}}\sim \left({\mathit{\rho }}_{\mathrm{w}}g{L}_{x}a\right)/\left(\mathit{\mu }{N}_{\mathrm{c}}\right)$, then κ∼1. We also note that α₅≪α₄, i.e. freezing is negligible.

For κ→∞, we expect steady-states with constant hydraulic potential as detailed in Eq. (15). Practically, this limit is attained for K_d∼10 m s⁻¹ (corresponding to κ≈10⁵). As long as κ≫1, we expect subglacial hydrology to be efficient; we call this limit hydraulically controlled. In this case, the leading-order water balance is

$$\begin{array}{}\text{(B6)}& \mathit{\kappa }{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial x}}\right)\sim {\mathit{\alpha }}_{\mathrm{4}}{\displaystyle \frac{N{u}^{\mathrm{1}/n+\mathrm{1}}}{u^{\mathrm{1}/n}+N}}.\end{array}$$

For κ≪1 we expect subglacial hydrology to become increasingly inefficient and the interplay between local melt production and freezing to dominate the basal energy balance. This requires thinner ice streams; we rescale below to capture this asymptotic limit.

The limit κ∼1 describes the intermediate regime. In this limit, the relative magnitudes of subglacial water flux, geothermal heat flux, basal energy dissipation, and conductive cooling determine the energy balance.

B2 Thermally controlled ice streams

In this limit, we choose the following scales

$$\begin{array}{}\text{(B7a)}& {\displaystyle }\left[x\right]& {\displaystyle }=L_x=\mathrm{1000}\text{km},\text{(B7b)}& {\displaystyle }\left[H\right]& {\displaystyle }={\displaystyle \frac{k\mathrm{\Delta }T}{q_{\text{geo}}}}\approx \mathrm{600}\text{m},\text{(B7c)}& {\displaystyle }\left[u\right]& {\displaystyle }={\displaystyle \frac{a\left[x\right]}{\left[H\right]}}\approx \mathrm{800}{\mathrm{m}\mathrm{yr}}^{-\mathrm{1}},\text{(B7d)}& {\displaystyle }\left[t\right]& {\displaystyle }={\displaystyle \frac{\left[H\right]}{a}}\approx \mathrm{1200}\text{years},\text{(B7e)}& {\displaystyle }\left[N\right]& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathit{\mu }}}{\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{[H{]}^{\mathrm{2}}}{\left[x\right]}}\approx \mathrm{7}\text{kPa}.\end{array}$$

We set $x=\left[x\right]x^{†}$, $H=\left[H\right]H^{†}$, $u=\left[u\right]u^{†}$, $t=\left[t\right]t^{†}$, $N=\left[N\right]N^{†}$, $h_{\mathrm{s}}=h_{\mathrm{0}}{h}_{\mathrm{s}}^{*}$, and $\mathrm{\Phi }={\mathit{\rho }}_{\mathrm{i}}g\left[H\right]{\mathrm{\Phi }}^{†}$. Dropping daggers immediately, we obtain the non-dimensional equations for mass balance (Eq. 1), momentum balance (Eq. 2), and water balance (Eq. 4),

$$\begin{array}{}\text{(B8a)}& {\displaystyle }& {\displaystyle }\begin{array}{r}{\displaystyle }\mathrm{2}{\mathit{\gamma }}_{\mathrm{1}}{\displaystyle \frac{\partial }{\partial x}}\left(H{\left|{\displaystyle \frac{\partial u}{\partial x}}\right|}^{\mathrm{1}/n-\mathrm{1}}{\displaystyle \frac{\partial u}{\partial x}}\right)-{\mathit{\gamma }}_{\mathrm{2}}H{u}^{\mathrm{1}/n}\\ {\displaystyle }-{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{3}}N{u}^{\mathrm{1}/n}}{{\mathit{\gamma }}_{\mathrm{3}}{u}^{\mathrm{1}/n}+N}}{\displaystyle \frac{u}{\left|u\right|}}=H{\displaystyle \frac{\partial H}{\partial x}},\end{array}\text{(B8b)}& {\displaystyle }& {\displaystyle \frac{\partial H}{\partial t}}+{\displaystyle \frac{\partial \left(uH\right)}{\partial x}}=\mathrm{1},\text{(B8c)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \left(e{h}_{\mathrm{s}}\right)}{\partial t}}& {\displaystyle }+\mathit{\kappa }{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial x}}\right)\\ {\displaystyle }& {\displaystyle }={\mathit{\gamma }}_{\mathrm{4}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{H}}\right)+{\mathit{\gamma }}_{\mathrm{5}}{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{3}}N{u}^{\mathrm{1}/n+\mathrm{1}}}{{\mathit{\gamma }}_{\mathrm{3}}{u}^{\mathrm{1}/n}+N}},\end{array}\end{array}$$

with non-dimensional groups

$$\begin{array}{}\text{(B9a)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{a^{\mathrm{1}/n}}{{\mathit{\rho }}_{\mathrm{i}}g{A}^{\mathrm{1}/n}[H{]}^{\mathrm{1}/n+\mathrm{1}}}}\approx {\mathrm{10}}^{-\mathrm{2}},\text{(B9b)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{2}}& {\displaystyle }=C_{\mathrm{w}}{\left({\displaystyle \frac{\left[x\right]}{W}}\right)}^{\mathrm{1}/n+\mathrm{1}}{\mathit{\gamma }}_{\mathrm{1}}\approx \mathrm{1.7},\text{(B9c)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{3}}& {\displaystyle }={\displaystyle \frac{C{a}^{\mathrm{1}/n}[x{]}^{\mathrm{1}/n+\mathrm{1}}}{\mathit{\mu }{\mathit{\rho }}_{\mathrm{i}}g[H{]}^{\mathrm{2}+\mathrm{1}/n}}}\approx \mathrm{140},\text{(B9d)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{4}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}{L}_h}}{\displaystyle \frac{k\mathrm{\Delta }T}{a{h}_{\mathrm{0}}}}\approx \mathrm{8},\text{(B9e)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{5}}& {\displaystyle }={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}g}{{\mathit{\rho }}_{\mathrm{w}}{L}_h{h}_{\mathrm{0}}}}{\left({\displaystyle \frac{k\mathrm{\Delta }T}{q_{\text{geo}}}}\right)}^{\mathrm{2}}\approx \mathrm{10}.\end{array}$$

As before, the dimensional group determining the role of subglacial hydrology is κ. Subglacial water flux is negligible for κ≪1, which describes thermally controlled oscillations. In this limit the leading-order water balance (neglecting terms of $O(\mathrm{1}/{\mathit{\gamma }}_{\mathrm{3}})$) is

$$\begin{array}{}\text{(B10)}& {\displaystyle \frac{\partial \left(e{h}_{\mathrm{s}}\right)}{\partial t}}\sim {\mathit{\gamma }}_{\mathrm{4}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{H}}\right)+{\mathit{\gamma }}_{\mathrm{5}}Nu.\end{array}$$

This limit is attained for $K_{\mathrm{d}}⪅{\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹.