ceaice_split_version/ceaice_part1/ceaice_appendix.tex

\section{Dynamics\label{app:dynamics}}

% \newcommand{\vek}[1]{\ensuremath{\vec{\mathbf{#1}}}}
% \newcommand{\vtau}{\vek{\mathbf{\tau}}}
For completeness we provide more details on the ice dynamics of the
sea-ice model. The momentum equation are
\begin{equation}
  \label{eq:momseaice}
  m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} +
  \vtau_{ocean} - m \nabla{\phi(0)} + \vek{F},
\end{equation}
where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area;
$\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector;
$\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$
directions, respectively;
$f$ is the Coriolis parameter;
$\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses,
respectively;
$g$ is the gravity accelation;
$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;
$\phi(0) = g\eta + p_{a}/\rho_{0} + mg/\rho_{0}$ is the sea surface
height potential in response to ocean dynamics ($g\eta$), to
atmospheric pressure loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a
reference density) and a term due to snow and ice loading \citep{campin08};
and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice
stress tensor $\sigma_{ij}$. %
Advection of sea-ice momentum is neglected. The wind and ice-ocean stress
terms are given by
\begin{align*}
  \vtau_{air}   = & \rho_{air}  C_{air}   |\vek{U}_{air}  -\vek{u}|
                   R_{air}  (\vek{U}_{air}  -\vek{u}), \\ 
  \vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}| 
                   R_{ocean}(\vek{U}_{ocean}-\vek{u}),
\end{align*}
where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere
and surface currents of the ocean, respectively; $C_{air/ocean}$ are
air and ocean drag coefficients; $\rho_{air/ocean}$ are reference
densities; and $R_{air/ocean}$ are rotation matrices that act on the
wind/current vectors.

For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can
be related to the ice strain rate and strength by a nonlinear
viscous-plastic (VP) constitutive law \citep{hibler79, zhang97}:
\begin{equation}
  \label{eq:vpequation}
  \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij} 
  + \left[\zeta(\dot{\epsilon}_{ij},P) -
    \eta(\dot{\epsilon}_{ij},P)\right]\dot{\epsilon}_{kk}\delta_{ij}  
  - \frac{P}{2}\delta_{ij}.
\end{equation}
The ice strain rate is given by
\begin{equation*}
  \dot{\epsilon}_{ij} = \frac{1}{2}\left( 
    \frac{\partial{u_{i}}}{\partial{x_{j}}} +
    \frac{\partial{u_{j}}}{\partial{x_{i}}}\right).
\end{equation*}
The maximum ice pressure $P_{\max}$, a measure of ice strength, depends on
both thickness $h$ and compactness (concentration) $c$:
\begin{equation}
  P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]},
\label{eq:icestrength}
\end{equation}
with the constants $P^{*}$ and $C^{*}=20$. The nonlinear bulk and shear 
viscosities $\eta$ and $\zeta$ are functions of ice strain rate
invariants and ice strength such that the principal components of the
stress lie on an elliptical yield curve with the ratio of major to
minor axis $e$ equal to $2$; they are given by:
\begin{align*}
  \zeta =& \min\left(\frac{P_{\max}}{2\max(\Delta,\Delta_{\min})},
   \zeta_{\max}\right) \\
  \eta =& \frac{\zeta}{e^2} \\
  \intertext{with the abbreviation}
  \Delta = & \left[
    \left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right)
    (1+e^{-2}) +  4e^{-2}\dot{\epsilon}_{12}^2 + 
    2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2})
  \right]^{\frac{1}{2}}.
\end{align*}
In the simulations of this paper, the bulk viscosities are bounded
above by imposing both a minimum
$\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ and a maximum $\zeta_{\max} =
P_{\max}/\Delta^*$, where
$\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress
tensor computation the replacement pressure $P = 2\,\Delta\zeta$
\citep{hibler95} is used so that the stress state always lies on the
elliptic yield curve by definition.

In the so-called truncated ellipse method (experiment TEM) the shear
viscosity $\eta$ is capped to suppress any tensile stress
\citep{hibler97, geiger98}:
\begin{equation}
  \label{eq:etatem}
  \eta = \min\left(\frac{\zeta}{e^2},
  \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
  {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2
      +4\dot{\epsilon}_{12}^2}}\right).
