manus_continuum_granular1

manuscript files for first continuum-till paper
git clone git://src.adamsgaard.dk/manus_continuum_granular1
Log | Files | Refs Back to index

commit 69590bb60c6d73629ec7ff89a29d3ec8567a1624
parent 66f60b572268222dc38647e01554fe839187c4f9
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Tue,  3 Dec 2019 09:53:27 +0100

Attempt to fix issues with AGU templates

Diffstat:
MBIBnew.bib | 24+++++++++++++++---------
Mcontinuum-granular-manuscript1.tex | 140+++++++++++++++++++++++++++++++++++++++----------------------------------------
2 files changed, 84 insertions(+), 80 deletions(-)

diff --git a/BIBnew.bib b/BIBnew.bib @@ -1215,7 +1215,7 @@ @Article{Iverson2010, Title = {Shear resistance and continuity of subglacial till: hydrology rules}, - Author = {Iverson, N. R.}, + Author = {N. R. Iverson}, Journal = {J. Glaciol.}, Year = {2010}, Number = {200}, @@ -1225,7 +1225,7 @@ @article{Iverson2010-2, title={Effects of soil aggregates on debris-flow mobilization: results from ring-shear experiments}, - author={Iverson, N. R. and Mann, J. E. and Iverson, R. M.}, + author={N. R. Iverson and J. E. Mann and R. M. Iverson}, journal={Engineering Geology}, volume={114}, number={1}, @@ -1251,7 +1251,7 @@ @Article{Iverson1999b, Title = {Coupling between a glacier and a soft bed: {II}. {M}odel results}, - Author = {Iverson, N. R.}, + Author = {N. R. Iverson}, Journal = {J. Glaciol.}, Year = {1999}, Number = {149}, @@ -1261,7 +1261,7 @@ @Article{Iverson1997, Title = {A ring-shear device for the study of till deformation: tests on tills with contrasting clay contents}, - Author = {Iverson, N. R. and Baker, R. W. and Hooyer, T. S.}, + Author = {N. R. Iverson and R. W. Baker and T. S. Hooyer}, Journal = {Quat. Sci. Rev.}, Year = {1997}, Number = {9}, @@ -1272,7 +1272,7 @@ @Article{Iverson1998, Title = {Ring-shear studies of till deformation: {C}oulomb-plastic behavior and distributed strain in glacier beds}, - Author = {Iverson, N. R. and Hooyer, T. S. and Baker, R. W.}, + Author = {N.R. Iverson and T. S. Hooyer and R. W. Baker}, Journal = {J. Glaciol.}, Year = {1998}, Pages = {634--642}, @@ -4499,7 +4499,7 @@ @article{Iverson2007, title={{Soft-bed experiments beneath Engabreen, Norway: regelation infiltration, basal slip and bed deformation}}, - author={Iverson, N. R. and Hooyer, T. S. and Fischer, U. H. and Cohen, D. and Moore, P. L. and Jackson, M. and Lappegard, G. and Kohler, J.}, + author={N. R. Iverson and T. S. Hooyer and U. H. Fischer and D. Cohen and P. L. Moore and M. Jackson and G. Lappegard and J. Kohler}, journal={J. Glaciol.}, volume={53}, number={182}, @@ -8746,7 +8746,7 @@ Winton and A. T. Wittenberg and F. Zeng and R. Zhang and J. P. Dunne}, volume = {122}, number = {12}, pages = {2302--2323}, - author = {N. R. Iverson and R. G. McCracken and L. K. Zoet and {\'{I}}. Ö. Benediktsson and A. Schomacker and M. D. Johnson and J. Woodard}, + author = {N. R. Iverson and R. G. McCracken and L. K. Zoet and {\'{I}}. \"{O} Benediktsson and A. Schomacker and M. D. Johnson and J. Woodard}, title = {A Theoretical Model of Drumlin Formation Based on Observations at M{\'{u}}lajökull, Iceland}, journal = {J. Geophys. Res.: Earth Surf.} } @@ -8926,7 +8926,7 @@ Winton and A. T. Wittenberg and F. Zeng and R. Zhang and J. P. Dunne}, @article{Iverson2000, doi={10.1126/science.290.5491.513}, year={2000}, - author={Iverson, R. M. and Reid, M. E. and Iverson, N. R. and LaHusen, R. G. and Lo- gan, M. and Mann, J. E. and Brien, D. L.}, + author={R. M. Iverson and M. E. Reid and N. R. Iverson and R. G. LaHusen and M. Logan and J. E. Mann and D. L. Brien}, volume={290}, number=5491, pages={513-–516}, @@ -9007,7 +9007,7 @@ Winton and A. T. Wittenberg and F. Zeng and R. Zhang and J. P. Dunne}, pages = {221--230}, author = {R. C. A. Hindmarsh}, title = {Coupled ice-till dynamics and the seeding of drumlins and bedrock forms}, - journals = {Ann. Glaciol.} + journal = {Ann. Glaciol.} } @article{Morey1984, doi = {10.1016/0165-232x(84)90048-x}, @@ -9211,3 +9211,9 @@ https://doi.org/10.1680/geot.1996.46.2.197 } https://doi.org/ https://doi.org/ +https://doi.org/10.7717/peerj.929/supp-1 + +@misc{1, + doi = {10.7717/peerj.929/supp-1}, + publisher = {{PeerJ}} +} diff --git a/continuum-granular-manuscript1.tex b/continuum-granular-manuscript1.tex @@ -1,7 +1,6 @@ \documentclass[draft]{agujournal2019} % substitute draft for final \usepackage[utf8]{inputenc} -\documentclass[draft]{agujournal2019} \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. @@ -20,6 +19,8 @@ % complete documentation is here: http://trackchanges.sourceforge.net/ %%%%%%% +\usepackage{amsmath} + \draftfalse \journalname{Geophysical Research Letters} @@ -35,7 +36,6 @@ \correspondingauthor{Anders Damsgaard}{anders@adamsgaard.dk} - \begin{keypoints} % each keypoint max 140 chars \item Non-local continuum model is consistent with subglacial till rheology \item Sediment advection depends on the current stress and the history of water-pressure forcings @@ -57,39 +57,37 @@ Deep deformation is most likely in tills with relatively high hydraulic permeabi TODO -\section{Introduction}% -\label{sec:introduction} - +\section{Introduction} Fast glacier and ice-sheet flow often ocurrs over weak sedimentary deposits, where basal slip accounts for nearly all movement (e.g., \citeA{Cuffey2010}). Basal sediments, called subglacial till, are diamictons commonly consisting of reworked sediments and erosional products (e.g., \citeA{Evans2006}). Meltwater fully saturates the pore space, and variations in subglacial water pressure are common and can be caused by internal variability (e.g., \citeA{Kavanaugh2009}) or external water input (e.g., \citeA{Andrews2014, Christoffersen2018}). In-situ field observations demonstrate that deformation of this layer can contribute significantly to the glacier movement (e.g., \citeA{Boulton1979, Humphrey1993, Truffer2000}). -\citet{Boulton1987} argued that a viscous rheological model with mild stress non-linearity appropriately describes subglacial till deformation. +\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, \citet{Kamb1991}, \citet{Iverson1998}, and \citet{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. +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. -\citet{Iverson2010} reviewed possible viscous contributions during till-water deformation, but deemed them to be of minor importance. +\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 (e.g., \citeA{Jenson1995}), landform formation (e.g., \citeA{Hindmarsh1999, Fowler2000}), localization of water drainage (e.g., \citeA{Walder1994, Ng2000b}), and ice-sheet behavior in a changing climate (e.g., \citeA{Pollard2009}). -Meanwhile, the Mohr-Coulomb plastic model continued to gain further empirical support from laboratory testing (e.g., \citeA{Rathbun2008, Iverson2015}), as well as field observations on mountain glaciers (e.g., \citeA{Hooke1997, Truffer2006, Iverson2007}), mountain glaciers (e.g., \citeA{Kavanaugh2006}), and ice sheets \citep[e.g.,][]{Tulaczyk2006, Gillet-Chaulet2016, Minchew2016}. +Meanwhile, the Mohr-Coulomb plastic model continued to gain further empirical support from laboratory testing (e.g., \citeA{Rathbun2008, Iverson2015}), as well as field observations on mountain glaciers (e.g., \citeA{Hooke1997, Truffer2006, Iverson2007}), mountain glaciers (e.g., \citeA{Kavanaugh2006}), and ice sheets (e.g., \citeA{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 \citet{Schoof2006} and \citet{Bueler2009} showed that Mohr-Coulomb friction can be included in ice-sheet models through mathematical reguralization. +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, \citet{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. +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. -\citet{Damsgaard2013} and \citet{Damsgaard2016} demonstrated that strain distribution and plasticity can be modeled in water-saturated tills by explicitly considering each sediment grain. +\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 \citet{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 (e.g., \citeA{Iverson1997}), we extend the \citet{Henann2013} model by including the notion of pore-pressure controlled effective stress, pore-pressure diffusion and add strength contributions from cohesion. +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 (e.g., \citeA{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 \citet{Henann2013} model and our modifications for modeling water-saturated and cohesive subglacial till. +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 \citep{Truffer2000, Kjaer2006}. +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 (e.g., \citeA{Fischer1997}) and ice streams (e.g., \citeA{Bindschadler2003}). Remnants of pressure deviations within the glacier bed cause hysteresis in stress and strain. 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. @@ -99,32 +97,32 @@ The model source code is constructed with minimal external dependencies, is free \label{sec:methods} Soils, tills, and other sediments are granular materials, consisting of discrete grains that interact with frictional losses. -\citet{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 \citet{daCruz2005} and \citet{Jop2006}, where stress and porosity depend on inertia in a non-linear manner. +\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 the spatially local state determines the local strain-rate response alone. 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 (e.g., \citeA{Nedderman1980, Forterre2008, Kamrin2018}). \subsection{Non-local granular fluidity (NGF) model}% \label{sub:granular_flow} -\citet{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 \citet{daCruz2005} and \citet{Jop2006}, and with non-local effects accurately describes strain distribution in a variety of experimental settings. -In the \citet{Henann2013} model, all material is assumed to have a uniform porosity and be in the critical state. +\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 (e.g., \citeA{Zhang2017}). The modeled sediment deforms as a highly nonlinear Bingham material with yield beyond the Mohr-Coulomb failure limit (e.g., \citeA{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}^\text{p}$) contributes all of the resultant deformation ($\dot{\gamma}$): +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}^\text{p} = g(\mu, \sigma_\text{n}') \mu, + \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_\text{n}'$ is the dimensionless ratio between shear stress ($\tau$ [Pa]) and effective normal stress ($\sigma_\text{n}' = \sigma_\text{n} - p_f$ [Pa]). -Water pressure is $p_\text{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 \citep{Zhang2017}, and consists of local and non-local components: +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_\text{local}), + \nabla^2 g = \frac{1}{\xi^2(\mu)} (g - g_\mathrm{local}), \label{eq:g} \end{equation} \end{linenomath*} @@ -132,38 +130,38 @@ The above Poisson-type equation acts to distribute strain in space according to The degree of non-locality is scaled by the cooperativity length $\xi$, wich, in turn, scales with nonlocal amplitude $A$ [-]: \begin{linenomath*} \begin{equation} - \xi(\mu) = \frac{Ad}{\sqrt{|(\mu - C/\sigma_\text{n}') - \mu_\text{s}|}}, + \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_\text{s}$ [-] is the static Coulomb yield coefficient, and $C$ [Pa] is the material cohesion. +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_\text{local}(\mu, \sigma_\text{n}') = + g_\mathrm{local}(\mu, \sigma_\mathrm{n}') = \begin{cases} - \sqrt{d^2 \sigma_\text{n}' / \rho_\text{s}} ((\mu - C/\sigma_\text{n}') - \mu_\text{s})/(b\mu) &\text{if } \mu - C/\sigma_\text{n}' > \mu_\text{s} \text{, and}\\ - 0 &\text{if } \mu - C/\sigma_\text{n}' \leq \mu_\text{s}. + \sqrt{d^2 \sigma_\mathrm{n}' / \rho_\mathrm{s}} ((\mu - C/\sigma_\mathrm{n}') - \mu_\mathrm{s})/(b\mu) &\mathrm{if } \mu - C/\sigma_\mathrm{n}' > \mu_\mathrm{s} \mathrm{, and}\\ + 0 & \mathrm{if } \mu - C/\sigma_\mathrm{n}' \leq \mu_\mathrm{s}. \end{cases} \label{eq:g_local} \end{equation} \end{linenomath*} -where $\rho_\text{s}$ is grain mineral density, and $b$ [-] controls the non-linear rate dependence beyond yield. +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. This characteristic also strengthens the material if the shear zone size is restricted by bed geometry. \subsection{Fluid-pressure evolution}% \label{sub:fluid_pressure_evolution} -We prescribe the transient evolution of pore-fluid pressure ($p_\text{f}$) by Darcian pressure diffusion (e.g., \citeA{Goren2010, Goren2011, Damsgaard2017b}): +We prescribe the transient evolution of pore-fluid pressure ($p_\mathrm{f}$) by Darcian pressure diffusion (e.g., \citeA{Goren2010, Goren2011, Damsgaard2017b}): \begin{linenomath*} \begin{equation} - \frac{\partial p_\text{f}}{\partial t} = \frac{1}{\phi\eta_\text{f}\beta_\text{f}} \nabla \cdot (k \nabla p_\text{f}), + \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_\text{f}$ denotes dynamic fluid viscosity [Pa s], $\beta_\text{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 \citet{Henann2013}. +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}). @@ -171,14 +169,14 @@ The local fluid pressure is used to determine the effective normal stress used i \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 \citep[e.g.,][]{Iverson2000, Iverson2010-2, Damsgaard2015} or compaction and dilation (e.g., \citeA{Dewhurst1996}). +For that reason the model is not able to simulate uniaxial compaction or shear-induced volume changes (e.g., \citeA{Iverson2000, Iverson2010-2, Damsgaard2015}) or compaction and dilation (e.g., \citeA{Dewhurst1996}). We currently have a transient granular continuum model with state-dependent porosity under development. -However, \citet{Iverson2010} argued that the majority of actively deforming subglacial sediment may be in the critical state. +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 \citep{Hooke1995}, the relation to the fluidity is even more obscure. +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 (e.g., \citeA{Tulaczyk1999}). @@ -191,37 +189,37 @@ We neglect the minuscule contribution to material shear strength by water viscos The effective normal stress within the layer is found by adding the lithostatic contribution that increases with depth to the normal stress applied from the top: \begin{linenomath*} \begin{equation} - \sigma_\text{n}(z) - = \sigma_\text{n,top} + (1 - \phi)\rho_\text{s} G (L_z - z), + \sigma_\mathrm{n}(z) + = \sigma_\mathrm{n,top} + (1 - \phi)\rho_\mathrm{s} G (L_z - z), \label{eq:sigma_n} \end{equation} \end{linenomath*} where $G$ [m s$^{-2}$] is gravitational acceleration, and \begin{linenomath*} \begin{equation} - \sigma_\text{n}'(z) = \sigma_\text{n}(z) - p_\text{f}(z). + \sigma_\mathrm{n}'(z) = \sigma_\mathrm{n}(z) - p_\mathrm{f}(z). \label{eq:sigma_n_eff} \end{equation} \end{linenomath*} -Normal stress $\sigma_\text{n}(z=L_z)$ and fluid pressure $p_\text{f}(z=L_z)$ at the top are described by the boundary condition as constant or time-variable values. +Normal stress $\sigma_\mathrm{n}(z=L_z)$ and fluid pressure $p_\mathrm{f}(z=L_z)$ at the top are described by the boundary condition as constant or time-variable values. The apparent friction coefficient $\mu$ is found as: \begin{linenomath*} \begin{equation} \mu(z) - = \mu_\text{top} \frac{\sigma_\text{n,top}'}{\sigma_\text{n}'(z)}. + = \mu_\mathrm{top} \frac{\sigma_\mathrm{n,top}'}{\sigma_\mathrm{n}'(z)}. \label{eq:tau} \end{equation} \end{linenomath*} -where $\mu_\text{top}$ is constant for stress-controlled experiments and dynamic for rate-controlled experiments, with the exact numerical procedure described later. +where $\mu_\mathrm{top}$ is constant for stress-controlled experiments and dynamic for rate-controlled experiments, with the exact numerical procedure described later. -We assign depth coordinates $z_i$, granular fluidity $g_i$, and fluid pressure $p_{\text{f},i}$ to a regular grid with ghost nodes and cell spacing $\Delta z$. -The fluidity field $g$ is solved for a set of mechanical forcings ($\mu$, $\sigma_\text{n}'$, boundary conditions for $g$), and material parameters ($A$, $b$, $d$). +We assign depth coordinates $z_i$, granular fluidity $g_i$, and fluid pressure $p_{\mathrm{f},i}$ to a regular grid with ghost nodes and cell spacing $\Delta z$. +The fluidity field $g$ is solved for a set of mechanical forcings ($\mu$, $\sigma_\mathrm{n}'$, boundary conditions for $g$), and material parameters ($A$, $b$, $d$). We rearrange Eq.~\ref{eq:g} and split the Laplace operator ($\nabla^2$) into a 1D central finite difference 3-point stencil. An iterative scheme is applied to relax the following equation at each grid node $i$: \begin{linenomath*} \begin{equation} g_i = {\left(1 + \alpha_i\right)}^{-1} - \left(\alpha_i g_\text{local}(\sigma_{\text{n},i}', \mu_i) + \left(\alpha_i g_\mathrm{local}(\sigma_{\mathrm{n},i}', \mu_i) + \frac{g_{i+1} + g_{i-1}}{2} \right), \label{eq:g_i} @@ -238,12 +236,12 @@ We apply fixed-value (Dirichlet) boundary conditions for the fluidity field ($g( This condition causes the velocity field transition towards a constant value at the domain edges. Neumann boundary conditions, which are not used here, create a velocity profile resembling a free surface flow. -The pore-pressure solution (Eq.~\ref{eq:p_f}) is constrained by a hydrostatic pressure gradient at the bottom ($dp_\text{f}/dz (z=0) = \rho_\text{f}G$), and a pressure forcing at the top, for example sinusoidal: $p_\text{f}(z = L_z) = A_\text{f} \sin(2\pi f t) + p_{\text{f},0}$. -Here, $A_\text{f}$ is the forcing amplitude [Pa], $f$ is the forcing frequency [1/s], and $p_{\text{f},0}$ is the mean pore pressure over time [Pa]. +The pore-pressure solution (Eq.~\ref{eq:p_f}) is constrained by a hydrostatic pressure gradient at the bottom ($dp_\mathrm{f}/dz (z=0) = \rho_\mathrm{f}G$), and a pressure forcing at the top, for example sinusoidal: $p_\mathrm{f}(z = L_z) = A_\mathrm{f} \sin(2\pi f t) + p_{\mathrm{f},0}$. +Here, $A_\mathrm{f}$ is the forcing amplitude [Pa], $f$ is the forcing frequency [1/s], and $p_{\mathrm{f},0}$ is the mean pore pressure over time [Pa]. As for the granular flow solution, we also use operator splitting and finite differences to solve the equation for pore-pressure diffusion (Eq.~\ref{eq:p_f}): \begin{linenomath*} \begin{equation} - \Delta p_{\text{f},i} = \frac{1}{\phi_i \eta_\text{f} \beta_\text{f}} + \Delta p_{\mathrm{f},i} = \frac{1}{\phi_i \eta_\mathrm{f} \beta_\mathrm{f}} \frac{\Delta t}{\Delta z} \left( \frac{2 k_{i+1} k_i}{k_{i+1} + k_i} \frac{p_{i+1} - p_i}{\Delta z} - @@ -263,10 +261,10 @@ However, certain experiments are best approached in a rate-controlled manner whe For our system of equations, this is an inverse problem that can be tackled by adjusting the applied friction at the top until the resultant shear velocity matches the desired value. We implement an automatic iterative procedure that can be set to match a shear velocity, or limit the velocity beneath a given value. First, an initial top velocity value $v_x^*$ is calculated in a forward manner on the base of an arbitrary value for friction $\mu^*$. -The difference between $v_x^*$ and the desired velocity $v_x^\text{d}$ is used in the calculation of a normalized residual $r$: +The difference between $v_x^*$ and the desired velocity $v_x^\mathrm{d}$ is used in the calculation of a normalized residual $r$: \begin{linenomath*} \begin{equation} - r = \frac{v_x^\text{d} - v_x^*}{v_x^* + 10^{-12}}. + r = \frac{v_x^\mathrm{d} - v_x^*}{v_x^* + 10^{-12}}. \label{eq:residual} \end{equation} \end{linenomath*} @@ -275,7 +273,7 @@ If the residual value $r$ is negative, the current applied friction produces a s If the absolute value of the residual exceeds the tolerance criteria ($|r| > 10^{-3}$), the applied friction is adjusted: \begin{linenomath*} \begin{equation} - \mu^*_\text{new} = \mu^* (1.0 + \theta r), + \mu^*_\mathrm{new} = \mu^* (1.0 + \theta r), \label{eq:stressiter} \end{equation} \end{linenomath*} @@ -292,7 +290,7 @@ For the first experiment with variable water pressure, we apply a water-pressure 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 \citep[e.g.][]{Bindschadler2003}. +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 (e.g., \citeA{Bindschadler2003}). A rate-controlled setup is the opposite end member, where changes in bed friction do not influence ice flow velocity. @@ -306,10 +304,10 @@ A rate-controlled setup is the opposite end member, where changes in bed frictio 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 \citep[after][]{Iverson2010}, and the NGF model with material friction $\mu_\text{s} = 0.5$ and effective normal stress $\sigma_\text{n}' = 100$ kPa. + 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}). c) Mohr-Coulomb analysis of till samples and NGF model. - d) Modeled strain distribution under varying effective normal stress ($\sigma_\text{n}'$) with the discrete-element method \citep[DEM,][]{Damsgaard2013} and the 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*} @@ -319,26 +317,26 @@ A rate-controlled setup is the opposite end member, where changes in bed frictio \section{Results}% \label{sec:results} -We first compare the modeled mechanical behavior to various tills tested in laboratory settings (Fig.~\ref{fig:rate_dependence} and~\ref{fig:mohr_coulomb}, after \citet{Iverson2010}). +We first compare the modeled mechanical behavior to various tills tested in laboratory settings (Fig.~\ref{fig:rate_dependence} and~\ref{fig:mohr_coulomb}, after \citeA{Iverson2010}). Over five orders of strain-rate magnitude, some tills show slight rate weakening and others are slightly rate strengthening (Fig.~\ref{fig:rate_dependence}a). -Shear-strain rates up to 5.000 a$^{-1}$ are realistic for natural glacier systems \citep{Cuffey2010}. +Shear-strain rates up to 5.000 a$^{-1}$ are realistic for natural glacier systems \cite{Cuffey2010}. Our model is effectively rate-independent over most of the range, but higher $b$ values provide larger frictional resistance at extreme shear-strain rates (Fig.~\ref{fig:rate_dependence}b), making the model under these conditions rate strengthening. -The modeled friction value can be linearly scaled by adjusting $\mu_\text{s}$ in Eqs.~\ref{eq:g_local} and~\ref{eq:cooperativity}. +The modeled friction value can be linearly scaled by adjusting $\mu_\mathrm{s}$ in Eqs.~\ref{eq:g_local} and~\ref{eq:cooperativity}. Our model can simulate any combination of effective friction (or friction angle $\varphi = \tan^{-1}(\mu_s)$) and cohesion (Fig.~\ref{fig:mohr_coulomb}), which is useful as these parameters are often constrained from Mohr-Coulomb tests on till samples. The NGF model contains parameter $A$ for adjusting the degree of material non-locality (Eq.~\ref{eq:cooperativity}). Unfortunately, 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 \citet{Damsgaard2013} which allow us to calibrate $A$. +Instead, we compare the modeled strain distribution with discrete-element results from \citeA{Damsgaard2013} which allow us to calibrate $A$. By inserting relevant material parameters for grain size, friction, stress, and shear velocity (DEM, Table~S1), the NGF model model approximates the strain distribution well (Fig.~\ref{fig:strain_distribution}). 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 adjustment towards the critical state. Figure~\ref{fig:parameter_test} is a systematic analysis of parameter influence under a constant shear rate. -All experiments are at a shear rate of 300 m a$^{-1}$ and a normal stress of $\sigma_\text{n}'$ = 100 kPa. +All experiments are at a shear rate of 300 m a$^{-1}$ and a normal stress of $\sigma_\mathrm{n}'$ = 100 kPa. Several observations emerge from this parameter sensitivity analysis. The representative grain size $d$ has a major influence on the strain distribution, where finer materials show deeper deformation. The material is slightly weaker with larger grain sizes. -The shear zone is more narrow with higher material static friction coefficients ($\mu_\text{s}$), as the material is less willing to fail. +The shear zone is more narrow with higher material static friction coefficients ($\mu_\mathrm{s}$), as the material is less willing to fail. Our implementation of cohesion does not influence strain after yield. Static friction and cohesion both linearly scale the bulk friction, as expected with Mohr-Coulomb materials (see also Fig.~\ref{fig:mohr_coulomb}). The non-local amplitude $A$ slightly changes the curvature of the shear strain profile, but does not affect the overall friction. @@ -403,7 +401,7 @@ The depth of maximum shear-strain rate corresponds to the depth of minimum in ef Figure~\ref{fig:stick_slip_depth} shows a time-stacked series of simulation state with depth. The experimental setup is rate-controlled and identical to Fig.~\ref{fig:stick_slip}b and~\ref{fig:hysteresis}b. -The water pressure perturbations decay exponentially with depth with a phase shift.% \citep[p.\ 271 in][]{Turcotte2002}. +The water pressure perturbations decay exponentially with depth with a phase shift. Deep deformation occurs when the effective normal stress is smaller at depth than at the top. We next perturb the top water pressure with pulses of triangular and square shape (Fig.~\ref{fig:pulse}). @@ -416,18 +414,18 @@ In this study it is assumed that there is a strong coupling between ice and bed. However, overpressurization and slip at the ice-bed interface may cause episodic decoupling at the interface and reduce bed deformation, as observed under Whillans Ice Stream, West Antarctica (e.g., \citeA{Engelhardt1998}), and in deposits from Pleistocene glaciations (e.g., \citeA{Piotrowski2001}). We see the presented framework as a significant improvement of treating sediment advection in ice-flow models, but acknowledge that a complete understanding of the sediment mass budget requires improved models of ice-bed interface physics. -The stress-dependent sediment advection without variations in the pore pressure observed in Fig.~\ref{fig:strain_distribution} is relevant for instability theories of subglacial landform development \citep{Hindmarsh1999, Fowler2000, Schoof2007, Fowler2018}. +The stress-dependent sediment advection without variations in the pore pressure observed in Fig.~\ref{fig:strain_distribution} is relevant for instability theories of subglacial landform development \cite{Hindmarsh1999, Fowler2000, Schoof2007, Fowler2018}. From geometrical considerations, it is likely that bed-normal stresses on the stoss side of subglacial topography are higher than on the lee side. With all other physical conditions being equal, our results indicate that shear-driven sediment advection would be larger on the stoss side of bed perturbations than behind them. Topography of non-planar ice-bed interfaces (proto-drumlins, ribbed moraines, etc.) may be transported and modulated through this variable transport capacity, unless stress differences are overprinted by spatial variations in water pressure (e.g., \citeA{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_\text{s}$ [m] in our system follows the form \citep[after Eq.~4.90 in][]{Turcotte2002}, +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 (after Eq.~4.90 in \citeA{Turcotte2002}), \begin{linenomath*} \begin{equation} - d_\text{s} + d_\mathrm{s} = \sqrt{ \frac{D P}{\pi} } - = \sqrt{ \frac{k}{\phi\eta_\text{f}\beta_\text{f}\pi f} }, + = \sqrt{ \frac{k}{\phi\eta_\mathrm{f}\beta_\mathrm{f}\pi f} }, \label{eq:skin_depth} \end{equation} \end{linenomath*} @@ -449,7 +447,7 @@ Lateral water input at depth may be a viable alternate mechanism for creating oc \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 \citep{Henann2013}, but is extended with cohesion and pore-pressure diffusion. +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. In agreement with laboratory results, the material is effectively rate-independent at glacial shear velocities. A simple shear experimental setup is adapted for analyzing the mechanical response under different stresses and water-pressure variations. @@ -458,7 +456,7 @@ Deep deformation may be common in coarse-grained subglacial tills with strong an Similarly, sudden water-pressure pulses are powerful drivers for single events of deep deformation. -\acknowledgements +\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.