seaice-experiments

main.tex (20080B)

---

 %\documentclass[11pt,a4paper]{article}
 \documentclass[11pt]{article}
 
 %\usepackage{a4wide}
 \usepackage[margin=1.4in]{geometry}
 
 \usepackage{graphicx}
 %\usepackage[german, english]{babel}
 %\usepackage{tabularx}
 %\usepackage{cancel}
 %\usepackage{multirow}
 %\usepackage{supertabular}
 %\usepackage{algorithmic}
 %\usepackage{algorithm}
 %\usepackage{amsthm}
 %\usepackage{float}
 %\usepackage{subfig}
 \usepackage{rotating}
 \usepackage{amsmath}
 
 \usepackage[T1]{fontenc} % Font encoding
 \usepackage{charter} % Serif body font
 \usepackage[charter]{mathdesign} % Math font
 \usepackage[scale=0.9]{sourcecodepro} % Monospaced fontenc
 \usepackage[lf]{FiraSans} % Sans-serif font
 \usepackage{textgreek}
 
 \usepackage{listings}
 
 \usepackage{hyperref}
 
 \usepackage{soul} % for st strikethrough command
 
 %\usepackage[round]{natbib}
 \usepackage[natbib=true, style=authoryear, bibstyle=authoryear-comp, 
 maxbibnames=10,
 maxcitenames=2, backend=bibtex8]{biblatex}
 \bibliography{/Users/ad/articles/own/BIBnew.bib}
 
 
 \begin{document}
 
 \title{Discrete-element simulation of sea-ice mechanics:\\
     Parameterization of contact mechanics}
 
 \author{\small Anders Damsgaard}
 \date{\small \today}
 
 \maketitle
 
 \section{Rationale}
 Lagrangian parameterizations of sea ice offer several advantages to Eulerian 
 continuum methods.  The continuum assumption of sea-ice mechanics becomes 
 increasingly questionable as spatial resolution and cell size approaches the 
 size of ice floes.  Ice-floe scale discretizations are natural for Lagrangian 
 models, which additionally offer the convenience of being able to handle 
 arbitrary sea-ice concentrations.  This is likely to improve model performance 
 in ice-marginal zones with strong advection.  Furthermore, phase transitions in 
 granular rheology around the jamming limit, such as observed when sea ice moves 
 through geometric confinements, includes sharp thresholds in effective 
 viscosity, which are typically ignored in Eulerian models.
 Granular jamming is a stochastic process dependent on having the right grains in 
 the right place at the right time, and the jamming likelihood over time can be 
 described by a probabilistic model
 \citep[e.g.][]{Tang2009}.  Difficult to parameterize in continuum formulations 
 \citep[e.g.][]{Rallabandi2017}, jamming comes natural to dense granular systems 
 simulated in a Lagrangian framework, and is a very relevant process complicating 
 sea-ice transport through narrow straits \citep[e.g.][]{Kwok2010}.
 However, computational performance is a central concern for coupled 
 ocean-atmosphere models, with traditional discrete-element method (DEM) 
 formulations being too costly for global-scale models.  Here we investigate the 
 behavior of a Lagrangian sea-ice model and assess how mechanical 
 parameterizations of differing complexity affect collective behavior at scales 
 much larger than the individual ice floes.
 
 
 \section{Contact models}
 The presented experiments compare jamming behavior between two differing 
 ice-floe contact models.  Common to both models, the resistive force 
 $\boldsymbol{f}_\text{n}$ to axial compressive strain between to cylindrical ice 
 floes $i$ and $j$ is provided by (Hookean) linear-elasticity based on the 
 overlap distance $\boldsymbol{\delta}_\text{n}$:
 \begin{equation}
     \boldsymbol{f}_\text{n}^{ij} =
     A^{ij} E^{ij} \boldsymbol{\delta}_\text{n}^{ij}
     \quad \text{when} \quad
     0 > |\boldsymbol{\delta}_\text{n}^{ij}| \equiv |\boldsymbol{x}^i - 
     \boldsymbol{x}^j| - (r^i + r^j)
     \label{eq:f_n}
 \end{equation}
 where the contact cross-sectional area $A^{ij} = R^{ij} \min(T^i, T^j)$ is 
 determined by the harmonic mean $R^{ij} = 2 r^i r^j/(r^i + r^j)$ of the ice-floe 
 radii $r^i$ and $r^j$, as well as the smallest of the involved ice-floe 
 thickness $T^i$ and $T^j$.  The linear elastic force resulting from axial strain 
 of a distance $|\boldsymbol{\delta}_\text{n}^{ij}|$ of the contact is scaled by 
 the harmonic mean of Young's modulus $E^{ij}$.  The stiffness scaling based on 
 Young's modulus is scale invariant \citep[e.g.][]{Obermayr2013}, and is based on 
 the principle that the elastic stiffness of the ice itself is constant, 
 regardless of ice-floe size.
 
 Nonlinear elasticity models based on Hertzian contact mechanics may 
 alternatively be applied to determine the stresses resulting from contact 
 compression \citep[e.g.][]{Herman2013}.  However, with nonlinear models the 
 numerical stability of the explicit temporal integration scheme will depend on 
 the stress and packing state of the granular assemblage, and will under 
 localized compressive extremes require very small discretization in time.
 
 As demonstrated later, models based compressive strength alone are not 
 sufficient to cause granular jamming.  We explore two different additions to the 
 contact model presented in Eq.~\ref{eq:f_n}.  The first approach is typical to 
 DEM models and is based on resolving including shear resistance through 
 tangential (contact parallel) elasticity and associated Coulomb frictional 
 limits.  An alternative approach, fundamentally complementary to compressive 
 elasticity and shear friction, is tensile strength of ice-floe contacts which 
 leads to a cohesive bulk granular rheology.
 
 \subsection{Contact-parallel elasticity with Coulomb friction}
 The contact-parallel (tangential) elasticity is resolved by determining the 
 contact travel distance $\boldsymbol{\delta}_\text{t}$ on the contact plane for 
 the duration of the contact $t_\text{c}$:
 \begin{equation}
     \boldsymbol{\delta}_\text{t} = \int_0^{t_\text{c}}
     \left[
         (\boldsymbol{v}^i - \boldsymbol{v}^j) \cdot \hat{\boldsymbol{t}}^{ij} -
         R^{ij} \left(\omega^i + \omega^j\right)
         \right]
     \label{eq:d_t}
 \end{equation}
 where $\boldsymbol{v}$ and $\omega$ denotes linear and rotational velocity, 
 respectively.
 The deformation distance $\boldsymbol{\delta}_\text{t}$ is corrected for contact 
 rotation over the lifetime of the interaction, and is used to determine the 
 contact-tangential elastic force:
 \begin{equation}
     \boldsymbol{f}_\text{t} = \frac{E^{ij} A^{ij}}{R^{ij}}
     \frac{2 (1 - \nu^2)}{(2 - \nu)(1 + \nu)} \boldsymbol{\delta}_\text{t}
 \end{equation}
 with $\nu$ being the dimensionless Poisson's ratio characteristic for the 
 material.  However, the frictional force above is restricted to the 
 Coulomb-frictional limit, relating the magnitude of the normal force to the 
 magnitude of the tangential force:
 \begin{equation}
     |\boldsymbol{f}_\text{t}| \leq \mu |\boldsymbol{f}_\text{n}|
     \label{eq:coulomb_friction}
 \end{equation}
 The Coulomb-frictional coefficient $\mu$ may be parameterized as state 
 dependent, for example with different values for sliding and stationary 
 contacts.  However, for the following a single and constant value is used for 
 the frictional coefficient $\mu$.
 
 In the case of slip ($|\boldsymbol{f}_\text{t}| > \mu 
 |\boldsymbol{f}_\text{n}|$) the length of the contact-parallel travel path 
 $\boldsymbol{\delta}_\text{t}$ is reduced to be consistent to be consistent with 
 the Coulomb limit.  This accounts for the tangential contact plasticity and 
 irreversible work associated with sliding.
 
 Since the above model of tangential shear resistance is based on deformation 
 distance on the inter-floe contact plane, it requires solving for ice-floe 
 rotational kinematics of each ice floe and a bookkeeping algorithm for storing 
 contact histories.
 
 \subsection{Tensile contact strength}
 Cohesion is introduced by parameterizing resistance to extension beyond the 
 overlap distance between a pair of ice floes (i.e.\ $\delta_\text{n}^{ij} > 0$).  
 In real settings tensile strength may be provided by refreezing processes at the 
 ice-floe interface or mechanical ridging.  An ideal formulation of bond 
 deformation includes resistance to bond tension, shear, twist, and rolling 
 \citep[e.g.][]{Potyondy2004, Obermayr2013, Herman2016}.  However, for this study 
 we explore the possibility of using bond \emph{tension} alone as a mechanical 
 component contributing to bulk granular shear strength.
 
 Tensile strength is parameterized by applying Eq.~\ref{eq:f_n} for the extensive 
 regime ($\delta_\text{n} > 0$).  Eq.~\ref{eq:f_n} is applied until the tensile 
 stress exceeds the tensile strength $\sigma_\text{c}$ defined for the bonds.  
 The implementation allows for time-dependent tensile strengthening based on 
 contact age, but for the following we parameterize the bonds to obtain full 
 tensile strength as soon as a pair of ice floes first undergo compression 
 ($\delta_\text{n} < 0$).
 
 
 \subsection{Implementation}
 The intent of this work is to develop and apply an offline sea-ice dynamical 
 code in order to explore strengths and limitations of different methods related 
 to sea-ice mechanics.  Once the appropriate parameterization has been identified 
 it will be built in to the GFDL model of coupled ocean/atmosphere dynamics.
 
 The presented experiments are performed with the 
 \texttt{SeaIce.jl}\footnote{\url{https://github.com/anders-dc/SeaIce.jl}} Julia 
 package for offline simulations of sea-ice mechanics, which includes the above 
 parameterizations.  Furthermore, the model includes the option to include linear 
 (Newtonian) contact viscosity in the normal and tangential components, but an 
 analysis of these is not included here.  Simulation-specific scripts are 
 provided in the repository 
 \texttt{SeaIce-experiments}\footnote{\url{https://github.com/anders-dc/SeaIce-experiments}}.
 
 The offline model implementation identifies ice-floe contacts and placement in 
 the ocean/atmosphere grid through a cell-based spatial discretization, which 
 reduces the computational overhead significantly (Fig.~\ref{fig:performance}, 
 and additionally mirrors the current Lagrangian ice-berg implementation in the 
 GFDL ocean model.
 
 \begin{figure}
     \begin{center}
         \includegraphics[width=0.7\textwidth]{graphics/profiling-cpu.pdf}
     \end{center}
     \caption{\label{fig:performance} Single-threaded performance of the DEM 
     algorithm with all-to-all contact search, or a contact search using spatial 
     decomposition based on the ocean/atmosphere grid (Profiling script: 
     \url{https://github.com/anders-dc/SeaIce.jl/blob/master/test/profiling.jl}).}
 \end{figure}
 
 The experiments rely on pseudo-random number generation for generating ice-floe 
 particle size distributions (PSDs).  In order to assess the statistical 
 probability of granular jamming we seed the system with different values and 
 repeat each experiment ten times.
 
 For the presented runs the ice floes are forced with a uniform and constant wind 
 field oriented from north to south, with a simplified channel-like constriction 
 in the middle.  The ocean also flows from north to south but is constrained by 
 the basin geometry.  The spatial velocity pattern is determined by a stream 
 function under mass conservation, meaning that currents increase as
 the channel narrows.  Ice floes are initially placed in a pseudo-random 
 arrangement north of the channel.  All parameters are kept constant between the 
 experiment sets unless explicitly noted.
 
 \section{Results}
 
 \subsection{Experiment set 1}
 
 Increasing frictional coefficients, increasing tensile strength.  Constant 
 uniform size distribution.
 Generated by running \texttt{`make`} from 
 \url{https://github.com/anders-dc/SeaIce-experiments/tree/master/cohesion-friction-comparison}
 Animation with \texttt{seed=1} can be found at:
 \url{https://youtu.be/JmWvoi9Z6Zk}
 
 \begin{sidewaysfigure}
     \begin{center}
         \includegraphics[width=0.99\textheight]%
         {graphics/cohesion-friction-comparison.png}
         \includegraphics[width=0.99\textheight]%
         {{graphics/cohesion-friction-comparison-jam_fraction.png}}
     \end{center}
     \caption{\label{fig:set1}%
     \textbf{Experiment set 1}:
     Comparison between the frictional and cohesive model with different 
     Coulomb-frictional coefficient values and tensile strengths.
     First row: Ice fluxes over time with increasing friction ($\mu = [0.0, 0.1, 
     0.2, 0.3, 0.4]$).
     Second row: Ice fluxes over time with increasing tensile strength 
     ($\sigma_\text{c} = [10, 100, 200, 400, 800]$ kPa).
     Third row: Fraction of runs jammed over time with increasing friction.
     Fourth row: Fraction of runs jammed over time with increasing tensile 
     strength.
     }
 \end{sidewaysfigure}
 
 While the frictionless runs rarely jam, we observe that increasing either 
 friction or cohesion greatly increases the likelihood of jamming 
 (Fig.~\ref{fig:set1}).  Under cohesive runs with large tensile strengths the ice 
 floes north of the assemblage behave like a single rigid block, while the 
 corresponding frictional runs show more spread in the jam-fraction probability.  
 At intermediate to low frictional coefficients ($\mu \leq 0.3$ or 
 $\sigma_\text{c} \leq$ 200 kPa) there is good agreement between frictional and 
 cohesive models, where the runs with $\mu = 0.3$ most closely resemble the 
 jam-fraction distribution of $\sigma_\text{c} = 200$ kPa.  We base further runs 
 on comparing the differences between these two parameter sets under varying 
 circumstances.
 
 \subsection{Experiment set 2}
 Constant friction coefficient ($\mu = 0.3$) for frictional runs, constant 
 tensile strength ($\sigma_\text{c} = 200$ kPa) for cohesive runs.  Variable 
 channel width.  Monodisperse ice-floe size distribution.
 Generated by running \texttt{`make`} from 
 \url{https://github.com/anders-dc/SeaIce-experiments/tree/master/cohesion-friction-width-monodisperse-comparison}.  
 Animation with \texttt{seed=1} can be found at:
 \url{https://youtu.be/tPMSxdw1UZ8}.
 
 \begin{sidewaysfigure}
     \begin{center}
         \includegraphics[width=0.99\textheight]%
         {graphics/cohesion-friction-width-monodisperse-comparison.png}
         \includegraphics[width=0.99\textheight]%
         {{graphics/cohesion-friction-width-monodisperse-comparison-jam_fraction.png}}
     \end{center}
     \caption{\label{fig:set2}%
     \textbf{Experiment set 2}:
     Comparison between the frictional ($\mu = 0.3$ and $\sigma_\text{c}$ = 0 
     kPa) and cohesive ($\mu = 0$ and $\sigma_\text{c}$ = 200 kPa) models with a 
     monodisperse (single) ice-floe size and increasing channel width.
     First row: Ice fluxes over time with friction ($\mu = 0.3$) and increasing 
     channel width.
     Second row: Ice fluxes over time with cohesion ($\sigma_\text{c} = 200$ kPa) 
     and increasing channel width.
     Third row: Fraction of frictional runs jammed over time.
     Fourth row: Fraction of cohesive runs jammed over time.
     }
 \end{sidewaysfigure}
 
 The jamming probability of two-dimensional granular assemblages under gravity is 
 well described in the literature, especially for monodisperse (same-sized) 
 granular materials \citep[e.g.][]{Beverloo1961, To2001, Tang2009}.
 While the ocean/atmosphere forcing is different from the gravity drag in the 
 silo flows in these publications, we observe the same relationship between 
 conduit width and grain size in monodisperse granular assemblages as previously 
 reported.  The consistency between the frictional and cohesive model 
 (Fig.~\ref{fig:set2}) indicates that the simpler cohesive model is performing 
 well in this setting.
 
 \subsection{Experiment set 3}
 Constant friction coefficient ($\mu = 0.3$) for frictional runs, constant 
 tensile strength ($\sigma_\text{c} = 200$ kPa) for cohesive runs.  Variable 
 channel width.  Constant and uniformly-distributed ice-floe size distribution.
 Generated by running \texttt{`make`} from 
 \url{https://github.com/anders-dc/SeaIce-experiments/tree/master/cohesion-friction-width-comparison}.
 
 \begin{sidewaysfigure}
     \begin{center}
         \includegraphics[width=0.99\textheight]%
         {graphics/cohesion-friction-width-comparison.png}
         \includegraphics[width=0.99\textheight]%
         {{graphics/cohesion-friction-width-comparison-jam_fraction.png}}
     \end{center}
     \caption{\label{fig:set3}%
     \textbf{Experiment set 3}:
     Comparison between the frictional ($\mu = 0.3$ and $\sigma_\text{c}$ = 0 
     kPa) and cohesive ($\mu = 0$ and $\sigma_\text{c}$ = 200 kPa) models with a 
     uniform ice-floe size distribution and increasing channel width.
     First row: Ice fluxes over time with friction ($\mu = 0.3$) and increasing 
     channel width.
     Second row: Ice fluxes over time with cohesion ($\sigma_\text{c} = 200$ kPa) 
     and increasing channel width.
     Third row: Fraction of frictional runs jammed over time.
     Fourth row: Fraction of cohesive runs jammed over time.
     }
 \end{sidewaysfigure}
 
 While the monodisperse ice-floe size used in Set 2 provides a comparative basis 
 to the literature on granular jamming, we have for Set 3 used a wide grain size 
 distribution and varied the channel width as before.  We observe that the width 
 of ice-floe size distribution decreases the likelihood of jamming, and that 
 these runs are less likely to jam than the corresponding monodisperse ensemble 
 runs.  Again we see a good correspondence between the frictional and cohesive 
 model runs (Fig.~\ref{fig:set3}).
 
 \subsection{Experiment set 4}
 Constant friction coefficient ($\mu = 0.3$) for frictional runs, constant 
 tensile strength ($\sigma_\text{c} = 200$ kPa) for cohesive runs.  Constant 
 channel width.  Different widths in power-law distributed ice-floe size 
 distribution.
 Generated by running \texttt{`make`} from 
 \url{https://github.com/anders-dc/SeaIce-experiments/tree/master/cohesion-friction-psd-comparison}.
 Animation with \texttt{seed=1} can be found at: 
 \url{https://youtu.be/I9P8XsA1S0I}
 
 \begin{sidewaysfigure}
     \begin{center}
         \includegraphics[width=0.99\textheight]%
         {graphics/cohesion-friction-psd-comparison.png}
         \includegraphics[width=0.99\textheight]%
         {{graphics/cohesion-friction-psd-comparison-jam_fraction.png}}
     \end{center}
     \caption{\label{fig:set4}%
     \textbf{Experiment set 4}:
     Comparison between the frictional ($\mu = 0.3$ and $\sigma_\text{c}$ = 0 
     kPa) and cohesive ($\mu = 0$ and $\sigma_\text{c}$ = 200 kPa) models with a 
     increasingly wide power-law ice-floe size distribution and constant channel 
     width.
     First row: Ice fluxes over time with friction ($\mu = 0.3$) and increasing 
     ice-floe size span.
     Second row: Ice fluxes over time with cohesion ($\sigma_\text{c} = 200$ kPa) 
     and increasing ice-floe size span.
     Third row: Fraction of frictional runs jammed over time.
     Fourth row: Fraction of cohesive runs jammed over time.
     }
 \end{sidewaysfigure}
 
 Ice-floe distributions in natural sea-ice settings have been observed to follow 
 a power-law probability distribution with an exponent of $\alpha \approx -1.8$ 
 \citep[e.g.][]{Steer2008, Herman2013}.  With constant channel size and an 
 increasing width of the ice-floe size distribution, we observe that the 
 probability of jamming decreases, and observe good agreement between the 
 frictional and cohesive model runs.
 
 \section{Summary}
 We construct a flexible discrete-element framework for simulating Lagrangian 
 sea-ice dynamics at the ice-floe scale, forced by ocean and atmosphere velocity 
 fields.  While frictionless contact models based on tensile stiffness alone are 
 very unlikely to jam, we describe two different approaches based on friction and 
 tensile strength, both resulting in increased bulk shear strength of the 
 granular assemblage.  We demonstrate that the discrete-element approach is able 
 to undergo dynamic jamming when forced through an idealized confinement, where 
 the probability of jamming is determined by the channel width, ice-floe size, 
 and ice-floe size variability.  The frictionless but cohesive contact model can 
 with certain tensile strength values display jamming behavior which on the large 
 scale is highly similar to a model with contact friction and ice-floe rotation.  
 The frictionless and cohesive model is likely an ideal candidate for coupled 
 computations due to its reduced algorithmic complexity.
 
 \printbibliography{}
 
 \end{document}