the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Surging of a Hudson Straitscale ice stream: subglacial hydrology matters but the process details mostly do not
Lev Tarasov
While subglacial hydrology is known to play a role in glacial dynamics on subannual to decadal scales, it remains unclear whether subglacial hydrology plays a critical role in ice sheet evolution on centennial or longer timescales. Furthermore, several drainage systems have been inferred, but it is unclear which is most applicable at the continental/glacial scale. More fundamentally, it is even unclear if the structural choice of subglacial hydrology truly matters for this context.
Here we compare the contribution to the surge behaviour of an idealized Hudson Straitlike ice stream from three subglacial hydrology systems. We use the newly updated BAsal Hydrology Model – BrAHMs2.0 – and provide model verification tests. BrAHMs2.0 incorporates two processbased representations of inefficient drainage dominant in the literature (linked cavity and poroelastic) and a nonmassconserving zerodimensional form (herein termed leaky bucket) coupled to an ice sheet systems model (the Glacial Systems Model, GSM). The linkedcavity and poroelastic configurations include an efficient drainage scheme while the leaky bucket does not. All three systems have a positive feedback on ice velocity, whereby faster basal velocities increase melt supply. The poroelastic and leakybucket systems have diagnostic effective pressure relationships – only the linkedcavity system has an additional negative feedback, whereby faster basal ice velocities increase the dynamical effective pressure due to higher cavity opening rates. We examine the contribution of mass transport, efficient drainage, and the linkedcavity negative feedback to surging. We also assess the likely bounds on poorly constrained subglacial hydrology parameters and adopt an ensemble approach to study their impact and interactions within those bounds.
We find that subglacial hydrology is an important system inductance for realistic ice stream surging but that the three formulations all exhibit similar surge behaviour within parametric uncertainties. Even a detail as fundamental as massconserving transport of subglacial water is not necessary for simulating a full range of surge frequency and amplitude. However, one difference is apparent: the combined positive and negative feedbacks of the linkedcavity system yield longer duration surges and a broader range of effective pressures than its poroelastic and leakybucket counterparts.
 Article
(5618 KB)  Fulltext XML
 BibTeX
 EndNote