\end{equation}

In the current implementation, the VP-model is integrated with the
semi-implicit line successive over relaxation (LSOR)-solver of
\citet{zhang97}, which allows for long time steps that, in our case,
are limited by the explicit treatment of the Coriolis term. The
explicit treatment of the Coriolis term does not represent a severe
limitation because it restricts the time step to approximately the
same length as in the ocean model where the Coriolis term is also
treated explicitly.

\citet{hunke97}'s introduced an elastic contribution to the strain
rate in order to regularize Eq.~\ref{eq:vpequation} in such a way that
the resulting elastic-viscous-plastic (EVP) and VP models are
identical at steady state,
\begin{equation}
  \label{eq:evpequation}
  \frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} +
  \frac{1}{2\eta}\sigma_{ij} 
  + \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij}  
  + \frac{P}{4\zeta}\delta_{ij}
  = \dot{\epsilon}_{ij}. 
\end{equation}
The EVP-model uses an explicit time stepping scheme with a short
timestep. According to the recommendation of \citet{hunke97}, the
EVP-model is stepped forward in time O(120) times within the physical
ocean model time step, to
allow for elastic waves to disappear.  Because the scheme does not
require a matrix inversion it is fast in spite of the small internal
timestep and simple to implement on parallel computers
\citep{hunke97}. For completeness, we repeat the equations for the
components of the stress tensor $\sigma_{1} =
\sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and
$\sigma_{12}$. Introducing the divergence $D_D =
\dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension
and shearing strain rates, $D_T =
\dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =
2\dot{\epsilon}_{12}$, respectively, and using the above
abbreviations, the equations~\ref{eq:evpequation} can be written as:
\begin{align}
  \label{eq:evpstresstensor1}
  \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +
  \frac{P}{2T} &= \frac{P}{2T\Delta} D_D \\
  \label{eq:evpstresstensor2}
  \frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T}
  &= \frac{P}{2T\Delta} D_T \\
  \label{eq:evpstresstensor12}
  \frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T}
  &= \frac{P}{4T\Delta} D_S 
\end{align}
Here, the elastic parameter $E$ is redefined in terms of a damping
timescale $T$ for elastic waves \[E=\frac{\zeta}{T}.\]
$T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and the external
(long) timestep $\Delta{t}$. We use $E_{0} = \frac{1}{3}$ which is
close to value of 0.36 used by \citet{hunke01}.

\section{Finite-volume discretization of the stress tensor
  divergence}
\label{app:discretization}
On an Arakawa C~grid, ice thickness and concentration and thus ice
strength $P$ and bulk and shear viscosities $\zeta$ and $\eta$ are
naturally defined at C-points in the center of the grid
cell. Discretization requires only averaging of $\zeta$ and $\eta$ to
vorticity or Z-points at the bottom left corner of the cell to give
$\overline{\zeta}^{Z}$ and $\overline{\eta}^{Z}$. In the following,
the superscripts indicate location at Z or C points, distance across
the cell (F), along the cell edge (G), between $u$-points (U),
$v$-points (V), and C-points (C). The control volumes of the $u$- and
$v$-equations in the grid cell at indices $(i,j)$ are $A_{i,j}^{w}$
and $A_{i,j}^{s}$, respectively. With these definitions (which follow
the model code documentation at \url{http://mitgcm.org} except that
vorticity or $\zeta$-points have been renamed to Z-points in order to
avoid confusion with the bulk viscosity $\zeta$), the strain
rates are discretized as:
\begin{linenomath*}\begin{align}
  \dot{\epsilon}_{11} &= \partial_{1}{u}_{1} + k_{2}u_{2} \\ \notag
  => (\epsilon_{11})_{i,j}^C &= \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} 
   + k_{2,i,j}^{C}\frac{v_{i,j+1}+v_{i,j}}{2} \\ 
  \dot{\epsilon}_{22} &= \partial_{2}{u}_{2} + k_{1}u_{1} \\\notag
  => (\epsilon_{22})_{i,j}^C &= \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} 
   + k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ 
   \dot{\epsilon}_{12} = \dot{\epsilon}_{21} &= \frac{1}{2}\biggl(
   \partial_{1}{u}_{2} + \partial_{2}{u}_{1} - k_{1}u_{2} - k_{2}u_{1}
   \biggr) \\ \notag
  => (\epsilon_{12})_{i,j}^Z &= \frac{1}{2}
  \biggl( \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^V} 
   + \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^U} \\\notag
  &\phantom{=\frac{1}{2}\biggl(}
   - k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2}
   - k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2}
   \biggr),
