A new methodology to simulate subglacial deformation of water-saturated granular material

. The dynamics of glaciers are to a large degree governed by processes operating at the ice–bed interface, and one of the primary mechanisms of glacier ﬂow over soft unconsolidated sediments is subglacial deformation. However, it has proven difﬁcult to constrain the mechanical response of subglacial sediment to the shear stress of an overriding glacier. In this study, we present a new methodology de-signed to simulate subglacial deformation using a coupled numerical model for computational experiments on grain-ﬂuid mixtures. The granular phase is simulated on a per-grain basis by the discrete element method. The pore water is mod-eled as a compressible Newtonian ﬂuid without inertia. The numerical approach allows close monitoring of the internal behavior under a range of conditions. Our computational experiments support the ﬁndings of previous studies where the rheology of a slowly deforming water-saturated granular bed in the steady state gener-ally conforms to the rate-independent plastic rheology. Before this so-called critical state, deformation is in many cases accompanied by volumetric changes as grain rearrangement in active shear zones changes the local porosity. For previously consolidated beds porosity increases can cause local pore-pressure decline, dependent on till permeability and shear rate. We observe that the pore-water pressure reduction strengthens inter-granular contacts, which results in increased shear strength of the granular material. In contrast, weakening takes place when shear deformation causes consolidation of dilated sediments or during rapid fabric development. Both processes of strengthening and weakening de-pend inversely on the sediment permeability and are transient phenomena tied to the porosity changes during the early stages of shear.

ice shelves leaves lateral (Whillans and van der Veen, 1997;Tulaczyk et al., 2000b;Price et al., 2002) and basal friction (Alley, 1993;MacAyeal et al., 1995;Stokes et al., 2007;Sergienko and Hindmarsh, 2013) as the main mechanical components resisting the flow. A fundamental understanding of subglacial dynamics is a requirement for our ability to predict future response of the ice sheets to climate change.
The pressure and flow of pore water in the subglacial bed can influence subglacial deformation in a number of ways. Assuming a Mohr-Coulomb constitutive relation of the basal till strength, an increase in pore-water pressure weakens the bed by reducing the effective stress, and this may facilitate basal movement if the driving shear stresses become sufficient to overcome the sediment yield strength (Kamb, 1991;Iverson et al., 1998;Tulaczyk et al., 2000a;Fischer and Clarke, 2001;Kavanaugh and Clarke, 2006).
If the hydraulic diffusivity of the bed is sufficiently low relative to the rate of deformation, a modulation of the porewater pressure at the ice-bed interface is, over time, carried into the subglacial bed, resulting in variable internal yield strength and ultimately variable shear strain rates with depth (Tulaczyk, 1999;Tulaczyk et al., 2000a;Kavanaugh and Clarke, 2006). Owing to local volumetric changes, variations from the hydrostatic fluid pressure distribution can be created inside the sediment by the onset and halt of granular deformation. This influences the local effective pressure and, in turn, the sediment yield strength (e.g., Iverson et al., 1998).
In the case of non-planar ice-bed geometry, excess porewater pressures can develop on the stoss side of objects ploughing through a subglacial bed (Iverson et al., 1994;Iverson, 1999;Thomason and Iverson, 2008). The elevated pore-water pressure weakens the sediment by lowering the effective stress, resulting in a net strain-rate weakening rheology (e.g., Iverson et al., 1998;Clark et al., 2003;Iverson, 2010), which has been associated with the stick-slip behavior of Whillans Ice Stream (Winberry et al., 2009).
Early field studies suggested a strain-rate strengthening Bingham or slightly non-linear viscous rheology of till (Boulton and Hindmarsh, 1987), which has been used to simplify analytical and numerical modeling of till mechanics (e.g., Alley et al., 1987;Hindmarsh, 1998;Fowler, 2000). Laboratory studies have, however, strongly favored the notion of till having a plastic, Mohr-Coulomb rheology, with a very small rate dependence in the case of critical-state deformation (Kamb, 1991;Iverson et al., 1998;Tulaczyk et al., 2000a;Rathbun et al., 2008). This Mohr-Coulomb rheology is also supported by field investigations (Truffer et al., 2000;Kavanaugh and Clarke, 2006). However, a rate-weakening rheology is expected in the case where obstacles plough through a soft and deformable bed (Iverson et al., 1994;Iverson, 1999; Thomason and Iverson, 2008).
The low viscosity of water facilitates flow at low stress and can only be expected to influence the overall strength