The role of subglacial hydrology at timescales longer than multiple decades and at ice sheet spatial scales is unclear. Previous studies have inferred that subglacial hydrology plays a strong role in internally (e.g. Siegfried et al., 2016) and externally (e.g. Joughin et al., 1996; Cook et al., 2021) driven ice sheet variability on subannual to multidecadal timescales (Retzlaff and Bentley, 1993; Alley et al., 1994; Ou, 2021; Bennett, 2003). Observations beyond these timescales do not exist.
Several subglacial hydrologic systems have been conceptualized (Flowers, 2015). Constraint of the role of hydrological systems is further challenged by the large parametric uncertainties for all choices of drainage system. For example, the bounds of hydraulic conductivity vary over several orders of magnitude and according to the particular system (Werder et al., 2013). These uncertainties hinder widespread adoption of subglacial hydrology models in Earth systems models in general and glacialcyclescale ice sheet modelling in particular (Flowers, 2018). As such, what is needed to adequately incorporate the subglacial hydrologic system into glacial cycle simulations is not understood.
We ask a basic question: does subglacial hydrology matter on longer than decadal timescales? And if so, to what extent are the structural details of the hydrological system important for this context, especially given the rest of the system uncertainties? Taking a modelling approach, we focus these broad questions on the following: is subglacial hydrology needed to capture Hudson Straitscale ice stream cyclicity? If so, should effective pressure be dynamically determined – based on fully massconserving lateral drainage? Or does a zerodimensional meltwater volume balance with a diagnostic pressure closure suffice? Turning to the parametric uncertainties, which are most important?
Previous modelbased tests of Hudson Strait ice stream surging (e.g. Calov et al., 2010; Payne et al., 2000; Payne and Dongelmans, 1997; MacAyeal, 1993) have focused on thermomechanical feedbacks but omitted the contribution from the subglacial hydrological system. While these studies capture surges in their simulations based on these limited feedbacks, all models except one (model (d), Calov et al., 2010) implemented an abrupt transition at the frozen–temperate thermal boundary, suddenly initiating largescale sliding within a grid cell. This abrupt thermal transition is physically unrealistic at the scales of ice sheet modelling: a region equivalent to a glacialcyclescale ice sheet model grid cell (100–2500 km^{2}) does not become instantly warm based with wholesale transition to basal sliding. Instead, the streaming portions of ice sheets transition to fastersliding velocities as their coupling to the bed (effective pressure) decreases. Subglacial hydrology is a potentially critical piece of the binge–purge conceptual model of internal oscillations (MacAyeal, 1993) as heat production from sliding and deformation work generates meltwater.
Here we examine the contribution to ice sheet internal oscillations from the three most dominant forms of distributed subglacial hydrology – linked cavity (Schoof, 2010), poroelastic (Flowers, 2000), and nonmasstransporting leaky bucket (Gandy et al., 2019) – relative to each other and to no hydrology at all. In each case, the frozentotemperate transition is smoothed following the work of Hank et al. (2023). We couple these processes to an ice sheet systems model, the Glacial Systems Model (GSM, Tarasov et al., 2012).
Simple configurations make system behaviours more interpretable (e.g. Calov et al., 2010; Payne et al., 2000). With a realistic bed and actual climate, spatiotemporal variations in model solutions are largely due to the variation in boundary conditions. We therefore model these coupled systems for a simplified North American analogue setup which implements a square bed and flat topography with soft beds in the southern latitudes and in the Hudson Strait and Bay area. The ice sheet is forced with a steady climate and firstorder feedbacks: the northwardcooling temperature trend, vertical lapse rate, and thermodynamic moisture control. The numerical model retains important processes while still being feasible to run large ensembles over a glacial cycle on continental scales to probe parametric uncertainties.
Below, we first test the BAsal Hydrology Model (BrAHMs). This includes a demonstration of the mass conservation, convergence, and symmetry of BrAHMs2.0 and verification of its solutions against another prominent model, the Glacier Drainage System model (GlaDS, (Werder et al., 2013). Next we show the sensitivity of ice sheet geometry to subglacial hydrologic parameters in comparison with climate and ice sheet parameters. Finally, we compare results from a set of four large ensembles (between 11 000 and 20 000 members each) using no hydrology and linkedcavity, poroelastic, and leakybucket hydrology.
In the context of continentalscale ice sheet modelling, resolving individual drainage elements and multiple topologies present within the domain is not computationally feasible. In this section we provide a brief overview of some structural choices made by others and present the options compared in this study, beginning with the current understanding of subglacial hydrology and progressing to increasingly approximate representations of it.
Water in the subglacial system flows either through inefficient drainage systems (pressure ∝ flux) or efficient drainage systems (pressure ∝ flux^{−1}, Flowers, 2015). Inefficient distributed networks are widespread under temperate areas of ice sheets, whereas efficient channel networks are discrete, localized elements. Each class evolves to the other and the change is controlled by system throughput, i.e. water flux. When the flux in an inefficient system rises above a threshold, the system transitions to efficient drainage. When the efficient system flux falls through different lower flux threshold, the system transitions to inefficient drainage, resulting in hysteresis (Schoof, 2010). Any masstransporting hydrology model should have three main components: mass conservation describing transport, a flow law describing flux as a function of hydraulic gradient and subglacial water thickness, and a pressure closure relationship.
2.1 Inefficient flow
In the inefficient drainage regime, flux and water pressure rise together. Several inefficient drainage systems have been theorized: thin film, poroelastic media, and linked cavities. Of these, poroelastic and linkedcavity systems (e.g. Flowers, 2000; Walder, 1986) dominate recently published models (Flowers, 2015; de Fleurian et al., 2018), and As such, these are the two systems we model and contrast herein.
In the poroelastic formulation, water can drain through the pore space of some permeable surficial material (e.g. till). Increasing subglacial water pressures expand the pore space and modify the permeability of the porous medium to flowing water. The conceptual basis for this system is examined in greater detail by Flowers and Clarke (2002). The pressure closure has no theoretical basis and is based on a power law with empirically constrained parameters (Flowers, 2000).
In the linkedcavity system, cavities within the base of the ice open up as basal ice flows over and around bed protrusions – fast flow and larger objects beget larger cavities (Kamb, 1987). As cavities grow larger and numerous they form a connected network linked through smaller orifices giving a tortuous drainage network.
The substrate type that controls which inefficient system dominates – i.e. till cover and roughness – is variable (Pelletier et al., 2016; Brubaker et al., 2013). Conceivably, while poroelastic drainage requires a porous ice sheet substrate, the cavities can form in any environment with bed protrusions which are less mobile than ice flow. A softbedded cavity has been seen at the base of a borehole in ice stream C (Carsey et al., 2002), and the theoretical basis for these cavities (Schoof, 2007) is motivated by drumlin formation (Fowler, 2009). However, cavities can only drain water once they grow enough to join and form a connected network (Rada and Schoof, 2018).
The contrast between kilometreorder or larger model scales and the metreorder or smaller process scales permits inefficient flow to be described as a continuum at the macro scale. On the macro scale, flux is related to water thickness and hydraulic gradient as follows (Flowers, 2015):
with flux Q, hydraulic conductivity k, and subglacial (basal) water thickness h_{wb}. The gradient of the hydraulic potential is given by
with subglacial water pressure P_{water}, density of freshwater ρ_{w}, gravitational acceleration g, and basal topographic elevation z_{b}. The exponents in Eq. (1) set laminar or turbulent flow. α=1 and β=2 give Darcy's law for laminar flow through porous media (Darcy, 1856; Muskat, 1934). $\mathit{\alpha}=\mathrm{5}/\mathrm{4}$ and $\mathit{\beta}=\mathrm{3}/\mathrm{2}$ give the Darcy–Weisbach relation for turbulent flow through conduits (Clarke, 1996; Weisbach, 1855). Equations (2) and (1) are combined with a water pressure closure relationship given by the underlying physical system to get the formulations in Sect. 2.1.1 and 2.1.2.
Water sheet thickness is a continuum property used to describe the average amount of water in a grid cell. Changes in water thickness are given by the fluxes and the aggregate of sources and sinks, m, in the water transport Eq. (3):
where Q is the subglacial water flux, and m is the aggregate of sources and sinks.
2.1.1 Poroelastic system
Pressurized subglacial water flows through the pore space of a layer between ice and bedrock, conceptualized as the interstitial space between till grains. As water pressure increases, the permeability of the porous medium rises. Water pressure is related to subglacial water thickness by a nonlinear function using porespace saturation (Eq. 4). This poroelastic drainage formulation is laid out in Flowers (2000). The flow law is Darcy's law describing laminar flux as a function of hydraulic gradient and subglacial water thickness. The pressure closure is an empirical relationship between the water column height in the elastic porespace and subglacial water pressure.
The Darcy flow law is Eq. (1) with α=1 and β=2. Water pressure in the elastic pore space is set by Eq. (4) (Flowers, 2000):
where P_{ice} is the pressure due to the weight of overbearing ice and h_{c} is the water thickness scalar interpreted as the thickness of the porespace accommodating water.
2.1.2 Linkedcavity system
As ice flows over protrusions in the bed, cavities open in the lee side. The faster ice flows and the higher the protrusion, the greater the opening rate. The weight of the overbearing ice acts to close the void through viscous creep. The tradeoff between these two rates determines the net cavity size change rate. These cavities are linked through smaller connections and form a drainage network whose throughput is controlled by orifice size and system tortuosity. As water flows more quickly in the drainage network, wall melting due to frictional heating at the ice–water interface further opens cavities and the interconnecting orifices, forming a more efficient system. The Darcy–Weisbach flow law for turbulent flux depends on the hydraulic gradient and subglacial water thickness. The pressure closure is based on cavity opening and closing velocities and mass balance. The Darcy–Weisbach flow law is Eq. (1) with $\mathit{\alpha}=\mathrm{5}/\mathrm{4}$ and $\mathit{\beta}=\mathrm{3}/\mathrm{2}$ with ${P}_{\mathrm{water}}={P}_{\mathrm{ice}}{N}_{\mathrm{eff}}$. For the linkedcavity system, k=k_{lc} in Eq. (1) aggregates quantities such as tortuosity, hydraulic gradient across the orifice, and cavity density. Completing the set of equations, the effective pressure, N_{eff}, is given by the opening–closing relationship for the cavity crosssectional area with respect to time in Eq. (5). This has three parts:

wall melting term ($\propto \mathit{Q}\cdot \mathit{\psi}$),

opening from sliding over bed protrusions (∝u_{b}h_{r}),

closing due to overburden pressure (creep) ($\propto {N}_{\mathrm{eff}}^{n}S$).
where S is the cavity size, c_{1} and c_{2} are constants, Q is flux, u_{b} is basal sliding velocity, h_{r} is bed protrusion height, and n is the exponent from Glen's flow law (Schoof, 2010; Werder et al., 2013). These three terms act to increase or decrease cavity area.
2.2 Efficient flow
In the efficient drainage regime, flux and water pressure are inversely related. Flux in the efficient system occurs in subglacial tunnels incised into overbearing ice (Röthlisberger, 1972), down into subglacial sediments (Walder and Fowler, 1994), or into hard bedrock (Alley, 1989). Channels eroded into bedrock remain in the same place through time, while those formed in ice or sediment can open, move, and close depending on overbearing ice and hydrologic conditions. The most commonly modelled efficient system is the Röthlisberger channel (R channel) carved up into the overbearing ice (de Fleurian et al., 2018). Dendritic subglacial tunnels open up into the ice from the base by wall melting due to frictional heat from the contact between ice and flowing water (Röthlisberger, 1972) – the faster the water, the larger the channel. Counter to the inefficient regime, water pressure and flux are inversely proportional (Schoof, 2010; Flowers, 2015). As water percolating through the inefficient system flows quickly enough to give significant wall melting, the system becomes unstable and quickly transitions to a channelized system (Schoof, 2010). Schoof (2010) showed that Eq. (5) bifurcates into the inefficient linkedcavity system and the efficient Rchannel system, the switch between the two controlled by flux in the subglacial system. At high fluxes, frictional melting of the tunnel ice wall from fastflowing water becomes a runaway effect opening an R channel into the ice. Canals likely open due to high flux as well in the subglacial system, where energetic water mobilizes sediment along its path (Alley, 1992; Walder and Fowler, 1994).
The conceptual basis for the efficient flow model herein is the R channel, which evolves out of the inefficient system based on high fluxes.
The model used here is a fully coupled system of hybrid SIA/SSA (shallow ice approximation, shallow shelf approximation) ice physics (Pollard and DeConto, 2012) and 3D ice thermodynamics and 1D bed thermodynamics (Tarasov and Peltier, 2007). The climate forcing imposes a background surface temperature trend and elevation dependencies for temperature and precipitation. The subglacial hydrology model includes a choice of linkedcavity, poroelastic, or leakybucket inefficient drainage, of which the linkedcavity and poroelastic drainage can be coupled to the efficient drainage tunnel solver. The transition from frozen to temperate is smoothed to more realistically capture the transition to sliding (as in model (d) of Calov et al., 2010) following the work of Hank et al. (2023). A more detailed description of the GSM is forthcoming (Tarasov et al., 2023).
3.1 Subglacial hydrology model
The subglacial hydrology model – BrAHMs2.0 – is an extensive update to version 1.0 (Kavanagh and Tarasov, 2018). The update includes the addition of linkedcavity and leakybucket systems, an updated generalized grid, modified convergence criteria, a modified flux limiter, and code restructuring. This model uses a finitevolume discretization with a staggered Arakawa C grid (fluxes at interfaces; Arakawa and Lamb, 1977). In the case of the 2D masstransporting hydrology setups (poroelastic and linked cavity), we implement the generalized flux calculation in Eq. (1) with a choice of either the pressuredetermining closure of Flowers (2000) or a modified version of Schoof (2010) as in that of Werder et al. (2013) and Bueler and van Pelt (2015), whereby the cavity opening rate is proportional to the difference in bed roughness and subglacial water sheet thickness. Schoof (2010) shows that the wall melting term in Eq. (5) is unimportant until a critical value is reached and the runaway effect opens tunnels (see assumption 1 below). As such, the wall melting term is assumed to be 0 until tunnelling is triggered.
In this model the cavities are described as a continuum: height of a cavity averaged over protrusion spacing (l_{r}) is given as ${h}_{\mathrm{cav}}=\frac{S}{{l}_{\mathrm{r}}}$:
The opening term is modified to drop as average cavity thickness rises over the bed protrusion u_{b}(h_{r}−h_{cav}) as in (e.g.) Werder et al. (2013), and cavities are assumed filled by subglacial water (h_{cav}=h_{wb}; see assumption 3). This leads to the relationship for water pressure evolution:
where ϕ_{eng} is the englacial porosity and m_{t} Eq. (8) is derived in Appendix E1 following Bueler (2014) and is similar to that used in Werder et al. (2013), Hewitt (2013), and Bueler and van Pelt (2015).
Bueler and van Pelt (e.g. 2015)(e.g. Flowers, 2000)(e.g. Flowers, 2000)While in the linkedcavity model the hydraulic conductivity is a single parameter, in the poroelastic model Flowers (2000) uses a meltwaterthicknessdependent arctan function for hydraulic conductivity to capture a transition from low to high permeability during the expansion of the pore space:
where k is the poroelastic hydraulic conductivity, k_{max} is the maximum conductivity, ${k}_{\mathrm{min}}={k}_{\mathrm{max}}/{k}_{\mathrm{ratio}}$ is the minimum conductivity, h_{c} is the critical water thickness in the pore space (h_{c} in Table 1), and k_{a} and k_{b} control the transition between the maximum and minimum conductivity.
Numerically, hydraulic conductivity in both the poroelastic and linkedcavity formulations is defined at the cell centres and is a function of cell temperature relative to the pressure melt point (T_{bp}). To account for the transition from fully cold based (frozen) to fully warm based (thawed), the bed is assumed to be fully frozen below T_{froz} and the hydraulic conductivity is given the following dependence on basal temperature relative to the pressure melt point (T_{bp}):
where k_{f} is the hydraulic conductivity of frozen till (effectively zero). As the flux should be a function of the potential difference across the interface, the harmonic mean of the adjacent cellcentred conductivities gives the most appropriate interface conductivity (Patankar, 1980).
In order to assess the importance of transport vs. pressure determination in surging, we implement a nonmassconserving zerothorder leakybucket scheme: a constant drainage rate (s_{drain}) counters the melt rate (s_{melt}) to give basal water thickness in that cell, following Gandy et al. (Eq. 3, 2019):
The leakybucket scheme uses the empirical pressuredetermining closure of Flowers (2000) shown in Eq. (4), with basal water thickness limited between 0 and the critical thickness of the pressure closure (h_{c} in Eq. 4).
Fully modelling the process of efficient drainage of water through the channel system would require very short time steps due to Courant–Friedrichs–Lewy (CFL) (Courant et al., 1928) restrictions and consequently prohibitively long run times. Given the disparity in timescale between efficient drainage (subannual) and the dynamical behaviour examined here (centennial to millennialscale surging), it is unlikely that dynamically versus diagnostically modelling the efficient drainage will have an effect on the longer timescale surging, though as the model is nonlinear there is potential for propagation across timescales. As such, an alternate scheme is used under the assumption that drainage happens far quicker than in the inefficient system (assumption 2), which should especially hold for ice sheet modelling contexts. If flux at a cell face exceeds the bifurcation threshold or “critical discharge” of Schoof (2010), water is routed down the background hydraulic gradient (i.e. from topography and ice sheet overburden), filling in potential lows along the way until routed water is depleted or exits the ice sheet. This subglacial meltwater routing scheme is a slight modification of the downslope surface meltwater routing scheme of Tarasov and Peltier (2006) (i.e. with a modified hydraulic gradient). This routing scheme is further discussed in Kavanagh and Tarasov (2018).
For details on the numerical solver used here, readers are invited to read Appendix B. Assumptions used in the design of this model are examined in Appendix C. The verification of the model implementation presented in Appendix D shows that the model

gives symmetric solutions given symmetric boundary conditions,

converges under increasing spatial and temporal resolution at a rate commensurate with the discretization schemes,

conserves mass,

gives similar solutions to another model using similar physics (GlaDS, Werder et al., 2013).
3.2 Basal drag coupling
The basal velocity is from either a hard or softbed sliding rule. For the hard bed the basal sliding rule is a fourthpower Weertman sliding law (Cuffey and Paterson, 2010):
where c_{hard} is a parameterized sliding coefficient which includes a parameterization for basal roughness, ${F}_{{N}_{\mathrm{eff}}}$ is the effective pressure normalization factor, N_{eff} is the effective pressure given by the subglacial hydrology model, and τ_{b} is the basal drag. The basal velocity for softbedded sliding is similarly a Weertmantype sliding law with integer exponent values between 1 and 7.
with separately parameterized softsliding coefficient c_{soft} (which also includes a parameterization for basal roughness).
Using a simple setup without externally driven variability from topography, a complex land–sea mask, and an unsteady climate, system behaviour is due to the initial transient response and internal feedbacks. Our Laurentide Ice Sheet square (LISsq) setup includes broad features of the North American bed (Fig. 1) and computationally cheap firstorder diagnostic climate imposed as a steady forcing with ice sheet thickness feedbacks. The simple climate allows a free southern margin determined by the background temperature and feedbacks, giving a dynamically determined ice sheet geometry at 50 km horizontal resolution. Next we present the design choices of this setup in three categories: bed, climate, and glacial systems.
4.1 Bed
LISsq aims to probe the effect of largescale hard to softbed transitions characteristic of North America. This simplified setup allows separating out the internal feedbacks from the externally forced elements (e.g. variability from real topography and a land–sea mask and unsteady, spatially varying climate). The shorter run times of this setup also allow larger ensembles, giving a better probe of the parameter space. The simplicity helps with model verification as any variability in the model stems purely from the encoded physical processes.
The majority of the inferred Late Pleistocene Laurentide substrate has been hard bedded (Clark et al., 2006), with unconsolidated sediment cover at the south and in the Hudson Bay and Strait. The Heinrich Event Intercomparison (HEINO) experiments were conducted over similar lengthscale hard beds with the same softbedded Hudson Bay and Strait at the centre of the hard bed (Calov et al., 2010). HEINO differed in that it included a circular continental configuration bounded by a highly ablating ocean – the ice sheet geometry was largely set. Here we wish to examine surge behaviour for a variety of ice sheet geometries within the roughly approximate range of the Laurentide length scales and bed. As such, a rectangular bed geometry is set with the boundary of the softbedded south at a constant latitude and an equilibrium line which is free to evolve with a changing ice sheet.
4.2 Climate
The LISsq climate prescribes a linear background temperature trend with lapse rate feedback. The annually averaged surface temperature, T_{surf}, is
where T_{north} is the ground level temperature at the northern end (^{∘}C), T_{grad} is the latitudinal warming rate (^{∘}C km^{−1}), L is the slope temperature lapse rate (^{∘}C km^{−1}), and H is ice sheet thickness (m; recall that the bed is at constant elevation and glacial isostatic adjustment is not included). The brackets, $\phantom{\rule{0.25em}{0ex}}\textcolor{red}{}$, denote the maximum.
These temperatures are then used together with a positive degree day scheme (PDD) to simulate net seasonal contribution to accumulation and ablation for an annual average temperature. The positive degree day sum assumes a 100 d melt season length with temperatures 10 ^{∘}C warmer than the annual mean, T_{surf}, and a melt coefficient in m per PDD. Ablation is then
where b_{melt} is ablation (m yr^{−1}) and T_{surf} is surface temperature (^{∘}C). Accumulation incorporates the thermodynamic effect on atmospheric moisture content using the August–Roche–Magnus approximation for the Clausius–Clapeyron relationship (Lawrence, 2005) with parameter ranges adjusted for undersaturated air. Accumulation, b_{accum}, is 0 where T_{surf}≥0 ^{∘}C:
where the reference precipitation rate, p_{ref}, and precipitation preexponential factor, h_{pre}, are ensemble parameters.
4.3 Glacial systems
We use a subset of the fullfeatured GSM for this setup. Here we omit glacial isostatic adjustment, surface meltwater drainage, sediment transport and production, and ice shelves with groundingline flux and calving model. This is in order to clearly show the effect of hydrology feedbacks on ice flow and ice thermomechanics.
4.4 Parameter range estimation
In this section we justify chosen parameter ranges based on physical and heuristic arguments and current understanding in the literature.
4.5 Hydraulic conductivity parameterization
The range of values appropriate for hydraulic conductivity varies according to whether the drainage system is assumed to be poroelastic or linked cavity or whether the flux is assumed to be laminar (Darcian) or turbulent (Darcy–Weisbach). Hydraulic conductivity is not truly known at the continuumlevel macroscale. Here we use a range based on bounding subglacial hydrologic flow velocities, typical hydraulic gradients, and subglacial water thicknesses.
The velocity of water flow in the subglacial channel endmember imposes an upper bound on the linkedcavity endmember flow velocity in the bifurcated channel–linkedcavity system (Schoof, 2010). Chandler et al. (2013) used dye tracing experiments at a landterminating west Greenland catchment to measure maximum velocities between moulin injection site and the margin. Their slowest first arrival time gave 1.00 m s^{−1} in the efficient drainage regime.
Fast ice velocities (e.g. ≈ 1 km yr^{−1}) give a loose lower bound on water flow speeds. Whereas the viscosities differ by many orders of magnitude (10^{14}Pa s for ice versus 10^{−3}Pa s for water; Cuffey and Paterson, 2010), the pressure gradient forces are less dissimilar. Assuming the hydraulic gradient is approximately equivalent to that imposed by ice sheet and bed topography (i.e. no contribution from basal water pressure) – around 1000 m per 56 km (Chandler et al., 2013) ice sheet surface gradient contribution and 500 m per 56 km for bed contribution (Morlighem et al., 2013) – gives a hydraulic gradient of ψ≈240 Pa m^{−1}. Assuming further ranges of 1 mm to 10 m of basal water thickness and Darcy–Weisbach flow speeds between $\mathrm{3}\times {\mathrm{10}}^{\mathrm{4}}$ and 1×10^{0} m s^{−1} gives a range of linkedcavity hydraulic conductivity (k_{cond}, Table 1) between $\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$ and $\mathrm{1}\times {\mathrm{10}}^{\mathrm{1}}$ m${}^{\mathrm{5}/\mathrm{4}}$ Pa${}^{\mathrm{1}/\mathrm{2}}$ s. To ensure complete bounding, we probe a wider range of $\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$ and 1×10^{1} m${}^{\mathrm{5}/\mathrm{4}}$ Pa${}^{\mathrm{1}/\mathrm{2}}$ s. This range encapsulates values suggested by Hager et al. (2022) and Werder et al. (2013). Flowers (2000) assessed the range of hydraulic conductivities to be k_{max}=1 m s^{−1} and ${k}_{\mathrm{min}}=\mathrm{10}\times {\mathrm{10}}^{\mathrm{7}}$ m s^{−1}. The hydraulic conductivity transitions from k_{max} to k_{min} according to Eq. (9).
4.6 Basal roughness
The height of bedrock protrusions relevant to subglacial cavity formation and its spatial variation lacks assessment in the literature and justified values are difficult to come by. The height of these protrusions, or terrain roughness, affects several basal processes in glaciated regions, including heat generation in basal ice, sliding, subglacial cavity opening, and bedrock quarrying. Length scales relevant to subglacial cavity formation have been estimated from chemical alteration of bedrock (the deposition of calcium carbonate precipitates) (Walder and Hallet, 1979). These cavity outlines form during slidingassociated regelation when water refreezes at the glacier substrate on the lee side of bedrock highs, precipitating dissolved carbonates. The deposits in this study indicate cavities 0.1–0.15 m high. Several studies then use a value in this range (e.g. Werder et al., 2013). Kingslake and Ng (2013) refers to Walder (1986) for this value, but Walder (1986) does not provide any justification for it in their Table 2 and does not refer explicitly to the earlier work of Walder and Hallet (1979).
In deglaciated areas with bed access, quantifying roughness at the ice sheet scale is a nonunique problem and measures abound. For example: standard deviation of elevation, power spectral density of elevation, and local bed slope. These are relative measures which do not identify the typical prominence of roughness features in a domain. What is needed for modelling linked cavities is the average height of bedrock protrusions relevant to the cavity scale (itself uncertain) at given wavelengths. How these heights vary spatially for previously glaciated regions has not been assessed. Identifying this as a gap in the current glaciological literature, we adopt similar scale values and probe a wide range in order to capture ice sheet sensitivity to the scale of cavityforming bump height. As stated above, Werder et al. (2013) and Kingslake and Ng (2013) both use h_{r}=0.1 with the latter referring to Walder (1986), who gives a range of 0.01 to 0.5 m for the relevant bump height. Iverson (2012) show cavities and quarrying are intrinsically linked. As such, the step size of quarried surfaces may indicate a scale for cavity growth. Anderson et al. (1982) mapped cavities forming along 1 m high steps at the base of Grinnell Glacier in Montana, United States. Following the same reasoning, the size of quarried boulders also gives an estimate of the upper bound for length scales. Though less common, 20 m boulders can be found (though if transported debris were comminuted in transit, the original size distribution would have been larger). As such, we use a range of ${h}_{\mathrm{r}}\in \left[\mathrm{0.01},\mathrm{20.0}\right]$ m and a range for the roughness wavelength as a function of roughness of ${l}_{\mathrm{r}}\in \left[\mathrm{1.0},\mathrm{20.0}\right]\times {h}_{\mathrm{r}}$.
4.7 Hydrology temperate transition
This parameter is used to interpolate between a conducting (at 0 ^{∘}C) and nonconducting (at a lowerbound temperature) hydrologic system with a logic similar to the temperature ramp reasoning. Thus, the range is based on Hank et al. (2023) and the lower bound of the interpolation is probed in the range of $\left[\mathrm{1.0},\mathrm{0.0}\right]$.
4.8 Tunnelswitching scalar
The flux threshold switch from inefficient to efficient drainage is given by the ratio of cavity opening due to sliding versus wall melting from viscous heating (Schoof, 2010):
where u_{b} is the velocity, ${h}_{\mathrm{r}}/{l}_{\mathrm{r}}$ the basal roughness ratio, c_{1} a scalar, α the Darcy–Weisbach water thickness exponent, ψ the hydraulic gradient. Q_{scale} is a scale factor adjusting for subgrid uncertainty – smallscale fluctuations in flux may trigger a runaway tunnelling positive feedback affecting the larger scale.
4.9 Effective pressure normalization
This is the value used to normalize the effective pressure in the basal sliding velocity calculation and is set based on typical effective pressures. Effective pressures greater than this parameter values should slow sliding and less than that should hasten sliding. We set this range to [10 kPa,1 MPa] based on the typical effective pressure values seen in Fig. 10. The effective pressure and normalization (${F}_{{N}_{\mathrm{eff}}}$) is incorporated into the hard and soft basal sliding velocities in Eqs. (13) and (14).
4.10 Basal sliding parameters
The soft and hard sliding factors used in Eqs. (13) and (14) were set to wide bounds somewhat outside the recommended range for the GSM (Tarasov et al., 2023); the power for softbedded sliding was kept within the typical range. These ranges were ${c}_{\mathrm{soft}}\in \left[\mathrm{0.01},\mathrm{4.0}\right]$ (set such that sliding speed is ≈3 km yr^{−1} for ${\left(\mathrm{30}\phantom{\rule{0.125em}{0ex}}\mathrm{kPa}\right)}^{{b}_{\mathrm{till}}}$ kPa of basal drag), ${c}_{\mathrm{hard}}\in \left[\mathrm{0.0},\mathrm{5.0}\right]$ (set such that sliding speed is ≈200 m yr^{−1} for 100 kPa of basal drag), and ${b}_{\mathrm{till}}\in \left[\mathrm{1},\mathrm{7}\right]$.
4.11 Climate parameters
A range of [5,10] ^{∘}C km^{−1} is used for slope lapse rate on the basis of PMIP2 Greenland model simulations in Erokhina et al. (2017). The range for T_{north} was obtained from the PMIP4 ensemble mean distribution of northern ($>\mathrm{75}{}^{\circ}$) latitude temperatures at the Last Glacial Maximum (LGM) in Kageyama et al. (2021) shown in Fig. 3. The precipitation parameter ranges in Eq. (17) were adjusted to bound the range of precipitation and temperatures below freezing in Kageyama et al. (2021), as shown in Fig. 3.
4.12 Ensemble design and parameter sensitivity
To understand the effect of hydrology, ensembles for different model configurations are compared (Table 2): linkedcavity (LC), poroelastic (PE), leakybucket (LB), and nohydrology (NH) ensembles – 18 816, 19 992, 15 288, and 11 760 runs in each ensemble respectively. Each ensemble varied the hydrology, ice sheet, and climate parameters simultaneously in order to capture parameter interactions and the number of runs was scaled with the number of parameters in each setup (15 in LC, 16 in PE, 12 in LB, and 9 in NH, shown in Table 1). The parameter space is sampled with the quasirandom lowdiscrepancy Saltelli extension of the Sobol sequence (Saltelli, 2002) as implemented in SALib (Herman and Usher, 2017) with secondorder terms enabled. Parameters are sampled with a loguniform distribution for parameter values which vary over orders of magnitude. Each run proceeded for 100 kyr with the first 50 kyr taken as spin up (from no ice, initial accumulation given by the background temperature from T_{north} and T_{grad}). Ensembles were run on a heterogenous Linux cluster with 24–32 Gb RAM and 8–24 Xeon or Opteron cores per node, clock speeds ranging from 2.4 to 2.7 GHz and a total of 652 cores. Run times averaged about 3 h.
Ice sheet geometries vary widely among runs for all model configurations. Maximum ice thickness ranges from 0 to ∼ 6000 m, while the maximum north–south extent ranges from 0 to 4500 km. Here we study surge behaviour at a scale similar to the Laurentide ice sheet by sieving (discarding runs outside the target metric ranges) the ensembles according to maximum ice sheet thickness and north–south extent. At the Last Glacial Maximum (LGM, 22 ka), the maximum ice thickness was plausibly around 4000 m (Tarasov et al., 2012). We use this estimate with a lower bound of 3000 m for the sieve in the main study and examined additional sieves with bounds [2500,3500] m and [3500,4500] m in Appendix A. LGM north–south extent was ≈ 4000 km, while the last margin to fully encircle the Hudson Bay and Strait (11.50 ka) extended ≈ 2500 km north to south (Dalton et al., 2020).
The importance of hydrology parameters for determining ice sheet geometry can be probed with sensitivity analysis. Local sensitivity analysis methods neglect interaction terms important for studying feedbacks in coupled models and so are not applicable here (Saltelli et al., 2008). Meanwhile variancebased (Sobol, 2001) methods require assumptions about the sampling structure of the underlying inputs. The trouble with coupled models is that they can be unstable; as such, there are incomplete runs which render sampling structure assumptions moot. There are other nonparametric sensitivity methods which do not require assumptions about the input sample distribution, but these require sample sizes even larger than those presented here in order to converge (e.g. Borgonovo, 2007; Pianosi and Wagener, 2015).
We develop a novel nonparametric method to measure sensitivity: we assess ice sheet geometry sensitivity to parameters by comparing the original uniform input parameter distribution with the parameter distribution corresponding to the sieved geometries (limiting the ensemble to those within geometric bounds). The nonparametric nature alleviates the need to make assumptions about the underlying parametric distribution class (e.g. variance is a normal distribution parameter). Using the impact of a sieve on parameter distribution to measure sensitivity means that assumptions about the sampling methodology are not required and that successive sieves can be applied to the ensembles to measure different aspects of model sensitivities. For example, in Sect. 5.1 we measure the sensitivity of surge frequency for those ensemble members which pass the geometry sieve by further sieving on surge frequency. Parameters which do not control the ice sheet geometry will have a similar distribution after selecting for that geometry range as the original input sample distribution. The more modified the distribution, the more sensitive the parameter. More precisely, each distribution is approximated with a kernel density function (KDF) normalized to unit area under the KDF. The sensitivity metric is then the integral of the absolute difference between the sieved and unsieved KDFs, i.e. the measure of how much the sieve modifies each parameter's KDF. For example, the maximum KDF difference would stem from a narrow spike on the sieved distribution, which would mean that parameter strongly controls the model output, e.g. the more limited range indicated for h_{pre} in Fig. 4a. We add a uniformly sampled dummy parameter not used by the model to set a threshold of accuracy of the sensitivity metric in each case. This dummy parameter has a very similar input and sieved distribution (with minor differences due to the essential random sampling from sieving), for example that for the LC geometry sieving in Fig. 4.
The sensitivity metrics for all parameters in Fig. 5 rise above the baseline significance level set by the dummy parameter in each ensemble. The temperature coefficient in the August–Roche–Magnus relation (h_{pre}), north–south temperature gradient and intercept (T_{grad} and T_{north}, Eq. 15), and lapse rate are the top four geometrycontrolling parameters in all cases except LB (though lapse rate is close for this ensemble as well). In the hydrology enabled setups, hydrology parameters rank in the top five.
The ranked parameter sensitivity for each model in Fig. 5 exhibits an inflection point in parametric sensitivity which we use to determine the number of controlling parameters. This inflection point is an approximate indication of the diminishing sensitivities in the model setup. As such, parameters to the right of this point are taken as sensitive and those to the left are considered insensitive and could be fixed for the purposes of geometry. Around twothirds of parameters fall on the righthand side of the inflection in each ensemble. For those hydrologybearing model configurations, half or more of the hydrology parameters lie in the sensitive zone. This shows that subglacial hydrology is even important at the scale of whole ice sheet geometry.
The most influential hydrology parameter in the LC setup (Table 2) is hydraulic conductivity, which controls the dynamic effective pressure, while in LB and PE the geometry is quite sensitive to the normalization of the effective pressure in basal drag (though PE is more sensitive to the tunnelling tendency, Q_{scale}). In the LC case, the ice sheet geometry is most sensitive to those parameters which control the dynamics of effective pressure themselves (k_{cond} and h_{r}). In the PE case, the parameters controlling the transition to efficient drainage (Q_{scale}) and effective pressure normalization are most important hydrology parameters. These parameters are both diagnostic controls on subglacial water balance and sliding velocity. Similar to PE, the most important subglacial hydrology parameter for LB is the effective pressure normalization.
4.13 Surge metric definition
The two most obvious measures of internal oscillation are amplitude and period. This highly nonlinear system does not exhibit sinusoidal behaviour, but we can pick surge metrics which approximate these measures. To this end, each surge type was evaluated in two ways – the number of surge events (an indication of periodicity, the number of red dots in Fig. 6) and strength (or speed increase, the height of vertical purple bars in Fig. 6) of surge events (i.e. amplitude).
The background sliding speed of the actual Hudson Strait Ice Stream (HSIS) in the nonsurging state is unknown. While this study does not aim to replicate the actual Hudson Strait (HS), we are studying the behaviour of an ice stream and sheet with similar dimensions to the HS and Laurentide ice sheets. As such, labelling and measuring the strength of a surge event needs to be agnostic of quiescentphase conditions between events. Ice stream acceleration at scales comparable to the HS has not been observed in the modern period. Though significantly smaller than the HSIS and its catchment, the Vavilov ice cap did accelerate from 12 to 75 m yr^{−1} between 1998 and 2011 CE (Willis et al., 2018). Satellite observations of the North East Greenland Ice Stream (NEGIS) combine with modelling to show acceleration greater than 1 m yr^{−2} in places between 1985 to 2018 CE (Grinsted et al., 2022). Therefore we define a surge event in this setup as a large increase in spatially averaged HS basal sliding speed (>1000 m yr^{−1} over a given 25–100year acceleration period) over the background, quiescentphase speed. Velocity can also change during a surge as portions of ice within the HS accelerate over others. Ice stream shear margins can be regions of the fastest velocity changes and ice stream geometry can change over time (Grinsted et al., 2022). As such, we do not define adjacent shortlived changes in velocity as separate surge events.
A typical run with surge events which passes the ${\mathcal{S}}_{\mathrm{geom}}^{\mathrm{0}}$ sieve is shown in Fig. 6. In order to label surge events (red dots in Fig. 6), we use peak prominence (Virtanen et al., 2020) – drawn from the concept of topographic prominence (height of local maxima above adjacent local minima) – to estimate surge events from the basal velocity time series (1year sample rate) for each run. This allowed surged metrics to be agnostic of any background value. In order to minimize spurious peaks picked on variations in velocity during a single event, a 401year median filter was applied. This means that abrupt velocity changes lasting ∼ 200 years or less will not get picked as events. This is less than the lower bound on HSIS surge duration inferred from ice rafted debris (IRD) by Dowdeswell et al. (1995), who estimate that those surges most likely lasted between 250 and 1250 years on the basis of Heinrich events interpreted in 50 North Atlantic drill cores. A comprehensive review of Heinrich events and IRD age intervals available in the literature by Hemming (Table 3, 2004) infers a mean duration of 495 years, where the lowest estimate is 208 years. The duration for these modelled surge events is calculated as full width at 80 % maximum prominence (height above adjacent local minima).
The sieves used for sensitivity analysis are shown in Table 3. Sieving the data by ice sheet geometry (${\mathcal{S}}_{\mathrm{geom}}^{\mathrm{0}}$) cuts the ensemble size to $\approx \mathrm{1}/\mathrm{6}$ to $\mathrm{1}/\mathrm{9}$: the poroelastic (PE) (Table 2) system has $\mathrm{3154}/\mathrm{19}\phantom{\rule{0.125em}{0ex}}\mathrm{992}$ (15.8 %), the linkedcavity (LC) system $\mathrm{3566}/\mathrm{18}\phantom{\rule{0.125em}{0ex}}\mathrm{816}$ (19.0 %), the leakybucket (LB) system $\mathrm{2721}/\mathrm{15}\phantom{\rule{0.125em}{0ex}}\mathrm{288}$ (17.8 %), and the nohydrology (NH) ensemble $\mathrm{1382}/\mathrm{11}\phantom{\rule{0.125em}{0ex}}\mathrm{760}$ (11.8 %) runs. The histograms in Fig. 7 show the frequency of surge events and strength of speedup of those events in the last 50 kyr of each simulation. The lower bound of HSIS surge frequency inferred from the Heinrich Event record (Hemming, 2004; Naafs et al., 2013) is three in 50 kyr. The rate of runs with 3 to 12 surge events in the sieved results is $\mathrm{423}/\mathrm{3154}$ (13.4 %) for PE, $\mathrm{504}/\mathrm{3566}$ (14.1 %) for LC, $\mathrm{836}/\mathrm{2721}$ (30.7 %) for LB, and $\mathrm{75}/\mathrm{1382}$ (5.4 %) for NH. The distribution of the frequency of surge events stemming from each hydrology setup is not significantly different from the others (though LB does have more surges in the 4–7$/\mathrm{50}$ kyr frequency range) nor is the magnitude of the ice stream speedup. The nohydrology case, however, does differ from those three: the rate of runs with surge events is significantly lower and the frequency and strength of events per run are also lower.
The duration of HS surge events highlights a difference between the three hydrologies: the linkedcavity system yields longer duration events and the trend in duration with increasing event frequency diverges between the linked cavity and the other two hydrology systems. As the duration of surge events necessarily depends on the frequency of those events (having more events in a time period decreases the maximum possible duration of those events), we examine surge duration as a function of the number of events (as shown by the horizontal purple lines in Fig. 6). In Fig. 8 we extract the median surge duration by selecting runs with a given number of events and comparing the duration–frequency trends between the four setups. Frequency levels with 10 or fewer runs passing the sieve are omitted as trends degrade around this level of membership.
As the frequency of surge events in each run increases, the median duration of surges in those runs stays largely flat, perhaps decreasing slightly for both the poroelastic and leakybucket hydrologies. This is not so for the linkedcavity system: the duration of surges increases up to seven surge levels, where it roughly doubles that of the poroelastic and leakybucket hydrologies. This relationship is stronger still when selecting thinner ice sheets with a mean maximum thickness between [2500,3500] m as shown in Fig. A2. In this geometry range, the surge duration decreases at first, reaching a minimum at three surges before steadily increasing in duration until it more than doubles the leakybucket surge duration (PE run counts are below the significance threshold). For thicker geometries no differences between the three hydrologies are apparent (Fig. A4).
5.1 Sensitivity of surge frequency
Applying a sieve on surge frequency in addition to the geometry sieve (𝒮_{surge} in Table 3) highlights the system sensitivity to subglacial hydrology. Figure 8 shows the result from selecting those runs with between 3 and 12 surge events, which is consistent with the minimum number of Hudson Strait surges inferred from the Heinrich Event record and the maximum number of events in the figure. The sensitivity ranking in Fig. 9 is insensitive to whether the sieve upper bound is 8 or 40 events, likely due to the fact that most runs have eight or fewer surge events. For all of the hydrology ensembles, the effective pressure normalization exerts the most control on surge frequency (Fig. 9). In the case of the PE and LB ensembles, hydrology parameters give the first and thirdhighest sensitivities. ${F}_{{N}_{\mathrm{eff}}}$ is highest in both cases, and k_{max} is third for PE and h_{c} is third for LB. For LC, the next hydrology parameters do not appear until seventh and eighth place. This may be due to the dual role ${F}_{{N}_{\mathrm{eff}}}$ plays in the linkedcavity system: it exerts influence on the sliding velocity, which in turn controls the cavity opening rate which is proportional to effective pressure. In the NH case, softbed sliding parameters c_{soft} (softbed sliding coefficient) and POWb_{till} (softbed sliding law power) are the most important for surge frequency. POWbtill is also the secondmost important parameter in both the LC and PE cases.
5.2 Relationship between effective pressure and sliding velocity
In Fig. 10, all warmbased points in the ensemble (across the parameter–space–time domain) for each hydrology configuration were crossplotted in log (N_{eff})−log (u_{b}) space in order to check for any systematic differences in velocity between the four configurations. If the configurations with subglacial hydrology had increased basal velocities at the ensemble level relative to the nohydrology case, then the conclusion that subglacial hydrology produces a wider distribution of surge characteristics would provide much less confidence.
The increased incidence of surge behaviour in the hydrology cases is not due to increased sliding – the nohydrology ensemble exhibits higher basal velocities than the three hydrology ensembles in Fig. 10. This check allowed for an interesting overall comparison between the hydrology configurations. The three hydrology formulations do exhibit differences in $\mathrm{log}\left({N}_{\mathrm{eff}}\right)\mathrm{log}\left({u}_{\mathrm{b}}\right)$ space (Fig. 10). Linkedcavity hydrology produces a bimodal clustering at lower velocities and higher effective pressures and higher velocities and lower effective pressures. This is a stark difference from the other two hydrologies, whose effective pressure distribution simply decays toward lower values. This bifurcation of the effective pressures from a linkedcavity system shows that it can sustain lower effective pressures than its poroelastic and leakybucket counterparts.
As we show above through sensitivity analysis and ensemble comparison of surge frequency and amplitude, subglacial hydrology is an important process that contributes to the feedbacks which govern Hudson Straitscale ice stream surging. While the process as a whole matters, the details matter less so – though it does depend on the aspects of ice stream surging under scrutiny. Across the three hydrology setups, the same range of HS basal velocity increase occurs: the magnitude of ice stream speedup is not dependent on the form of the subglacial hydrologic system, and the three models can attain the same velocities within parametric uncertainty. This means that for model experiments looking to realistically capture ice stream surges, a leakybucket hydrology (the computationally cheapest of the three) is sufficient. Additionally, the range of frequency of HS surge occurrences is quite similar across the three hydrologies. However, the nohydrology case falls short of covering the range inferred for actual Heinrich events attributed to HS surging (Naafs et al., 2013). This indicates that the inclusion of some form of coupled subglacial hydrology is important for modelling largescale surge periodicities on geologic timescales. Once again, however, the exact form of the subglacial hydrology does not matter for the periodicity of the surge onsets.
Plausibly, one might expect that simply increasing the sliding coefficient in the nohydrology case would generate more surges. We therefore compared the basal velocity distributions between the configurations (Sect. 5.2). The velocity distributions (minimum, mode, maximum) in the hydrology ensembles were slower relative to the nohydrology configuration in Fig. 10. The range of softbed sliding coefficient covered in each ensemble approaches the bounds of plausibility: ${c}_{\mathrm{soft}}\in \left[\mathrm{0.01},\mathrm{4.0}\right]$, where c_{soft} is scaled to give a 3 km yr^{−1} sliding velocity for 30 kPa basal drag. HS surge behaviour cannot be captured by increasing the sliding coefficient.
Increasing the lapse rate to nonphysical bounds can increase the incidence of HS surge events in the nohydrology case. In the main experiments, the lapse rate is limited to the range [5,10] ^{∘}C km^{−1}. However, increasing the lapse rates to [10,20] ^{∘}C km^{−1} increases the rate of surge events. This is because decreasing the surface temperature of ice in the Hudson Bay and Strait both increases the vertical heat diffusion and decreases the temperature of ice advected to the base during a surge event. This enables a stronger thermomechanical surge termination mechanism.
Surge initiation at peak velocity for Hudson Straitscale ice streams as soon as the pressure melt point is reached is physically implausible. Basal velocity increases after ice becomes warm based and the effective pressure decreases. Inclusion of subglacial hydrology in the coupled system accomplishes this. The accommodation of increasing amounts of basal meltwater and pressurization (in the case that channelization does not occur) acts as a system inductance, and the ice stream continues to speed up after becoming warm based. This inductance does not require the lateral transport of meltwater, only the balance of meltwater and a pressure closure dependence on subglacial water thickness.
Though periodicity and strength of surges are similar between the three hydrologybearing experiments, an interesting distinction occurs when examining the duration of events at varied frequencies. The stabilizing negative feedback of increasing effective pressure at higher basal velocities in the linkedcavity pressure closure gives surge durations longer (up to double the time, depending on frequency) than those of the diagnostic pressure closure of the poroelastic and leakybucket hydrologies. This feedback also results in a bimodal effective pressure distribution (i.e. Fig. 10). When studying ice stream surge behaviour, any of the hydrologies may give the same surge response in terms of frequency and strength of surges. If the study requires a more granular understanding of how long the surge was active, for example when studying the surge timing of multiple ice streams in a catchment (e.g. Payne, 1998; Anandakrishnan and Alley, 1997) or the lifespan of palaeoice streams, our results suggest that accounting for the appropriate hydrology system is required.
It is not possible to simulate fully dynamic channelized drainage at the scale studied here; the CFL criterion (Courant et al., 1928) would impose prohibitively long run times for our context. For illustration, Chandler et al. (2013) measured a lowerbound water velocity of 1 m s^{−1} in the channel system. At this speed with our (coarse) resolution of 50 km, a time step of 0.00158 model years is required. The downgradient routing scheme representing the efficient drainage system is not restricted by CFL, and so the time step depends only on the inefficient system, which is typically in the range of 0.5 to 0.25 model years. A dynamic model of the efficient system would increase BrAHMs run time anywhere from a 150 to a > 300fold rendering simulation of millennialscale variability infeasible.
Dynamical changes in flow through the efficient system occur on diurnal to seasonal timescales, while the timescales of system features examined here are centennial to millennial. This separation in scale by several orders of magnitude makes it unlikely that dynamical changes in the efficient system (requiring a dynamic model) would be a significant control on the longerscale variability. However, in a nonlinear system, such a control across scales cannot be fully ruled out.
While the treatment of efficient drainage in the model makes it more difficult to closely examine its role in the overall surging system, it is possible to evaluate its role at the ensemble level. At this level it is apparent that efficient drainage does not play a significant role in surging at this scale. Three points bring this to light:

The impact on effective pressure from the downgradient tunnel routing scheme is exaggerated as its modification of the basal water distribution is immediate instead of smooth. This modifies the effective pressure field in both the poroelastic and linkedcavity systems through the ${h}_{\mathrm{wb}}/{h}_{\mathrm{c}}$ term in Eq. (4) and the cavity opening rate term in Eq. (8) respectively.

The tunnelswitching criterion is well established from a physical mechanistic standpoint (Schoof, 2010). Here we have included a tunnelswitching subgrid uncertainty factor (Sect. 4.8). Sensitivity analysis shows that this parameter, which varies by 3 orders of magnitude, plays little role in surge generation (Fig. 9) where the parameter ranks last in the linkedcavity model and sixth from last in the poroelastic model. The role of efficient drainage may be greater in the poroelastic system than the linkedcavity system, suggested by its higher ranking in both the surge sensitivity (Fig. 9) and geometry sensitivity (Fig. 5). This difference in sensitivity may result from the fact that, whereas in the linkedcavity system the rate of change in effective pressure is proportional to the basal water thickness, in the poroelastic system the effective pressure is directly proportional to the basal water thickness.

Though the efficient system is not included in the leakybucket configuration, there is little difference in the range of surge frequency and amplitude with respect to the other two systems. The distinction in surge duration stems from the dynamic pressure closure of the linkedcavity system and its direct twoway feedback with sliding velocity.
The model presented herein passes multiple verification tests and as such is dependable for comparing the effects of structural choices of subglacial hydrology. The sensitivity analysis and ensemble comparison shows that subglacial hydrology is an important control on both ice sheet geometry and on the surging of major ice streams similar in scale to the Hudson Strait Ice Stream. However, depending on the characteristics of interest, the process details do not matter within current parametric uncertainties. The details do not matter for surge periodicity nor strength, but when studying the surge duration the hydrologic details are essential.
Surge behaviours can be produced in the absence of modelling a subglacial hydrology system, but this requires unrealistic assumptions: pushing lapse rates to unrealistic ranges or implementing an unphysical sudden thaw in a large grid cell when the temperature reaches the pressure melt point. Subglacial hydrology provides a system inductance necessary for realistic ice speedup at the temperate transition. The critical components are the accommodation of meltwater and a meltwater pressure closure, not the massconserving meltwater transport itself.
BrAHMs2.0 solves the conservative transport equation for the distribution of subglacial water (Eq. 3) and the effective pressure evolution equation (Eq. 8) using combined explicit and semiimplicit methods. Time integration is done first with Heun's method for the initial time step followed by a leapfrog trapezoidal predictor–corrector method (Kavanagh and Tarasov, 2018). To avoid time splitting, Heun's method is called after every 10 leapfrog steps (varying the number of leapfrog steps had little effect on the solution in tests).
The verification of this scheme and its implementation is presented in Appendix D with a fourpronged approach. The model is shown to give spatially symmetric solutions given symmetric boundary conditions. The convergence is examined for the spatial and temporal discretizations and found to approximately match the expected rate for each scheme: the firstorder upstream finitevolume implementation spatially converges at a linear rate, while the secondorder leapfrog trapezoidal implementation temporally converges nearly quadratically. The associated partial differential equations, however, are nonlinear, coupled, and likely to have nonlocal responses. As such, assessing the expected convergence rate of this system is not straightforward (Tadmor, 2012). In Appendix D the model is shown to also conserve mass and match the solution of another numerical model with similar physics (Werder et al., 2013).
The physics of the linkedcavity system are highly nonlinear. As such, a set of simplifying assumptions is required to make numerical modelling of this framework feasible.

Wall melting is not a control on cavity size until tunnels are opened. Drainage systems switch from inefficient to efficient for a given value of flux. Schoof (2010) showed that the evolution of the subglacial drainage system (described in Eq. 5) gives a bifurcation between cavity style and tunnel style drainage networks. Given effective pressure, the cavity opening speed is dominated by basal sliding below a certain flux and by runaway wall melting above it.

At timescales of continentalscale ice sheets, tunnels drain water instantaneously. The timescale of drainage through subglacial tunnels is less than a single melt season, much shorter than the centennial to millennialscale changes this model is applied to. This assumption alleviates CFL violations from fast tunnel flux, which would render modelling on the long timescales of glacial cycles infeasible.

Cavities are filled with water. Consider the timescale for the closure of a recently drained cavity given various combinations of ice sheet overburden (thickness, m) and sliding velocities. This timescale for closure (from Eq. 6) is given by
$$\begin{array}{}\text{(C1)}& T={\displaystyle \frac{S}{{u}_{\mathrm{b}}{h}_{\mathrm{r}}{c}_{\mathrm{2}}{N}^{\mathrm{3}}S}}.\end{array}$$The range of timescales, assuming speed in the range of 1–1000 m yr^{−1} and ice overburden thickness greater than 200 m, is shown in Fig. C1, where the maximum time for closure is around 2 weeks, less than the minimum time step of 0.125 years in the hydrology model.
Oreskes et al. (1994) describe model verification in general as the task of demonstrating model veracity, correctly asserting that no model can ever be proven – only disproven. However, this problem is not unique to computational model testing, this is a more philosophical epistemological problem. As Sornette et al. (2007) identifies, we do not prove models, we simply build our trust in them through a series of failed attempts to disprove them. In this section, we document performance on some simple tests which every model should pass before any amount of confidence can be conferred.
Following others (e.g. Sornette et al., 2007), we take model verification to be more pedestrian than validation: a test that the computational model actually solves the model equations as intended – or, as Roache (1997) defines it, “solving the equations right.” Meanwhile, we take validation as the converse of Roache (1997), i.e. “solving the right equations.” Validationwise, in this work we are showing not that the right equations were solved but that this seems to be of low consequence.
The results presented in this section were reached in an effort to expose errors in the models, the lowesthanging fruit in gaining confidence in the model solutions. The verification strategy in this section is to satisfy the following:

model solutions are symmetric given symmetric input,

model solutions converge under increasing spatial and temporal resolution,

mass is conserved, and

models using similar physics should have similar solutions.
Using simplified setups, expected behaviours are straightforward and in some cases may be calculated by hand (though hand calculations are not shown here). By using a progression of most simple to increasingly complex model setups for testing, model behaviour can be verified against expected behaviour and shown capable of simulating increasingly realistic environments. Here we demonstrate that the model correctly solves the equations. A progression of forcings and couplings were used, of which the transient, twoway coupled solutions from the least stable parameters (while still physical) are shown.
Parabolic surface topographies haven been used to approximate nonstreaming ice sheet topographies (e.g. Mathews, 1974). The Subglacial Hydrology Model Intercomparison Project (SHMIP) (de Fleurian et al., 2018) uses such an ice sheet surface (depicted in Fig. D3) and provides solutions to models using similar physics as the model herein. This therefore provides an appropriate test bed. This SQRT_TOPO surface is given by
and the flat base z_{b}=0.
Testing of the linkedcavity system with a Darcy–Weisbach flux model configuration (Eqs. 8 and 1) is presented here as this is the most nonlinear form and a new addition to the model.
The basal sliding velocity is determined by the effective pressure from Eq. (8):
where ${k}_{\mathrm{slide}}=\mathrm{5.0}\times {\mathrm{10}}^{\mathrm{1}}$ m s^{−1} is a scaling constant, an effective pressure regularization 10 Pa is applied for numerical stability, and basal shear stress (τ_{b}) is calculated from the constant driving stress (τ_{d}):
D1 Symmetry test
Spatial symmetry at each spatial resolution was calculated as the sum of the difference between the two ice sheet halves across the divide. This difference is 0 for all fields showing perfect symmetry.
D2 Temporal resolution test
Here we test the effect of changing the length of the time step in the basal hydrology on model solution using the SHMIP SQRT_TOPO setup (depicted in Fig. D3). The perrun discrepancy with respect to the shortest time step shown in Fig. D1 is calculated as
As a first test of convergence under increasing temporal resolution (decreasing time step length), the hydrology model was run to steady state under SHMIP scenario A (constant 2.5 mm yr^{−1}). Seeing convergence at shorter time steps for the steady forcing, an unsteady sinusoidal meltwater forcing was applied (50year period, 3.5 mm yr^{−1} amplitude). The convergence for the unsteady case is shown in Fig. D1 with the error metric of Eq. (D4). The rate of convergence is approximately quadratic as expected for the 𝒪(Δt^{2}) leapfrog trapezoidal scheme.
D3 Spatial resolution test
Here we show the effect of varying spatial resolution on the model solution. The model was run to steady state with prescribed melt and basal velocity (1.75 m yr^{−1} and 2.0 m yr^{−1} respectively). For this test, the SHMIP setup was used as shown in Fig. D3. The SQRT_TOPO flowline length from divide to toe was set to 2500 km, and the number of grid cells was adjusted, so that the highest resolution was Δx_{i}=22.66 km: $\left\{{n}_{i}=\mathrm{2}i+\mathrm{1}\right\}\text{for}i\in \left[\mathrm{11},\mathrm{121}\right],\phantom{\rule{0.33em}{0ex}}i\in \mathbb{N}$ and $\mathrm{\Delta}{x}_{i}=\mathrm{2500}\phantom{\rule{0.125em}{0ex}}\mathrm{km}/{n}_{i}$. The model solution at each resolution was linearly interpolated to the highestresolution grid, and the sum of the elementwise difference with the highest resolution was used for the (L1 norm) error metric, in keeping with Eq. (D4):
Figure D2 shows the convergence of model solutions (same set as Sect. D2) at increasing spatial resolution (shorter cell width). The numerical order of the upwind and finitevolume schemes used here is 𝒪(Δx). The approximately linear rate of convergence in Fig. D2 matches this numerical order.
D4 Mass conservation
Mass conservation is demonstrated by comparing flux at the margin to source rates of water or sediment within the ice sheet: the integral of the melt rate over the ice sheet less the total flux through the margin will give the change in basal water volume over time. Integrating this change up to each time step will give the basal water volume at each time step, which can be compared to modelcalculated basal water volume in order to assess mass conservation.
To test mass conservation with unsteady input, we applied a sinusoidal meltwater forcing,
(where T=50 kyr is the longest and highestamplitude period and melt =3.5 mm yr^{−1}), to the SQRT_TOPO setup (depicted in Fig. D3) and calculated basal sliding velocity dynamically as in Eq. (D2). Here we assume incompressibility of water such that volume is scaled mass.
A net volume of basal water time series was calculated by timeintegrating the net of input and output, net${}_{\mathrm{hyd}}^{{t}_{i}}$, up to each time step t_{i}:
where 𝒜 is the area covered by ice, m_{t} is the melt at the ice sheet base (Eq. D6), 𝒮 is the ice margin (interface beyond which ice thickness is 0), and $\mathit{Q}\cdot \widehat{\mathit{n}}$ is the flux through the margin, The net${}_{\mathrm{hyd}}^{{t}_{i}}$ time series was then compared against modelled total water volume (${V}_{\mathrm{hyd}}^{{t}_{i}}$) to calculate mass conservation error (ERR${}^{{t}_{i}}$):
where V_{hyd} is the volume of water under the ice sheet.
The dynamic model outputs from this test are summarized in Fig. D4. This mass conservation test shows a maximum error of 0.052 % between the model output and the calculation in Eq. (D7) (given in Eq. D8).
D5 Comparison with Werder et al. (2013) results for SHMIP
Results for this model are compared with output of the Glacier Drainage System model (GlaDS, Werder et al., 2013) employing the same physics: a continuum representation of a linkedcavity system with Darcy–Weisbach flux shown in Fig. D5. While the model of Werder et al. (2013) is similar to this one, there are noteworthy differences. Werder et al. (2013) use an unstructured mesh and finiteelement discretization, and the channel elements are always active (with water exchanged between the channels at the edges and the distributed system at the subdomains). This is in contrast to BrAHMs2.0, in which the channel system switches on in a particular cell given a flux criterion (and uses finitevolume discretization with a regular Cartesian grid). We therefore use the SHMIP scenario in which the least amount of channelized flux is active in order to get the most structurally consistent comparison between the two models. BrAHMs2.0 closely reproduces the flux and effective pressure solutions for this scenario, concluding our verification that we solved the equations “right”.
E1 Pressure closure of Bueler and van Pelt (2015)
Here we use the timevarying water pressure calculation of Bueler and van Pelt (2015). The rationale summarized here is shown by Bueler (2014). Here the subglacial and englacial hydrologic systems are assumed to be in perfect communication and their coevolution is described. The englacial hydrologic system is analogized to a rigid “pore space” (comprised of crevasses, moulins, englacial channels, and intergranular porosity). The total volume of water is the sum of englacial and subglacial water:
and the mass balance for incompressible water is
from total flux in and out of a control section of the system plus any sources (volume water, $\frac{m}{{\mathit{\rho}}_{\mathrm{w}}}$) within that section. This section is of the area Δx by Δy, and pressure in the connected englacialsubglacial system is given by the hydrostatic head in the englacial part:
The effective englacial porosity (ice volume relative proportion of connected englacial void space) is ϕ_{eng}. Cavity volume within an area of bed with roughness wavelength l_{r} (cavitygenerating obstacle spacing) is
where n_{cav} is the number of cavities in the given bed section and V_{cav} is their average volume. Differentiating this gives the change in pressure with time:
The $\left(\mathrm{1}/{l}_{\mathrm{r}}^{\mathrm{2}}\right)\frac{\partial {V}_{\mathrm{cav}}}{\partial t}=\frac{\partial {h}_{\mathrm{cav}}}{\partial t}$ derivative is given by the opening and closing balance in Eq. (7):
Here opening due to wall melting has been omitted (see assumption 1) relative to what is shown by Bueler (2014). As Δx→0 and Δy→0, the difference between the fluxes in versus out of the control section goes to the divergence of the fluxes within it.
with m_{t} the source of water in thickness per unit time. We assume that water only travels laterally through the subglacial system, and so all fluxes are through the linked cavities.
Analysis and data are available at Zenodo (https://doi.org/10.5281/zenodo.10018301, Drew, 2023).
MD performed subglacial hydrology model development, experimental design and execution, analysis, and writing. LT maintains the GSM, advised on experimental design, and edited the paper.
The contact author has declared that neither of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
The figures in this work were generated with Matplotlib (Hunter, 2007). Computational resources were provided by the Canadian Foundation for Innovation (CFI) and ACENET – a regional partner of the Digital Research Alliance of Canada. We thank the two anonymous reviewers for their helpful suggestions.
This work was supported by NSERC Discovery (grant no. RGPIN201806658), the Canadian Foundation for Innovation, and the German Federal Ministry of Education and Research initiative for Research for Sustainability through the PalMod project.
This paper was edited by Alexander Robinson and reviewed by two anonymous referees.
Alley, R.: WaterPressure Coupling of Sliding and Bed Deformation: I. Water System, J. Glaciol., 35, 108–118, https://doi.org/10.3189/002214389793701527, 1989. a
Alley, R. B.: How can lowpressure channels and deforming tills coexist subglacially?, J. Glaciol., 38, 200–207, https://doi.org/10.3189/s0022143000009734, 1992. a
Alley, R. B., Anandakrishnan, S., Bentley, C. R., and Lord, N.: A waterpiracy hypothesis for the stagnation of Ice Stream C, Antarctica, Ann. Glaciol., 20, 187–194, https://doi.org/10.3189/1994aog201187194, 1994. a
Anandakrishnan, S. and Alley, R. B.: Stagnation of Ice Stream C, West Antarctica by water piracy, Geophys. Res. Lett., 24, 265–268, https://doi.org/10.1029/96gl04016, 1997. a
Anderson, R. S., Hallet, B., Walder, J., and Aubry, B. F.: Observations in a cavity beneath grinnell glacier, Earth Surf. Proc. Land., 7, 63–70, https://doi.org/10.1002/esp.3290070108, 1982. a
Arakawa, A. and Lamb, V. R.: Computational design of the basic dynamical processes of the UCLA general circulation model, General circulation models of the atmosphere, 17, 173–265, https://doi.org/10.1016/B9780124608177.500094, 1977. a
Bennett, M. R.: Ice streams as the arteries of an ice sheet: their mechanics, stability and significance, EarthSci. Rev., 61, 309–339, https://doi.org/10.1016/s00128252(02)001307, 2003. a
Borgonovo, E.: A new uncertainty importance measure, Reliab. Eng. Syst. Safe., 92, 771–784, https://doi.org/10.1016/j.ress.2006.04.015, 2007. a
Brubaker, K. M., Myers, W. L., Drohan, P. J., Miller, D. A., and Boyer, E. W.: The Use of LiDAR Terrain Data in Characterizing Surface Roughness and Microtopography, Appl. Environ. Soil Sci., 2013, 1–13, https://doi.org/10.1155/2013/891534, 2013. a
Bueler, E.: Extending the lumped subglacial–englacial hydrology model of Bartholomaus and others (2011), J. Glaciol., 60, 808–810, https://doi.org/10.3189/2014jog14j075, 2014. a, b, c
Bueler, E. and van Pelt, W.: Massconserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635, https://doi.org/10.5194/gmd816132015, 2015. a, b, c, d, e
Calov, R., Greve, R., AbeOuchi, A., Bueler, E., Huybrechts, P., Johnson, J. V., Pattyn, F., Pollard, D., Ritz, C., Saito, F., and Tarasov, L.: Results from the IceSheet Model Intercomparison Project–Heinrich Event Intercomparison (ISMIP HEINO), J. Glaciol., 56, 371–383, https://doi.org/10.3189/002214310792447789, 2010. a, b, c, d, e
Carsey, F., Behar, A., Lonne Lane, A., Realmuto, V., and Engelhardt, H.: A borehole camera system for imaging the deep interior of ice sheets, J. Glaciol., 48, 622–628, https://doi.org/10.3189/172756502781831124, 2002. a
Chandler, D. M., Wadham, J. L., Lis, G. P., Cowton, T., Sole, A., Bartholomew, I., Telling, J., Nienow, P., Bagshaw, E. B., Mair, D., Vinen, S., and Hubbard, A.: Evolution of the subglacial drainage system beneath the Greenland Ice Sheet revealed by tracers, Nat. Geosci., 6, 195–198, https://doi.org/10.1038/ngeo1737, 2013. a, b, c
Clark, P. U., Archer, D., Pollard, D., Blum, J. D., Rial, J. A., Brovkin, V., Mix, A. C., Pisias, N. G., and Roy, M.: The middle Pleistocene transition: characteristics, mechanisms, and implications for longterm changes in atmospheric pCO_{2}, Quaternary Sci. Rev., 25, 3150–3184, https://doi.org/10.1016/j.quascirev.2006.07.008, 2006. a
Clarke, G. K. C.: Lumpedelement analysis of subglacial hydraulic circuits, J. Geophys. Res.Sol. Ea., 101, 17547–17559, https://doi.org/10.1029/96jb01508, 1996. a
Cook, S. J., Christoffersen, P., and Todd, J.: A fullycoupled 3D model of a large Greenlandic outlet glacier with evolving subglacial hydrology, frontal plume melting and calving, J. Glaciol., 68, 486–502, https://doi.org/10.1017/jog.2021.109, 2021. a
Courant, R., Friedrichs, K., and Lewy, H.: Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann., 100, 32–74, https://doi.org/10.1007/bf01448839, 1928. a, b
Cuffey, K. and Paterson, W. S. B.: The Physics of Glaciers, Academic Press, 4 Edn., ISBN 9780123694614, 2010. a, b
Dalton, A. S., Margold, M., Stokes, C. R., Tarasov, L., Dyke, A. S., Adams, R. S., Allard, S., Arends, H. E., Atkinson, N., Attig, J. W., Barnett, P. J., Barnett, R. L., Batterson, M., Bernatchez, P., Borns, H. W., Breckenridge, A., Briner, J. P., Brouard, E., Campbell, J. E., Carlson, A. E., Clague, J. J., Curry, B. B., Daigneault, R.A., DubéLoubert, H., Easterbrook, D. J., Franzi, D. A., Friedrich, H. G., Funder, S., Gauthier, M. S., Gowan, A. S., Harris, K. L., Hétu, B., Hooyer, T. S., Jennings, C. E., Johnson, M. D., Kehew, A. E., Kelley, S. E., Kerr, D., King, E. L., Kjeldsen, K. K., Knaeble, A. R., Lajeunesse, P., Lakeman, T. R., Lamothe, M., Larson, P., Lavoie, M., Loope, H. M., Lowell, T. V., Lusardi, B. A., Manz, L., McMartin, I., Nixon, F. C., Occhietti, S., Parkhill, M. A., Piper, D. J., Pronk, A. G., Richard, P. J., Ridge, J. C., Ross, M., Roy, M., Seaman, A., Shaw, J., Stea, R. R., Teller, J. T., Thompson, W. B., Thorleifson, L. H., Utting, D. J., Veillette, J. J., Ward, B. C., Weddle, T. K., and Wright, H. E.: An updated radiocarbonbased ice margin chronology for the last deglaciation of the North American Ice Sheet Complex, Quaternary Sci. Rev., 234, 106223, https://doi.org/10.1016/j.quascirev.2020.106223, 2020. a
Darcy, H.: Les Fontaines Publiques De La Ville De Dijon: Exposition Et Application Des Principes A Suivre Et Des Formules A Employer Dans Les Questions De Distribution D'eau, Libraire des corps impériaux des ponts et chaussées et des mines, 49 Quai des Augustins, Paris, France, 1st Edn., 1856. a
de Fleurian, B., Werder, M. A., Beyer, S., Brinkerhoff, D. J., Delaney, I., Dow, C. F., Downs, J., Gagliardini, O., Hoffman, M. J., Hooke, R. L., Seguinot, J., and Sommers, A. N.: SHMIP The subglacial hydrology model intercomparison Project, J. Glaciol., 64, 897–916, https://doi.org/10.1017/jog.2018.78, 2018. a, b, c
Dowdeswell, J. A., Maslin, M. A., Andrews, J. T., and McCave, I. N.: Iceberg production, debris rafting, and the extent and thickness of Heinrich layers (H1, H2) in North Atlantic sediments, Geology, 23, 301, https://doi.org/10.1130/00917613(1995)023<0297:ipdrat>2.3.co;2, 1995. a
Drew, M.: Analysis scripts and data for paper “Surging of a Hudson Strait Scale Ice Stream: Subglacial hydrology matters but the process details mostly don't”, Zenodo [code and data set], https://doi.org/10.5281/zenodo.10018301, 2023. a
Erokhina, O., Rogozhina, I., Prange, M., Bakker, P., Bernales, J., Paul, A., and Schulz, M.: Dependence of slope lapse rate over the Greenland ice sheet on background climate, J. Glaciol., 63, 568–572, https://doi.org/10.1017/jog.2017.10, 2017. a
Flowers, G.: A Multicomponent Coupled Model of Glacier Hydrology, Ph.D. thesis, University of British Columbia, Department of Earth, Ocean and Atmospheric Sciences Faculty of Science 2020–2207 Main Mall Vancouver, BC Canada V6T 1Z4, 2000. a, b, c, d, e, f, g, h, i, j, k
Flowers, G. E.: Modelling water flow under glaciers and ice sheets, Proceedings of the Royal Society A: Mathematical, Phys. Eng. Sci., 471, 1–41, https://doi.org/10.1098/rspa.2014.0907, 2015. a, b, c, d, e
Flowers, G. E.: Hydrology and the future of the Greenland Ice Sheet, Nat. Commun., 9, 2729, https://doi.org/10.1038/s41467018050020, 2018. a
Flowers, G. E. and Clarke, G. K. C.: A multicomponent coupled model of glacier hydrology 1. Theory and synthetic examples, J. Geophys. Res.Sol. Ea., 107, ECV 9–1–ECV 9–17, https://doi.org/10.1029/2001jb001122, 2002. a
Fowler, A. C.: Instability modelling of drumlin formation incorporating leeside cavity growth, P. Roy. Soc. A, 465, 2681–2702, https://doi.org/10.1098/rspa.2008.0490, 2009. a
Gandy, N., Gregoire, L. J., Ely, J. C., Cornford, S. L., Clark, C. D., and Hodgson, D. M.: Exploring the ingredients required to successfully model the placement, generation, and evolution of ice streams in the BritishIrish Ice Sheet, Quaternary Sci. Rev., 223, 105915, https://doi.org/10.1016/j.quascirev.2019.105915, 2019. a, b
Grinsted, A., Hvidberg, C. S., Lilien, D. A., Rathmann, N. M., Karlsson, N. B., Gerber, T., Kjær, H. A., Vallelonga, P., and DahlJensen, D.: Accelerating ice flow at the onset of the Northeast Greenland Ice Stream, Nat. Commun., 13, 5589, https://doi.org/10.1038/s41467022329992, 2022. a, b
Hager, A. O., Hoffman, M. J., Price, S. F., and Schroeder, D. M.: Persistent, extensive channelized drainage modeled beneath Thwaites Glacier, West Antarctica, The Cryosphere, 16, 3575–3599, https://doi.org/10.5194/tc1635752022, 2022. a
Hank, K., Tarasov, L., and Mantelli, E.: Modeling sensitivities of thermally and hydraulically driven ice stream surge cycling, Geosci. Model Dev., 16, 5627–5652, https://doi.org/10.5194/gmd1656272023, 2023. a, b, c
Hemming, S. R.: Heinrich events: Massive late Pleistocene detritus layers of the North Atlantic and their global climate imprint, Rev. Geophys., 42, RG1005, https://doi.org/10.1029/2003rg000128, 2004. a, b
Herman, J. and Usher, W.: SALib: An opensource Python library for Sensitivity Analysis, J. Open Source Softw., 2, 97, https://doi.org/10.21105/joss.00097, 2017. a
Hewitt, I.: Seasonal changes in ice sheet motion due to melt water lubrication, Earth Planet. Sc. Lett., 371–372, 16–25, https://doi.org/10.1016/j.epsl.2013.04.022, 2013. a
Hunter, J. D.: Matplotlib: A 2D graphics environment, Comput. Sci. Eng., 9, 90–95, https://doi.org/10.1109/MCSE.2007.55, 2007. a
Iverson, N. R.: A theory of glacial quarrying for landscape evolution models, Geology, 40, 679–682, https://doi.org/10.1130/g33079.1, 2012. a
Joughin, I., Tulaczyk, S., Fahnestock, M., and Kwok, R.: A MiniSurge on the Ryder Glacier, Greenland, Observed by Satellite Radar Interferometry, Science, 274, 228–230, https://doi.org/10.1126/science.274.5285.228, 1996. a
Kageyama, M., Harrison, S. P., Kapsch, M.L., Lofverstrom, M., Lora, J. M., Mikolajewicz, U., SherriffTadano, S., Vadsaria, T., AbeOuchi, A., Bouttes, N., Chandan, D., Gregoire, L. J., Ivanovic, R. F., Izumi, K., LeGrande, A. N., Lhardy, F., Lohmann, G., Morozova, P. A., Ohgaito, R., Paul, A., Peltier, W. R., Poulsen, C. J., Quiquet, A., Roche, D. M., Shi, X., Tierney, J. E., Valdes, P. J., Volodin, E., and Zhu, J.: The PMIP4 Last Glacial Maximum experiments: preliminary results and comparison with the PMIP3 simulations, Clim. Past, 17, 1065–1089, https://doi.org/10.5194/cp1710652021, 2021. a, b, c
Kamb, B.: Glacier surge mechanism based on linked cavity configuration of the basal water conduit system, J. Geophys. Res., 92, 9083–9100, https://doi.org/10.1029/jb092ib09p09083, 1987. a
Kavanagh, M. and Tarasov, L.: BrAHMs V1.0: a fast, physically based subglacial hydrology model for continentalscale application, Geosci. Model Dev., 11, 3497–3513, https://doi.org/10.5194/gmd1134972018, 2018. a, b, c
Kingslake, J. and Ng, F.: Modelling the coupling of flood discharge with glacier flow during jökulhlaups, Ann. Glaciol., 54, 25–31, https://doi.org/10.3189/2013aog63a331, 2013. a, b
Lawrence, M. G.: The Relationship between Relative Humidity and the Dewpoint Temperature in Moist Air: A Simple Conversion and Applications, B. Am. Meteorol. Soc., 86, 225–234, https://doi.org/10.1175/bams862225, 2005. a
MacAyeal, D. R.: Binge/purge oscillations of the Laurentide Ice Sheet as a cause of the North Atlantic’s Heinrich events, Paleoceanography, 8, 775–784, https://doi.org/10.1029/93pa02200, 1993. a, b
Mathews, W. H.: Surface profiles of the laurentide ice sheet in its marginal areas, J. Glaciol., 13, 37–43, https://doi.org/10.1017/s0022143000023352, 1974. a
Morlighem, M., Rignot, E., Mouginot, J., Wu, X., Seroussi, H., Larour, E., and Paden, J.: Highresolution bed topography mapping of Russell Glacier, Greenland, inferred from Operation IceBridge data, J. Glaciol., 59, 1015–1023, https://doi.org/10.3189/2013jog12j235, 2013. a
Muskat, M.: The Flow of Compressible Fluids Through Porous Media and Some Problems in Heat Conduction, Physics, 5, 71–94, https://doi.org/10.1063/1.1745233, 1934. a
Naafs, B., Hefter, J., and Stein, R.: Millennialscale ice rafting events and Hudson Strait Heinrich(like) Events during the late Pliocene and Pleistocene: a review, Quaternary Sci. Rev., 80, 1–28, https://doi.org/10.1016/j.quascirev.2013.08.014, 2013. a, b
Oreskes, N., ShraderFrechette, K., and Belitz, K.: Verification, Validation, and Confirmation of Numerical Models in the Earth Sciences, Science, 263, 641–646, https://doi.org/10.1126/science.263.5147.641, 1994. a
Ou, H.W.: A theory of glacier dynamics and instabilities Part 1: Topographically confined glaciers, J. Glaciol., 68, https://doi.org/10.1017/jog.2021.20, 2021. a
Patankar, S. V.: Numerical Heat Transfer and Fluid Flow, McGrawHill Book Company, 1221 Avenue of the Americas, New York, New York, USA, 1st Edn., https://doi.org/10.1201/9781482234213, 1980. a
Payne, A. J.: Dynamics of the Siple Coast ice streams, west Antarctica: Results from a thermomechanical ice sheet model, Geophys. Res. Lett., 25, 3173–3176, https://doi.org/10.1029/98gl52327, 1998. a
Payne, A. J. and Dongelmans, P. W.: Selforganization in the thermomechanical flow of ice sheets, J. Geophys. Res.Sol. Ea., 102, 12219–12233, https://doi.org/10.1029/97jb00513, 1997. a
Payne, A. J., Huybrechts, P., AbeOuchi, A., Calov, R., Fastook, J. L., Greve, R., Marshall, S. J., Marsiat, I., Ritz, C., Tarasov, L., and Thomassen, M. P. A.: Results from the EISMINT model intercomparison: the effects of thermomechanical coupling, J. Glaciol., 46, 227–238, https://doi.org/10.3189/172756500781832891, 2000. a, b
Pelletier, J. D., Broxton, P. D., Hazenberg, P., Zeng, X., Troch, P. A., Niu, G.Y., Williams, Z., Brunke, M. A., and Gochis, D.: A gridded global data set of soil, intact regolith, and sedimentary deposit thicknesses for regional and global land surface modeling, J. Adv. Model. Earth Sy., 8, 41–65, https://doi.org/10.1002/2015ms000526, 2016. a
Pianosi, F. and Wagener, T.: A simple and efficient method for global sensitivity analysis based on cumulative distribution functions, Environ. Modell. Softw., 67, 1–11, https://doi.org/10.1016/j.envsoft.2015.01.004, 2015. a
Pollard, D. and DeConto, R. M.: A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica, The Cryosphere, 6, 953–971, https://doi.org/10.5194/tc69532012, 2012. a
Rada, C. and Schoof, C.: Channelized, distributed, and disconnected: subglacial drainage under a valley glacier in the Yukon, The Cryosphere, 12, 2609–2636, https://doi.org/10.5194/tc1226092018, 2018. a
Retzlaff, R. and Bentley, C. R.: Timing of stagnation of Ice Stream C, West Antarctica, from shortpulse radar studies of buried surface crevasses, J. Glaciol., 39, 553–561, https://doi.org/10.3189/s0022143000016440, 1993. a
Roache, P. J.: Quantification of uncertainty in computational fluid dynamics, Annu. Rev. Fluid Mech., 29, 123–160, https://doi.org/10.1146/annurev.fluid.29.1.123, 1997. a, b
Röthlisberger, H.: Water Pressure in Intra and Subglacial Channels, J. Glaciol., 11, 177–203, https://doi.org/10.3189/s0022143000022188, 1972. a, b
Saltelli, A.: Making best use of model evaluations to compute sensitivity indices, Comput. Phys. Commun., 145, 280–297, https://doi.org/10.1016/s00104655(02)002801, 2002. a
Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., and Tarantola, S.: Global Sensitivity Analysis, John Wiley & Sons, Ltd., The Atrium, Southern Gate, Chichester, West Sussex PO19 8SQ, England, 1st Edn., ISBN 9780470059975, 2008. a
Schoof, C.: Cavitation on Deformable Glacier Beds, SIAM J. Appl. Math., 67, 1633–1653, https://doi.org/10.1137/050646470, 2007. a
Schoof, C.: Icesheet acceleration driven by melt supply variability, Nature, 468, 803–806, https://doi.org/10.1038/nature09618, 2010. a, b, c, d, e, f, g, h, i, j, k, l, m
Siegfried, M. R., Fricker, H. A., Carter, S. P., and Tulaczyk, S.: Episodic ice velocity fluctuations triggered by a subglacial flood in West Antarctica, Geophys. Res. Lett., 43, 2640–2648, https://doi.org/10.1002/2016gl067758, 2016. a
Sobol, I.: Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Math. Comput. Sim., 55, 271–280, https://doi.org/10.1016/s03784754(00)002706, 2001. a
Sornette, D., Davis, A. B., Ide, K., Vixie, K. R., Pisarenko, V., and Kamm, J. R.: Algorithm for model validation: Theory and applications, P. Natl. Acad. Sci. USA, 104, 6562–6567, https://doi.org/10.1073/pnas.0611677104, 2007. a, b
Tadmor, E.: A review of numerical methods for nonlinear partial differential equations, B. Am. Math. Soc., 49, 507–554, https://doi.org/10.1090/s027309792012013794, 2012. a
Tarasov, L. and Peltier, W.: A calibrated deglacial drainage chronology for the North American continent: evidence of an Arctic trigger for the Younger Dryas, Quaternary Sci. Rev., 25, 659–688, https://doi.org/10.1016/j.quascirev.2005.12.006, 2006. a
Tarasov, L. and Peltier, W. R.: Coevolution of continental ice cover and permafrost extent over the last glacialinterglacial cycle in North America, J. Geophys. Res., 112, F02S08, https://doi.org/10.1029/2006jf000661, 2007. a
Tarasov, L., Dyke, A. S., Neal, R. M., and Peltier, W.: A datacalibrated distribution of deglacial chronologies for the North American ice complex from glaciological modeling, Earth Planet. Sc. Lett., 315–316, 30–40, https://doi.org/10.1016/j.epsl.2011.09.010, 2012. a, b
Tarasov, L., Hank, K., and Lecavalier, B.: The glacial systems model (GSM) Version 23C, in preparation, 2023. a, b
Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, İ., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nat. Methods, 17, 261–272, https://doi.org/10.1038/s4159201906862, 2020. a
Walder, J. and Hallet, B.: Geometry of Former Subglacial Water Channels and Cavities, J. Glaciol., 23, 335–346, https://doi.org/10.1017/s0022143000029944, 1979. a, b
Walder, J. S.: Hydraulics of Subglacial Cavities, J. Glaciol., 32, 439–445, https://doi.org/10.3189/s0022143000012156, 1986. a, b, c, d
Walder, J. S. and Fowler, A.: Channelized subglacial drainage over a deformable bed, J. Glaciol.y, 40, 3–15, https://doi.org/10.3189/s0022143000003750, 1994. a, b
Weisbach, J.: Lehrbuch der ingenieur und maschinenmechanik, Vieweg, Braunschweig, Germany, 1st Edn., https://babel.hathitrust.org/cgi/pt?id=wu.89088908009&view=1up&seq=8 (last access: 12 December 2023), 1855. a
Werder, M. A., Hewitt, I. J., Schoof, C. G., and Flowers, G. E.: Modeling channelized and distributed subglacial drainage in two dimensions, J. Geophys. Res.Earth, 118, 2140–2158, https://doi.org/10.1002/jgrf.20146, 2013. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p
Willis, M. J., Zheng, W., Durkin, W. J., Pritchard, M. E., Ramage, J. M., Dowdeswell, J. A., Benham, T. J., Bassford, R. P., Stearns, L. A., Glazovsky, A. F., Macheret, Y. Y., and Porter, C. C.: Massive destabilization of an Arctic ice cap, Earth Planet. Sc. Lett., 502, 146–155, https://doi.org/10.1016/j.epsl.2018.08.049, 2018. a
 Abstract
 Introduction
 Subglacial hydrology
 Model description
 LISsq experimental design
 HS surging results
 Discussion of surge contribution
 Conclusions
 Appendix A: Surging with thinner and thicker ice sheets
 Appendix B: Subglacial hydrology model solver
 Appendix C: Subglacial hydrology model assumptions
 Appendix D: Subglacial hydrology model verification
 Appendix E: Discretization
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Subglacial hydrology
 Model description
 LISsq experimental design
 HS surging results
 Discussion of surge contribution
 Conclusions
 Appendix A: Surging with thinner and thicker ice sheets
 Appendix B: Subglacial hydrology model solver
 Appendix C: Subglacial hydrology model assumptions
 Appendix D: Subglacial hydrology model verification
 Appendix E: Discretization
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References