commit 29380028bf88e3afe1d8a6935a9a3006bc3281e1
parent cb25136d1a2087b21154c01d5c7576d2432ce30e
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Tue, 12 Nov 2019 16:49:57 +0100
Incorporate more of Liran's feedback
Diffstat:
2 files changed, 63 insertions(+), 25 deletions(-)
diff --git a/BIBnew.bib b/BIBnew.bib
@@ -9184,3 +9184,40 @@ Winton and A. T. Wittenberg and F. Zeng and R. Zhang and J. P. Dunne},
title={Impact of the ice strength formulation on the performance of sa sea ice thickness distribution model in the Arctic},
journal = {J. Geophys. Res.: Oceans}
}
+
+@article{Kavanaugh2009,
+ doi = {10.1029/2008jf001036},
+ url = {https://doi.org/10.1029%2F2008jf001036},
+ year = 2009,
+ month = {feb},
+ publisher = {American Geophysical Union ({AGU})},
+ volume = {114},
+ number = {F1},
+ author = {J. L. Kavanaugh},
+ title = {Exploring glacier dynamics with subglacial water pressure pulses: Evidence for self-organized criticality?},
+ journal = {J. Geophys. Res.}
+}
+@article{Christoffersen2018,
+ doi = {10.1038/s41467-018-03420-8},
+ url = {https://doi.org/10.1038%2Fs41467-018-03420-8},
+ year = 2018,
+ month = {mar},
+ publisher = {Springer Science and Business Media {LLC}},
+ volume = {9},
+ number = {1},
+ author = {P. Christoffersen and M. Bougamont and A. Hubbard and S. H. Doyle and S. Grigsby and R. Pettersson},
+ title = {Cascading lake drainage on the Greenland Ice Sheet triggered by tensile shock and fracture},
+ journal = {Nature Commun.}
+}
+@article{Palmer2015,
+ doi = {10.1038/ncomms9408},
+ url = {https://doi.org/10.1038%2Fncomms9408},
+ year = 2015,
+ month = {oct},
+ publisher = {Springer Science and Business Media {LLC}},
+ volume = {6},
+ number = {1},
+ author = {S. Palmer and M. McMillan and M. Morlighem},
+ title = {Subglacial lake drainage detected beneath the Greenland ice sheet},
+ journal = {Nature Commun.}
+}+
\ No newline at end of file
diff --git a/continuum-granular-manuscript1.tex b/continuum-granular-manuscript1.tex
@@ -162,39 +162,39 @@ The local contribution to fluidity is defined as:
\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_\text{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 the above framework, the material strengthens when the shear zone size is restricted by thickness of the granular bed.
+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 \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}),
+ \frac{\partial p_\text{f}}{\partial t} = \frac{1}{\phi\eta_\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$].
+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}.
-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}).
+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 \citep[e.g.,][]{Iverson2000, Iverson2010-2, Damsgaard2015} or compaction \citep[e.g.,][]{Dewhurst1996}.
+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 \citep[e.g.,][]{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.
-For that reason, we see this contribution as a valuable first step.
+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.
-However, it is awkward to describe grain size distributions of diamictons with fractal grain size distribution with a single length scale \citep{Hooke1995}.
-We expect that a volumetrically dominant grain size dominates the strain distribution, outside of effects of ploughing by large clasts \citep[e.g.,][]{Tulaczyk1999}.
-Future research will investigate how wide grain-size distributions affect strain distribution, and will benchmark against specifically designed laboratory experiments on tills.
+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.
+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 \citep[e.g.,][]{Tulaczyk1999}.
\subsection{Numerical solution procedure}%
\label{sub:numerical_solution_procedure}
@@ -202,7 +202,7 @@ We apply the model in a one-dimensional setup where simple shear occurs along a
The spatial domain is $L_z = 8$ m long and is discretized into cells with equal size to the representative grain size $d$.
The upper boundary, i.e.\ the ``ice-bed interface'', exerts effective normal stress and shear stress on the granular assemblage.
We neglect the minuscule 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:
+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)
@@ -210,16 +210,15 @@ The effective normal stress is found by adding the lithostatic contribution that
\label{eq:sigma_n}
\end{equation}
\end{linenomath*}
+where $G$ [m s$^{-2}$] is gravitational acceleration, and
\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:
+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.
+The apparent friction coefficient $\mu$ is found as:
\begin{linenomath*}
\begin{equation}
\mu(z)
@@ -227,6 +226,7 @@ The shear friction is through the depth of the model found as:
\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.
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$).
@@ -252,12 +252,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 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].
+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].
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}}
+ \Delta p_{\text{f},i} = \frac{1}{\phi_i \eta_\text{f} \beta_\text{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} -
@@ -266,10 +266,10 @@ As for the granular flow solution, we also use operator splitting and finite dif
\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}.
+For each time step $\Delta t$, a solution to Eq.~\ref{eq:p_f_solution} is found by the Crank-Nicholson (CN) method \citep[e.g.,][]{Patankar1980, Ferziger2002, Press2007}.
+In the procedure the pressure field at $t + \Delta t$ is found by mixing explicit and implicit solutions with equal weight.
The method is unconditionally stable and second-order accurate in time and space.
+Our implementation of grain and fluid dynamics is highly efficient, and for the presented experiments each time step completes in less than 1 ms on a single CPU core.
\subsubsection{Rate-controlled experiments} % quick edit, needs rewrite. perhaps also move somewhere else
The continuum model is in its presented form suited for resolving strain rate and shear velocity from a given stress forcing, i.e. in a stress-controlled setup.
@@ -505,7 +505,7 @@ As long as fluid and diffusion properties are constant, an analytical solution t
\begin{equation}
d_\text{s}
= \sqrt{ \frac{D P}{\pi} }
- = \sqrt{ \frac{k}{\phi\mu_\text{f}\beta_\text{f}\pi f} },
+ = \sqrt{ \frac{k}{\phi\eta_\text{f}\beta_\text{f}\pi f} },
\label{eq:skin_depth}
\end{equation}
\end{linenomath*}