The fluid model
The inter-particle fluid is handled by conventional continuum computation dynamics (CFD). The implementation follows the compressible Darcian flow presented by Goren et al. (2011). This approach was favored over a full Stokes solution of fluid flow (Gidaspow, 1994;Zhu et al., 2007;Zhou et al. Kloss et al., 2012) since it allows for convenient parameterization of the hy ical permeabilities. The model assumes insignificant fluid inertia, which is priate for the subglacial setting. The volumetric fraction of the fluid phase (the porosity, ) is incorpor the Eulerian formulations of the compressible continuity equation and mom equation using the local average method (Anderson and Jackson, 1967;Yu, 1997). The Darcy constitutive equation is used for conserving mom (Eq. 5.5) (McNamara et al., 2000;Goren et al., 2011): where v f is the fluid velocity, v p is the particle velocity, k is the hydraulic ability, f is the adiabatic fluid compressibility and µ f is the dynamic fluid vi The continuity equation (Eq. 5.4) is in the form of a transient diffusion eq with the forcing term acting as a source/sink for the fluid pressure. The pr p f , is the pressure deviation from the hydrostatic pressure distribution. Th sure deviation is sometimes referred to as the excess pressure. We refrai using this term, as it may be misleading for pressures that are smaller th of subglacial materials in a few select scenarios, e.g., when pore-water flow due to porosity change in early stages of shear deformation is limited by low hydraulic diffusivities (e.g., Iverson, 2010). The mechanics of coupled granularfluid mixtures have previously been numerically investigated for studies of fluidized beds (e.g., Anderson and Jackson, 1967;Gidaspow et al., 1992;Hoomans et al., 1996;Xu and Yu, 1997;McNamara et al., 2000;Feng and Yu, 2004;Jajcevic et al., 2013), the stability of inclined, fluid-immersed granular materials (e.g., Topin et al., 2011;Mutabaruka et al., 2014), mechanics during confined deformation (e.g., Goren et al., 2011;Catalano et al., 2014), debris flow (e.g., Hutter et al., 1994;Mangeney et al., 2007;Goren et al., 2011), and for the design of industrial components, e.g., hydrocyclones (e.g., Wang et al., 2007;Zhou et al., 2010), or silos and hoppers (Kloss et al., 2012).
In this study, we explore the interaction between the fluid and granular phases in water-saturated consolidated particle assemblages undergoing slow shear deformation. A dry granular assemblage deforms rate-independently in a pseudostatic manner when deformational rates are sufficiently low (GDR-MiDi, 2004;Damsgaard et al., 2013). Critical-state deformation of water-saturated granular materials is rate independent during slow shear, but transient changes in pore pressure may cause initiation of liquefaction at higher shear velocities (Goren et al., 2011). The particle-fluid mixture is in this study sheared at velocities and stresses comparable to those found in subglacial settings. The computational approach allows for investigating the internal granular mechanics and feedbacks during progressive shear deformation.
In the following section, we present the details of the numerical implementation of particle-fluid flow and describe The Cryosphere, 9, 2183-2200, 2015 www.the-cryosphere.net/9/2183/2015/ the experimental setup. We then present and discuss the modeled deformational properties of the particle-fluid mixture. Finally, we analyze how the fluid influences formation of shear zones and under which conditions deformation is rate dependent.

The granular model
We use the discrete element method (DEM) (Cundall and Strack, 1979) to simulate the granular deformation. With the DEM, particles are treated as separate, cohesionless entities, which interact by soft-body deformation defined by a prescribed contact law. The contact mechanics are micromechanically parameterized. The temporal evolution is handled by integration of the momentum equations of translation: and rotation: (2) i and j are particle indexes, m is the particle mass, I is the particle rotational inertia, and x and are linear and rotational particle positions, respectively. The model considers the transient kinematics during time t. Gravitational acceleration is included as the vector g. The inter-particle contact force vectors are denoted f n and f t in the normal and tangential direction relative to the contact interface. The particle radius is included as r. The fluid-particle interaction force is denoted f i (Fig. 1), and the inter-particle normal vector is termed n. Here, δ n is the virtual inter-particle overlap at the contact, scaling the compressive strain of the particle. The inter-particle contact forces are determined by a linear-elastic contact model causing repulsion between particle pairs and providing frictional contact strength. The elastic stiffness is k n in the normal direction to the contact and k t parallel (tangential) to the contact interface. The magnitude of the tangential force f t is limited by the Coulomb frictional coefficient µ (Cundall and Strack, 1979;Luding, 2008;Radjaï and Dubois, 2011;Damsgaard et al., 2013). The bulk Coulomb friction coefficient of a granular assemblage depends on the elastic properties, the inter-particle frictional coefficient, and the granular packing arrangement (e.g., Booth et al., 2014).
The vector δ t is the tangential displacement on the interparticle interface when corrected for contact rotation. In the case of slip, the length of δ t is adjusted to a length consistent with Coulomb's condition (||δ t || = −µ||f n ||/k t ) (Luding, 2008;Radjaï and Dubois, 2011). The linear elasticity allows temporal integration with a constant time step length t.