\end{align}\end{linenomath*}
so that the diagonal terms of the strain rate tensor are naturally
defined at C-points and the symmetric off-diagonal term at
Z-points. No-slip boundary conditions ($u_{i,j-1}+u_{i,j}=0$ and
$v_{i-1,j}+v_{i,j}=0$ across boundaries) are implemented via
``ghost-points''; for free slip boundary conditions
$(\epsilon_{12})^Z=0$ on boundaries.

For a spherical polar grid, the coefficients of the metric terms are
$k_{1}=0$ and $k_{2}=-\tan\phi/a$, with the spherical radius $a$ and
the latitude $\phi$; $\Delta{x}_1 = \Delta{x} = a\cos\phi
\Delta\lambda$, and $\Delta{x}_2 = \Delta{y}=a\Delta\phi$. For a
general orthogonal curvilinear grid as used in this paper, $k_{1}$ and
$k_{2}$ can be approximated by finite differences of the cell widths:
%\enlargethispage{1cm}
\begin{linenomath*}\begin{align}
  k_{1,i,j}^{C} &= \frac{1}{\Delta{y}_{i,j}^{F}}
  \frac{\Delta{y}_{i+1,j}^{G}-\Delta{y}_{i,j}^{G}}{\Delta{x}_{i,j}^{F}} \\
  k_{2,i,j}^{C} &= \frac{1}{\Delta{x}_{i,j}^{F}}
  \frac{\Delta{x}_{i,j+1}^{G}-\Delta{x}_{i,j}^{G}}{\Delta{y}_{i,j}^{F}} \\
  k_{1,i,j}^{Z} &= \frac{1}{\Delta{y}_{i,j}^{U}}
  \frac{\Delta{y}_{i,j}^{C}-\Delta{y}_{i-1,j}^{C}}{\Delta{x}_{i,j}^{V}} \\
  k_{2,i,j}^{Z} &= \frac{1}{\Delta{x}_{i,j}^{V}}
  \frac{\Delta{x}_{i,j}^{C}-\Delta{x}_{i,j-1}^{C}}{\Delta{y}_{i,j}^{U}}
\end{align}\end{linenomath*}

The stress tensor is given by the constitutive viscous-plastic
relation $\sigma_{\alpha\beta} = 2\eta\dot{\epsilon}_{\alpha\beta} +
[(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2
]\delta_{\alpha\beta}$ \citep{hibler79}. The stress tensor divergence
$(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is
discretized in finite volumes. This conveniently avoids dealing with
further metric terms, as these are ``hidden'' in the differential cell
widths. For the $u$-equation ($\alpha=1$) we have:
\begin{linenomath*}\begin{align}
  (\nabla\sigma)_{1}: \phantom{=}&
  \frac{1}{A_{i,j}^w}
  \int_{\mathrm{cell}}(\partial_1\sigma_{11}+\partial_2\sigma_{21})\,dx_1\,dx_2
  \\\notag
  =& \frac{1}{A_{i,j}^w} \biggl\{
  \int_{x_2}^{x_2+\Delta{x}_2}\sigma_{11}dx_2\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
  + \int_{x_1}^{x_1+\Delta{x}_1}\sigma_{21}dx_1\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
  \biggr\} \\ \notag
  \approx& \frac{1}{A_{i,j}^w} \biggl\{
  \Delta{x}_2\sigma_{11}\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
  + \Delta{x}_1\sigma_{21}\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
  \biggr\} \\ \notag
  =& \frac{1}{A_{i,j}^w} \biggl\{
  (\Delta{x}_2\sigma_{11})_{i,j}^C - (\Delta{x}_2\sigma_{11})_{i-1,j}^C \\\notag
  \phantom{=}& \phantom{\frac{1}{A_{i,j}^w} \biggl\{}
  + (\Delta{x}_1\sigma_{21})_{i,j+1}^Z - (\Delta{x}_1\sigma_{21})_{i,j}^Z
  \biggr\}
  \intertext{with}
  (\Delta{x}_2\sigma_{11})_{i,j}^C =& \phantom{+}
  \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
  \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag
  &+ \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
  k_{2,i,j}^C \frac{v_{i,j+1}+v_{i,j}}{2} \\ \notag
  \phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
  \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} \\ \notag
  \phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
  k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag
  \phantom{=}& - \Delta{y}_{i,j}^{F} \frac{P}{2} \\
  % 
  (\Delta{x}_1\sigma_{21})_{i,j}^Z =& \phantom{+}
  \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
  \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\ \notag
  & + \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
  \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\ \notag
  & - \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j} 
  k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2} \\ \notag
  & - \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j} 
  k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2}
