commit 976b9545d218a275fbe2d38a9c749f1adbf24acc
parent 4584e4efb0732a405ee22afeccf02c412273ae8a
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Thu, 10 Oct 2019 14:53:37 +0200
Change to review format, add some references in discussion
Diffstat:
2 files changed, 73 insertions(+), 6 deletions(-)
diff --git a/BIBnew.bib b/BIBnew.bib
@@ -9160,3 +9160,28 @@ Winton and A. T. Wittenberg and F. Zeng and R. Zhang and J. P. Dunne},
title = {Quantitative Rheological Model for Granular Materials: The Importance of Particle Size},
booktitle = {Handbook of Materials Modeling}
}
+@article{Fowler2018,
+ doi = {10.1080/11035897.2018.1444671},
+ year = 2018,
+ month = {apr},
+ publisher = {Informa {UK} Limited},
+ volume = {140},
+ number = {2},
+ pages = {93--105},
+ author = {A. C. Fowler},
+ title = {The philosopher in the kitchen: the role of mathematical modelling in explaining drumlin formation},
+ journal = {{GFF}}
+}
+
+@article{Hermanowski2019,
+ doi = {10.1029/2018jf004939},
+ year = 2019,
+ month = {jul},
+ publisher = {American Geophysical Union ({AGU})},
+ volume = {124},
+ number = {7},
+ pages = {1720--1741},
+ author = {P. Hermanowski and J. A. Piotrowski},
+ title = {Groundwater Flow Under a Paleo-Ice Stream of the Scandinavian Ice Sheet and Its Implications for the Formation of Stargard Drumlin Field, {NW} Poland},
+ journal = {J. Geophys. Res.: Earth Surf.}
+}
diff --git a/continuum-granular-manuscript1.tex b/continuum-granular-manuscript1.tex
@@ -14,6 +14,16 @@
\usepackage[lf]{FiraSans} % Sans-serif font
\usepackage{csquotes}
+%% Layout
+\usepackage[a4paper, margin=1.2in]{geometry}
+
+%% Review formatting
+\usepackage{lineno}
+\linenumbers % add empty line before equation environments
+\usepackage{setspace}
+%\onehalfspacing
+\doublespacing
+
%% Graphics
\usepackage{graphicx}
@@ -50,13 +60,15 @@ We show that past pulses in water pressure can transfer shear away from the ice-
\section{Introduction}%
\label{sec:introduction}
Subglacial sediment deformation is in many settings of primary importance to glacier flow \citep[e.g.][]{Boulton1974, Engelhardt1990, Fischer1994, Truffer2006}.
-Sediment mechanics influence glacier stability, sediment transport, and bedform genesis, which is why till rheology is long debated \citep[e.g.][]{Alley1986, Boulton1987, Kamb1991, Iverson1995, Hindmarsh1997, Hooke1997, Fowler2003, Kavanaugh2006, Iverson2010, Hart2011}.
+Sediment mechanics influence glacier stability, sediment transport, and bedform genesis, which is why till rheology is long debated \citep[e.g.][]{Alley1986, Boulton1987, Kamb1991, Iverson1995, Hindmarsh1997, Hooke1997, Fowler2003, Kavanaugh2006, Iverson2010, Hart2011, Fowler2018}.
Modeling of till transport requires that the strain distribution in the soft bed can be described by the stress field and material properties.
-The simplest invoked equation is of the form
+The simplest invoked equation is of the form:
+\begin{linenomath*}
\begin{equation}
\dot{\gamma} = a \frac{\tau^n}{(\sigma_\text{n}')^m},
\label{eq:rheology}
\end{equation}
+\end{linenomath*}
where $\dot{\gamma}$ is the shear-strain rate [s$^{-1}$], and $a$ is a material rate dependence [s$^{-1}$].
The shear- and effective normal stress [Pa] is denoted by $\tau$ and $\sigma_\text{n}'$, respectively.
The stress exponent ($n$ and $m$) values characterize the mechanical non-linearity.
@@ -101,10 +113,12 @@ Importantly, the Mohr-Coulomb relationship only describes the yield point of ine
\citet{Bagnold1954} realized that granular flows show a complex rate dependence.
It was later shown that a dimensionless inertia number $I$ summarizes the mechanical behavior of cohesionless, critical-state, granular flows \citep{GDR-MiDi2004}:
+\begin{linenomath*}
\begin{equation}
I = \frac{\dot{\gamma} d}{\sqrt{\sigma_\text{n}/\rho_\text{s}}},
\label{eq:inertia_number}
\end{equation}
+\end{linenomath*}
where $d$ is the representative grain size, $\sigma_\text{n}$ is the normal stress, and $\rho_\text{s}$ is the density of the grain material.
$I$ describes the relative importance of the microscopic and macroscopic time scales.
For example, if a granular material is sheared quickly, $I$ goes up as inertia increases and grains spend less time in a locked arrangement.
@@ -140,14 +154,17 @@ In this paper we expand the steady-state NGF continuum model for granular flow b
\label{sub:granular_flow}
In the NGF model, the sediment deforms as a highly nonlinear Bingham material with yield beyond the Mohr-Coulomb failure limit \citep[e.g.][]{Henann2013, Henann2016}.
We assume that elasticity is negligible and set the total shear rate $\dot{\gamma}$ to consist of a plastic contribution $\dot{\gamma}^\text{p}$:
+\begin{linenomath*}
\begin{equation}
\dot{\gamma} \approx \dot{\gamma}^\text{p} = g(\mu_\text{c}, \sigma_\text{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 is a kinematic variable governed by grain velocity fluctuations and packing fraction \citep{Zhang2017}, and consists of local and non-local components.
The local contribution to fluidity is defined as:
+\begin{linenomath*}
\begin{equation}
g_\text{local}(\mu, \sigma_\text{n}') =
\begin{cases}
@@ -156,29 +173,36 @@ The local contribution to fluidity is defined as:
\end{cases}
\label{eq:g_local}
\end{equation}
+\end{linenomath*}
where $d$ [m] is the representative grain diameter, $\mu_\text{s}$ [-] is the static Coulomb yield coefficient, $C$ [Pa] is the material cohesion, and $b$ [-] is the non-linear rate dependence beyond yield.
The failure point is determined by the Mohr-Coulomb constituent relation.
Beyond failure, the flow is governed by a Poisson-type equation that distributes strain in space according to material properties and stress state.
The degree of non-locality is scaled by the cooperativity length $\xi$:
+\begin{linenomath*}
\begin{equation}
\nabla^2 g = \frac{1}{\xi^2(\mu)} (g - g_\text{local}),
\label{eq:g}
\end{equation}
+\end{linenomath*}
where
+\begin{linenomath*}
\begin{equation}
\xi(\mu) = \frac{Ad}{\sqrt{|(\mu - C/\sigma_\text{n}') - \mu_\text{s}|}}.
\label{eq:cooperativity}
\end{equation}
+\end{linenomath*}
The non-locality scales with nonlocal amplitude $A$ [-] and grain size $d$.
In the above framework, the material strengthens when the shear zone size is restricted by thickness of the granular bed.
\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 \citep[e.g.][]{Goren2010, Goren2011, Damsgaard2017b}:
+\begin{linenomath*}
\begin{equation}
\frac{\partial p_\text{f}}{\partial t} = \frac{1}{\phi\mu_\text{f}\beta_\text{f}} \nabla \cdot (k \nabla p_\text{f}),
\label{eq:p_f}
\end{equation}
+\end{linenomath*}
where $\mu_\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}.
The 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}).
@@ -202,29 +226,36 @@ The spatial domain is $L_z = 8$ m long and is discretized into cells with equal
The upper boundary, i.e.\ the ``ice-bed interface'', exerts effective normal stress and shear stress on the granular assemblage.
We neglect the miniscule contribution to material shear strength by water viscosity.
The effective normal stress 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),
\label{eq:sigma_n}
\end{equation}
+\end{linenomath*}
+\begin{linenomath*}
\begin{equation}
\sigma_\text{n,top}' = \sigma_\text{n,top} - p_\text{f,top}
\quad \text{and} \quad
\sigma_\text{n}'(z) = \sigma_\text{n}(z) - p_\text{f}(z).
\label{eq:sigma_n_eff}
\end{equation}
+\end{linenomath*}
Normal stress ($\sigma_\text{n,top}$) and fluid pressure ($p_\text{f,top}$) at the top are given as constant or time-variable values.
The shear friction is through the depth of the model found as:
+\begin{linenomath*}
\begin{equation}
\mu(z)
= \mu_\text{top} \frac{\sigma_\text{n,top}'}{\sigma_\text{n}'(z)}.
\label{eq:tau}
\end{equation}
+\end{linenomath*}
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 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)
@@ -232,11 +263,14 @@ An iterative scheme is applied to relax the following equation at each grid node
\right),
\label{eq:g_i}
\end{equation}
+\end{linenomath*}
where
+\begin{linenomath*}
\begin{equation}
\alpha_i = \frac{\Delta z^2}{2\xi^2(\mu_i)}.
\label{eq:alpha}
\end{equation}
+\end{linenomath*}
We apply fixed-value (Dirichlet) boundary conditions for the fluidity field ($g(z=0) = g(z=L_z) = 0$).
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.
@@ -244,6 +278,7 @@ Neumann boundary conditions, which are not used here, create a velocity profile
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 sinusoidal pressure forcing at the top ($p_\text{f}(z = L_z) = A \sin(2\pi f t) + p_{\text{f},0}$).
Here, $A$ 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].
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 \mu_\text{f} \beta_\text{f}}
\frac{\Delta t}{\Delta z}
@@ -253,6 +288,7 @@ As for the granular flow solution, we also use operator splitting and finite dif
\right).
\label{eq:p_f_solution}
\end{equation}
+\end{linenomath*}
For each time step $\Delta t$, a solution to Eq.~\ref{eq:p_f_solution} is first found by explicit temporal integration.
We then use Jacobian iterations to find an implicit solution to the same equation using underrelaxation.
For the final pressure field at $t + \Delta t$ we mix the explicit and implicit solutions with equal weight, which is known as the Crank-Nicholson method \citep[e.g.][]{Patankar1980, Ferziger2002, Press2007}.
@@ -265,17 +301,21 @@ For our system of equations this is an inverse problem that can be tackled by ad
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 calculated as a normalized residual $r$:
+\begin{linenomath*}
\begin{equation}
r = \frac{v_x^\text{d} - v_x^*}{v_x^* + 10^{-12}}.
\label{eq:residual}
\end{equation}
+\end{linenomath*}
The value $10^{-12}$ is added to the denominator to avoid division by zero if the initial applied friction value $\mu^*$ does not cause yield.
If the residual value $r$ is negative, the current applied friction produces a shear velocity that exceeds the desired value, and vice versa.
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),
\label{eq:stressiter}
\end{equation}
+\end{linenomath*}
where $\theta = 10^{-2}$ is a chosen relaxation factor.
The computations are then rerun with the new applied stress until the tolerance criteria is met.
In rate-\emph{limited} experiments, the iterative procedure is only performed for negative residual values.
@@ -467,29 +507,31 @@ We next perturb the top water pressure with pulses of triangular and square shap
\section{Discussion}%
\label{sec:discussion}
-The stress-dependt sediment advection observed here is relevant for instability theories of subglacial landform development \citep{Hindmarsh1999, Fowler2000, Schoof2007}.
+The stress-dependt sediment advection observed here is relevant for instability theories of subglacial landform development \citep{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, Figure~\ref{fig:strain_distribution} indicates that sediment advection through shear is stress dependent.
-Topography of non-planar ice-bed interfaces (proto-drumlins, ribbed moraines, etc.) may be modulated or amplified through the variable transport capacity, unless normal stress variations are overprinted by spatial variations in water pressure \citep[e.g.][]{McCracken2016, Iverson2017b}.
+Topography of non-planar ice-bed interfaces (proto-drumlins, ribbed moraines, etc.) may be modulated or amplified through the variable transport capacity, unless normal stress variations are overprinted by spatial variations in water pressure \citep[e.g.][]{McCracken2016, Iverson2017b, Hermanowski2019}.
Previously, \citet{Iverson2001} modeled the subglacial bed as a series of parallel Coulomb-frictional slabs.
They demonstrated that random perturbations in effective stress at depth can distribute deformation away from the ice-bed interface.
\citet{Tulaczyk1999} and \citet{Tulaczyk2000} demonstrated that Darcian diffusion of pore-pressure variations into the bed can distribute strain away form the ice-bed interface, without a lengthscale controlling deformation.
We couple the water-pressure diffusion with a more complex sediment rheology.
-The slight rate dependence (Fig.~\ref{fig:rate-dependence}) makes it relatively trivial to couple to ice-flow models, while retaining realistic sediment physics.
+The slight rate dependence (Fig.~\ref{fig:rate_dependence}) makes it relatively trivial to couple to ice-flow models, while retaining realistic sediment physics.
The water pressure variations vary with the same periodocity as the forcing, but with exponential decay in amplitude and increasing lag at depth.
The skin depth is defined as the distance where the fluctuation amplitude decreases to $1/e$ of its surface value \citep[e.g.][]{Cuffey2010}.
As long as fluid and diffusion properties are constant, an analytical solution to skin depth $d_\text{s}$ [m] in our system follows the form \citep[after Eq.~4.90 in][]{Turcotte2002},
+\begin{linenomath*}
\begin{equation}
d_\text{s}
= \sqrt{ \frac{D P}{\pi} }
= \sqrt{ \frac{k}{\phi\mu_\text{f}\beta_\text{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.
The relation implies that the amplitude in water-pressure forcing does not influence the maximum depth of slip.
@@ -520,7 +562,7 @@ Practically all of the shear strain through a perturbation cycle occurs above th
\section*{Appendix}%
\label{sec:appendix}
The source code for the grain-water model is available at \url{https://src.adamsgaard.dk/1d_fd_simple_shear}.
-All results and figures can be reproduced by following the instructions in the experiment repository for this publication, available at \url{https://src.adamsgaard.dk/continuum_granular_exp_manus1}.
+All results and figures 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
\printbibliography{}