The fluid model
The inter-particle fluid is handled by using conventional continuum mechanics. The implementation follows the compressible Darcian flow model presented by Goren et al. (2011). This approach was favored over a full Navier-Stokes solution of fluid flow (Gidaspow, 1994;Zhu et al., 2007;Zhou et al., 2010;Kloss et al., 2012) because it allows for convenient parameterization of the hydrological permeabilities. The model assumes insignificant fluid inertia, which is appropriate for subglacial settings. Fluid compressibility is only important under undrained or very low-permeability conditions (e.g., Goren et al., 2011). We include fluid compressibility since the numerical method does not require us to impose the common simplification of incompressibility. The volumetric fraction of the fluid phase (the porosity, φ) is incorporated in the Eulerian formulations of the compressible continuity equation and momentum equation using the local average method (Anderson and Jackson, 1967;Xu and Yu, 1997). The Darcy constitutive equation is used for conserving momentum (Eq. 5) (McNamara et al., 2000;Goren et al., 2011): where v f is the fluid velocity, v p is the particle velocity, and k is the hydraulic permeability estimated from the local porosity. The adiabatic fluid compressibility is denoted β f , and µ f is the dynamic fluid viscosity. The continuity equation (Eq. 4) is in the form of a transient diffusion equation with the forcing term acting as a source/sink for the fluid pressure.
The pressure, p f , is the pressure deviation from the hydrostatic pressure distribution. This pressure deviation is somewww.the-cryosphere.net/9/2183/2015/ The Cryosphere, 9, 2183-2200, 2015 times referred to as the excess pressure. We refrain from using this term, as it may be misleading for pressures that are smaller than the hydrostatic value. The simulation domain is discretized in a regular rectilinear orthogonal grid. The pressure is found using the Crank-Nicolson method of mixed explicit and implicit temporal integration, which is unconditionally stable and second-order accurate in time and space (e.g., Patankar, 1980;Ferziger and Perić, 2002;Press et al., 2007). The implicit solution is obtained using the iterative Jacobi relaxation method (e.g., Ferziger and Perić, 2002;Press et al., 2007;Gerya, 2010), which is light on memory requirements and ideal in terms of parallelism for our graphics processing unit (GPU) implementation, although not optimal in terms of convergence. The numerical solution is continuously checked by the Courant-Friedrichs-Lewy condition (Courant et al., 1967). The partial derivatives are approximated by finite differences.

Porosity
The local porosity is determined at the fluid cell center. For a cell with a set of N grains in its vicinity, it is determined by inverse-distance weighting the grains using a bilinear interpolation scheme (McNamara et al., 2000;Goren et al., 2011). The weight function s has the value 1 at the cell center and decreases linearly to 0 at a distance equal to the cell width ( x, Fig. 2): V i g is the volume of grain i. x 3 is the fluid cell volume, and x f is the cell center position.
is the product operator. The average grain velocity at the cell center (v) is found using the same weighting function described above (Eq. 6). Additionally, large grains with velocity v i contribute to the cell-average velocity with a greater magnitude: v The change in porosity is the main forcing the particles exert onto the fluid (Eq. 5). At time step n it is estimated by central differences for second-order accuracy: The porosity at n + 1 is found by estimating the upcoming particle positions from temporal integration of their current positions and velocities.

Permeability
Relations between grain size and hydraulic properties of sediments supported by empirical evidence (e.g., Hazen, 1911;Kozeny, 1927;Carman, 1937;Krumbein and Monk, 1943;Harleman et al., 1963;Schwartz and Zhang, 2003). The Kozeny-Carman estimation of permeability k is commonly used and is of the form where d is the representative grain diameter. Due to constraints on the computational time step we are unable to simulate fine grain sizes with realistic elastic properties within a reasonable time frame. In order to give a first-order estimate of the deformational behavior of fine-grained sediments, we therefore use a modified version of the above relationship, where the permeability varies as a function of the porosity and a predefined permeability pre-factor k c : The Cryosphere, 9, 2183-2200, 2015 www.the-cryosphere.net/9/2183/2015/ Using this approach we can simulate large particles with the hydrological properties of fine-grained materials, while retaining the effect of local porosity variations on the intrinsic permeability. We do note, however, that the dilative magnitude during deformation is likely different for clay materials due to their plate-like shape. This is because sediments with a considerable amount of arbitrarily oriented clay minerals are likely to compact during deformation as the clay particles align to accommodate shear strain.