\end{align}\end{linenomath*}

Similarly, we have for the $v$-equation ($\alpha=2$):
\begin{linenomath*}\begin{align}
  (\nabla\sigma)_{2}: \phantom{=}&
  \frac{1}{A_{i,j}^s}
  \int_{\mathrm{cell}}(\partial_1\sigma_{12}+\partial_2\sigma_{22})\,dx_1\,dx_2 
  \\\notag
  =& \frac{1}{A_{i,j}^s} \biggl\{
  \int_{x_2}^{x_2+\Delta{x}_2}\sigma_{12}dx_2\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
  + \int_{x_1}^{x_1+\Delta{x}_1}\sigma_{22}dx_1\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
  \biggr\} \\ \notag
  \approx& \frac{1}{A_{i,j}^s} \biggl\{
  \Delta{x}_2\sigma_{12}\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
  + \Delta{x}_1\sigma_{22}\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
  \biggr\} \\ \notag
  =& \frac{1}{A_{i,j}^s} \biggl\{
  (\Delta{x}_2\sigma_{12})_{i+1,j}^Z - (\Delta{x}_2\sigma_{12})_{i,j}^Z
  \\ \notag
  \phantom{=}& \phantom{\frac{1}{A_{i,j}^s} \biggl\{}
  + (\Delta{x}_1\sigma_{22})_{i,j}^C - (\Delta{x}_1\sigma_{22})_{i,j-1}^C
  \biggr\} 
  \intertext{with}
  (\Delta{x}_1\sigma_{12})_{i,j}^Z =& \phantom{+}
  \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
  \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\\notag
  &+ \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
  \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\ \notag
  &- \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
  k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2} \\ \notag
  &- \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
  k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2} \\ \notag
  % 
  (\Delta{x}_2\sigma_{22})_{i,j}^C =& \phantom{+}
  \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
  \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag
  &+ \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
  k_{2,i,j}^{C} \frac{v_{i,j+1}+v_{i,j}}{2} \\ \notag
  & + \Delta{x}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
  \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} \\ \notag
  & + \Delta{x}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
  k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag
  & -\Delta{x}_{i,j}^{F} \frac{P}{2}
