manus_continuum_granular1

continuum-granular-manuscript1.tex (34083B)

---

 \documentclass[draft]{agujournal2019}
 % substitute draft for final \usepackage[utf8]{inputenc}
 
 \usepackage{url} %this package should fix any errors with URLs in refs.
 \usepackage{lineno}
 \usepackage[inline]{trackchanges} %for better track changes. finalnew option will compile document with changes incorporated.
 \usepackage{soul}
 \linenumbers
 %%%%%%%
 % As of 2018 we recommend use of the TrackChanges package to mark revisions.
 % The trackchanges package adds five new LaTeX commands:
 %
 %  \note[editor]{The note}
 %  \annote[editor]{Text to annotate}{The note}
 %  \add[editor]{Text to add}
 %  \remove[editor]{Text to remove}
 %  \change[editor]{Text to remove}{Text to add}
 %
 % complete documentation is here: http://trackchanges.sourceforge.net/
 %%%%%%%
 
 \usepackage{amsmath}
 
 \draftfalse
 
 \journalname{Geophysical Research Letters}
 
 \begin{document}
 
 %%% TITLE
 \title{Deep slip in subglacial till from water-pressure memory}
 \authors{Anders Damsgaard\affil{1}, Jenny Suckale\affil{1}, and Liran Goren\affil{2}}
 
 \affiliation{1}{Department of Geophysics, Stanford University}
 \affiliation{2}{Department of Geological and Environmental Sciences, Ben-Gurion University of the Negev}
 
 \correspondingauthor{Anders Damsgaard}{anders@adamsgaard.dk}
 
 \begin{keypoints}  % each keypoint max 140 chars
 \item Non-local continuum model of deforming grains is consistent with subglacial till rheology
 \item Sediment transport depends on the subglacial stress and sediment properties
 \item Diffusion of water-pressure perturbations generates deep shear zones within the till
 \end{keypoints}
 
 \begin{abstract}
 The dynamic interplay between ice flow, meltwater drainage, and basal sediment deformation is crucial for understanding glacier and ice-sheet behavior.
 Field observations indicate that subglacial sediment (till) layers often are mobile, and that a significant amount of the ice translation is contributed by till deformation.
 In turn, subglacial sediment transport constructs landforms that influence glacier stress balance, hydrology, and geomorphology.
 However, many numerical and theoretical models ignore till deformation and assume ice movement is contributed by interficial slip between ice and bed.
 Geotechnical analyses concluded that till deforms according to Mohr-Coulomb plasticity, which proved difficult to describe in models.
 Here we present a new continuum model that is both inspired by the granular nature of till layers and is consistent with Mohr-Coulomb mechanics, allowing modeling of the critical aspects of glacier-sediment-hydrology system.
 Our model results indicate that pulses in water pressure can shift deformation away from the ice-bed interface and far into the bed, causing episodes of significant till advection.
 Deep deformation is most likely in tills with relatively high hydraulic permeability, forced by long-lasting and large water-pressure variations.
 \end{abstract}
 
 \section*{Plain Language Summary}
 TODO
 
 
 \section{Introduction}
 Fast glacier and ice-sheet flow often ocurrs over weak sedimentary deposits. 
 While internal ice deformation is small, basal sediment deformation and slip at the ice-bed interface accounts for nearly all movement \cite <e.g.,>[]{Cuffey2010}.
 Basal sediments, called subglacial till, are diamictons commonly consisting of reworked sediments and erosional products \cite <e.g.,>[] {Evans2006}.
 Meltwater fully saturates the pore space, and variations in subglacial water pressure are common and can be caused by internal variability \cite<e.g.,>[] {Kavanaugh2009} or external water input \cite<e.g.,>[] {Andrews2014, Christoffersen2018}.
 In-situ field observations demonstrate that deformation of this layer can contribute significantly to the glacier movement \cite<e.g.,>[] {Boulton1979, Humphrey1993, Truffer2000}.
 \citeA{Boulton1987} argued that a viscous rheological model with mild stress non-linearity appropriately describes subglacial till deformation.
 A viscous rheology implies that the stress required to deform the till is strongly dependent on how fast it is deformed.
 However, \citeA{Kamb1991}, \citeA{Iverson1998}, and \citeA{Tulaczyk2000} demonstrated from laboratory shear tests that rate-independent Mohr-Coulomb plasticity, as common for sedimentary materials, is a far better rheological description for subglacial till.
 Mohr-Coulomb plastic materials have a yield strength that linearly scales with effective stress and is insensitive to strain rate.
 \citeA{Iverson2010} reviewed possible viscous contributions during till-water deformation, but deemed them to be of minor importance.
 In spite of a limited observational basis, viscous rheologies continued to be applied as they allow for mathematical modeling of till advection.
 Tills with viscous rheology were used to explain coupled ice-bed processes including subglacial sediment transport \cite<e.g.,>[] {Jenson1995}, landform formation \cite<e.g.,>[] {Hindmarsh1999, Fowler2000}, localization of water drainage \cite<e.g.,>[] {Walder1994, Ng2000b}, and ice-sheet behavior in a changing climate \cite<e.g.,>[] {Pollard2009}.
 Meanwhile, the Mohr-Coulomb plastic model continued to gain further empirical support from laboratory testing \cite<e.g.,>[] {Rathbun2008, Iverson2015}, as well as field observations on mountain glaciers \cite<e.g.,>[] {Hooke1997, Truffer2006, Kavanaugh2006, Iverson2007}, and ice sheets \cite<e.g.,>[] {Tulaczyk2006, Gillet-Chaulet2016, Minchew2016}.
 Inconveniently, the plastic rheology caused a deadlock for the typical continuum modeling of ice and till, as the Mohr-Coulomb constitutive model offers no direct relation between stress and strain rate.
 The deadlock was partially resolved when \citeA{Schoof2006} and \citeA{Bueler2009} showed that Mohr-Coulomb friction can be included in ice-sheet models through mathematical reguralization.
 While the methods describe the mechanical effect of the bed on the flowing ice, they offer no treatment of sediment erosion, transport, and deposition as strain in the sedimentary bed is not included.
 Yet, \citeA{Ritz2015} demonstrated that the description of basal friction is highly influential for future Antarctic ice-sheet contributions to global-mean sea level rise, where sliding over plastic beds is likely to increase future contributions to sea-level rise.
 
 \citeA{Damsgaard2013} and \citeA{Damsgaard2016} demonstrated that strain distribution and plasticity can be modeled in water-saturated tills by explicitly considering each sediment grain.
 Unfortunately, intense computational requirements associated with the grain-scale modeling entirely outrule model applicability for coupled ice-till simulations.
 Instead, simulation of landform to ice-sheet scale requires continuum models.
 In this study we build on continuum-modeling advances in granular mechanics and produce a model appropriate for water-saturated sediment deformation in the subglacial environment.
 We rely on the original model by \citeA{Henann2013} that was developed for critical state deformation of dry and cohesionless granular materials.
 To correctly account for subglacial tills that are water saturated and often contain a certain amount of cohesion that generally increases with clay content \cite<e.g.,>[] {Iverson1997}, we extend the \citeA{Henann2013} model by including the notion of pore-pressure controlled effective stress, pore-pressure diffusion and add strength contributions from cohesion.
 The resultant model contributes the methodological basis required for understanding the coupled dynamics of ice flow over deformable beds.
 Different from previous continuum models for till, our model remains true to rheological properties observed in laboratory and field settings.
 
 In the following, we present the \citeA{Henann2013} model and our modifications for modeling water-saturated and cohesive subglacial till.
 We discuss its applicability and technical limitations before comparing the simulated sediment behavior to published results from laboratory experiments on tills.
 The model produces rich dynamics that is consistent and could explain previously poorly-understood field observations.
 In particular, the model demonstrates its ability to produce deformation deep away from the ice-bed interface, as occasionally observed in field settings \cite{Truffer2000, Kjaer2006}.
 The model produces stick-slip dynamics under variable water pressures, as observed in mountain glaciers \cite<e.g.,>[] {Fischer1997} and ice streams \cite<e.g.,>[] {Bindschadler2003}. 
 The model source code is constructed with minimal external dependencies, is freely available, and is straight-forward to couple to models of ice-sheet dynamics and glacier hydrology.
 
 
 \section{Methods}%
 \label{sec:methods}
 
 Soils, tills, and other sediments are granular materials, consisting of discrete grains that interact with frictional losses.
 \citeA{GDR-MiDi2004} introduced a non-dimensional inertia number that summarizes the mechanical behavior of dry and dense granular deformation.
 This finding evolved into an empirical continuum rheology in \citeA{daCruz2005} and \citeA{Jop2006}, where stress and porosity depend on inertia in a non-linear manner.
 However, these continuum models are \emph{local}, meaning that local stresses determines the local strain-rate response alone.
 This means that material properties do not influence shear zone width.
 Granular deformation contains numerous non-local effects, where flow rates in neighboring areas influence the tendency of a sediment parcel to deform.
 Granular shear zones are an example of the non-locality as they have a minimum width dependent on grain characteristics \cite<e.g.,>[] {Nedderman1980, Forterre2008, Kamrin2018}.
 
 \subsection{Non-local granular fluidity (NGF) model}%
 \label{sub:granular_flow}
 \citeA{Henann2013} presented the non-local granular fluidity (NGF) model where a \emph{fluidity} field variable accounts for the non-local effects on deformation.
 The model builds on the previous continuum rheology for granular materials by \citeA{daCruz2005} and \citeA{Jop2006}, and with non-local effects accurately describes strain distribution in a variety of experimental settings.
 In the \citeA{Henann2013} model, all material is assumed to have a uniform porosity and be in the critical state.
 Fluidity acts as a state variable, describing the phase transition between non-deforming (jammed) and actively deforming (flowing) parts of the sediment \cite<e.g.,>[] {Zhang2017}.
 The modeled sediment deforms as a highly nonlinear Bingham material with yield beyond the Mohr-Coulomb failure limit \cite<e.g.,>[] {Henann2013, Henann2016}, but unlike classical plastic models it includes a closed form relation that predicts the stress-strain rate relation beyond yield.
 For the purposes in this paper, we assume that plastic shear strain ($\dot{\gamma}^\mathrm{p}$) contributes all of the resultant deformation ($\dot{\gamma}$):
 \begin{linenomath*}
 \begin{equation}
 	\dot{\gamma} \approx \dot{\gamma}^\mathrm{p} = g(\mu, \sigma_\mathrm{n}') \mu,
 	\label{eq:shear_strain_rate}
 \end{equation}
 \end{linenomath*}
 where $\mu = \tau/\sigma_\mathrm{n}'$ is the dimensionless ratio between shear stress ($\tau$ [Pa]) and effective normal stress ($\sigma_\mathrm{n}' = \sigma_\mathrm{n} - p_f$ [Pa]).
 Water pressure is $p_\mathrm{f}$ [Pa] and $g$ [s$^{-1}$] is the granular fluidity.
 The fluidity $g$ is a kinematic variable governed by grain velocity fluctuations and packing fraction \cite{Zhang2017}, and consists of local and non-local components:
 \begin{linenomath*}
 \begin{equation}
 	\nabla^2 g = \frac{1}{\xi^2(\mu)} (g - g_\mathrm{local}),
 	\label{eq:g}
 \end{equation}
 \end{linenomath*}
 The above Poisson-type equation acts to distribute strain in space according to material properties and non-local stress state.
 The degree of non-locality is scaled by the cooperativity length $\xi$, which, in turn, scales with nonlocal amplitude $A$ [-]:
 \begin{linenomath*}
 \begin{equation}
 	\xi(\mu) = \frac{Ad}{\sqrt{|(\mu - C/\sigma_\mathrm{n}') - \mu_\mathrm{s}|}},
 	\label{eq:cooperativity}
 \end{equation}
 \end{linenomath*}
 where $d$ [m] is the representative grain diameter, $\mu_\mathrm{s}$ [-] is the static Coulomb yield coefficient, and $C$ [Pa] is the material cohesion.
 The local contribution to fluidity is defined as:
 \begin{linenomath*}
 \begin{equation}
 	g_\mathrm{local}(\mu, \sigma_\mathrm{n}') =
 	\begin{cases}
 	    \sqrt{d^2 \sigma_\mathrm{n}' / \rho_\mathrm{s}} ((\mu - C/\sigma_\mathrm{n}') - \mu_\mathrm{s})/(b\mu) &\mathrm{if}\quad \mu - C/\sigma_\mathrm{n}' > \mu_\mathrm{s} \mathrm{,}\quad\mathrm{and}\\
 	    0 & \mathrm{if}\quad \mu - C/\sigma_\mathrm{n}' \leq \mu_\mathrm{s}.
 	\end{cases}
 	\label{eq:g_local}
 \end{equation}
 \end{linenomath*}
 where $\rho_\mathrm{s}$ is grain mineral density, and $b$ [-] controls the non-linear rate dependence beyond yield.
 The failure point is principally determined by the Mohr-Coulomb constituent relation in the conditional of Eq.~\ref{eq:g_local}.
 However, the non-locality in Eq.~\ref{eq:g} infers that deformation can occur in places that otherwise would not fail, in cases where the surrounding areas have a high local fluidity.
 
 \subsection{Fluid-pressure evolution}%
 \label{sub:fluid_pressure_evolution}
 We prescribe the transient evolution of pore-fluid pressure ($p_\mathrm{f}$) by Darcian pressure diffusion \cite<e.g.,>[] {Goren2010, Goren2011, Damsgaard2017b}:
 \begin{linenomath*}
 \begin{equation}
 	\frac{\partial p_\mathrm{f}}{\partial t} = \frac{1}{\phi\eta_\mathrm{f}\beta_\mathrm{f}} \nabla \cdot (k \nabla p_\mathrm{f}),
 	\label{eq:p_f}
 \end{equation}
 \end{linenomath*}
 where $\eta_\mathrm{f}$ denotes dynamic fluid viscosity [Pa s], $\beta_\mathrm{f}$ is adiabatic fluid compressibility [Pa$^{-1}$], and $k$ is intrinsic permeability [m$^2$].
 The sediment is assumed to be in the critical state throughout the domain, as in the original formulation by \citeA{Henann2013}.
 This means that porosity and fluid pressure does not change as a function of granular deformation.
 The local fluid pressure is used to determine the effective normal stress used in the granular flow calculations (Eq.~\ref{eq:shear_strain_rate} and~\ref{eq:g_local}).
 
 \subsection{Limitations of the continuum model}%
 \label{sub:limitations_of_the_continuum_model}
 The presented model considers the material to be in the critical (steady) state throughout the domain.
 Consequently, porosity is prescribed as a constant and material-specific parameter.
 For that reason the model is not able to simulate uniaxial compaction or shear-induced volume changes \cite<e.g.,>[] {Iverson2000, Iverson2010b, Damsgaard2015} or compaction and dilation \cite<e.g.,>[] {Dewhurst1996}.
 We currently have a transient granular continuum model with state-dependent porosity under development.
 However, \citeA{Iverson2010} argued that the majority of actively deforming subglacial sediment may be in the critical state.
 For that reason, and due to the emerging insights from the current simple model, we see this contribution as a valuable first step that allows us to isolate the dynamic that emerges in the critical state.
 
 In the NGF model, the representative grain size $d$ scales the non-locality and strain distribution.
 The relation between the physical grain size distribution and the effective grain size $d$ that controls the fluidity has so far not been sufficiently explored.
 For the case of till layers, that are characterized by a fractal grain size distribution \cite{Hooke1995}, the relation to the fluidity is even more obscure.
 The issue is left for future research and here we assume that $d$ corresponds to the volumetrically dominant grain size, if one exists.
 Specifically designed laboratory experiments with various tills should inform the treatment of length scale, outside of ploughing effects by large clasts protruding from the basal ice \cite<e.g.,>[] {Tulaczyk1999}.
 
 \subsection{Simulation setup}
 We apply the model in a one-dimensional setup with simple shear (Fig.~\ref{fig:validation}a).
 Parameter values are listed in Table~S1.
 The lower boundary condition for the granular phase is no slip ($v_x(L=0) = 0$).
 The upper boundary condition for the granular phase is fixed shear friction $\mu(z=L_z)$ under stress controlled settings, or fixed shear velocity $v_x(z=L_z)$ for rate-controlled experiments.
 The upper normal stress ($\sigma_\mathrm{n}(z=L_z)$) is constant, and normal stress linearly increases with depth due to material weight.
 Effective normal stress ($\sigma_\mathrm{n}' = \sigma_\mathrm{n} - p_\mathrm{f}$) varies if water pressure $p_\mathrm{f}$ changes.
 For the water-pressure solver, the top pressure ($p_\mathrm{f}(z=L_z)$) is either constant or varied through time.
 The water pressure is set to follow the hydrostatic gradient at the lower boundary ($dp_\mathrm{f}/dz(z=0) = \rho_\mathrm{f}G$).
 For the experiments with variable water pressure, we apply a water-pressure forcing amplitude of $A_\mathrm{f}$ = 80 kPa that modulates effective normal stress at the top around 100 kPa.
 Unless noted, we use the value 0.94 for the dimensionless rate-dependence parameter $b$, which is impirically constrained from laboratory experiments on glass beads \cite{Forterre2003, Henann2016}.
 Under stress-controlled configurations, the model formulation produces a shear strain rate from sediment properties and applied stress.
 Rate-controlled conditions are implemented by iteratively adjusting the applied stress until the output shear velocity matches the desired value (full details in Text S1.1).
 Many simulations are performed under both stress- and rate-controlled shear, which both idealize the driving glacier physics.
 Real glacier settings fall somewhere in between, depending on how important basal friction is to the overall stress balance.
 Stress-controlled conditions approximate a setting where ice flow directly responds to changes in subglacial strain rates.
 Whillans Ice Plain, West Antarctica is an example of this setting, where a low surface slope and low driving stress results in stick-slip movement \cite<e.g.,>[] {Bindschadler2003}.
 A rate-controlled setup is the opposite end member, where changes in bed friction do not influence ice flow velocity.
 
 \begin{figure*}[htbp]
 	\begin{center}
 		\includegraphics[width=0.49\textwidth]{experimental-setup.pdf}
 		\includegraphics[width=0.49\textwidth]{experiments/fig-rate_dependence.pdf}\\
 		\includegraphics[width=0.49\textwidth]{experiments/fig-mohr_coulomb.pdf}
 		\includegraphics[width=0.49\textwidth]{experiments/fig-strain_distribution.pdf}
 		\caption{\label{fig:validation}%
 			Comparison between shear experiments on subglacial till and the non-local granular fluidity (NGF) model applied in this study.
 			a) Experimental setup for the NGF model in one-dimensional shear.
 			The upper boundary is constant friction $\mu$ (stress controlled), or constant velocity $v_x$ (velocity controlled).
 			b) Rate dependence of critical-state friction in laboratory experiments on till (after \citeA{Iverson2010}), and the NGF model with material friction $\mu_\mathrm{s}$ = 0.5 and effective normal stress $\sigma_\mathrm{n}'$ = 100 kPa.
 			The dimensionles parameter $b$ controls frictional rate dependence (Eq.~\ref{eq:g_local}), where $b$ = 0.94 is appropriate for glass beads \cite{Henann2016}.
 			c) Mohr-Coulomb analysis of till samples and NGF model.
 			d) Modeled strain distribution under varying effective normal stress ($\sigma_\mathrm{n}'$) with the discrete-element method (DEM, \citeA{Damsgaard2013}) and the NGF model.
 		}
 	\end{center}
 \end{figure*}
 
 
 
 \section{Results}%
 \label{sec:results}
 We first compare the modeled mechanical behavior to various tills tested in laboratory settings.
 Over five orders of strain-rate magnitude, some tills show slight rate weakening and others are slightly rate strengthening (Fig.~\ref{fig:validation}b).
 Shear-strain rates up to $\sim5 \times 10^3$ a$^{-1}$ are realistic for natural glacier systems \cite{Cuffey2010}.
 The model is effectively rate-independent over most of the range, but higher values of $b$, which determines the degree of rate dependence, provide larger frictional resistance at extreme shear-strain rates (Fig.~\ref{fig:validation}b).
 The modeled friction value can be linearly scaled by adjusting $\mu_\mathrm{s}$ in Eqs.~\ref{eq:cooperativity} and~\ref{eq:g_local}, and the model can simulate any combination of effective friction (or friction angle $\varphi = \tan^{-1}(\mu_s)$) and cohesion (Fig.~\ref{fig:validation}c).
 
 The NGF model contains parameter $A$ for adjusting the degree of material non-locality (Eq.~\ref{eq:cooperativity}).
 However, at present no laboratory experiment exists in the literature where the effects of normal stress are analysed for changes in strain distribution in the till.
 Instead, we compare the modeled strain distribution with discrete-element results from \citeA{Damsgaard2013}.
 By inserting relevant material parameters for grain size, friction, stress, and shear velocity (DEM, Table~S1), the NGF model approximates the strain distribution well (Fig.~\ref{fig:validation}d).
 Both models show that sediment advection is pressure dependent, with low effective normal stresses producing shallow deformation, and high effective normal stresses deepening the material mobilization.
 The DEM results took more than two months of computational time, whereas the continuum model is completed in a fraction of a second, albeit without detail of individual particle kinematics and dynamical adjustment towards the critical state.
 
 Next, we vary the water pressure at the top in a sinusoidal manner and observe the resultant shear dynamics over a simulation time of seven days (Fig.~S1).
 The experiments are performed in both stress and rate-controlled configurations.
 Figure~\ref{fig:stick_slip_depth} shows depth variations in the rate-controlled configuration through a day of simulation time, which corresponds to a single wavelength of the water-pressure forcing.
 The effective normal stress generally increases with depth according to the difference in grain and fluid density (Fig.~\ref{fig:stick_slip_depth}a).
 However, water pressure variations at the ice-bed interface can reverse this depth trend, and create minima in effective normal stress away from the ice-bed interface.
 Due to diffusion, water pressure perturbations decay exponentially with depth and travel with a phase shift (Fig.~\ref{fig:stick_slip_depth}b).
 Maximum shear deformation moves downwards into the bed as effective normal stress at the top boundary increases (Fig.~\ref{fig:stick_slip_depth}c).
 Granular failure generally occurs where effective normal stress is at its minimum, resulting in a plug-like flow (Fig.~\ref{fig:stick_slip_depth}c,d) during reversal of the depth trend in effective normal stress (Fig.~\ref{fig:stick_slip_depth}b).
 
 Over the entire simulation period, the system shows stick-slip behavior under stress-controlled conditions where velocities range from 0 to $\sim$9 km/d with strong hysteresis (Fig.~\ref{fig:hysteresis}a).
 The response during the first cycle ($t<1$ d) is slightly different from later cycles ($t>1$ d) as the model is initialized with a hydrostatic water-pressure distribution.
 Under rate-controlled conditions, the shear stress varies linearly as predicted for Mohr-Coulomb materials but with hysteresis at high effective normal stresses (Fig.~\ref{fig:hysteresis}d).
 The shear stress is higher during deep shear as the lithostatic contribution increases effective normal stress at depth.
 Both driving modes produce deep deformation during periods where water-pressure at the ice-bed interface is at its lowest magnitude.
 Under stress-controlled conditions, the deep deformation produces an insignificant amount till flux (Fig.~\ref{fig:hysteresis}c), as the shear stress at this time as the shear stress is insufficient for slip (Fig.~\ref{fig:hysteresis}a).
 Instead, the majority of sediment transport occurs as shallow deformation during rapid slip events and high water pressures (Fig.~\ref{fig:hysteresis}b,c).
 On the contrary, the largest sediment transport under rate-controlled conditions occurs with low water pressures and high effective normal-stresses at the ice-bed interface (Fig.~\ref{fig:hysteresis}e,f).
 Pulse perturbations of various shape in water pressure are also able to cause deep deformation, and the maximum deformation depth increases with increasing perturbation amplitude (Fig.~S2).
 
 \begin{figure*}[htbp]
 	\begin{center}
 		\includegraphics[width=15.0cm]{experiments/fig-stick_slip_rate_depth.pdf}
 		\caption{\label{fig:stick_slip_depth}%
 			Pore-pressure diffusion and strain distribution with depth with a sinusoidal water-pressure forcing from the top.
 			The bed thickness is $L_z$ = 8 m, but here only the upper 4 m are shown.
 			The water-pressure forcing has a daily periodocity, and plot lines are one hour apart in simulation time.
 			The full horizontal line marks the skin depth from Eq.~\ref{eq:skin_depth}, and the dashed line marks the expected maximum deformation depth from Eq.~\ref{eq:max_depth}.
 		}
 	\end{center}
 \end{figure*}
 
 \begin{figure}[htbp]
 	\begin{center}
 		\includegraphics[width=0.49\textwidth]{experiments/fig-hysteresis_stress.pdf}
 		\includegraphics[width=0.49\textwidth]{experiments/fig-hysteresis_rate.pdf}
 		\caption{\label{fig:hysteresis}%
 			Stress, velocity, and sediment flux under sinusoidal forcing in water pressure and effective normal stress at the ice-bed interface.
 			a-c) is under stress-controlled shear, d-f) is under rate-controlled shear.
 			Black arrows denote the temporal evolution.
 		}
 	\end{center}
 \end{figure}
 
 
 \section{Discussion}%
 \label{sec:discussion}
 This study has the specific aim of quantifying advective sediment transport under shear.
 As granular deformation is associated with finite length scales, it is crucial to include non-local terms in the granular model equations instead of applying simpler \emph{local} sediment rheology models \cite<e.g.,>[] {daCruz2005, Jop2006, Forterre2008}.
 However, the modeled sediment flux presented here may present an upper bound since we assume that there is a strong mechanical coupling between ice and bed.
 Overpressurization and slip at the ice-bed interface may cause episodic decoupling and reduce bed deformation.
 Interface slip is observed both under contemporary ice streams \cite<e.g.,>[] {Engelhardt1998}, and in deposits from past glaciations \cite<e.g.,>[] {Piotrowski2001}.
 A complete understanding of subglacial sediment transport requires further research of ice-bed interface physics.
 Furthermore, we assume that the onset of deformation is not affected by initial hardening caused by shear-zone dilation \cite<e.g.,>[] {Iverson1998, Moore2002, Damsgaard2015}.
 
 In simulations without water-pressure variations, the sediment transport is dependent on the magnitude of the effective normal stress (Fig.~\ref{fig:validation}d).
 Such a dependence is a prerequisite for instability theories of subglacial landform development \cite{Hindmarsh1999, Fowler2000, Schoof2007, Fowler2018}.
 From geometrical considerations, it is likely that up-ice sloping stoss sides of subglacial topography experience higher bed-normal stress than down-ice sloping lee sides.
 With all other physical conditions being equal, our results indicate that shear-driven sediment advection would be larger on the stoss side than on the lee side.
 Topography of non-planar ice-bed interfaces (proto-drumlins, ribbed moraines, etc.) may be transported and modulated through this spatially variable transport capacity, unless stress differences are overprinted by variations in water pressure \cite<e.g.,>[] {Sergienko2013, McCracken2016, Iverson2017b, Hermanowski2019b}.
 
 At depth, the water pressure variations display exponential decay in amplitude and increasing lag.
 As long as fluid and diffusion properties are constant and the layer is sufficiently thick, an analytical solution to skin depth $d_\mathrm{s}$ [m] in our system follows the form \cite<e.g.,>[] {Turcotte2002},
 \begin{linenomath*}
 \begin{equation}
 	d_\mathrm{s}
 	= \sqrt{ \frac{D P}{\pi} }
 	= \sqrt{ \frac{k}{\phi\eta_\mathrm{f}\beta_\mathrm{f}\pi f} },
 	\label{eq:skin_depth}
 \end{equation}
 \end{linenomath*}
 where $D$ is the hydraulic diffusivity [m$^2$/s] and $P$ [s] is the period of the oscillations.
 The remaining terms were previously defined.
 Figure~\ref{fig:skin_depth}a shows the skin depth for water at 0$^\circ$C under a range of permeabilities and forcing frequencies.
 However, as the deformation pattern depends on both hydraulic properties and the pressure-perturbation amplitude (Fig.~S2,~S4), the skin depth alone is insufficient to judge the occurence of deep deformation.
 Therefore, we constrain an analytical solution for diffusive pressure perturbation to find the largest depth $z' = L_z - z$ that contains a minimum in effective normal stress over the cause of a pressure-perturbation cycle (see Supplementary Information Text S2 for full derivation):
 \begin{linenomath*}
 \begin{equation}
 	0 =
 	\sin \left( \frac{3\pi}{2} - \frac{z'}{d_\mathrm{s}} \right)
 	+ \cos \left( \frac{3\pi}{2} - \frac{z'}{d_\mathrm{s}} \right)
 	+ \frac{(\rho_\mathrm{s} - \rho_\mathrm{f}) G d_\mathrm{s}}{A_\mathrm{f}}
 	\exp \left( \frac{z'}{d_\mathrm{s}} \right).
 	\label{eq:max_depth}
 \end{equation}
 \end{linenomath*}
 When using the same stress, hydraulic parameters, and pressure-perturbation amplitude, the analytical solution matches our NGF simulations well (dashed horizontal line in Fig.~\ref{fig:stick_slip_depth}).
 
 \begin{figure}[htbp]
 	\begin{center}
 		\includegraphics[width=15cm]{experiments/fig-skin_depth.pdf}
 		\caption{\label{fig:skin_depth}%
 			a) Skin depth of pore-pressure fluctuations (Eq.~\ref{eq:skin_depth}) with forcing frequencies ranging from yearly to hourly periods.
 			The permeability spans values common for tills \cite{Schwartz2003}
 			b) Depth of maximum deformation with a water-pressure forcing of $A_\mathrm{f}$ = 10 kPa.
 			c) Same as (b) but with $A_\mathrm{f}$ = 100 kPa.
 		}
 	\end{center}
 \end{figure}
 
 Figure~\ref{fig:skin_depth}b and~c show the maximum expected deformation depth from solutions to Eq.~\ref{eq:max_depth} through pressure perturbations enforced from the ice-bed interface.
 Deep deformation occurs when the combination of forcing amplitude ($A_\mathrm{f}$) and skin depth (i.e.\ the $k/f$ ratio) is optimal.
 A combination of highly permeability and slow water-pressure forcing frequency results in a large skin depth (Fig.~\ref{fig:skin_depth}a), but the rapid pore-pressure diffusion is unlikely to overcome the lithostatic stress increase with depth. 
 Conversely, low skin depths (i.e.\ relatively impermeable materials and/or fast water-pressure forcing frequencies) are associated with limited penetration distance into the bed, limiting the maximum deformation depth.
 For that reason, coarse and relatively permeabile tills are more susceptible to deep deformation due to larger associated skin depths (Fig.~\ref{fig:skin_depth}a), as long as the perturbation amplitude is sufficiently large (Fig.~\ref{fig:skin_depth}b,c).
 Fine-grained and relatively impermeable tills are unlikely to undergo deep deformation as pressure perturbations propagate at slower velocities.
 Albeit at shallower depths, deformation is still expected to occasionally occur away from the ice-bed interface in poor hydraulic conductors (Fig.~\ref{fig:skin_depth}b,c).
 
 Other water-pressure forcings may be additional mechanisms for deep deformation if they cause minima in effective stress at depth.
 For example, lateral water input from lake drainage or hydrological rerouting may also create episodes of deep slip, in particular when horizontal bedding decreases vertical permeability \cite<e.g.,>[] {Kjaer2006}.
 Water-escape structures are commonly observed in till deposits \cite<e.g.,> {van_der_Meer2009, Knight2015, Salamon2016, Ingolfsson2016}, which indicates that overpressurization of water within till beds may be extremely common.
 Since tills fail at their weakest spot, we suggest that subglacial deformation may not only be a patchy mosaic of deforming spots in the horizontal plane \cite<e.g.,> {Piotrowski2004}, but also at depth.
 
 
 \section{Conclusion}%
 \label{sec:conclusion}
 We present a new model for coupled computation of subglacial till and water.
 The model is based on the concept of non-local granular fluidity \cite{Henann2013}, but is extended with cohesion and pore-pressure diffusion.
 The mechanics adhere to Mohr-Coulomb plasticity, with a weak and non-linear rate dependence governed by stress and sediment properties.
 The material is effectively rate-independent at glacial shear velocities, generally in agreement with laboratory results.
 We adapt a simple shear experimental setup for analyzing the mechanical response under different stresses and water-pressure variations.
 With cyclical water-pressure variations at the ice-bed interface, deep deformation occurs when remnant high water pressures at depth overcome the effective lithostatic gradient.
 Deep deformation may be common in coarse-grained subglacial tills with strong annual water-pressure differences.
 Similarly, sudden water-pressure pulses are powerful drivers for single events of deep deformation.
 
 
 \acknowledgments
 A.D. benefited from conversations with Dongzhuo Li, Indraneel Kasmalkar, Jason Amundson, Martin Truffer, Dougal Hansen, and Lucas Zoet during model development.
 Analysis and visualization of model output was performed with Gnuplot.
 
 \section*{Appendix}%
 \label{sec:appendix}
 The grain-water model is written in C and is available under free-software licensing at \url{https://src.adamsgaard.dk/1d_fd_simple_shear}.
 All simulation parameters can be specified as command-line arguments.
 Run \texttt{./1d\_fd\_simple\_shear -h} after compiling for usage information.
 The results and figures in this paper can be reproduced by following the instructions in the experiment repository for this publication, available at \url{https://src.adamsgaard.dk/manus_continuum_granular1_exp}.
 
 \bibliography{BIBnew}
 
 \end{document}