Particle-fluid interaction
The dynamic coupling from the pore fluid to the solid particles acts through the particle-fluid force (f i ) in Eq. (1). Our implementation of this coupling follows the procedure outlined by Xu and Yu (1997), Feng and Yu (2004), and Zhou et al. (2010) (scheme 3). In a complete formulation, the interaction force between particles and the pore fluid consists of several components (e.g., Zhu et al., 2007). The most well-known interaction force is the drag force, which is evident in Stokes' experiments of particle settling. This force is typically incorporated in particle-fluid models by semi-empirical relationships. It scales linearly with the relative velocity difference between particle and fluid and the particle-surface drag coefficient (Ergun, 1952;Wen and Yu, 1966;Gidaspow et al., 1992;Di Felice, 1994;Zhu et al., 2007).
The pressure-gradient force is another important fluidinteraction force (Anderson and Jackson, 1967;Zhu et al., 2007). It results when there is a difference in pressure across a surface and is well known from buoyancy. In liquids with hydrostatic pressure distributions the pressure grows linearly with depth as the weight of the overburden fluid increases. The resultant surface pressure on submerged body increases correspondingly with depth and causes a buoyant force upwards. In fluids without hydrostatic pressure distribution, the pressure-gradient force can be a highly variable vector field.
Initial tests using a full Navier-Stokes solution for the fluid phase showed us that the pressure-gradient force was by far the dominant interaction force in our pseudo-static shear experiments. The drag force was the second-most important force but 2 orders of magnitude weaker than the pressuregradient force due to very slow fluid and particle flow. Since we neglect fluid inertia in our fluid description, we included only the pressure-gradient force. The fluid pressure in our model records the pressure difference from the hydrostatic pressure. For this reason we add a term to the pressuregradient force which describes the buoyancy of a fully sub-

Meth
Normal stress 0 on wall with fixed p . The bottom boundary is impermeable (no flow, free slip Neumann). The fluid grid is extended upwards to allow for dilation and movement of the upper wall. The granular phase is compressed between a fixed wall at the bottom, and a dynamic top wall, which exerts a normal stress ( 0 ) downwards. The material is sheared by moving the topmost particles parallel to the x-axis.
The highest value (k c = 3.5 ⇥ 10 13 m 2 ) matches a permeability of 1.9 ⇥ 10 14 These two end-member permeabilities roughly correspond to what Iverson et (1997a) and Iverson et al. (1998) estimated for the clay-rich Two Rivers till the clay-poor Storgläciaren till, respectively.
The lower boundary is impermeable, and a fixed fluid pressure is specified the top boundary. These boundary conditions imply that the ice-bed interfac a relatively permeable zone with rapid diffusion of hydrological pressure, wh is likely for subglacial beds with low permeability (e.g. Alley, 1989;Creyts Schoof, 2009;Kyrke-Smith et al., 2014). In coarse-grained tills it is likely that subglacial till diffusivity exceeds the hydraulic diffusivity at the ice bed interfa The lateral boundaries are periodic (wrap-around). If a particle moves outside grid on the right side it immediately reappears on the left side. Likewise, part pairs can be in mechanical contact although placed on opposite sides of the g at the periodic boundaries.
The particle size distribution is narrow compared to that of subglacial t which often display a fractal size distribution (e.g. Hooke and Iverson, 1995). F tal size distributions minimize internal granular compressive stress heterogenei (Hooke and Iverson, 1995;Iverson et al., 1996), but, in the absence of grain cru . The bottom boundary is impermeable (no flow, free slip Neumann). The fluid grid is extended upwards to allow for dilation and movement of the upper wall. The granular phase is compressed between a fixed wall at the bottom and a dynamic top wall, which exerts a normal stress (σ 0 ) downwards. The material is sheared by moving the topmost particles parallel to the x axis. merged particle as the weight of the displaced fluid: where V g is the volume of the particle, ρ f is the fluid density, and g is the vector of gravitational acceleration. The particlefluid interaction force is added to the sum of linear forces per particle (Eq. 1). The particle force is not added to the fluid pressure equation (Eq. 4) since fluid inertia is ignored. The fluid is instead forced by variations in porosity.