\end{align}\end{linenomath*}


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "ceaice_part1"
%%% End: 
1	mlosch	1.7	\section{Dynamics\label{app:dynamics}}
2
3			% \newcommand{\vek}[1]{\ensuremath{\vec{\mathbf{#1}}}}
4			% \newcommand{\vtau}{\vek{\mathbf{\tau}}}
5			For completeness we provide more details on the ice dynamics of the
6			sea-ice model. The momentum equation are
7			\begin{equation}
8			\label{eq:momseaice}
9			m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} +
10			\vtau_{ocean} - m \nabla{\phi(0)} + \vek{F},
11			\end{equation}
12			where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area;
13			$\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector;
14			$\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$
15			directions, respectively;
16			$f$ is the Coriolis parameter;
17			$\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses,
18			respectively;
19			$g$ is the gravity accelation;
20			$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;
21			$\phi(0) = g\eta + p_{a}/\rho_{0} + mg/\rho_{0}$ is the sea surface
22			height potential in response to ocean dynamics ($g\eta$), to
23			atmospheric pressure loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a
24			reference density) and a term due to snow and ice loading \citep{campin08};
25			and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice
26			stress tensor $\sigma_{ij}$. %
27			Advection of sea-ice momentum is neglected. The wind and ice-ocean stress
28			terms are given by
29			\begin{align*}
30			\vtau_{air} = & \rho_{air} C_{air} \|\vek{U}_{air} -\vek{u}\|
31			R_{air} (\vek{U}_{air} -\vek{u}), \\
32			\vtau_{ocean} = & \rho_{ocean}C_{ocean} \|\vek{U}_{ocean}-\vek{u}\|
33			R_{ocean}(\vek{U}_{ocean}-\vek{u}),
34			\end{align*}
35			where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere
36			and surface currents of the ocean, respectively; $C_{air/ocean}$ are
37			air and ocean drag coefficients; $\rho_{air/ocean}$ are reference
38			densities; and $R_{air/ocean}$ are rotation matrices that act on the
39			wind/current vectors.
40
41			For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can
42			be related to the ice strain rate and strength by a nonlinear
43			viscous-plastic (VP) constitutive law \citep{hibler79, zhang97}:
44			\begin{equation}
45			\label{eq:vpequation}
46			\sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij}
47			+ \left[\zeta(\dot{\epsilon}_{ij},P) -
48			\eta(\dot{\epsilon}_{ij},P)\right]\dot{\epsilon}_{kk}\delta_{ij}
49			- \frac{P}{2}\delta_{ij}.
50			\end{equation}
51			The ice strain rate is given by
52			\begin{equation*}
53			\dot{\epsilon}_{ij} = \frac{1}{2}\left(
54			\frac{\partial{u_{i}}}{\partial{x_{j}}} +
55			\frac{\partial{u_{j}}}{\partial{x_{i}}}\right).
56			\end{equation*}
57			The maximum ice pressure $P_{\max}$, a measure of ice strength, depends on
58			both thickness $h$ and compactness (concentration) $c$:
59			\begin{equation}
60			P_{\max} = P^{}c\,h\,e^{[C^{}\cdot(1-c)]},
61			\label{eq:icestrength}
62			\end{equation}
63			with the constants $P^{}$ and $C^{}=20$. The nonlinear bulk and shear
64			viscosities $\eta$ and $\zeta$ are functions of ice strain rate
65			invariants and ice strength such that the principal components of the
66			stress lie on an elliptical yield curve with the ratio of major to
67			minor axis $e$ equal to $2$; they are given by:
68			\begin{align*}
69			\zeta =& \min\left(\frac{P_{\max}}{2\max(\Delta,\Delta_{\min})},
70			\zeta_{\max}\right) \\
71			\eta =& \frac{\zeta}{e^2} \\
72			\intertext{with the abbreviation}
73			\Delta = & \left[
74			\left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right)
75			(1+e^{-2}) + 4e^{-2}\dot{\epsilon}_{12}^2 +
76			2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2})
77			\right]^{\frac{1}{2}}.
78			\end{align*}
79			In the simulations of this paper, the bulk viscosities are bounded
80			above by imposing both a minimum
81			$\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ and a maximum $\zeta_{\max} =
82			P_{\max}/\Delta^*$, where
83			$\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress
84			tensor computation the replacement pressure $P = 2\,\Delta\zeta$
85			\citep{hibler95} is used so that the stress state always lies on the
86			elliptic yield curve by definition.
87
88			In the so-called truncated ellipse method (experiment TEM) the shear
89			viscosity $\eta$ is capped to suppress any tensile stress
90			\citep{hibler97, geiger98}:
91			\begin{equation}
92			\label{eq:etatem}
93			\eta = \min\left(\frac{\zeta}{e^2},
94			\frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
95			{\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2
96			+4\dot{\epsilon}_{12}^2}}\right).
97			\end{equation}
98
99			In the current implementation, the VP-model is integrated with the
100			semi-implicit line successive over relaxation (LSOR)-solver of
101			\citet{zhang97}, which allows for long time steps that, in our case,
102			are limited by the explicit treatment of the Coriolis term. The
103			explicit treatment of the Coriolis term does not represent a severe
104			limitation because it restricts the time step to approximately the
105			same length as in the ocean model where the Coriolis term is also
106			treated explicitly.
107
108			\citet{hunke97}'s introduced an elastic contribution to the strain
109			rate in order to regularize Eq.~\ref{eq:vpequation} in such a way that
110			the resulting elastic-viscous-plastic (EVP) and VP models are
111			identical at steady state,
112			\begin{equation}
113			\label{eq:evpequation}
114			\frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} +
115			\frac{1}{2\eta}\sigma_{ij}
116			+ \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij}
117			+ \frac{P}{4\zeta}\delta_{ij}
118			= \dot{\epsilon}_{ij}.
119			\end{equation}
120			The EVP-model uses an explicit time stepping scheme with a short
121			timestep. According to the recommendation of \citet{hunke97}, the
122			EVP-model is stepped forward in time O(120) times within the physical
123			ocean model time step, to
124			allow for elastic waves to disappear. Because the scheme does not
125			require a matrix inversion it is fast in spite of the small internal
126			timestep and simple to implement on parallel computers
127			\citep{hunke97}. For completeness, we repeat the equations for the
128			components of the stress tensor $\sigma_{1} =
129			\sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and
130			$\sigma_{12}$. Introducing the divergence $D_D =
131			\dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension
132			and shearing strain rates, $D_T =
133			\dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =
134			2\dot{\epsilon}_{12}$, respectively, and using the above
135			abbreviations, the equations~\ref{eq:evpequation} can be written as:
136			\begin{align}
137			\label{eq:evpstresstensor1}
138			\frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +
139			\frac{P}{2T} &= \frac{P}{2T\Delta} D_D \\
140			\label{eq:evpstresstensor2}
141			\frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T}
142			&= \frac{P}{2T\Delta} D_T \\
143			\label{eq:evpstresstensor12}
144			\frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T}
145			&= \frac{P}{4T\Delta} D_S
146			\end{align}
147			Here, the elastic parameter $E$ is redefined in terms of a damping
148			timescale $T$ for elastic waves \[E=\frac{\zeta}{T}.\]
149			$T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and the external
150			(long) timestep $\Delta{t}$. We use $E_{0} = \frac{1}{3}$ which is
151			close to value of 0.36 used by \citet{hunke01}.
152
153	mlosch	1.3	\section{Finite-volume discretization of the stress tensor
154			divergence}
155			\label{app:discretization}
156			On an Arakawa C~grid, ice thickness and concentration and thus ice
157			strength $P$ and bulk and shear viscosities $\zeta$ and $\eta$ are
158	mlosch	1.5	naturally defined at C-points in the center of the grid
159	mlosch	1.4	cell. Discretization requires only averaging of $\zeta$ and $\eta$ to
160			vorticity or Z-points at the bottom left corner of the cell to give
161			$\overline{\zeta}^{Z}$ and $\overline{\eta}^{Z}$. In the following,
162			the superscripts indicate location at Z or C points, distance across
163			the cell (F), along the cell edge (G), between $u$-points (U),
164			$v$-points (V), and C-points (C). The control volumes of the $u$- and
165			$v$-equations in the grid cell at indices $(i,j)$ are $A_{i,j}^{w}$
166			and $A_{i,j}^{s}$, respectively. With these definitions (which follow
167			the model code documentation at \url{http://mitgcm.org} except that
168	mlosch	1.6	vorticity or $\zeta$-points have been renamed to Z-points in order to
169			avoid confusion with the bulk viscosity $\zeta$), the strain
170	mlosch	1.5	rates are discretized as:
171	mlosch	1.4	\begin{linenomath*}\begin{align}
172	mlosch	1.3	\dot{\epsilon}_{11} &= \partial_{1}{u}_{1} + k_{2}u_{2} \\ \notag
173			=> (\epsilon_{11})_{i,j}^C &= \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}}
174			+ k_{2,i,j}^{C}\frac{v_{i,j+1}+v_{i,j}}{2} \\
175			\dot{\epsilon}_{22} &= \partial_{2}{u}_{2} + k_{1}u_{1} \\\notag
176			=> (\epsilon_{22})_{i,j}^C &= \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}}
177			+ k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\
178			\dot{\epsilon}_{12} = \dot{\epsilon}_{21} &= \frac{1}{2}\biggl(
179			\partial_{1}{u}_{2} + \partial_{2}{u}_{1} - k_{1}u_{2} - k_{2}u_{1}
180			\biggr) \\ \notag
181			=> (\epsilon_{12})_{i,j}^Z &= \frac{1}{2}
182			\biggl( \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^V}
183			+ \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^U} \\\notag
184			&\phantom{=\frac{1}{2}\biggl(}
185			- k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2}
186			- k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2}
187			\biggr),
188	mlosch	1.4	\end{align}\end{linenomath*}
189	mlosch	1.3	so that the diagonal terms of the strain rate tensor are naturally
190			defined at C-points and the symmetric off-diagonal term at
191			Z-points. No-slip boundary conditions ($u_{i,j-1}+u_{i,j}=0$ and
192			$v_{i-1,j}+v_{i,j}=0$ across boundaries) are implemented via
193			``ghost-points''; for free slip boundary conditions
194			$(\epsilon_{12})^Z=0$ on boundaries.
195
196			For a spherical polar grid, the coefficients of the metric terms are
197			$k_{1}=0$ and $k_{2}=-\tan\phi/a$, with the spherical radius $a$ and
198			the latitude $\phi$; $\Delta{x}_1 = \Delta{x} = a\cos\phi
199			\Delta\lambda$, and $\Delta{x}_2 = \Delta{y}=a\Delta\phi$. For a
200			general orthogonal curvilinear grid as used in this paper, $k_{1}$ and
201			$k_{2}$ can be approximated by finite differences of the cell widths:
202			%\enlargethispage{1cm}
203	mlosch	1.4	\begin{linenomath*}\begin{align}
204	mlosch	1.3	k_{1,i,j}^{C} &= \frac{1}{\Delta{y}_{i,j}^{F}}
205			\frac{\Delta{y}_{i+1,j}^{G}-\Delta{y}_{i,j}^{G}}{\Delta{x}_{i,j}^{F}} \\
206			k_{2,i,j}^{C} &= \frac{1}{\Delta{x}_{i,j}^{F}}
207			\frac{\Delta{x}_{i,j+1}^{G}-\Delta{x}_{i,j}^{G}}{\Delta{y}_{i,j}^{F}} \\
208			k_{1,i,j}^{Z} &= \frac{1}{\Delta{y}_{i,j}^{U}}
209			\frac{\Delta{y}_{i,j}^{C}-\Delta{y}_{i-1,j}^{C}}{\Delta{x}_{i,j}^{V}} \\
210			k_{2,i,j}^{Z} &= \frac{1}{\Delta{x}_{i,j}^{V}}
211			\frac{\Delta{x}_{i,j}^{C}-\Delta{x}_{i,j-1}^{C}}{\Delta{y}_{i,j}^{U}}
212	mlosch	1.4	\end{align}\end{linenomath*}
213	heimbach	1.1
214	mlosch	1.3	The stress tensor is given by the constitutive viscous-plastic
215			relation $\sigma_{\alpha\beta} = 2\eta\dot{\epsilon}_{\alpha\beta} +
216			[(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2
217			]\delta_{\alpha\beta}$ \citep{hibler79}. The stress tensor divergence
218			$(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is
219			discretized in finite volumes. This conveniently avoids dealing with
220			further metric terms, as these are ``hidden'' in the differential cell
221			widths. For the $u$-equation ($\alpha=1$) we have:
222	mlosch	1.4	\begin{linenomath*}\begin{align}
223	mlosch	1.3	(\nabla\sigma)_{1}: \phantom{=}&
224			\frac{1}{A_{i,j}^w}
225			\int_{\mathrm{cell}}(\partial_1\sigma_{11}+\partial_2\sigma_{21})\,dx_1\,dx_2
226			\\\notag
227			=& \frac{1}{A_{i,j}^w} \biggl\{
228			\int_{x_2}^{x_2+\Delta{x}_2}\sigma_{11}dx_2\biggl\|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
229			+ \int_{x_1}^{x_1+\Delta{x}_1}\sigma_{21}dx_1\biggl\|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
230			\biggr\} \\ \notag
231			\approx& \frac{1}{A_{i,j}^w} \biggl\{
232			\Delta{x}_2\sigma_{11}\biggl\|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
233			+ \Delta{x}_1\sigma_{21}\biggl\|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
234			\biggr\} \\ \notag
235			=& \frac{1}{A_{i,j}^w} \biggl\{
236			(\Delta{x}_2\sigma_{11})_{i,j}^C - (\Delta{x}_2\sigma_{11})_{i-1,j}^C \\\notag
237			\phantom{=}& \phantom{\frac{1}{A_{i,j}^w} \biggl\{}
238			+ (\Delta{x}_1\sigma_{21})_{i,j+1}^Z - (\Delta{x}_1\sigma_{21})_{i,j}^Z
239			\biggr\}
240			\intertext{with}
241			(\Delta{x}_2\sigma_{11})_{i,j}^C =& \phantom{+}
242			\Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
243			\frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag
244			&+ \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
245			k_{2,i,j}^C \frac{v_{i,j+1}+v_{i,j}}{2} \\ \notag
246			\phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
247			\frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} \\ \notag
248			\phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
249			k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag
250			\phantom{=}& - \Delta{y}_{i,j}^{F} \frac{P}{2} \\
251			%
252			(\Delta{x}_1\sigma_{21})_{i,j}^Z =& \phantom{+}
253			\Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
254			\frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\ \notag
255			& + \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
256			\frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\ \notag
257			& - \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
258			k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2} \\ \notag
259			& - \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
260			k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2}
261	mlosch	1.4	\end{align}\end{linenomath*}
262	heimbach	1.1
263	mlosch	1.3	Similarly, we have for the $v$-equation ($\alpha=2$):
264	mlosch	1.4	\begin{linenomath*}\begin{align}
265	mlosch	1.3	(\nabla\sigma)_{2}: \phantom{=}&
266			\frac{1}{A_{i,j}^s}
267			\int_{\mathrm{cell}}(\partial_1\sigma_{12}+\partial_2\sigma_{22})\,dx_1\,dx_2
268			\\\notag
269			=& \frac{1}{A_{i,j}^s} \biggl\{
270			\int_{x_2}^{x_2+\Delta{x}_2}\sigma_{12}dx_2\biggl\|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
271			+ \int_{x_1}^{x_1+\Delta{x}_1}\sigma_{22}dx_1\biggl\|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
272			\biggr\} \\ \notag
273			\approx& \frac{1}{A_{i,j}^s} \biggl\{
274			\Delta{x}_2\sigma_{12}\biggl\|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
275			+ \Delta{x}_1\sigma_{22}\biggl\|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
276			\biggr\} \\ \notag
277			=& \frac{1}{A_{i,j}^s} \biggl\{
278			(\Delta{x}_2\sigma_{12})_{i+1,j}^Z - (\Delta{x}_2\sigma_{12})_{i,j}^Z
279			\\ \notag
280			\phantom{=}& \phantom{\frac{1}{A_{i,j}^s} \biggl\{}
281			+ (\Delta{x}_1\sigma_{22})_{i,j}^C - (\Delta{x}_1\sigma_{22})_{i,j-1}^C
282			\biggr\}
283			\intertext{with}
284			(\Delta{x}_1\sigma_{12})_{i,j}^Z =& \phantom{+}
285			\Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
286			\frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\\notag
287			&+ \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
288			\frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\ \notag
289			&- \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
290			k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2} \\ \notag
291			&- \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
292			k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2} \\ \notag
293			%
294			(\Delta{x}_2\sigma_{22})_{i,j}^C =& \phantom{+}
295			\Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
296			\frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag
297			&+ \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
298			k_{2,i,j}^{C} \frac{v_{i,j+1}+v_{i,j}}{2} \\ \notag
299			& + \Delta{x}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
300			\frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} \\ \notag
301			& + \Delta{x}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
302			k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag
303			& -\Delta{x}_{i,j}^{F} \frac{P}{2}
304	mlosch	1.4	\end{align}\end{linenomath*}
305	heimbach	1.1
306
307			%%% Local Variables:
308			%%% mode: latex
309	mlosch	1.3	%%% TeX-master: "ceaice_part1"
310	heimbach	1.1	%%% End: