Hostname: page-component-7c8c6479df-nwzlb Total loading time: 0 Render date: 2024-03-17T14:53:57.674Z Has data issue: false hasContentIssue false

Sediment behavior controls equilibrium width of subglacial channels

Published online by Cambridge University Press:  27 November 2017

ANDERS DAMSGAARD*
Affiliation:
Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography, University of California, San Diego, CA, USA
JENNY SUCKALE
Affiliation:
Department of Geophysics, Stanford University, Stanford, CA, USA Institute of Computational and Mathematical Engineering, Stanford University, Stanford, CA, USA
JAN A. PIOTROWSKI
Affiliation:
Department of Geoscience, Aarhus University, Aarhus, Denmark
MORGANE HOUSSAIS
Affiliation:
Levich Institute, City College of CUNY, 140th Street and Convent Avenue, New York, NY, USA
MATTHEW R. SIEGFRIED
Affiliation:
Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography, University of California, San Diego, CA, USA Department of Geophysics, Stanford University, Stanford, CA, USA
HELEN A. FRICKER
Affiliation:
Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography, University of California, San Diego, CA, USA
*
Correspondence: Anders Damsgaard <andersd@princeton.edu>
Rights & Permissions [Opens in a new window]

Abstract

Flow-frictional resistance at the base of glaciers and ice sheets is strongly linked to subglacial water pressure. Understanding the physical mechanisms that govern meltwater fluxes in subglacial channels is hence critical for constraining variations in ice flow. Previous mathematical descriptions of soft-bed subglacial channels assume a viscous till rheology, which is inconsistent with laboratory data and the majority of field studies. Here, we use a grain-scale numerical formulation coupled to pore-water dynamics to analyze the structural stability of channels carved into soft beds. Contrary to the soft-bed channel models assuming viscous till rheology, we show that the flanks of till channels can support substantial ice loads without creep closure of the channel, because the sediment has finite frictional strength. Increased normal stress on the channel flanks causes plastic failure of the sediment, and the channel rapidly shrinks to increase the ice-bed contact area. We derive a new parameterization for subglacial channelized flow on soft beds and show that channel dynamics are dominated by fluvial erosion and deposition processes with thresholds linked to the plastic rheology of subglacial tills. We infer that the described limits to channel size may cause subglacial drainage to arrange in networks of multiple closely spaced channels.