Computational experiments
The computational fluid dynamics algorithm is implemented in CUDA C (NVIDIA, 2013) in order to allow a direct integration with the GPU-based particle solver. The coupled particle-fluid code is free software (source code available at https://github.com/anders-dc/sphere), licensed under the GNU Public License v.3. The simulations were performed on a GNU/Linux system with a pair of NVIDIA Tesla K20c GPUs. The experimental results are visualized using Par-aView (Henderson et al., 2007) and Matplotlib (Hunter, 2007). The experimental setup is a rectangular volume (Fig. 3) where a fluid-saturated particle assemblage deforms due to www.the-cryosphere.net/9/2183/2015/ The Cryosphere, 9, 2183-2200, 2015 To determine the effects of the pore water, we perform experiments with and without fluids, and for the experiments with fluids present, the permeability pre-factor k c is varied to constrain the effect of the hydraulic conductivity and diffusivity on the overall deformation style. The low value used for k c (3.5 × 10 −15 m 2 ) results in an intrinsic permeability of k = 1.9 × 10 −16 m 2 for a porosity of 0.3 (Eq. 11). The highest value (k c = 3.5 × 10 −13 m 2 ) matches a permeability of 1.9 × 10 −14 m 2 . These two end-member permeabilities roughly correspond to what Iverson et al. (1997a) and Iverson et al. (1998) estimated for the clay-rich Two Rivers till and the clay-poor Storglaciären till, respectively.
The lower boundary is impermeable, and a fixed fluid pressure is specified for the top boundary. These boundary conditions imply that the simulated ice-bed interface is a relatively permeable zone with rapid diffusion of hydrological pressure, which is likely for subglacial beds with low permeability (e.g., Alley, 1989;Creyts and Schoof, 2009;Kyrke-Smith et al., 2014). In coarse-grained tills the subglacial till diffusivity may exceed the hydraulic diffusivity at the ice-bed interface, unless a network of linked cavities causes rapid diffusion of pore-water pressures here (e.g., Kamb, 1987). The lateral boundaries are periodic (wrap-around). When a particle moves outside the grid on the right side it immediately reappears on the left side. Likewise, particle pairs can be in mechanical contact although placed on opposite sides of the grid.
The particle size distribution is equidimensional, a stark contrast relative to subglacial tills which often display a fractal size distribution (e.g., Hooke and Iverson, 1995). Fractal size distributions minimize internal stress heterogeneities (Hooke and Iverson, 1995;Iverson et al., 1996), but, in the absence of grain crushing, an assemblage with a wide particle size distribution dilates from a consolidated state with the same magnitude as assemblages with a narrow particle size distribution (Morgan, 1999) and displays the same frictional strength (Morgan, 1999;Mair et al., 2002;Mair and Hazzard, 2007). The comparable dilation magnitude justifies the computationally efficient narrow particle size distribution used here. However, as previously noted, shear zones in clayrich materials can compact during shear due to preferential parallel alignment, which is not possible to capture with the methodology presented here.
The simulated particle size falls in the gravel category of grain size (Table 1). The large size allows us to perform the temporal integration with larger time steps (Radjaï and Dubois, 2011;Damsgaard et al., 2013). The frictional force between two bodies is independent of their size (Amontons' second law) but is proportional to the normal force on the contact interface (Mitchell and Soga, 2005), as reflected in the contact law in the discrete element method (Eq. 3). We prescribe the normal forcing at the boundary as a normal stress, which implies that the normal force exerted onto a particle assemblage at the boundaries scales with domain size. For a number of total particles in a given packing configuration the ratio between particle size and inter-particle force is constant, which causes the shear strength to be independent of simulated particle size. This scale independence is The shear friction values (top) are smoothed with a moving Hanning window function to approximate the strength of larger particle assemblages. The material peak friction increases with strain rate due to reductions of internal fluid pressure. This strengthening is taking place when the dilation rate exceeds the dissipation rate of the fluid. The diagonal hatched pattern marks the earliest state of shear characterized by rapid dilation.
verified in laboratory experiments, where the granular shear strength of non-clay materials is known to be mainly governed by grain shape and surface roughness instead of grain size (Schellart, 2000;Mitchell and Soga, 2005).

Experiment preparation and procedure
The particles are initially positioned in a dry, tall volume, from where gravity allows them to settle into a dense state at the bottom of the bounding box. The particle assemblage is then consolidated by moving the top wall downwards until the desired level of consolidation stress is reached for an extended amount of time. The same top wall is thereafter used to shear the material in a fluid-saturated state (Fig. 3).
For the shear experiments, the uppermost particles are forced to move with the top wall at a prescribed horizontal velocity (Fig. 3). The particles just above the bottom wall are prescribed to be neither moving or rotating. The physical size of the domain is 0.52 m long in the direction of shear (x), 0.26 m in the transverse horizontal direction (y), and 0.55 m tall (z). We use the exact same material for all shear experiments, allowing for perfect reproducibility and effectively removing inter-experimental variability caused by differences in grain packing. The micromechanical properties and geometrical values used are listed in Table 1  Rate-independent Figure 6. Peak frictional strength before the critical state of the lowpermeability granular bed (k c = 3.5 × 10 −15 m 2 ) at different shear velocities. The frictional strength is constant and rate independent at velocities lower than 10 1 m a −1 as pore-pressure diffusion rates far exceed rates in volumetric change.

Scaling of the shear velocity
The heavy computational requirements of the discrete element method necessitates upscaling of the shearing velocity in order to reach a considerable shear strain within a manageable length of time. Temporal upscaling does not influence the mechanical behavior of dry granular materials, as long as the velocity is below a certain limiting velocity (GDR-MiDi, 2004;Damsgaard et al., 2013;Gu et al., 2014). The shearing velocity used here (2.32×10 −2 m s −1 ), although roughly 3 orders of magnitude greater than the velocities observed in subglacial environments (e.g., 316 m a −1 = 10 −5 m s −1 ), guarantees quasi-static, rate-independent deformation in the granular phase, identical to the behavior at lower strain rates. The particle inertia parameter, I , quantifies the influence of grain inertia in dry granular materials (GDR-MiDi, 2004). Values of I below 10 −3 correspond to pseudo-static and rateindependent shear deformation in dry granular materials. I has a value of 1.7 × 10 −4 in the shear experiments of this study.
The fluid phase needs separate treatment in order to correctly resolve slow shear behavior at faster shearing velocities. This behavior is achieved by linearly scaling the fluid dynamic viscosity with the relationship between actual shearing velocity and the reference glacial sliding velocity. By decreasing the viscosity, the fluid is allowed to more quickly adjust to external and internal forcings. The velocity scaling adjusts the time-dependent parameters of hydraulic conductivity and diffusivity correctly. The intrinsic permeability k is time independent, and the values produced here are directly comparable with real geological materials. The fluid viscosity is scaled to a lower value of 1.797×10 −6 Pa s, consistent with the scaling factor used for the shearing velocity. We test the influence of shearing rate by varying this parameter.

Results
First we investigate the strain-rate dependence of the sediment strength and dilation by shearing a relatively impermeable sediment (k c = 3.5 × 10 −15 m 2 ) at different shear velocities. The shear velocity directly influences the magnitude of the peak shear strength, dilation, and internal fluid pressure (Figs. 4 and 5). At the reference shearing velocity the peak shear frictional strength is 0.71, which corresponds to a shear stress of 14 kPa at an effective stress of 20 kPa (Figs. 4,top left,and 6). When sheared 100 times slower, the peak shear friction has decreased to 0.62, corresponding to a shear stress of 12 kPa (Figs. 4, top right, and 6). The peak values are measured during the transition from the dense and consolidated pre-failure state to the critical state where a shear zone is fully established. This transition is characterized by rapid dilation due to porosity increases in the shear zone ( Fig. 4, middle). During fast shearing velocities the volumetric change outpaces the diffusion of fluid pressure, causing the internal pore-water pressure in the shear zone to decline (Figs. 4,bottom,and 5). Dilatant hardening causes the peak shear strength to increase at large shear velocities (Fig. 6), while the strength reduces to values corresponding to the rate-independent granular strength for lower velocities.
In this model framework, adjusting the hydraulic permeability of the same coarse sediment leads to similar conditional strengthening as shearing the sediment at different rates (Fig. 7). Without fluids (the dry experiment), the peak shear friction (Fig. 7, left) is relatively low and is domi-nated by high-frequency fluctuations. The fluid-saturated experiment with the relatively high permeability (k c = 3.5 × 10 −13 m 2 ) has similar shear strength, but the high-frequency oscillations in shear friction are reduced by the fluid presence. The dilation is similar to the dry experiment but with slightly decreased magnitude. The mean fluid pressure deviation from hydrostatic values (Fig. 7, bottom left) is close to 0. The low-permeability experiment (Fig. 7, right) is characterized by the largest initial peak strength and lowest magnitude of dilation. Compared to the other experiments, the dilation reaches its maximum values at lower shear strain. The fluid pressure decreases almost instantaneously at first, whereafter it equilibrates towards the hydrostatic value (0 Pa).
At constant shearing rate with different permeabilities (Fig. 8, top) or at variable shearing rates with constant permeability (Fig. 8, bottom), we observe that pore-water dynamics have a significant effect on the distribution of strain. The effects of the fluid are visible at different depths within the deforming material (Figs. 9 and 10). The deformation is pervasive with depth for the relatively permeable experiment ( Fig. 9 top), and the fluid pressures deviate only slightly from the hydrostatic values. The relatively small pressure gradients cause only weak fluid forces on the particles in this experiment. Contrasting these results, deformation in the impermeable experiment is primarily governed by shallow deformation beneath the top wall ( Figs. 9 and 10, bottom).
Differences in hydraulic permeability influence the dynamics of the fluid over time, as illustrated in Fig. 11. The fluid pressures in the permeable material (top) are initially predominantly negative, reflecting the increasing dilation (Fig. 7, middle). In the critical state (after a shear strain value of 0.1), the fluid pressures fluctuate around the hydrostatic value (0 Pa).

Strain-rate dependency
Several studies have highlighted the importance of feedbacks between the solid and fluid phases during granular deformation (e.g., Iverson et al., 1994Iverson et al., , 1997bIverson et al., , 1998Pailha et al., 2008;Iverson, 2010;Rondon et al., 2011;Mutabaruka et al., 2014). A shear-rate dependence in a grain-fluid mixture can only originate from the fluid phase, since dry granular materials deform rate-independently under pseudo-static shear deformation (GDR-MiDi, 2004;Damsgaard et al., 2013). Rate dependence emerges, however, as soon as the flow of viscous pore fluids starts to influence the solid phase.
Water has a relatively low viscosity, which implies that the shear stress required to deform the fluid phase alone is extremely low. However, the fluid phase influences the bulk rheology if diffusion of fluid pressures is limited relative to volumetric forcing rates, as in a rapidly deforming but relatively impermeable porous material. The coupled particle-The Cryosphere, 9, 2183-2200, 2015 www.the-cryosphere.net/9/2183/2015/ fluid interactions cause the material to respond as a low-pass filter when forced with changes in volume and porosity. The reequilibration of pressure anomalies depends on the volumetric strain rate, water viscosity, and material permeability. Any forcing that affects local porosity causes the material to respond in part like a viscous dashpot due to internal fluid flow. The pore-fluid viscosity likely increases if suspended clayparticles or ice crystals are present. The exact rheology of such mixtures may very well be non-Newtonian, and we have therefore not included such effects in this study.

Dilatant hardening: effects on sediment strength and deformation depth
Transient porosity changes take place when granular materials deform before attaining the critical state (e.g., Reynolds, 1885;Mead, 1925;Schofield and Wroth, 1968;Iverson et al., 2000). These changes can cause deviation from the hydrostatic pressure distribution, which affects the magnitude of the local effective normal stress. Considering the Mohr-Coulomb constitutive relation for till rheology, a reduction in pore-water pressure increases the effective stress, which in turn strengthens the material in the shear zone ( Fig. 12) (e.g., Iverson et al., 1997bIverson et al., , 1998Moore and Iverson, 2002;Iverson, 2010). In our results we obtain insight into the micromechanical mechanisms of strengthening without explicitly specifying Mohr-Coulomb constitutive behavior. Deviations from the hydrostatic pressure distribution cause pressuregradient forces, which during dilation resist the development by pushing the grains together towards the low pressure of the shear zone (Fig. 13). The tangential strength of interparticle contacts in the DEM is determined by Coulomb friction (Eq. 3), which implies a linear correlation between con- igure 5.11. Temporal evolution (x-axis) of horizontally averaged uid pressures (y-axis). The permeable material (top) is able to uickly respond to internal volumetric changes, which are shortived and of small magnitude. The low-permeable material (botom) is dominated by large pressure reductions and relatively slow elaxation.
69 Figure 11. Temporal evolution (x axis) of horizontally averaged fluid pressures (y axis). The permeable material (top) is able to quickly respond to internal volumetric changes, which are shortlived and of small magnitude. The low-permeability material (bottom) is dominated by large pressure reductions and relatively slow relaxation.
tact normal force and tangential contact strength. Heavily loaded particle contacts are thus less likely to slip, and chains of particles with strong contacts cause increased resistance to deformation (Damsgaard et al., 2013). The fluid force on the particles strengthens the inter-particle contacts and increases the shear friction until hydrostatic pressure conditions are reestablished.
The dilatant hardening causes a shallower deformational profile in these experiments (Figs. 12,9 left,and 8 After Iverson et al. 1998, JoG Mohr-Coulomb plastic: |· | = C + ‡ 0 tan " where ‡ 0 = ‡ i ≠ Figure 12. Particle-fluid interaction during deformation of a consolidated sediment (after Iverson et al., 1998). Figure 13. Micromechanical cause of dilatant hardening. A consolidated sediment (top) is deformed with a vertical gradient in velocity. The grains are forced past each other in order to accommodate the shear strain. The deformation causes dilation, which increases porosity locally and decreases fluid pressure (bottom). The gradient in fluid pressure pulls particles together (Eq. 12), which increases the load on inter-particle contacts. The larger inter-particle normal stress increases the shear strength of the contact (Eq. 3) resulting in a stronger sediment. sure distribution is recovered at all depths of the material, and the deformational profiles progressively deepen until they match the dry experiment at higher shear-strain values ( 2). This progressive widening of the shear zone causes dilation to slowly increase while shear strength only displays very slight decrease (Figs. 4 and 7). The shallow deformation at low strains is consistent with the laboratory results by Iverson et al. (1997a), where the shear zone in the coarse-grained Storglaciären till in all cases was wider than the shear zone of the fine-grained Two Rivers till. The velocity profile of www.the-cryosphere.net/9/2183/2015/ The Cryosphere, 9, 2183-2200, 2015 the shear zone determines the material flux. A shallower deformation depth and a transiently lower subglacial sediment transport rate is thus to be expected from subglacial shearing of compacted, low-permeable sediments, relative to permeable counterparts. These results are consistent with observations of very shallow deformation in subglacial tills with a relatively low permeability (Engelhardt and Kamb, 1998;Piotrowski et al., 2004), although shallow deformation is also consistent with fine grain sizes (Tulaczyk, 1999), low effective stresses (Damsgaard et al., 2013), and rate-weakening ploughing at the ice-till interface (Thomason and Iverson, 2008).
Owing to the granularity of the material, the dilation displays small fluctuations around the critical-state value. The small volumetric oscillations create new fluid-pressure deviations from the hydrostatic value, which alternately slightly weaken or strengthen the sediment (Fig. 7, top). In cases where the shear stress is close to the sediment shear strength, the hardening may be sufficient to stabilize patches of the bed (Piotrowski, 1987), or the weakening may trigger slip (Goren et al., 2011).
Shear deformation can cause compaction instead of dilation if the material porosity exceeds the critical-state porosity, if elongated grains rotate into a more compact arrangement with high fabric strength, or if relatively soft granular components such as aggregates disintegrate. Compactioninduced weakening is a primary mechanism for debris-flow mobilization in landslides, as demonstrated experimentally at small scales  and at field scales (Iverson et al., 2000). The granular model applied here is not able to reproduce the shear-induced compaction that, e.g., clayrich materials can display during early shear due to clayfabric development (e.g., Dewhurst et al., 1996). Compaction causes increased pore-water pressure in the shear zone in cases where the volumetric change exceeds the timescale of pore-water pressure diffusion. Some of the effective normal stress on the shear zone is consequentially reduced, which in turn decreases the material strength. The reduction of strength due to compaction is rate dependent like the dilative hardening.
Our results confirm that the interplay between the solid and fluid phases can influence the sediment strength in a transient manner (e.g., Moore and Iverson, 2002;Iverson et al., 1998;Iverson, 2010). Pore-water pressures decrease during deformation, and shear strength increases until deformation ceases or the critical state is reached. Once the local and regional hydraulic system recovers from the pore-pressure reduction, the sediment strength is once again reduced and a new deformation phase may be initiated. The magnitude of strengthening is dictated by consolidation at times between slip events as well as the ability of the subglacial hydrological system to accommodate reductions in pressure at the icebed interface.
A variable shear strength of the till influences ice flow when the basal shear stress is in the range of till strength variability. Since surface slopes of ice streams are low, driving stresses tend to be low as well. Inferred values of driving stresses at the Northeast Greenland ice stream (Joughin et al., 2001), Whillans Ice Stream and ice plain (Engelhardt and Kamb, 1998), and Pine Island Glacier  lie within the range of 2 to 23 kPa (Alley and Whillans, 1991;Cuffey and Paterson, 2010) and are thus potentially sensitive to the variability in till strength. If the glacier moves with variable velocities in a stick-slip or surging manner, periods of stagnation may consolidate and strengthen the sediment, in effect delaying the following slip event (Iverson, 2010). Future studies will focus on investigating the mechanical consequences of transiently variable stresses, caused by changes in pore-water pressure or shear stress.

Conclusions
We numerically simulate a two-way coupled particle-fluid mixture under pseudo-static shear deformation and confirm results on transient strengthening from previous physical experiments (e.g., Moore and Iverson, 2002) and calculations (e.g., Iverson et al., 1998).
The grains are simulated individually by the discrete element method, while the fluid phase is treated as a compressible and slowly flowing fluid adhering to Darcy's law. The fluid influences the particles through local deviations from the hydrostatic pressure distribution. Due to the extremely low viscosity of water, the deformational behavior of dense granular material is governed by inter-grain contact mechanics. The porosity of a granular packing evolves asymptotically towards a critical-state value with increasing shear strain. Changes in porosity cause deviations from the hydrostatic pressure if the rate of porosity change exceeds the rate of pressure diffusion. The rate of pressure diffusion is governed by the fluid viscosity, the local porosity, and the hydraulic permeability. Low fluid pressures developing due to sediment dilation increase the frictional strength of intergrain contacts by increasing the loading between grain pairs. The magnitude of the strengthening effect is rate dependent and increases with shear velocity and decreases with increasing hydraulic permeability. The rheology is plastic but ratedependent dilative strengthening can contribute to the material strength during early stages of fast deformation. If the till is very porous or the deformation is accompanied by strong fabric development or grain crushing, compaction in the shear zone is expected to weaken the sediment, causing a rate weakening with increased shear rate until the excess pressures are reduced by hydraulic diffusion.
We furthermore show that for fast shear velocities permeable sediments are only weakly influenced by the fluid phase, resulting in little shear strengthening and a deep decimeter-scale deformation dictated by the normal stress and grain sizes. Impermeable and consolidated sediments display slight dilatant strengthening at high shear velocity.
The Cryosphere, 9, 2183-2200, 2015 www.the-cryosphere.net/9/2183/2015/ The strengthening causes deformation to focus at the ice-bed interface if pore-water pressures are higher there. The resultant depth of deformation is on the millimeter-to-centimeter scale. Actively deforming patches in the subglacial mosaic of deforming and stable spots act as sinks for meltwater. Additionally, sediment dilation can cause substantial thinning of a water film at the ice-bed interface. If the subglacial shearing movement halts, the sediment gradually weakens as the fluid pressure readjusts to the hydrostatic value. These temporal changes in sediment strength may explain observed variability in glacier movement.