Type
Papers
Creative Commons
Creative Common License - CCCreative Common License - BYCreative Common License - NCCreative Common License - SA
This is an Open Access article, distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike licence (http://creativecommons.org/licenses/by-nc-sa/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the same Creative Commons licence is included and the original work is properly cited. The written permission of Cambridge University Press must be obtained for commercial re-use.
Copyright
Copyright © The Author(s) 2017

1. INTRODUCTION

Channelization of subglacial meltwater is well documented from the sedimentary and geomorphological record (e.g., Piotrowski, Reference Piotrowski1999; Piotrowski and others, Reference Piotrowski, Geletneky and Vater1999; Jørgensen and Sandersen, Reference Jørgensen and Sandersen2006; Greenwood and others, Reference Greenwood, Clason, Nyberg, Jakobsson and Holmlund2016; Bjarnadóttir and others, Reference Bjarnadóttir, Winsborrow and Andreassen2017; Simkins and others, Reference Simkins2017). Several geophysical investigations support that channelization of subglacial meltwater is an important process under contemporary glaciers and ice sheets (e.g., Hubbard and Nienow, Reference Hubbard and Nienow1997; Winberry and others, Reference Winberry, Anandakrishnan and Alley2009; Horgan and others, Reference Horgan2013; Le Brocq and others, Reference Le Brocq2013; Schroeder and others, Reference Schroeder, Blankenship and Young2013; Gimbert and others, Reference Gimbert, Tsai, Amundson, Bartholomaus and Walter2016; Drews and others, Reference Drews2017). Numerical models of Antarctic ice-stream flow have shown that channelized water flow is necessary for producing sufficient variations in water pressure and subglacial friction, to match observed variations in ice-surface velocity (e.g., Thompson and others, Reference Thompson, Simons and Tsai2014; Rosier and others, Reference Rosier, Gudmundsson and Green2015), and meltwater channelization appears to be necessary for stabilizing ice-stream shear margins (Suckale and others, Reference Suckale, Platt, Perol and Rice2014; Perol and Rice, Reference Perol and Rice2015; Perol and others, Reference Perol, Rice, Platt and Suckale2015; Elsworth and Suckale, Reference Elsworth and Suckale2016). Channelized meltwater flow is especially important for water and ice dynamics during subglacial lake drainage events (e.g., Bartholomew and others, Reference Bartholomew2010; Stevens and others, Reference Stevens2015; Brinkerhoff and others, Reference Brinkerhoff, Meyer, Bueler, Truffer and Bartholomaus2016; Dow and others, Reference Dow, Werder, Nowicki and Walker2016; Fricker and others, Reference Fricker, Siegfried, Carter and Scambos2016; Carter and others, Reference Carter, Fricker and Siegfried2017; Simkins and others, Reference Simkins2017). Recent studies have detailed that the influence of channel dynamics extends downstream of the grounding line, since point release of water through channels at grounding lines greatly enhances basal melt rates through the formation of turbulent buoyant plumes (e.g., Jenkins, Reference Jenkins2011; Xu and others, Reference Xu, Rignot, Menemenlis and Koppes2012; Le Brocq and others, Reference Le Brocq2013; Marsh and others, Reference Marsh2016; Alley and others, Reference Alley, Scambos, Siegfried and Fricker2016; Drews and others, Reference Drews2017).

Channelization is not only relevant for projecting ice fluxes, but also for understanding subglacial landforms as water flow in the channels may be a primary mechanism for eroding and delivering sediment from the ice-sheet or glacier interior to its margin (e.g., Swift and others, Reference Swift, Nienow, Spedding and Hoey2002; Alley and others, Reference Alley, Lawson, Larson, Evenson and Baker2003; Kehew and others, Reference Kehew, Piotrowski and Jørgensen2012; Bjarnadóttir and others, Reference Bjarnadóttir, Winsborrow and Andreassen2017). Recent discoveries of curvilinear and erosive landforms near the last-glacial maximum terminus (Fig. 1) have raised questions on how subglacial drainage systems arrange on soft beds (Lesemann and others, Reference Lesemann, Piotrowski and Wysota2010, Reference Lesemann, Piotrowski and Wysota2014).

Fig. 1. Glacial curvilineations incised into sedimentary plateaus by subglacial meltwater erosion during the last glaciation near Zbójno, Poland, from Lesemann and others (Reference Lesemann, Piotrowski and Wysota2010). Water and ice-sheet flow was generally from the top left (NW) to the lower right (SE).

Glacier beds provide resistive friction to glaciers and ice sheets and limit how fast they flow, and are typically subdivided into hard (rigid and impermeable) or soft (deformable and typically sedimentary) types. The amount of resistance provided by the bed is strongly dependent on subglacial water pressure (e.g., Weertman, Reference Weertman1957; Lliboutry, Reference Lliboutry1968; Budd and others, Reference Budd, Keage and Blundy1979; Iken, Reference Iken1981; Bindschadler, Reference Bindschadler1983; Boulton and Hindmarsh, Reference Boulton and Hindmarsh1987; Fowler, Reference Fowler1987; Kamb, Reference Kamb1991; Hooke and others, Reference Hooke, Hanson, Iverson, Jansson and Fischer1997), and subglacial channels are the most efficient drainage type and water-pressure altering component at the ice–bed interface (e.g., Flowers, Reference Flowers2015). Meltwater evacuation through channels generally decreases pore-water pressure and increases basal strength in the adjacent areas under the ice (e.g., Shoemaker, Reference Shoemaker1986; Boulton and Hindmarsh, Reference Boulton and Hindmarsh1987; Rempel, Reference Rempel2009). Since changes in subglacial water-flow patterns directly affect the glacier stress balance, changes in subglacial hydrology can cause large-scale rearrangement in the flow patterns of ice sheets (e.g., Raymond, Reference Raymond2000; Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000a,Reference Tulaczyk, Kamb and Engelhardtb; Boulton and others, Reference Boulton, Hagdorn, Maillot and Zatsepin2009; Piotrowski and others, Reference Piotrowski, Hermanowski and Piechota2009; Bougamont and others, Reference Bougamont2015; Elsworth and Suckale, Reference Elsworth and Suckale2016).

The importance of changes in subglacial hydrology on ice flow has sparked efforts to develop and improve mathematical descriptions of subglacial water-pressure evolution in ice-flow models (e.g., Bougamont and others, Reference Bougamont, Price, Christoffersen and Payne2011; van der Wel and others, Reference van der Wel, Christoffersen and Bougamont2013; Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013; de Fleurian and others, Reference de Fleurian2014; Kyrke-Smith and others, Reference Kyrke-Smith, Katz and Fowler2015). Building on the pioneering work of Röthlisberger (Reference Röthlisberger1972), Shreve (Reference Shreve1972), and Nye (Reference Nye1976), which established the mathematical framework of dynamic water flow through channels melted into basal ice (R-channels), recent studies have improved our understanding of subglacial water drainage as a coupled system on hard beds. Differences in drainage properties between distributed and discrete drainage modes hint towards a complex interplay during non-steady drainage (e.g., Schoof, Reference Schoof2010; Hewitt, Reference Hewitt2011, Reference Hewitt2013; Kingslake and Ng, Reference Kingslake and Ng2013; Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013). Schoof (Reference Schoof2010) and Hewitt (Reference Hewitt2013) demonstrated that variability in meltwater input has the potential to cause ice-flow acceleration, because transient increases in meltwater supply overwhelm the water-transport capacity of the subglacial hydrological system. However, the majority of ice in contemporary and past ice sheets moves over soft beds (e.g., Alley and others, Reference Alley, Blankenship, Bentley and Rooney1986; Engelhardt and others, Reference Engelhardt, Humphrey, Kamb and Fahnestock1990; Anandakrishnan and others, Reference Anandakrishnan, Blankenship, Alley and Stoffa1998; Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000a; Kamb, Reference Kamb2001; Stokes and Clark, Reference Stokes and Clark2001; Piotrowski and others, Reference Piotrowski, Larsen and Junge2004; King and others, Reference King, Hindmarsh and Stokes2009; Christianson and others, Reference Christianson2014; Greenwood and others, Reference Greenwood, Clason, Nyberg, Jakobsson and Holmlund2016; Bjarnadóttir and others, Reference Bjarnadóttir, Winsborrow and Andreassen2017; Kulessa and others, Reference Kulessa2017; Simkins and others, Reference Simkins2017), but soft-bed channelized drainage is not included in current ice-sheet models despite efforts to parameterize the governing characteristics.

At present, our understanding of soft-bed hydrological processes is limited. Shoemaker (Reference Shoemaker1986) presented the first mathematical analysis of subglacial hydrology and sediment stability for a soft-bedded ice sheet drained by a combination of porous (groundwater) flow and channelized drainage. He assumed that erosion of sediment at the channel floor would incise channels into the soft bed, but argued that the increase in effective stress on the channel flanks would increase sediment strength sufficiently to make channels stable at any size.

Alley (Reference Alley1992) investigated the mechanical stability of a simplified channel geometry, and applied analytical relationships for till mechanics, including perfect plasticity and Bingham rheology, in part based on parameter fits by Boulton and Hindmarsh (Reference Boulton and Hindmarsh1987). While offering convenience by uniquely linking stress and strain in mathematical models, the weakly non-linear viscous till rheologies proposed by Boulton and Hindmarsh (Reference Boulton and Hindmarsh1987) disagree with the nearly plastic sediment rheology known from fundamental granular and soil mechanics (e.g., Schofield and Wroth, Reference Schofield and Wroth1968; Nedderman, Reference Nedderman1992; Terzaghi and others, Reference Terzaghi, Peck and Mesri1996; Mitchell and Soga, Reference Mitchell and Soga2005), field measurements on subglacial till deformation (e.g., Hooke and others, Reference Hooke, Hanson, Iverson, Jansson and Fischer1997; Kavanaugh and Clarke, Reference Kavanaugh and Clarke2006), laboratory deformation experiments on till (e.g., Kamb, Reference Kamb1991; Iverson and others, Reference Iverson, Hooyer and Baker1998; Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000a; Iverson and others, Reference Iverson2007), inversion of subglacial mechanics from ice-surface velocities (e.g., Tulaczyk, Reference Tulaczyk2006; Walker and others, Reference Walker, Christianson, Parizek, Anandakrishnan and Alley2012; Goldberg and others, Reference Goldberg, Schoof and Sergienko2014; Minchew and others, Reference Minchew2016; Gillet-Chaulet and others, Reference Gillet-Chaulet2016), and numerical experiments (e.g., Iverson and Iverson, Reference Iverson and Iverson2001; Kavanaugh and Clarke, Reference Kavanaugh and Clarke2006; Damsgaard and others, Reference Damsgaard2013, Reference Damsgaard2016).

Walder and Fowler (Reference Walder and Fowler1994) combined an earlier mathematical formulation of R-channel closure (Nye, Reference Nye1976) with the mildly non-linear viscous closure relationship for conduits in till (Boulton and Hindmarsh, Reference Boulton and Hindmarsh1987; Fowler and Walder, Reference Fowler and Walder1993) and derived a new mathematical model for subglacial channels in soft beds. Contrary to the Bingham visco-plastic relationship applied by Alley (Reference Alley1992), Walder and Fowler (Reference Walder and Fowler1994) ignored sediment yield strength and parameterized till to continuously creep toward the channel and counteract fluvial erosion. In their mathematical framework, they demonstrated that the type of channelized drainage is a function of surface slope. R-channels are likely to form under steep-sloped mountain glaciers, while fluvial incision into the soft bed is hypothesized dominate under relatively flat parts of ice sheets. Ng (Reference Ng2000a) further developed the mathematical theory by Walder and Fowler (Reference Walder and Fowler1994), and demonstrated that sediment erosion and deposition by fluvial transport is more important than creep closure of the idealized viscous sediment. A similar model, also with viscous till rheology and fluvial channel incision into soft beds, was shown to effectively approximate drainage histories of Antarctic subglacial lakes better than a R-channel model with erosion into the ice base (Carter and others, Reference Carter, Fricker and Siegfried2017).

However, the models assuming viscous till rheology provide a continuous flux of sediment towards the channel, and disregard sediment yield strength and associated plastic failure limits. Yet, subglacial till is known to behave like other clastic sediments with a nearly perfect plastic and rate-independent rheology with a yield strength τ y governed by the Mohr–Coulomb constitutive relation, τ y = C + N tan(ϕ u), where C is the effective sediment cohesion, N is the effective stress normal to the shear plane, and ϕ u is the angle of internal friction (e.g., Terzaghi, Reference Terzaghi1943). As a matter of fact, the plastic behavior of sediment beds in general, with or without cohesive forces, is related to the general behavior of granular material (e.g., GDR-MiDi, 2004; Houssais and others, Reference Houssais, Ortiz, Durian and Jerolmack2015; Houssais and Jerolmack, Reference Houssais and Jerolmack2017). In this study, we use numerical simulations of purely granular material (without cohesion) to move one step further in understanding the impact of sediment bed plasticity on the subglacial channel shape dynamics and the subglacial hydrology. In particular, we will study the effect of the effective normal stress at the ice–bed interface, and investigate how strong horizontal pore-pressure gradients toward the channel impact the shape and stability of the channel.

In the next section, we present the physical background of subglacial channel modeling, the simulation method we used, and how simulations of granular material were initialized and performed. Afterwards we present and discuss our results and combine our findings with established approaches for sediment and water in a new continuum formulation for subglacial soft-bed channels.

2. BACKGROUND AND METHODS

2.1. Continuum modeling of subglacial channels

Mathematical models of subglacial hydrology contain several balance equations, related to conservation of water mass, transients in hydraulic properties related to conduit opening and closure, and conservation of water momentum and energy (e.g., Nye, Reference Nye1976; Walder and Fowler, Reference Walder and Fowler1994; Clarke, Reference Clarke2005; Schoof, Reference Schoof2010; Hewitt, Reference Hewitt2011; Kingslake and Ng, Reference Kingslake and Ng2013; Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013; Flowers, Reference Flowers2015). The cross-sectional size of a soft-bed subglacial channel evolves by the combined effect of melting and creep closure of the channel roof, and fluvial sediment erosion and deposition at the channel base. Here, we assume that the channel is governed by sediment flux alone, implying that the ice interface remains planar. Simultaneous incision of both the ice roof and sedimentary bed, may occur under more energetic conditions (e.g., Alley, Reference Alley1989; Walder and Fowler, Reference Walder and Fowler1994; Carter and others, Reference Carter, Fricker and Siegfried2017). Changes in the channel cross-sectional area over time ∂S/∂t are balanced by the along-channel gradient of fluvial sediment flux Q s through the Exner equation with a representative bed porosity Φ (e.g., Ng, Reference Ng2000a):

(1) $$ {\partial S \over \partial t} = {1 \over 1 - \Phi} {\partial Q_{\rm s} \over \partial s} $$

where s is the channel streamwise dimension. Significant effort has been devoted to constraining the relationships behind sediment transport in fluvial settings, and this is an ongoing topic of research (e.g., Lajeunesse and others, Reference Lajeunesse, Malverti and Charru2010). As a result there are numerous propositions for parameterizing the sediment transport, Q s (e.g., Shields, Reference Shields1936; Meyer-Peter and Müller, Reference Meyer-Peter and Müller1948; Einstein, Reference Einstein1950; Parker, Reference Parker1978; Stock and Montgomery, Reference Stock and Montgomery1999; Whipple and Tucker, Reference Whipple and Tucker1999; Ng, Reference Ng2000a; Davy and Lague, Reference Davy and Lague2009; Lajeunesse and others, Reference Lajeunesse, Malverti and Charru2010). The commonly used empirical relationships for sediment transport are a function of the shear stress generated by fluid flow near the bed, τ, generally reported as the dimensionless Shields number, τ* = τ/[(ρ g − ρ w)gD] where ρ g, ρ w , g, and D are the particle density, fluid density, gravity, and particle diameter, respectively. The onset of sediment transport is classically associated with a critical Shields number, τ*c, that is required to overcome grain friction (Shields, Reference Shields1936). Meyer-Peter and Müller (Reference Meyer-Peter and Müller1948) is a well-supported sediment transport relationship describing bed-load transport in a turbulent flow:

(2) $$ Q_{\rm s} = 8 \max \lpar 0\comma \tau^{\ast} - \tau_{\rm c}^{\ast} \rpar ^{3/2} W \sqrt{{\rho_{\rm g} - \rho_{\rm w} \over \rho_{\rm w}}gD^3} $$

where W is the channel width. Importantly, τ*c is highly material dependent, and, in particular, increases as particle sizes become small enough that cohesion forces become significant. We note that the sediment–flux relationship presented above may fall short for very clay-rich or multimodal beds (e.g., Wilcock, Reference Wilcock1998; Houssais and Lajeunesse, Reference Houssais and Lajeunesse2012). The fluid shear stress τ along the hydraulic perimeter can be inferred through the Darcy–Weisbach formula (e.g., Henderson, Reference Henderson1966; Walder and Fowler, Reference Walder and Fowler1994; Ng, Reference Ng2000a; Creyts and Clarke, Reference Creyts and Clarke2010; Carter and others, Reference Carter, Fricker and Siegfried2017), $\tau = 0.125 f^{\prime} \rho _{\rm w} \left(QS^{-1} \right)^2$ , where f′ is a dimensionless friction factor.

As for hard-bed channel models the water flux Q along the channel length axis s can be described by a turbulent flow law (Hewitt, Reference Hewitt2011):

(3) $$ Q = \sqrt{F^{-1} S^{8/3} \left(\psi - {\partial P_{\rm c} \over \partial s}\right)}\comma $$

where P c is the averaged effective pressure in the channel, and ψ is the topographically constrained pressure gradient:

(4) $$ \psi = -\rho_{\rm i} g {\partial \lpar b + H \rpar \over \partial s} - \lpar \rho_{\rm w} - \rho_{\rm i} \rpar g {\partial b\over \partial s}. $$

Here, b is the bed topography and H is the ice thickness (e.g., Hewitt, Reference Hewitt2011), and F is a function of conduit geometry and the Manning friction coefficient n′: F = ρ w gn′[2(π + 2)2/π]2/3. The water-flow law (Eqn 3) is typically rearranged to solve for along-channel change in effective pressure P c (e.g., Ng, Reference Ng2000b; Kingslake and Ng, Reference Kingslake and Ng2013; Carter and others, Reference Carter, Fricker and Siegfried2017):

(5) $$ {\partial P_{\rm c} \over \partial s} = \psi - {FQ^2 \over S^{8/3}}. $$

The water flux Q in the above equation is usually found from an equation of water conservation, often by assuming water incompressibility, negligible change in water storage along the flow path, and negligible sediment fluxes relative to the water discharge. The change in water flux downstream is given by a local source term, ${\dot m}$ (e.g., Schoof, Reference Schoof2010):

(6) $$ {\partial Q \over \partial s} = {\dot m}. $$

The source term ${\dot m}$ can be comprised several contributions, including influx from the surrounding bed through groundwater seepage, inflow along the ice–bed interface, and input from englacial storage. In the following, we apply a numerical model of granular deformation to explore the conditional sediment stability a channel, and tweak the above formulation for channel size (Eqn 1) to capture the essential behavior associated with Mohr–Coulomb plasticity of the subglacial sediment.

2.2. Grain-scale sediment modeling

2.2.1. Model description

We use a discrete element method (DEM, e.g. Cundall and Strack, Reference Cundall and Strack1979; Damsgaard and others, Reference Damsgaard2013) in order to resolve the sediment mechanics in the bed surrounding an idealized subglacial channel. The DEM is a Lagrangian-type numerical approach of multi-body classical mechanics. Newton's Second Law of motion is explicitly integrated to find translational and rotational acceleration, velocity, and position for each grain through time. For a grain i in contact with jN other grains, the sum of forces consists of gravitational pull ( f g), grain-to-grain elastic-frictional contact forces ( f n and f t), and the fluid-pressure force ( f f),

(7) $$ {\partial^2 {\bf x}^i \over \partial t^2} m^i = {\bi f}_{\rm g}^i + \sum_{\,j \in N} \left({\bi f}_{\rm n}^{i\comma j} + {\bi f}_{{\rm t}}^{i\comma j} \right)+ {\bi f}_{{\rm f}}^i\comma $$

where x is the grain center position, m is the grain mass, and t is the time. A similar equation conserves angular momentum:

(8) $$ {\partial^2 {\bf \Omega}^i\over \partial t^2} I^i = \sum_{\,j \in N} \left(-\left(r^i {\bi n}^{i\comma j} \times \ {\bi f}_{\rm t}^{i\comma j} \right)\right)$$

where n is the grain-contact normal vector, Ω is the angular particle position, and I is the grain rotational inertia.

The forces from grain interactions ( f n and f t) are determined by a Hookean (linear elastic) rheology. The contact-normal interaction force is found from the contact-normal component of the inter-grain overlap distance vector δ :

(9) $$ {\bi f}_{\rm n}^{ij} = -k_{\rm n} {\bi \delta}_{\rm n}^{i\comma j}. $$

The tangential (contact-parallel) interaction force is similarly found from the contact-tangential component of the overlap distance vector, but is in magnitude limited by the Coulomb frictional coefficient μ:

(10) $$ {\bi f}_{\rm t}^{ij} = - \min \left[k_{\rm t} {\bi \delta}_{\rm t}^{i\comma j}\comma \,\,\mu \Vert {\bi f}_{\rm n}^{i\comma j}\Vert \right]{{\bi \delta}_{\rm t}^{i\comma j}\over \Vert {\bi \delta}_{\rm t}^{i\comma j}\Vert }. $$

The contact-normal and tangential travel vectors ( δ n and δ t) are incrementally adjusted over the duration of a grain-to-grain interaction, and continuously corrected in the case of contact rotation. If Coulomb failure occurs at the contact, the tangential travel is readjusted to correspond to a length consistent with Coulomb's condition (|| δ t|| = μ|| f n||k −1 t, e.g. Luding, Reference Luding2008; Radjaï and Dubois, Reference Radjaï and Dubois2011). The contact stiffnesses for the elastic interactions (k n and k t) are, contrary to our previous study using this model (e.g., Damsgaard and others, Reference Damsgaard2015), determined by scaling against a macroscopic elasticity:

(11) $$ k_{\rm n,t} = {{E \pi \lpar r_i + r_j \rpar} \over 2}\comma $$

where E is the Young's modulus, and r i and r j are the radii of two interacting grains. This approach makes the bulk elastoplastic behavior independent of the chosen grain size (Ergenzinger and others, Reference Ergenzinger, Seifried and Eberhard2011; Obermayr and others, Reference Obermayr, Dressler, Vrettos and Eberhard2013).

The fluid-pressure force on each grain is determined by the local gradient of the water-pressure field (p f), as well as buoyant uplift from the weight of displaced fluid (Goren and others, Reference Goren, Aharonov, Sparks and Toussaint2011; Damsgaard and others, Reference Damsgaard2015):

(12) $$ {\bi f}_{\rm f}^i = -V^i \nabla p_{\rm f} - \rho_{\rm f} V^i {\bi g} $$

where V is the grain volume, ρ f is the fluid density, and g is the gravitational acceleration vector. We ignore other and weaker fluid-interaction forces (Stokes drag, Saffman force, Magnus force, virtual mass force) (e.g., Zhu and others, Reference Zhu, Zhou, Yang and Yu2007; Zhou and others, Reference Zhou, Kuang, Chu and Yu2010), which become important with faster fluid flow and associated larger Reynolds numbers.

Once the sum of forces (right-hand side of Eqn 7) and torques (right-hand side of Eqn 8) have been determined, we perform explicit and third-order temporal integrations with Taylor expansions in order to determine the new kinematic state (e.g., Kruggel-Emden and others, Reference Kruggel-Emden, Sturm, Wirtz and Scherer2008). The maximum admissible time step is tied to the propagation of seismic (elastic) waves through the granular assemblage, and is determined by the density, size, and elastic stiffnesses in the granular system (e.g., Radjaï and Dubois, Reference Radjaï and Dubois2011):

(13) $$ \Delta t = {\epsilon \over \sqrt{\lpar {\max\lpar k_{\rm n}\comma k_{\rm t}\rpar }/{\min\lpar m\rpar }\rpar }} $$

with a constant safety factor (ε = 0.07).

2.2.2. Meltwater in sediment and channel

We describe pore-water pressure in the sediment by a time-dependent diffusion equation with a forcing term related to porosity change. The rate of pressure diffusion scales according to Darcy's law in heterogeneous materials:

(14) $$ {\partial p_{\rm f} \over \partial t} = {1 \over \beta_{\rm f} \phi \mu_{\rm f}} \nabla \cdot \left(k \nabla p_{\rm f} \right)- {1 \over \beta_{\rm f} \lpar 1-\phi\rpar } \left({\partial \phi \over \partial t} + {\bar{\bi v}} \cdot \nabla \phi \right)$$

where p f is the fluid-pressure deviation from the hydrostatic pressure, β f is the adiabatic fluid compressibility, and μ f is the fluid dynamic viscosity. The locally averaged grain velocity is denoted ${\bar {\bi v}}$ . The local permeability is denoted k and is determined by empirical relations as a function of local porosity ϕ (Goren and others, Reference Goren, Aharonov, Sparks and Toussaint2011; Damsgaard and others, Reference Damsgaard2015). The second term on the right-hand side forces a response in water pressure as local porosity changes, and is corrected for advection of porosity. For the sake of simplicity, the fluid-flow model considers the whole domain as porous, with a very high porosity and permeability in the channel (see Fig. 2). The above equation is not adequate for simulating turbulent water flow and sediment in bedload inside the channel.

Fig. 2. Overview of the granular assemblage of 58 000 particles and the discretization of the fluid grid. We assume symmetry around the channel center (+x) and limit our simulation domain to one of the sides. The model domains of the two phases are superimposed during the simulations.

We calculate the fluid-phase dynamics (Eqn 14) on an Eulerian and regular grid, superimposed over the granular assemblage. The fluid is fully two-way coupled to the granular phase. The fluid forces the sediment through spatial pressure gradients different from the hydrostatic pressure distribution, and the sediment forces the fluid phase by local changes in porosity (e.g., McNamara and others, Reference McNamara, Flekkøy and Måløy2000; Goren and others, Reference Goren, Aharonov, Sparks and Toussaint2011; Damsgaard and others, Reference Damsgaard2015). Contrary to our previous two-phase modeling studies (Damsgaard and others, Reference Damsgaard2015, Reference Damsgaard2016), the fluid grid is adjusted in size to the spatial domain of the granular phase during the simulations in order to correctly resolve the dynamics during volumetric changes, without requiring a constant-pressure boundary condition at the top.

2.3. Design of numerical experiments

The purpose of our numerical experiments is to test the hypothesis that sediment dynamics plays an important role in shaping subglacial channels. We design our experiments to simulate the sediment surrounding a channel cavity, which is stressed by a virtual ice–bed interface modeled as a horizontal wall with a pre-defined normal stress. The model is three dimensional in order to correctly resolve grain rotation and interlocking, but the simulation domain is shortened in the along-flow dimension (y) in order to reduce the computational overhead.

We perform dry (water-free) and wet (water-saturated) simulations, with the hypothesis that the presence of water without large pressure gradients does not influence sediment stability, as long as the effective normal stress is equal. The dry simulations (Fig. 3) are designed to inform about steady-state channel geometries under a constant normal stress on the top boundary, which spans the ice–sediment interface and the channel conduit. We perform these simulations without water in order to illustrate pure granular mechanics without transient hardening or softening associated with a pore-fluid pressurization (e.g., Iverson and others, Reference Iverson, Hooyer and Baker1998; Moore and Iverson, Reference Moore and Iverson2002; Goren and others, Reference Goren, Aharonov, Sparks and Toussaint2011; Damsgaard and others, Reference Damsgaard2015). The omission of fluid in these simulations additionally serves to reduce the required computational time.

Fig. 3. Boundary conditions and initial state in the presented simulations for the granular phase in the dry experiments. The frictionless lateral boundaries imply no boundary-parallel movement.

For the second set of experiments, we employ wet (water-saturated) simulations (Fig. 4) to investigate sediment stability under various defined water-pressure gradients. The water pressures are kept constant at the channel and at the far-away lateral boundary (x = 0 m in Fig. 2), and water pressures are inside the sediment initialized according to a linear interpolation between these Dirichlet boundary conditions. After the initialization procedure, described below, the porosity at the top of the sediment is naturally larger because of the flat ice–bed interface. Additionally, the porosity generally tends to slightly decrease with depth as lithostatic pressure and packing density increases. Local deviations in water pressure from the hydrostatic pressure distribution directly affect the stress balance inside the sediment and at the upper boundary.

Fig. 4. Boundary conditions and initial state in the presented simulations for the granular and fluid phases in the wet experiments.

The granular assemblage is prepared by letting 1000 uniformly distributed and spherical particles settle in a small cubic volume under the influence of gravity. Afterwards the assemblage is duplicated and repositioned 58 times in space to construct the desired geometry for a total of 58 000 grains. We assume that the channel is horizontally symmetrical around its center and simulate half of the channel width. Before proceeding with our main experiments we perform a relaxation step where we allow the grains to settle and adjust their arrangement to the new geometry. Figures 3, 4 show the boundary conditions, as well as the geometry of the simulated sediment at the beginning of the experiments. Table 1 lists the values of the relevant geometrical, mechanical, and temporal parameters.

Table 1. Simulation parameters and their values for the granular experiments

* Simulation time and time-dependent parameters are listed in scaled model time, see Damsgaard and others (Reference Damsgaard2015, Reference Damsgaard2016).

The governing equation for fluid dynamics (Eqn 14) is in its standard form singular for areas that do not contain grains (ϕ = 1), such as the channel cavity. We therefore impose an upper limit on porosity of ϕ = 0.9 in Eqn (14) and the other fluid-related equations in order to allow our fluid-phase formulation to work for the channel conduit. Nonetheless, the strong non-linearity of the permeability relation (k = 3.5 × 10−13m2 ϕ 3(1 − ϕ)−2) causes much faster diffusion of fluid pressures in the conduit than internally in the sediment, consistent with our expectations of how the hydraulic system should behave. However, we are for the purposes of this study mainly interested in the mechanical behavior of the sediment surrounding the channel. Similarly, we do not resolve water flow along the channel length axis, but focus our experimental analysis of grain-fluid interaction inside of the sediment with the channel cavity acting more as a boundary condition.

We choose a relatively low value for Young's modulus in Eqn (11) in order to increase the computational time step (Eqn 13 and Table 1), but the resultant elastic softening has negligible effect on porosity and mechanical behavior. Compressive strain on grain contacts is in our model due to the softening expected to be a hundred times greater than theoretical values for quartz, but the actual strains remain very small (2 × 10−3% of the grain radius under a compressive stress of 10 kPa). Furthermore, we chose a relatively large grain size and narrow grain-size distribution to increase computational efficiency during the grain-to-grain contact searches. At the same time, the hydraulic properties are decoupled from sediment size to allow for low-permeability simulations with relatively large grains (Goren and others, Reference Goren, Aharonov, Sparks and Toussaint2011; Damsgaard and others, Reference Damsgaard2015). We do note that the simplification in grain size, shape, and morphology will influence granular deformation patterns and sediment frictional properties (e.g., Iverson and others, Reference Iverson, Hooyer and Baker1998; Damsgaard and others, Reference Damsgaard2013). We also omit tensile interactions between grains giving cohesionless behavior. However, apparent cohesion observed in laboratory shear tests on real geological materials often vanishes at very low normal stresses (e.g., Schellart, Reference Schellart2000), making grain-tension less important at free sediment surfaces.

We assume that channels in subglacial beds inherently form due to flow instabilities by differential erosion, similar to the processes driving distributed flow to R-channel drainage on hard beds. Intrinsic channelization during fluid flow has been demonstrated in physical dam-breach experiments (e.g., Walder and others, Reference Walder, Iverson, Godt, Logan and Solovitz2015) and in smaller experimental studies (e.g., Catania and Paola, Reference Catania and Paola2001; Mahadevan and others, Reference Mahadevan, Orpe, Kudrolli and Mahadevan2012; Kudrolli and Clotet, Reference Kudrolli and Clotet2016; Métivier and others, Reference Métivier, Lajeunesse and Devauchelle2017). Channelization occurs when fluid flux is able to erode the porous medium through which it is flowing. In turn, this process reduces local friction against the fluid flux and further enhances flow localization. In this study, we investigate the mechanical state and deformation patterns around an already established channel conduit incised into the sedimentary bed, and are not able to include the process of channel formation itself.

3. RESULTS

We do not observe significant differences in steady-state channel geometry between the dry and wet simulations at equal effective normal stresses, which is expected as the fluid only affects the granular force balance in the presence of water-pressure gradients (Eqn 12). However, the wet simulations take three times longer to complete.

In the dry experiments, we observe that the effective normal stress on the channel flanks results in rapid failure, occurring over a few seconds (Fig. 5). Larger normal stresses result in more deformation. Consistent with sediment plasticity, the deformation stops when a new stress balance has been established. The final and steady-state geometries are shown in Figure 6, where larger normal stresses result in smaller channel cavities. We do expect slow creep in the sediment, especially if there are large fluctuations in water pressure and granular stresses (Pons and others, Reference Pons, Amon, Darnige, Crassous and Clément2015; Damsgaard and others, Reference Damsgaard2016) or strong water flow (Houssais and others, Reference Houssais, Ortiz, Durian and Jerolmack2015, Reference Houssais, Ortiz, Durian and Jerolmack2016). However, the associated creep is expected to occur with a significantly more non-linear rheology than previously used for soft-bed channels (e.g., Alley, Reference Alley1992; Walder and Fowler, Reference Walder and Fowler1994; Ng, Reference Ng2000b; Carter and others, Reference Carter, Fricker and Siegfried2017), and we do not detect measurable creep on the timescales considered here (Fig. 5).

Fig. 5. Averaged grain displacements through time for a range of effective normal stresses, N in the dry granular experiments.

Fig. 6. Channel geometries in the steady state at different effective normal stresses N in dry experiments. Grains are colored according to their cumulative vertical displacement and rendered in an orthogonal projection. The absence of water does not influence the steady-state geometry unless the sediment is flushed by strong water-pressure gradients.

We can describe the relationship between the imposed effective normal stress and observed maximum channel width reasonably well by a linear fit (Fig. 7):

(15) $$ W_{\rm max} = \lpar a N + b\rpar $$

with fitting parameter values a = −0.118 m kPa−1 and b = 4.60 m (associated standard deviations (std dev.) reported in Fig. 7). By assuming a simple geometry where the slope of the channel flanks is given by the sediment angle of repose θ (Fig. 8), we can approximate the limit to channel the cross-sectional area, S max as a function of the channel width:

(16) $$ S_{\rm max} = {1\over 4} W_{\rm max}^2 \tan \theta. $$

We constrain the channel size in the continuum channel model (Eqns 16) in order to capture the stress-based limits. We do not view extrapolation of the channel width below the tasted range (N = [2.5;40] kPa) as physically meaningful. Geometries with slopes beneath the angle of repose should be stable without load from the ice overburden under events with overpressurized water (N ≤ 0 kPa), allowing for the existence of larger subglacial cavities.

Fig. 7. Observed maximum channel width, W max, under a constant imposed value for effective normal stress, N. Data points are fitted using a non-linear least-squares fit, with fitted parameter values and their corresponding std dev. noted.

Fig. 8. Idealized schematic of effective stress magnitude and till yield strength in a cross-section around a subglacial channel incised in the sedimentary bed. Sediments beneath the channel floor are not stressed by ice weight on the flanks, and comprise a weak channel floor wedge (CFW).

The grain-scale stress balance (Fig. 9) shows that stresses from the ice load on channel flanks are transmitted predominantly downwards, and do not affect the stress balance in the sediment beneath the channel conduit. The sediment beneath the channel conduit is loaded exclusively by its own weight (lithostatic pressure), which increases linearly with depth. This behavior contrasts with viscous till models, that assume that the stress on the channel flanks contributes to till uplift at the bottom of the channel.

Fig. 9. Contact pressure on the grains $\left(\sum \vert\vert{\bi f}_{\rm n}\vert\vert 4 \pi r^2 \right)$ in the steady state with an effective stress of N = 7 kPa.

In our wet simulations, forced with water-pressure differences between channel and bed, the pressure-gradients create drag forces on the grains oriented towards the channel (Fig. 10). The presence of smaller water-pressure gradients cause minor sediment rearrangement, but deformation rapidly decays and the channel continues to be stable (Fig. 11, left and center). The channel geometry is not significantly different from the dry experiment under small water-pressure gradients (compare the simulation with N = 10 kPa in Fig. 6 with Fig. 11, left and center). However, when larger pressure gradients are prescribed the sediment flushes into the channel, causing rapid channel infill and closure (Fig. 11, right).

Fig. 10. Fluid-pressure forces (Eqn 12) on the individual sediment grains visualized by arrows. Top: Under low gradients in pore-water pressure (here dp f/dx = −0.1 kPa m−1), the fluid-pressure force primarily contributes weak buoyant uplift on the grains. Bottom: Larger pressure gradients (dp f /dx =  − 10 kPa m−1) destabilize the sediment through liquefaction at the channel floor and cause collapse of the sediment with associated chaotic interactions (e.g., large internal forces). Individual force vectors are not representative of the grain ensemble and should not be interpreted directly. Both simulation snapshots are from t = 1.5 s.

Fig. 11. Total per-grain spatial displacements at different times for three wet (water-saturated) simulations with different imposed water-pressure gradients and an effective stress of N = 10 kPa. The water-pressure gradients cause flow and drag forces toward the channel, and destabilize the conduit at higher values.

Figure 12 shows the results of an example implementation of the continuum channel model (Eqns 16, with limits to channel size imposed by Eqns 15 and 16). The solution is found iteratively with a Dirichlet boundary condition of P c = 0 Pa at the terminus. Gradients are estimated with upwind finite difference approximations. We assume that channel effective pressure (P c) equals the magnitude of the effective normal stress on the channel flanks (N). In real settings channel effective pressure will most commonly be slightly higher than normal stress on the flanks (P c > || N ||), as water pressure differences cause dominant water flow towards the channel.

Fig. 12. Example run of the soft-bed subglacial channel model outlined in Eqns (16, 15 and 16), with a square-root ice geometry and linearly increasing meltwater influx ( ${\dot m}$ ) towards +s. Top: Channel effective pressure (P c), ice-overburden pressure (P i), water-pressure (P w) and water flux (Q). Middle: Sediment flux increases non-linearly with water flux (Eqn 2). Bottom: Maximum (S max) and actual channel cross-sectional size (S), together with channel growth rate (dS/dt). This example is with a constant forcing and t = 2 days. The minimum channel size is set to S = 0.01 m2 for numerical considerations. In this example the sediment yield strength prevents channel existence except near the terminus where effective pressure is relatively low.

In this example, the water fluxes are increasing downstream due to the water supply ( ${\dot m}$ ). The increasing water flux generally increases the sediment-bedload flux (Q s) downstream. At s = 0 to 0.8 km the large effective pressures (P c) inhibit channel development due to sediment failure. Close to the margin (s = 0.8–1.0 km) the effective pressure is sufficiently low to allow for channel existence. The increase in channel cross-sectional size near the margin decreases the fluvial shear stress and causes a rapid drop in fluvial sediment flux Q s.

4. DISCUSSION

We have demonstrated that sediment yield strength can prevent collapse of channels. Because of sediment plasticity, channel erosion into subglacial sediments is not balanced by a continuous flux of viscous till deformation into the channel cavity, which was assumed in prior parameterizations (e.g., Fowler and Walder, Reference Fowler and Walder1993; Walder and Fowler, Reference Walder and Fowler1994; Ng, Reference Ng2000b). Instead, channel cross-sectional geometry is governed by fluvial erosive and depositional dynamics until channel size is limited by sediment yield strength.

Infilled fossil subglacial channels seen in field sections along the southern margin of the last Scandinavian Ice Sheet (Fig. 13) show similar geometry to the stable channel conduits observed in our granular experiments (e.g., Figs 6, 11), where the sediment angle of repose is a principal control on channel cross-sectional geometry. The channel sizes observed in the field (Fig. 13) are within the limits observed in our granular experiments (Fig. 7) and the simplified ice-marginal area of our continuum channel model (Fig. 12).

Fig. 13. Examples of infilled subglacial channels in the southern, marginal part of the Scandinavian Ice Sheet from the last glaciation (Piotrowski et al., in prep.). (a) Subglacial channel at Ebeltoft, Djursland, Denmark. The channel is found within a single till unit about 10 km inside the ice margin. It is flat-topped and filled with parallel-bedded outwash sand and gravel intercalated with layers of silt. Single outsized stones, possibly dropped from the channel roof are randomly dispersed in the outwash deposit. Along the channel bottom and on its left-hand side the infill material is deformed into attenuated folds, irregular detached sediment pods and flame structures. Since the flanks of the channel are below the angle of repose of sand, the deformation structures suggest syndepositional sediment intrusion into the channel driven by a pressure gradient oriented toward the channel axis. (b). Subglacial channel at Glaznoty, north-central Poland, about 25 km within the ice limit. The channel occurs at the interface between proglacial outwash deposits (below) and till (above). The channel is distinctly lens-shaped with an upward-convex top suggesting an R-channel incised upward into the ice at a late stage of formation. It is infilled with massive coarse-grained sand and gravel.

While Shoemaker (Reference Shoemaker1986) assumed that the bed would strengthen against failure under high effective stresses, we show that the magnitude of effective stress on the bed surrounding the channel puts an upper limit on channel dimensions. The exact relationship between effective normal stress and channel size (Fig. 7) will be material dependent, and the granular model applied here includes several simplifying assumptions related to grain size and grain shape. The grain angularity and size distribution in subglacial tills might result in slightly larger yield strengths and allow for the existence of larger channels. Overconsolidation of subglacial tills additionally contributes to shear strength (e.g. Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000a), which will increase stability and limits to channel size. However, the stress–size relationship presented here should be a reasonable approximation for a simple model. While our results are based on numerical experiments on simplistic granular materials, our interpretation stands that it is unlikely for soft-bed channels to exist at effective subglacial normal stresses more than ~100 kPa in magnitude (Fig. 7). We propose that the constrained relationship between channel geometry and effective normal stress makes it theoretically possible to asses water pressure and hydrological conditions under past ice sheets, e.g. when meltwater channels are identified from geomorphological interpretation of subaerial or subaqueous topography (e.g., Lesemann and others, Reference Lesemann, Piotrowski and Wysota2010; Greenwood and others, Reference Greenwood, Clason, Nyberg, Jakobsson and Holmlund2016; Bjarnadóttir and others, Reference Bjarnadóttir, Winsborrow and Andreassen2017), or from channel-sediment geometries in the glaciogeological sedimentary record (e.g., Piotrowski, Reference Piotrowski1999; Tylmann and others, Reference Tylmann, Piotrowski and Wysota2013).

Ice overburden stress does not impact the sediments beneath the channel floor directly; the compressive stress acting on channel-floor sediments is a result of their own weight as lithostatic pressure increases with depth (Fig. 9). Sediment dynamics at the floor of subglacial channels are then similar to those in subaerial rivers, although the drainage system arrangement may differ (Catania and Paola, Reference Catania and Paola2001; Métivier and others, Reference Métivier, Lajeunesse and Devauchelle2017). Due to the plastic rheology of sediment, erosion of the channel floor will be counteracted by an immediate sediment response when the channel is at its size limit, acting to reestablish stress balance in the bed surrounding the channel. If differential sediment advection toward the cavity bends the ice–bed interface downwards progressively, the resultant channel landform of erosive origin is likely to appear wider than the channel cavity itself (Boulton and Hindmarsh, Reference Boulton and Hindmarsh1987; Ng, Reference Ng2000b).

The evolution of subaerial river size is governed by feedbacks and stability limits (e.g. Métivier and others, Reference Métivier, Lajeunesse and Devauchelle2017). We conclude that there are two distinct feedbacks stabilizing subglacial channel size on soft beds: (1) If a channel becomes sufficiently efficient in evacuating subglacial meltwater, the pore-water pressure decreases and the effective normal stress on the surrounding areas of the bed increases. If this stress increase causes the channel size to exceed the stability limit, the channel spontaneously reduces in size which decreases the hydrological transport capacity which, in turn, can increase the water pressure and decrease the effective normal stress. (2) For a given water flux, an increase in channel cross-sectional size due to fluvial erosion decreases the water-driven shear stress on the channel bottom. The decrease in shear stress decreases sediment transport (Eqn 2), which decreases channel growth (decreasing ∂Q s/∂s in Eqn 1). These feedbacks may ultimately lead to stabilization of subglacial sliding and hydrology, tied to the plastic failure limit of the bed and the basal hydraulic transmissivity.

For the subglacial drainage model presented here, we consider channelized flow alone although it could be significantly improved by coupling it to a diffusive flow model accounting for sheet flow at the ice–bed interface, groundwater flow, and/or R-channel incision (e.g., Creyts and Schoof, Reference Creyts and Schoof2009; Hewitt, Reference Hewitt2011; Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013; Flowers, Reference Flowers2015). Furthermore, the parameterization for water flux in the channel can be improved by including both laminar and turbulent descriptions, dependent on the Reynolds number of the water flow (e.g., van der Wel and others, Reference van der Wel, Christoffersen and Bougamont2013).

With increasing subglacial discharge, we expect a hydraulic transition from distributed drainage to sediment-incised channels (e.g., Mahadevan and others, Reference Mahadevan, Orpe, Kudrolli and Mahadevan2012; Kudrolli and Clotet, Reference Kudrolli and Clotet2016). However, yield failure of the channel flanks imposes a limit to their cross-sectional area and hydraulic transport capacity. If the hydraulic transmissivity becomes insufficient against the water fluxes, we either expect R-channel incision into the ice base, or the formation of multiple parallel drainage channels incised into the sediment. If subglacial soft-bed channels were able to grow without size constraints, it would be thermodynamically more efficient to gather drainage in fewer and larger channels.

An arrangement of closely spaced channel drainage is observed at Glaznoty, north-central Poland, where the channel in Figure 13 is one of a series of small, similar-sized channels occurring at the paleo-ice–bed interface. Parallel sediment landforms (glacial curvilineations) have been observed from under the palao-Scandinavian Ice Sheet (Fig. 1, Lesemann and others, Reference Lesemann, Piotrowski and Wysota2010, Reference Lesemann, Piotrowski and Wysota2014), which may be governed by the channel-size limits described here. Furthermore, radar-echo soundings of the sedimentary bed of Thwaites Glacier, West Antarctica have been interpreted to reflect hydraulic transitions between few to many parallel and closely spaced channelized drainage elements incised into the sedimentary bed (Schroeder and others, Reference Schroeder, Blankenship and Young2013). We do note that the landforms and channel-drainage elements observed in the glacial curvilineations are generally larger than what is predicted from our limits related to plastic failure (Fig. 7). We believe that this discrepancy is mainly related to the fact that the numerical material with spherical and smooth grains is mechanically weaker than subglacial tills with elongated and angular grains.

Liquefaction at the earth surface and in subaqueous environments is known to be initiated by overpressurization in the pore space, effectively reducing the compressive stress to low or even negative values (e.g., Zhang and Campbell, Reference Zhang and Campbell1992; Terzaghi and others, Reference Terzaghi, Peck and Mesri1996; Xu and Yu, Reference Xu and Yu1997; Mitchell and Soga, Reference Mitchell and Soga2005). In situ measurements of subglacial water pressure indicate that water pressures are highly variable through time (e.g., Hooke, Reference Hooke1984; Engelhardt and Kamb, Reference Engelhardt and Kamb1997; Hooke and others, Reference Hooke, Hanson, Iverson, Jansson and Fischer1997; Bartholomaus and others, Reference Bartholomaus, Anderson and Anderson2008; Andrews and others, Reference Andrews2014; Schoof and others, Reference Schoof, Rada, Wilson, Flowers and Haseloff2014). We propose that events of liquefaction in subglacial channels may be common when water-pressure in the channel rapidly decreases, and this process may be able to make significant volumes of weak sediment from the channel-floor wedge available for fluvial transport, especially in sedimentary beds with low permeability that require long timescales to respond to changes in interfacial water pressure.

As previously discussed, we are not able to include the effect of subglacial deformation due to ice movement along the channel length, which might be important for long-term channel stability. We do note that subglacial shear could be included by setting the y boundaries to be periodic (e.g., Damsgaard and others, Reference Damsgaard2013), but the spatial domain length along y in the current setup is too short to allow for shearing without geometrical instabilities. Sediment advection associated with shear deformation causes frictional heating and granular diffusion (e.g., Hooyer and Iverson, Reference Hooyer and Iverson2000; Utter and Behringer, Reference Utter and Behringer2004); these process are likely to drive channel closure with a rate proportional to the shear-strain rate in the sediment. We also assume that the ice–bed interface remains flat and rigid over time, while Ng (Reference Ng2000b) demonstrated that differential ice and till advection toward the channel conduit bends the interface over longer timescales. Our model approach could be improved by simulating a dynamically evolving ice-bed interficial geometry. However, we believe that inclinations in the ice–bed interface are unlikely to fundamentally alter the principal stresses in the surrounding sediment, and assume that the ice is responding elastically over the timescales investigated in the experiments (<1 min).

5. CONCLUSIONS

Current relationships for subglacial channel dynamics incised into sedimentary beds assume linear to mildly non-linear viscous relationships for till rheology, which results in continuous sediment flux toward the channel balancing erosion by water flow. However, sediments are known to be nearly perfect plastic with a yield strength dependent on the confining stress.

We have coupled two separate models to gain a multi-scale understanding between sediment deformation and subglacial channel stability. Our granular model informs about sediment stability under different effective normal stresses and water-pressure forcings. We observe that the channel conduit size is strongly limited by the magnitude of effective normal stress on the channel flanks, and that creep closure is negligible. The compressive stresses from the ice–bed interface on the channel flanks are oriented subvertically instead of being directed towards the channel floor. The channel-flooring sediments are only compacted by their own weight. Strong water-pressure differences between the channel and its surrounding parts can cause horizontal infilling by sediment movement along the ice–bed interface.

We use the results from our granular simulations to include effects of sediment plasticity in a continuum model of soft-bed subglacial channels. The channel size is limited by the yield strength of the sediment, which in turn depends on effective normal stress on the channel flanks. The size limit implies that multiple closely spaced channels are needed for transporting large amounts of water, which corresponds to geophysical observations under contemporary ice sheets and geomorphological signatures from previously glaciated areas. The presented continuum model for channelized drainage, derived from our suite of numerical simulations, increases the realism of hydrology models for ice sheets and glaciers residing on soft beds.

ACKNOWLEDGMENTS

A. Damsgaard was supported by a postdoctoral fellowship by the Cecil H. and Ida M. Green Foundation and NSF-PLR-1543396. A. Damsgaard acknowledges the NVIDIA Corporation for a hardware grant and the Extreme Science and Engineering Discovery Environment (XSEDE) for a high-performance computing startup allocation. In preparation for this paper, A. Damsgaard has benefited from discussions with C. W. Elsworth, L. Goren, I. J. Hewitt, S. P. Carter, A. Cabrales-Vargas, and V. Tsai. We thank Jonathan Kingslake and anonymous reviewer for constructive comments, which improved the manuscript and Neil Glasser for handling the editorial process.

APPENDIX A. SOURCE CODE AVAILABILITY

The source code for the grain-fluid model is available at https://github.com/anders-dc/sphere, where the script channel-wet.py can be used as a template for model runs. An example implementation of the subglacial hydrology model built on channelization dynamics is available from https://github.com/anders-dc/granular-channel-hydro/blob/master/1d-channel.py.

Footnotes

*

Present address: Geophysical Fluid Dynamics Laboratory, Princeton University, NJ, USA.

References

REFERENCES

Alley, KE, Scambos, TA, Siegfried, MR and Fricker, HA (2016) Impacts of warm water on Antarctic ice shelf stability through basal channel formation. Nat. Geosci., 9(4), 290293 Google Scholar
Alley, RB (1989) Water-pressure coupling of sliding and bed deformation: I. Water system. J. Glaciol., 35(119), 108118 Google Scholar
Alley, RB (1992) How can low-pressure channels and deforming tills coexist subglacially? J. Glaciol., 38, 200207 CrossRefGoogle Scholar
Alley, RB, Blankenship, DD, Bentley, CR and Rooney, ST (1986) Deformation of till beneath ice stream B, West Antarctica. Nature, 322, 5759 Google Scholar
Alley, RB, Lawson, DE, Larson, GJ, Evenson, EB and Baker, GS (2003) Stabilizing feedbacks in glacier-bed erosion. Nature, 424(6950), 758760 Google Scholar
Anandakrishnan, S, Blankenship, DD, Alley, RB and Stoffa, PL (1998) Influence of subglacial geology on the position of a West Antarctic ice stream from seismic observations. Nature, 394(6688), 6265 Google Scholar
Andrews, LC and 7 others (2014) Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet. Nature, 514, 8083 Google Scholar
Bartholomaus, TC, Anderson, RS and Anderson, SP (2008) Response of glacier basal motion to transient water storage. Nat. Geosci., 1(1), 3337 CrossRefGoogle Scholar
Bartholomew, I and 5 others (2010) Seasonal evolution of subglacial drainage and acceleration in a Greenland outlet glacier. Nat. Geosci., 3(6), 408411 Google Scholar
Bindschadler, R (1983) The importance of pressurized subglacial water in separation and sliding at the glacier bed. J. Glaciol., 29(101), 319 Google Scholar
Bjarnadóttir, LR, Winsborrow, MCM and Andreassen, K (2017) Large subglacial meltwater features in the central Barents Sea. Geology, 45(2), 159162 CrossRefGoogle Scholar
Bougamont, M, Price, S, Christoffersen, P and Payne, AJ (2011) Dynamic patterns of ice stream flow in a 3-D higher-order ice sheet model with plastic bed and simplified hydrology. J. Geophys. Res. Earth Surf., 116(F04018)Google Scholar
Bougamont, M and 5 others (2015) Reactivation of Kamb Ice Stream tributaries triggers century-scale reorganization of Siple Coast ice flow in West Antarctica. Geophys. Res. Lett., 42(20), 84718480 Google Scholar
Boulton, GS and Hindmarsh, RCA (1987) Sediment deformation beneath glaciers: rheology and geological consequences. J. Geophys. Res., 92(B9), 90599082 Google Scholar
Boulton, GS, Hagdorn, M, Maillot, PB and Zatsepin, S (2009) Drainage beneath ice sheets: groundwater–channel coupling, and the origin of esker systems from former ice sheets. Quat. Sci. Rev., 28(7–8), 621638 Google Scholar
Brinkerhoff, DJ, Meyer, CR, Bueler, E, Truffer, M and Bartholomaus, TC (2016) Inversion of a glacier hydrology model. Ann. Glaciol., 57, 112 Google Scholar
Budd, WF, Keage, PL and Blundy, NA (1979) Empirical studies of ice sliding. J. Glaciol., 23, 157170 Google Scholar
Carter, SP, Fricker, HA and Siegfried, MR (2017) Antarctic subglacial lakes drain through sediment-floored canals: theory and model testing on real and idealized domains. Cryosphere, 11(1), 381405 CrossRefGoogle Scholar
Catania, G and Paola, C (2001) Braiding under glass. Geology, 29(3), 259 2.0.CO;2>CrossRefGoogle Scholar
Christianson, K and 7 others (2014) Dilatant till facilitates ice-stream flow in northeast Greenland. Earth Planet. Sci. Lett., 401, 5769 CrossRefGoogle Scholar
Clarke, GKC (2005) Subglacial processes. Annu. Rev. Earth Planet. Sci., 33, 247276 CrossRefGoogle Scholar
Creyts, TT and Clarke, GKC (2010) Hydraulics of subglacial supercooling: theory and simulations for clear water flows. J. Geophys. Res., 115(F3)Google Scholar
Creyts, TT and Schoof, CG (2009) Drainage through subglacial water sheets. J. Geophys. Res. Earth Surf., 114(F4)Google Scholar
Cundall, PA and Strack, ODL (1979) A discrete numerical model for granular assemblies. Géotechnique, 29, 4765 Google Scholar
Damsgaard, A and 5 others (2013) Discrete element modeling of subglacial sediment deformation. J. Geophys. Res. Earth Surf., 118, 22302242 Google Scholar
Damsgaard, A and 5 others (2015) A new methodology to simulate subglacial deformation of water-saturated granular material. Cryosphere, 9, 21832200 CrossRefGoogle Scholar
Damsgaard, A and 6 others (2016) Ice flow dynamics forced by water pressure variations in subglacial granular beds. Geophys. Res. Lett., 43 Google Scholar
Davy, P and Lague, D (2009) Fluvial erosion/transport equation of landscape evolution models revisited. J. Geophys. Res., 114(F3)Google Scholar
de Fleurian, B and 6 others (2014) A double continuum hydrological model for glacier applications. Cryosphere, 8(1), 137153 CrossRefGoogle Scholar
Dow, CF, Werder, MA, Nowicki, S and Walker, RT (2016) Modeling Antarctic subglacial lake filling and drainage cycles. Cryosphere, 10(4), 13811393 CrossRefGoogle Scholar
Drews, R and 9 others (2017) Actively evolving subglacial conduits and eskers initiate ice shelf channels at an Antarctic grounding line. Nat. Commun., 8, 15228 Google Scholar
Einstein, HA (1950) The bed-load function for sediment transportation in open channel flows, vol. 1026. U.S. Department of Agricure, Washington, DC Google Scholar
Elsworth, CW and Suckale, J (2016) Rapid ice flow rearrangement induced by subglacial drainage in West Antarctica. Geophys. Res. Lett., 43(22), 1169711707 CrossRefGoogle Scholar
Engelhardt, H and Kamb, B (1997) Basal hydraulic system of a West Antarctic ice stream: Constraints from borehole observations. J. Glaciol., 43(144), 207230 Google Scholar
Engelhardt, H, Humphrey, N, Kamb, B and Fahnestock, M (1990) Physical conditions at the base of a fast moving Antarctic ice stream. Science, 248(4951), 5759 Google Scholar
Ergenzinger, C, Seifried, R and Eberhard, P (2011) A discrete element model to describe failure of strong rock in uniaxial compression. Granul. Matter, 13, 341364 CrossRefGoogle Scholar
Flowers, GE (2015) Modelling water flow under glaciers and ice sheets. Proc. R. Soc. Lond. A, R. Soc., 471, 20140907Google Scholar
Fowler, A and Walder, J (1993) Creep closure of channels in deforming subglacial till. Proc. R. Soc. Lond. A, R. Soc., 441, 1731 Google Scholar
Fowler, AC (1987) Sliding with cavity formation. J. Glaciol., 33(115), 255267 CrossRefGoogle Scholar
Fricker, HA, Siegfried, MR, Carter, SP and Scambos, TA (2016) A decade of progress in observing and modelling Antarctic subglacial water systems. Phil. Trans. R. Soc. A, 374(2059), 20140294Google Scholar
GDR-MiDi (2004) On dense granular flows. Euro. Phys. J. E, 14, 341365 Google Scholar
Gillet-Chaulet, F and 6 others (2016) Assimilation of surface velocities acquired between 1996 and 2010 to constrain the form of the basal friction law under Pine Island Glacier. Geophys. Res. Lett., 43(19), 1031110321 CrossRefGoogle Scholar
Gimbert, F, Tsai, VC, Amundson, JM, Bartholomaus, TC and Walter, JI (2016) Subseasonal changes observed in subglacial channel pressure, size, and sediment transport. Geophys. Res. Lett., 43(8), 37863794 CrossRefGoogle Scholar
Goldberg, DN, Schoof, C and Sergienko, OV (2014) Stick-slip motion of an Antarctic Ice Stream: The effects of viscoelasticity. J. Geophys. Res. Earth Surf., 119(7), 15641580 CrossRefGoogle Scholar
Goren, L, Aharonov, E, Sparks, D and Toussaint, R (2011) The mechanical coupling of fluid-filled granular material under shear. Pure Appl. Geophys., 168(12), 22892323 Google Scholar
Greenwood, SL, Clason, CC, Nyberg, J, Jakobsson, M and Holmlund, P (2016) The Bothnian Sea ice stream: early Holocene retreat dynamics of the southcentral Fennoscandian Ice Sheet. Boreas, 46(2), 346362 Google Scholar
Henderson, FM (1966) Open Channel Flow. MacMillan, New York Google Scholar
Hewitt, IJ (2011) Modelling distributed and channelized subglacial drainage: the spacing of channels. J. Glaciol., 57(202), 302314 CrossRefGoogle Scholar
Hewitt, IJ (2013) Seasonal changes in ice sheet motion due to melt water lubrication. Earth Planet. Sci. Lett., 371-372, 1625 Google Scholar
Hooke, RLeB (1984) On the role of mechanical energy in maintaining subglacial water conduits at atmospheric pressure. J. Glaciol., 30, 180187 Google Scholar
Hooke, RLeB, Hanson, B, Iverson, NR, Jansson, P and Fischer, UH (1997) Rheology of till beneath Storglaciären, Sweden. J. Glaciol., 43(143), 172179 Google Scholar
Hooyer, TS and Iverson, NR (2000) Diffusive mixing between shearing granular layers: constraints on bed deformation from till contacts. J. Glaciol., 46(155), 641651 Google Scholar
Horgan, HJ and 7 others (2013) Estuaries beneath ice sheets. Geology, 41(11), 11591162 Google Scholar
Houssais, M and Jerolmack, DJ (2017) Toward a unifying constitutive relation for sediment transport across environments. Geomorphology, 277, 251264 Google Scholar
Houssais, M and Lajeunesse, E (2012) Bedload transport of a bimodal sediment bed. J. Geophys. Res.: Earth Surf., 117(F4), F04015 Google Scholar
Houssais, M, Ortiz, CP, Durian, DJ and Jerolmack, DJ (2015) Onset of sediment transport is a continuous transition driven by fluid shear and granular creep. Nat. Commun., 6 Google Scholar
Houssais, M, Ortiz, CP, Durian, DJ and Jerolmack, DJ (2016) Rheology of sediment transported by a laminar flow. Phys. Rev. E, 94(6)Google Scholar
Hubbard, B and Nienow, P (1997) Alpine subglacial hydrology. Quat. Sci. Rev., 16, 939955 Google Scholar
Iken, A (1981) The effect of the subglacial water pressure on the sliding velocity of a glacier in an idealized numerical model. J. Glaciol., 27(97), 407421 Google Scholar
Iverson, NR and Iverson, RM (2001) Distributed shear of subglacial till due to Coulomb slip. J. Glaciol., 47(158), 481488 Google Scholar
Iverson, NR, Hooyer, TS and Baker, RW (1998) Ring-shear studies of till deformation: Coulomb-plastic behavior and distributed strain in glacier beds. J. Glaciol., 148, 634642 Google Scholar
Iverson, NR and 7 others (2007) Soft-bed experiments beneath Engabreen, Norway: regelation infiltration, basal slip and bed deformation. J. Glaciol., 53(182), 323340 Google Scholar
Jenkins, A (2011) Convection-Driven Melting near the Grounding Lines of Ice Shelves and Tidewater Glaciers. J. Phys. Oceanogr., 41(12), 22792294 Google Scholar
Jørgensen, F and Sandersen, PBE (2006) Buried and open tunnel valleys in Denmark—erosion beneath multiple ice sheets. Quat. Sci. Rev., 25(11–12), 13391363 Google Scholar
Kamb, B (1991) Rheological nonlinearity and flow instability in the deforming bed mechanism of ice stream motion. J. Geophys. Res., 96(B10), 1658516595 CrossRefGoogle Scholar
Kamb, B (2001) The West Antarctic Ice Sheet: Behavior and Environment, in Alley RB and Bindshadler RA (eds) Basal Zone of the West Antarctic Ice Streams and its Role in Lubrication of Their Rapid Motion, pp. 157–199. American Geophysical Union, Washington DCGoogle Scholar
Kavanaugh, JL and Clarke, GKC (2006) Discrimination of the flow law for subglacial sediment using in situ measurements and an interpretation model. J. Geophys. Res. Earth Surf., 111(F1)Google Scholar
Kehew, AE, Piotrowski, JA and Jørgensen, F (2012) Tunnel valleys: Concepts and controversies — A review. Earth-Sci. Rev., 113, 3358 Google Scholar
King, EC, Hindmarsh, RCA and Stokes, CR (2009) Formation of mega-scale glacial lineations observed beneath a West Antarctic ice stream. Nat. Geosci., 2(8), 585588 Google Scholar
Kingslake, J and Ng, F (2013) Modelling the coupling of flood discharge with glacier flow during jkulhlaups. Ann. Glaciol., 54(63), 2531 Google Scholar
Kruggel-Emden, H, Sturm, M, Wirtz, S and Scherer, V (2008) Selection of an appropriate time integration scheme for the discrete element method (DEM). Comput Chem. Eng., 32(10), 22632279 Google Scholar
Kudrolli, A and Clotet, X (2016) Evolution of porosity and channelization of an erosive medium driven by fluid flow. Phys. Rev. Lett., 117, 028001Google Scholar
Kulessa, B and 10 others (2017) Seismic evidence for complex sedimentary control of Greenland Ice Sheet flow. Sci. Adv., 3(8), e1603071 Google Scholar
Kyrke-Smith, TM, Katz, RF and Fowler, AC (2015) Subglacial hydrology as a control on emergence, scale, and spacing of ice streams. J. Geophys. Res. Earth Surf., 120, 15011514 Google Scholar
Lajeunesse, E, Malverti, L and Charru, F (2010) Bed-load transport in turbulent flow at the grain scale: experiments and modeling. J. Geophys. Res., 115(F4)Google Scholar
Le Brocq, AM and 10 others, 2013. Evidence from ice shelves for channelized meltwater flow beneath the Antarctic Ice Sheet, Nat. Geosci., 6(11), 945948 CrossRefGoogle Scholar
Lesemann, J-E, Piotrowski, JA and Wysota, W (2010) “Glacial curvilineations”: new glacial landforms produced by longitudinal vortices in subglacial meltwater flows. Geomorphology, 120(3–4), 153161 Google Scholar
Lesemann, J-E, Piotrowski, JA and Wysota, W (2014) Genesis of the “glacial curvilineation” landscape by meltwater processes under the former Scandinavian Ice Sheet, Poland. Sediment. Geol., 312, 118 Google Scholar
Lliboutry, L (1968) General theory of subglacial cavitation and sliding of temperate glaciers. J. Glaciol., 7, 2158 Google Scholar
Luding, S (2008) Introduction to DEMs: basic of contact force models and how to perform the micro-macro transition to continuum theory. Euro. J. Environ. Civil Eng., 12(7–8), 785826 Google Scholar
Mahadevan, A, Orpe, AV, Kudrolli, A and Mahadevan, L (2012) Flow-induced channelization in a porous medium. Eur. Phys. Lett., 98(5), 58003 Google Scholar
Marsh, OJ and 6 others (2016) High basal melting forming a channel at the grounding line of Ross Ice Shelf, Antarctica. Geophys. Res. Lett., 43(1), 250255 Google Scholar
McNamara, S, Flekkøy, EG and Måløy, KJ (2000) Grains and gas flow: Molecular dynamics with hydrodynamic interactions. Phys. Rev. E, 61, 40544059 Google Scholar
Métivier, F, Lajeunesse, E and Devauchelle, O (2017) Laboratory rivers: Lacey's law, threshold theory, and channel stability. Earth Surf. Dyn., 5(1), 187198 Google Scholar
Meyer-Peter, E and Müller, R (1948) Formulas for bed-load transport. Int. Assoc. for Hydraul. Res., Stockholm Google Scholar
Minchew, B and 7 others (2016) Plastic bed beneath Hofsjökull Ice Cap, central Iceland, and the sensitivity of ice flow to surface meltwater flux. J. Glaciol., 62, 147158 Google Scholar
Mitchell, JK and Soga, K (2005) Fundamentals of soil behavior. Wiley, New York Google Scholar
Moore, PL and Iverson, NR (2002) Slow episodic shear of granular materials regulated by dilatant strengthening. Geology, 30(9), 843846 2.0.CO;2>CrossRefGoogle Scholar
Nedderman, RM (1992) Statics and kinematics of granular materials. Cambridge University Press, Cambridge Google Scholar
Ng, FSL (2000a) Canals under sediment-based ice sheets. Ann. Glaciol., 30, 146152 Google Scholar
Ng, FSL (2000b) Coupled ice-till deformation near subglacial channels and cavities. J. Glaciol., 46(155), 580598 Google Scholar
Nye, JF (1976) Water flow in glaciers: jökulhlaups, tunnels and veins. J. Glaciol., 17, 181207 Google Scholar
Obermayr, M, Dressler, K, Vrettos, C and Eberhard, P (2013) A bonded-particle model for cemented sand. Comput. Geotech., 49, 299313 Google Scholar
Parker, G (1978) Self-formed straight rivers with equilibrium banks and mobile bed. Part 1. The sand-silt river. J. Fluid Mech., 89(1), 109 Google Scholar
Perol, T and Rice, JR (2015) Shear heating and weakening of the margins of West Antarctic ice streams. Geophys. Res. Lett., 42(9), 34063413 Google Scholar
Perol, T, Rice, JR, Platt, JD and Suckale, J (2015) Subglacial hydrology and ice stream margin locations. J. Geophys. Res. Earth Surf., 120, 13521368 Google Scholar
Piotrowski, JA (1999) Channelized subglacial drainage under soft-bedded ice sheets: evidence from small N-channels in Central European Lowland. Geol. Q., 43(2), 153162 Google Scholar
Piotrowski, JA, Geletneky, J and Vater, R (1999) Soft-bedded subglacial meltwater channel from the Welzow-Süd open-cast lignite mine, Lower Lusatia, eastern Germany. Boreas, 28(3), 363374 Google Scholar
Piotrowski, JA, Larsen, NK and Junge, FW (2004) Reflections on soft subglacial beds as a mosaic of deforming and stable spots. Quat. Sci. Rev., 23(9), 9931000 Google Scholar
Piotrowski, JA, Hermanowski, P and Piechota, AM (2009) Meltwater discharge through the subglacial bed and its land-forming consequences from numerical experiments in the Polish lowland during the last glaciation. Earth Surf. Process. Landforms, 34(4), 481492 CrossRefGoogle Scholar
Pons, A, Amon, A, Darnige, T, Crassous, J and Clément, E (2015) Mechanical fluctuations suppress the threshold soft-glassy solids: The secular drift scenario. Phys. Rev. E, 92(2)Google Scholar
Radjaï, F and Dubois, F (2011) Discrete-element modeling of granular materials. Wiley-Iste, London, UK Google Scholar
Raymond, CF (2000) Energy balance of ice streams. J. Glaciol., 46(155), 665674 Google Scholar
Rempel, AW (2009) Effective stress profiles and seepage flows beneath glaciers and ice sheets. J. Glaciol., 55(191), 431443 Google Scholar
Rosier, SHR, Gudmundsson, GH and Green, JAM (2015) Temporal variations in the flow of a large Antarctic ice-stream controlled by tidally induced changes in the subglacial water system. Cryosphere Discuss., 9, 23972429 Google Scholar
Röthlisberger, H (1972) Water pressure in intra-and subglacial channels. J. Glaciol., 11(62), 177203 Google Scholar
Schellart, WP (2000) Shear test results for cohesion and friction coefficients for different granular materials: scaling implications for their usage in analogue modelling. Tectonophysics, 324(1), 116 Google Scholar
Schofield, AN and Wroth, P (1968) Critical state soil mechanics. McGraw-Hill, London Google Scholar
Schoof, C (2010) Ice-sheet acceleration driven by melt supply variability. Nature, 468(7325), 803806 Google Scholar
Schoof, C, Rada, CA, Wilson, NJ, Flowers, GE and Haseloff, M (2014) Oscillatory subglacial drainage in the absence of surface melt. Cryosphere, 7(6), 959976 Google Scholar
Schroeder, DM, Blankenship, DD and Young, DA (2013) Evidence for a water system transition beneath Thwaites Glacier, West Antarctica. Proc. Natl. Acad. Sci. USA, 110(30), 1222512228 Google Scholar
Shields, IA (1936) Anwendung der Aehnlichkeitsmechanik und der Turbulenzforschung auf die Geschiebebewegung. Mitt. Preuss Ver.-Anst, 26 Google Scholar
Shoemaker, EM (1986) Subglacial hydrology for an ice sheet resting on a deformable aquifer. J. Glaciol., 32(110), 2030 CrossRefGoogle Scholar
Shreve, RL (1972) Movement of water in glaciers. J. Glaciol., 11(62), 205214 Google Scholar
Simkins, LM and 8 others (2017) Anatomy of a meltwater drainage system beneath the ancestral East Antarctic ice sheet. Nat. Geosci. 10, 691697 Google Scholar
Stevens, LA and 7 others (2015) Greenland supraglacial lake drainages triggered by hydrologically induced basal slip. Nature, 522(7554), 7376 CrossRefGoogle ScholarPubMed
Stock, JD and Montgomery, DR (1999) Geologic constraints on bedrock river incision using the stream power law. J. Geophys. Res. Solid Earth, 104(B3), 49834993 Google Scholar
Stokes, CR and Clark, CD (2001) Palaeo-ice streams. Quat. Sci. Rev., 20(13), 14371457 Google Scholar
Suckale, J, Platt, JD, Perol, T and Rice, JR (2014) Deformation-induced melting in the margins of the West Antarctic ice streams. J. Geophys. Res. Earth Surf., 119(5), 10041025 Google Scholar
Swift, DA, Nienow, PW, Spedding, N and Hoey, TB (2002) Geomorphic implications of subglacial drainage configuration: rates of basal sediment evacuation controlled by seasonal drainage system evolution. Sediment. Geol., 149, 519 Google Scholar
Terzaghi, K (1943) Theoretical soil mechanics. Wiley. London, UK Google Scholar
Terzaghi, K, Peck, RB and Mesri, G (1996) Soil mechanics in engineering practice. John Wiley & Sons, London, UK Google Scholar
Thompson, J, Simons, M and Tsai, VC (2014) Modeling the elastic transmission of tidal stresses to great distances inland in channelized ice streams. Cryosphere, 8(6), 20072029 Google Scholar
Tulaczyk, S (2006) Scale independence of till rheology. J. Glaciol., 52(178), 377380 Google Scholar
Tulaczyk, S, Kamb, WB and Engelhardt, HF (2000a) Basal mechanics of ice stream B, West Antarctica I. Till mechanics. J. Geophys. Res., 105(B1), 463481 Google Scholar
Tulaczyk, S, Kamb, WB and Engelhardt, HF (2000b) Basal mechanics of ice stream B, West Antarctica II. Undrained plastic-bed model. J. Geophys. Res., 105(B1), 483494 Google Scholar
Tylmann, K, Piotrowski, JA and Wysota, W (2013) The ice/bed interface mosaic: deforming spots intervening with stable areas under the fringe of the Scandinavian Ice Sheet at Sampława, Poland. Boreas, 42(2), 428441 Google Scholar
Utter, B, Behringer, RP (2004) Self-diffusion in dense granular shear flows. Phys. Rev. E, 69(3)Google Scholar
van der Wel, N, Christoffersen, P and Bougamont, M (2013) The influence of subglacial hydrology on the flow of Kamb Ice Stream, West Antarctica, J. Geophys. Res. Earth Surf., 118(1), 97110 Google Scholar
Walder, JS and Fowler, A (1994) Channelized subglacial drainage over a deformable bed. J. Glaciol., 40(134)Google Scholar
Walder, JS, Iverson, RM, Godt, JW, Logan, M and Solovitz, SA (2015) Controls on the breach geometry and flood hydrograph during overtopping of noncohesive earthen dams. Water Resour. Res., 51, 67016724 Google Scholar
Walker, RT, Christianson, K, Parizek, BR, Anandakrishnan, S and Alley, RB (2012) A viscoelastic flowline model applied to tidal forcing of Bindschadler Ice Stream, West Antarctica. Earth Planet. Sci. Lett., 319, 128132 Google Scholar
Weertman, J (1957) On the sliding of glaciers. J. Glaciol., 3(21), 3338 Google Scholar
Werder, MA, Hewitt, IJ, Schoof, CG and Flowers, GE (2013) Modeling channelized and distributed subglacial drainage in two dimensions. J. Geophys. Res. Earth Surf., 118(4), 21402158 Google Scholar
Whipple, KX and Tucker, GE (1999) Dynamics of the stream-power river incision model: implications for height limits of mountain ranges, landscape response timescales, and research needs. J. Geophys. Res. Solid Earth, 104(B8), 1766117674 CrossRefGoogle Scholar
Wilcock, PR (1998) Two-fraction model of initial sediment motion in gravel-bed rivers. Science, 280(5362), 410412 Google Scholar
Winberry, JP, Anandakrishnan, S and Alley, RB (2009) Seismic observations of transient subglacial water-flow beneath MacAyeal Ice Stream, West Antarctica. Geophys. Res. Lett., 36(11)Google Scholar
Xu, BH and Yu, AB (1997) Numerical simulation of the gas-solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics. Chem. Eng. Sci., 52(16), 27852809 Google Scholar
Xu, Y, Rignot, E, Menemenlis, D and Koppes, M (2012) Numerical experiments on subaqueous melting of Greenland tidewater glaciers in response to ocean warming and enhanced subglacial discharge. Ann. Glaciol., 53(60), 229234 Google Scholar
Zhang, Y and Campbell, CS (1992) The interface between fluid-like and solid-like behavior in two-dimensional granular flows. J. Fluid Mech., 237 541 Google Scholar
Zhou, ZY, Kuang, SB, Chu, KW and Yu, AB (2010) Discrete particle simulation of particle-fluid flow: model formulations and their applicability. J. Fluid Mech., 661, 482510 Google Scholar
Zhu, HP, Zhou, ZY, Yang, RY and Yu, AB (2007) Discrete particle simulation of particulate systems: Theoretical developments. Chem. Eng. Sci., 62, 33783396 Google Scholar
Figure 0

Fig. 1. Glacial curvilineations incised into sedimentary plateaus by subglacial meltwater erosion during the last glaciation near Zbójno, Poland, from Lesemann and others (2010). Water and ice-sheet flow was generally from the top left (NW) to the lower right (SE).

Figure 1

Fig. 2. Overview of the granular assemblage of 58 000 particles and the discretization of the fluid grid. We assume symmetry around the channel center (+x) and limit our simulation domain to one of the sides. The model domains of the two phases are superimposed during the simulations.

Figure 2

Fig. 3. Boundary conditions and initial state in the presented simulations for the granular phase in the dry experiments. The frictionless lateral boundaries imply no boundary-parallel movement.

Figure 3

Fig. 4. Boundary conditions and initial state in the presented simulations for the granular and fluid phases in the wet experiments.

Figure 4

Table 1. Simulation parameters and their values for the granular experiments

Figure 5

Fig. 5. Averaged grain displacements through time for a range of effective normal stresses, N in the dry granular experiments.

Figure 6

Fig. 6. Channel geometries in the steady state at different effective normal stresses N in dry experiments. Grains are colored according to their cumulative vertical displacement and rendered in an orthogonal projection. The absence of water does not influence the steady-state geometry unless the sediment is flushed by strong water-pressure gradients.

Figure 7

Fig. 7. Observed maximum channel width, Wmax, under a constant imposed value for effective normal stress, N. Data points are fitted using a non-linear least-squares fit, with fitted parameter values and their corresponding std dev. noted.

Figure 8

Fig. 8. Idealized schematic of effective stress magnitude and till yield strength in a cross-section around a subglacial channel incised in the sedimentary bed. Sediments beneath the channel floor are not stressed by ice weight on the flanks, and comprise a weak channel floor wedge (CFW).

Figure 9

Fig. 9. Contact pressure on the grains $\left(\sum \vert\vert{\bi f}_{\rm n}\vert\vert 4 \pi r^2 \right)$ in the steady state with an effective stress of N = 7 kPa.

Figure 10

Fig. 10. Fluid-pressure forces (Eqn 12) on the individual sediment grains visualized by arrows. Top: Under low gradients in pore-water pressure (here dpf/dx = −0.1 kPa m−1), the fluid-pressure force primarily contributes weak buoyant uplift on the grains. Bottom: Larger pressure gradients (dpf/dx =  − 10 kPa m−1) destabilize the sediment through liquefaction at the channel floor and cause collapse of the sediment with associated chaotic interactions (e.g., large internal forces). Individual force vectors are not representative of the grain ensemble and should not be interpreted directly. Both simulation snapshots are from t = 1.5 s.

Figure 11

Fig. 11. Total per-grain spatial displacements at different times for three wet (water-saturated) simulations with different imposed water-pressure gradients and an effective stress of N = 10 kPa. The water-pressure gradients cause flow and drag forces toward the channel, and destabilize the conduit at higher values.

Figure 12

Fig. 12. Example run of the soft-bed subglacial channel model outlined in Eqns (1–6, 15 and 16), with a square-root ice geometry and linearly increasing meltwater influx (${\dot m}$) towards +s. Top: Channel effective pressure (Pc), ice-overburden pressure (Pi), water-pressure (Pw) and water flux (Q). Middle: Sediment flux increases non-linearly with water flux (Eqn 2). Bottom: Maximum (Smax) and actual channel cross-sectional size (S), together with channel growth rate (dS/dt). This example is with a constant forcing and t = 2 days. The minimum channel size is set to S = 0.01 m2 for numerical considerations. In this example the sediment yield strength prevents channel existence except near the terminus where effective pressure is relatively low.

Figure 13

Fig. 13. Examples of infilled subglacial channels in the southern, marginal part of the Scandinavian Ice Sheet from the last glaciation (Piotrowski et al., in prep.). (a) Subglacial channel at Ebeltoft, Djursland, Denmark. The channel is found within a single till unit about 10 km inside the ice margin. It is flat-topped and filled with parallel-bedded outwash sand and gravel intercalated with layers of silt. Single outsized stones, possibly dropped from the channel roof are randomly dispersed in the outwash deposit. Along the channel bottom and on its left-hand side the infill material is deformed into attenuated folds, irregular detached sediment pods and flame structures. Since the flanks of the channel are below the angle of repose of sand, the deformation structures suggest syndepositional sediment intrusion into the channel driven by a pressure gradient oriented toward the channel axis. (b). Subglacial channel at Glaznoty, north-central Poland, about 25 km within the ice limit. The channel occurs at the interface between proglacial outwash deposits (below) and till (above). The channel is distinctly lens-shaped with an upward-convex top suggesting an R-channel incised upward into the ice at a late stage of formation. It is infilled with massive coarse-grained sand and gravel.