articles/ceaice/ceaice_model.tex

\section{Model Formulation}
\label{sec:model}

The MITgcm sea ice model (MITsim) is based on a variant of the
viscous-plastic (VP) dynamic-thermodynamic sea ice model
\citep{zhang97} first introduced by \citet{hibler79, hibler80}. In
order to adapt this model to the requirements of coupled
ice-ocean state estimation, many important aspects of the original code have
been modified and improved:
\begin{itemize}
\item the code has been rewritten for an Arakawa C-grid, both B- and
  C-grid variants are available; the C-grid code allows for no-slip
  and free-slip lateral boundary conditions;
\item two different solution methods for solving the nonlinear
  momentum equations have been adopted: LSOR \citep{zhang97}, and EVP
  \citep{hunke97};
\item ice-ocean stress can be formulated as in \citet{hibler87} or as in \citet{cam08}; 
\item ice variables \ml{can be} advected by sophisticated, \ml{conservative}
  advection schemes \ml{with flux limiting};
\item growth and melt parameterizations have been refined and extended
  in order to allow for more stable automatic differentiation of the code.
\end{itemize}

The sea ice model is tightly coupled to the ocean compontent of the MITgcm
\citep{mar97a}.  Heat, fresh water fluxes and surface stresses are computed
from the atmospheric state and modified by the ice model at every time step.


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

The momentum equation of the sea-ice model is
\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^{*}$. 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*}
The bulk viscosities are bounded above by imposing both a minimum
$\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) 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 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{zhang98}, 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.\refeq{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}
%In the EVP model, equations for the components of the stress tensor
%$\sigma_{ij}$ are solved explicitly. Both model formulations will be
%used and compared the present sea-ice model study.
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 120 times within the physical
ocean model time step (although this parameter is under debate), 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\refeq{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}$. \citet{hunke97} recommend
$E_{0} = \frac{1}{3}$.

Our discretization differs from \citet{zhang98, zhang03} in the
underlying grid, namely the Arakawa C-grid, but is otherwise
straightforward. The EVP model, in particular, is discretized
naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the
center points and $\sigma_{12}$ on the corner (or vorticity) points of
the grid. With this choice all derivatives are discretized as central
differences and averaging is only involved in computing $\Delta$ and
$P$ at vorticity points.

For a general curvilinear grid, one needs in principle to take metric
terms into account that arise in the transformation of a curvilinear
grid on the sphere. For now, however, we can neglect these metric
terms because they are very small on the \ml{[modify following
  section3:] cubed sphere grids used in this paper; in particular,
only near the edges of the cubed sphere grid, we expect them to be
non-zero, but these edges are at approximately 35\degS\ or 35\degN\ 
which are never covered by sea-ice in our simulations.  Everywhere
else the coordinate system is locally nearly cartesian.}  However, for
last-glacial-maximum or snowball-earth-like simulations the question
of metric terms needs to be reconsidered.  Either, one includes these
terms as in \citet{zhang03}, or one finds a vector-invariant
formulation for the sea-ice internal stress term that does not require
any metric terms, as it is done in the ocean dynamics of the MITgcm
\citep{adcroft04:_cubed_sphere}.

Lateral boundary conditions are naturally ``no-slip'' for B-grids, as
the tangential velocities points lie on the boundary. For C-grids, the
lateral boundary condition for tangential velocities is realized via
``ghost points'', allowing alternatively no-slip or free-slip
conditions. In ocean models free-slip boundary conditions in
conjunction with piecewise-constant (``castellated'') coastlines have
been shown to reduce in effect to no-slip boundary conditions
\citep{adcroft98:_slippery_coast}; for sea-ice models the effects of
lateral boundary conditions have not yet been studied.

Moving sea ice exerts a stress on the ocean which is the opposite of
the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is
applied directly to the surface layer of the ocean model. An
alternative ocean stress formulation is given by \citet{hibler87}.
Rather than applying $\vtau_{ocean}$ directly, the stress is derived
from integrating over the ice thickness to the bottom of the oceanic
surface layer. In the resulting equation for the \emph{combined}
ocean-ice momentum, the interfacial stress cancels and the total
stress appears as the sum of windstress and divergence of internal ice
stresses: $\delta(z) (\vtau_{air} + \vek{F})/\rho_0$, \citep[see also
Eq.\,2 of][]{hibler87}. The disadvantage of this formulation is that
now the velocity in the surface layer of the ocean that is used to
advect tracers, is really an average over the ocean surface
velocity and the ice velocity leading to an inconsistency as the ice
temperature and salinity are different from the oceanic variables.

%\subparagraph{boundary conditions: no-slip, free-slip, half-slip}
%\begin{itemize}
%\item transition from existing B-Grid to C-Grid
%\item boundary conditions: no-slip, free-slip, half-slip
%\item fancy (multi dimensional) advection schemes
%\item VP vs.\ EVP \citep{hunke97}
%\item ocean stress formulation \citep{hibler87}
%\end{itemize}

\subsection{Thermodynamics}
\label{app:thermodynamics}

In its original formulation the sea ice model \citep{menemenlis05}
uses simple thermodynamics following the appendix of
\citet{semtner76}. This formulation does not allow storage of heat,
that is, the heat capacity of ice is zero. Upward conductive heat flux
is parameterized assuming a linear temperature profile and together
with a constant ice conductivity. It is expressed as
$(K/h)(T_{w}-T_{0})$, where $K$ is the ice conductivity, $h$ the ice
thickness, and $T_{w}-T_{0}$ the difference between water and ice
surface temperatures. This type of model is often refered to as a
``zero-layer'' model. The surface heat flux is computed in a similar
way to that of \citet{parkinson79} and \citet{manabe79}. 

The conductive heat flux depends strongly on the ice thickness $h$.
However, the ice thickness in the model represents a mean over a
potentially very heterogeneous thickness distribution.  In order to
parameterize a sub-grid scale distribution for heat flux
computations, the mean ice thickness $h$ is split into seven thickness
categories $H_{n}$ that are equally distributed between $2h$ and a
minimum imposed ice thickness of $5\text{\,cm}$ by $H_n=
\frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each
thickness category is area-averaged to give the total heat flux
\citep{hibler84}.

\ml{[This is Ian Fenty's work and we may want to remove this paragraph
  from the paper]: %
The atmospheric heat flux is balanced by an oceanic heat flux from
below.  The oceanic flux is proportional to
$\rho\,c_{p}\left(T_{w}-T_{fr}\right)$ where $\rho$ and $c_{p}$ are
the density and heat capacity of sea water and $T_{fr}$ is the local
freezing point temperature that is a function of salinity. Contrary to
\citet{menemenlis05}, this flux is not assumed to instantaneously melt
or create ice, but a time scale of three days is used to relax $T_{w}$
to the freezing point.}
%
The parameterization of lateral and vertical growth of sea ice follows
that of \citet{hibler79, hibler80}.

On top of the ice there is a layer of snow that modifies the heat flux
and the albedo \citep{zhang98}. Snow modifies the effective
conductivity according to 
\[\frac{K}{h} \rightarrow \frac{1}{\frac{h_{s}}{K_{s}}+\frac{h}{K}},\]
where $K_s$ is the conductivity of snow and $h_s$ the snow thickness.
If enough snow accumulates so that its weight submerges the ice and
the snow is flooded, a simple mass conserving parameterization of
snowice formation (a flood-freeze algorithm following Archimedes'
principle) turns snow into ice until the ice surface is back at $z=0$
\citep{leppaeranta83}.

Effective ice thickness (ice volume per unit area,
$c\cdot{h}$), concentration $c$ and effective snow thickness
($c\cdot{h}_{s}$) are advected by ice velocities:
\begin{equation}
  \label{eq:advection}
  \frac{\partial{X}}{\partial{t}} = - \nabla\cdot\left(\vek{u}\,X\right) +
  \Gamma_{X} + D_{X}
\end{equation}
where $\Gamma_X$ are the thermodynamic source terms and $D_{X}$ the
diffusive terms for quantities $X=(c\cdot{h}), c, (c\cdot{h}_{s})$.
%
From the various advection scheme that are available in the MITgcm
\citep{mitgcm02}, we choose flux-limited schemes
\citep[multidimensional 2nd and 3rd-order advection scheme with flux
limiter][]{roe85, hundsdorfer94} to preserve sharp gradients and edges
that are typical of sea ice distributions and to rule out unphysical
over- and undershoots (negative thickness or concentration). These
scheme conserve volume and horizontal area and are unconditionally
stable, so that we can set $D_{X}=0$.  \ml{[do we need to proove that?
  can we proove that? citation?]}

There is considerable doubt about the reliability of such a simple
thermodynamic model---\citet{semtner84} found significant errors in
phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in
such models---, so that today many sea ice models employ more complex
thermodynamics. A popular thermodynamics model of \citet{winton00} is
based on the 3-layer model of \citet{semtner76} and treats brine
content by means of enthalphy conservation. This model requires in
addition to ice-thickness and compactness (fractional area) additional
state variables to be advected by ice velocities, namely enthalphy of
the two ice layers and the thickness of the overlying snow layer.
\ml{[Jean-Michel, your turf: ]Care must be taken in advecting these
  quantities in order to ensure conservation of enthalphy. Currently
  this can only be accomplished with a 2nd-order advection scheme with
  flux limiter \citep{roe85}.}

Further documentation and model code can be found at
\url{http://mitgcm.org}.


%\subsection{C-grid}
%\begin{itemize}
%\item no-slip vs. free-slip for both lsr and evp; 
%  "diagnostics" to look at and use for comparison
%  \begin{itemize}
%  \item ice transport through Fram Strait/Denmark Strait/Davis
%    Strait/Bering strait (these are general) 
%  \item ice transport through narrow passages, e.g.\ Nares-Strait
%  \end{itemize}
%\item compare different advection schemes (if lsr turns out to be more
%  effective, then with lsr otherwise I prefer evp), eg. 
%  \begin{itemize}
%  \item default 2nd-order with diff1=0.002
%  \item 1st-order upwind with diff1=0.
%  \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)
%  \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)
%  \end{itemize}
%  That should be enough. Here, total ice mass and location of ice edge
%  is interesting. However, this comparison can be done in an idealized
%  domain, may not require full Arctic Domain?
%\item
%Do a little study on the parameters of LSR and EVP
%\begin{enumerate}
%\item convergence of LSR, how many iterations do you need to get a
%  true elliptic yield curve 
%\item vary deltaTevp and the relaxation parameter for EVP and see when
%  the EVP solution breaks down (relative to the forcing time scale).
%  For this, it is essential that the evp solver gives use "stripeless"
%  solutions, that is your dtevp = 1sec solutions/or 10sec solutions
%  with SEAICE\_evpDampC = 615.
%\end{enumerate}

%\end{itemize}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "ceaice"
%%% End: 
1	mlosch	1.4	\section{Model Formulation}
2	dimitri	1.1	\label{sec:model}
3
4	mlosch	1.4	The MITgcm sea ice model (MITsim) is based on a variant of the
5			viscous-plastic (VP) dynamic-thermodynamic sea ice model
6			\citep{zhang97} first introduced by \citet{hibler79, hibler80}. In
7			order to adapt this model to the requirements of coupled
8	dimitri	1.11	ice-ocean state estimation, many important aspects of the original code have
9	mlosch	1.4	been modified and improved:
10			\begin{itemize}
11			\item the code has been rewritten for an Arakawa C-grid, both B- and
12			C-grid variants are available; the C-grid code allows for no-slip
13			and free-slip lateral boundary conditions;
14			\item two different solution methods for solving the nonlinear
15	dimitri	1.11	momentum equations have been adopted: LSOR \citep{zhang97}, and EVP
16	mlosch	1.4	\citep{hunke97};
17	dimitri	1.11	\item ice-ocean stress can be formulated as in \citet{hibler87} or as in \citet{cam08};
18	mlosch	1.10	\item ice variables \ml{can be} advected by sophisticated, \ml{conservative}
19			advection schemes \ml{with flux limiting};
20			\item growth and melt parameterizations have been refined and extended
21	dimitri	1.11	in order to allow for more stable automatic differentiation of the code.
22	mlosch	1.4	\end{itemize}
23	dimitri	1.1
24	dimitri	1.11	The sea ice model is tightly coupled to the ocean compontent of the MITgcm
25			\citep{mar97a}. Heat, fresh water fluxes and surface stresses are computed
26			from the atmospheric state and modified by the ice model at every time step.
27	dimitri	1.12
28
29			\subsection{Dynamics}
30			\label{app:dynamics}
31
32			The momentum equation of the sea-ice model is
33			\begin{equation}
34			\label{eq:momseaice}
35			m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} +
36			\vtau_{ocean} - m \nabla{\phi(0)} + \vek{F},
37			\end{equation}
38			where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area;
39			$\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector;
40			$\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$
41			directions, respectively;
42			$f$ is the Coriolis parameter;
43			$\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses,
44			respectively;
45			$g$ is the gravity accelation;
46			$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;
47			$\phi(0) = g\eta + p_{a}/\rho_{0} + mg/\rho_{0}$ is the sea surface
48			height potential in response to ocean dynamics ($g\eta$), to
49			atmospheric pressure loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a
50			reference density) and a term due to snow and ice loading \citep{campin08};
51			and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice
52			stress tensor $\sigma_{ij}$. %
53			Advection of sea-ice momentum is neglected. The wind and ice-ocean stress
54			terms are given by
55			\begin{align*}
56			\vtau_{air} = & \rho_{air} C_{air} \|\vek{U}_{air} -\vek{u}\|
57			R_{air} (\vek{U}_{air} -\vek{u}), \\
58			\vtau_{ocean} = & \rho_{ocean}C_{ocean} \|\vek{U}_{ocean}-\vek{u}\|
59			R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\
60			\end{align*}
61			where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere
62			and surface currents of the ocean, respectively; $C_{air/ocean}$ are
63			air and ocean drag coefficients; $\rho_{air/ocean}$ are reference
64			densities; and $R_{air/ocean}$ are rotation matrices that act on the
65			wind/current vectors.
66
67			For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can
68			be related to the ice strain rate and strength by a nonlinear
69			viscous-plastic (VP) constitutive law \citep{hibler79, zhang97}:
70			\begin{equation}
71			\label{eq:vpequation}
72			\sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij}
73			+ \left[\zeta(\dot{\epsilon}_{ij},P) -
74			\eta(\dot{\epsilon}_{ij},P)\right]\dot{\epsilon}_{kk}\delta_{ij}
75			- \frac{P}{2}\delta_{ij}.
76			\end{equation}
77			The ice strain rate is given by
78			\begin{equation*}
79			\dot{\epsilon}_{ij} = \frac{1}{2}\left(
80			\frac{\partial{u_{i}}}{\partial{x_{j}}} +
81			\frac{\partial{u_{j}}}{\partial{x_{i}}}\right).
82			\end{equation*}
83			The maximum ice pressure $P_{\max}$, a measure of ice strength, depends on
84			both thickness $h$ and compactness (concentration) $c$:
85			\begin{equation}
86			P_{\max} = P^{}c\,h\,e^{[C^{}\cdot(1-c)]},
87			\label{eq:icestrength}
88			\end{equation}
89			with the constants $P^{}$ and $C^{}$. The nonlinear bulk and shear
90			viscosities $\eta$ and $\zeta$ are functions of ice strain rate
91			invariants and ice strength such that the principal components of the
92			stress lie on an elliptical yield curve with the ratio of major to
93			minor axis $e$ equal to $2$; they are given by:
94			\begin{align*}
95			\zeta =& \min\left(\frac{P_{\max}}{2\max(\Delta,\Delta_{\min})},
96			\zeta_{\max}\right) \\
97			\eta =& \frac{\zeta}{e^2} \\
98			\intertext{with the abbreviation}
99			\Delta = & \left[
100			\left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right)
101			(1+e^{-2}) + 4e^{-2}\dot{\epsilon}_{12}^2 +
102			2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2})
103			\right]^{\frac{1}{2}}.
104			\end{align*}
105			The bulk viscosities are bounded above by imposing both a minimum
106			$\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a
107			maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where
108			$\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress
109			tensor computation the replacement pressure $P = 2\,\Delta\zeta$
110			\citep{hibler95} is used so that the stress state always lies on the
111			elliptic yield curve by definition.
112
113			In the so-called truncated ellipse method the shear viscosity $\eta$
114			is capped to suppress any tensile stress \citep{hibler97, geiger98}:
115			\begin{equation}
116			\label{eq:etatem}
117			\eta = \min\left(\frac{\zeta}{e^2},
118			\frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
119			{\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2
120			+4\dot{\epsilon}_{12}^2}}\right).
121			\end{equation}
122
123			In the current implementation, the VP-model is integrated with the
124			semi-implicit line successive over relaxation (LSOR)-solver of
125			\citet{zhang98}, which allows for long time steps that, in our case,
126			are limited by the explicit treatment of the Coriolis term. The
127			explicit treatment of the Coriolis term does not represent a severe
128			limitation because it restricts the time step to approximately the
129			same length as in the ocean model where the Coriolis term is also
130			treated explicitly.
131
132			\citet{hunke97}'s introduced an elastic contribution to the strain
133			rate in order to regularize Eq.\refeq{vpequation} in such a way that
134			the resulting elastic-viscous-plastic (EVP) and VP models are
135			identical at steady state,
136			\begin{equation}
137			\label{eq:evpequation}
138			\frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} +
139			\frac{1}{2\eta}\sigma_{ij}
140			+ \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij}
141			+ \frac{P}{4\zeta}\delta_{ij}
142			= \dot{\epsilon}_{ij}.
143			\end{equation}
144			%In the EVP model, equations for the components of the stress tensor
145			%$\sigma_{ij}$ are solved explicitly. Both model formulations will be
146			%used and compared the present sea-ice model study.
147			The EVP-model uses an explicit time stepping scheme with a short
148			timestep. According to the recommendation of \citet{hunke97}, the
149			EVP-model is stepped forward in time 120 times within the physical
150			ocean model time step (although this parameter is under debate), to
151			allow for elastic waves to disappear. Because the scheme does not
152			require a matrix inversion it is fast in spite of the small internal
153			timestep and simple to implement on parallel computers
154			\citep{hunke97}. For completeness, we repeat the equations for the
155			components of the stress tensor $\sigma_{1} =
156			\sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and
157			$\sigma_{12}$. Introducing the divergence $D_D =
158			\dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension
159			and shearing strain rates, $D_T =
160			\dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =
161			2\dot{\epsilon}_{12}$, respectively, and using the above
162			abbreviations, the equations\refeq{evpequation} can be written as:
163			\begin{align}
164			\label{eq:evpstresstensor1}
165			\frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +
166			\frac{P}{2T} &= \frac{P}{2T\Delta} D_D \\
167			\label{eq:evpstresstensor2}
168			\frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T}
169			&= \frac{P}{2T\Delta} D_T \\
170			\label{eq:evpstresstensor12}
171			\frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T}
172			&= \frac{P}{4T\Delta} D_S
173			\end{align}
174			Here, the elastic parameter $E$ is redefined in terms of a damping timescale
175			$T$ for elastic waves \[E=\frac{\zeta}{T}.\]
176			$T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and
177			the external (long) timestep $\Delta{t}$. \citet{hunke97} recommend
178			$E_{0} = \frac{1}{3}$.
179
180			Our discretization differs from \citet{zhang98, zhang03} in the
181			underlying grid, namely the Arakawa C-grid, but is otherwise
182			straightforward. The EVP model, in particular, is discretized
183			naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the
184			center points and $\sigma_{12}$ on the corner (or vorticity) points of
185			the grid. With this choice all derivatives are discretized as central
186			differences and averaging is only involved in computing $\Delta$ and
187			$P$ at vorticity points.
188
189			For a general curvilinear grid, one needs in principle to take metric
190			terms into account that arise in the transformation of a curvilinear
191			grid on the sphere. For now, however, we can neglect these metric
192			terms because they are very small on the \ml{[modify following
193			section3:] cubed sphere grids used in this paper; in particular,
194			only near the edges of the cubed sphere grid, we expect them to be
195			non-zero, but these edges are at approximately 35\degS\ or 35\degN\
196			which are never covered by sea-ice in our simulations. Everywhere
197			else the coordinate system is locally nearly cartesian.} However, for
198			last-glacial-maximum or snowball-earth-like simulations the question
199			of metric terms needs to be reconsidered. Either, one includes these
200			terms as in \citet{zhang03}, or one finds a vector-invariant
201			formulation for the sea-ice internal stress term that does not require
202			any metric terms, as it is done in the ocean dynamics of the MITgcm
203			\citep{adcroft04:_cubed_sphere}.
204
205			Lateral boundary conditions are naturally ``no-slip'' for B-grids, as
206			the tangential velocities points lie on the boundary. For C-grids, the
207			lateral boundary condition for tangential velocities is realized via
208			``ghost points'', allowing alternatively no-slip or free-slip
209			conditions. In ocean models free-slip boundary conditions in
210			conjunction with piecewise-constant (``castellated'') coastlines have
211			been shown to reduce in effect to no-slip boundary conditions
212			\citep{adcroft98:_slippery_coast}; for sea-ice models the effects of
213			lateral boundary conditions have not yet been studied.
214
215			Moving sea ice exerts a stress on the ocean which is the opposite of
216			the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is
217			applied directly to the surface layer of the ocean model. An
218			alternative ocean stress formulation is given by \citet{hibler87}.
219			Rather than applying $\vtau_{ocean}$ directly, the stress is derived
220			from integrating over the ice thickness to the bottom of the oceanic
221			surface layer. In the resulting equation for the \emph{combined}
222			ocean-ice momentum, the interfacial stress cancels and the total
223			stress appears as the sum of windstress and divergence of internal ice
224			stresses: $\delta(z) (\vtau_{air} + \vek{F})/\rho_0$, \citep[see also
225			Eq.\,2 of][]{hibler87}. The disadvantage of this formulation is that
226			now the velocity in the surface layer of the ocean that is used to
227			advect tracers, is really an average over the ocean surface
228			velocity and the ice velocity leading to an inconsistency as the ice
229			temperature and salinity are different from the oceanic variables.
230
231			%\subparagraph{boundary conditions: no-slip, free-slip, half-slip}
232			%\begin{itemize}
233			%\item transition from existing B-Grid to C-Grid
234			%\item boundary conditions: no-slip, free-slip, half-slip
235			%\item fancy (multi dimensional) advection schemes
236			%\item VP vs.\ EVP \citep{hunke97}
237			%\item ocean stress formulation \citep{hibler87}
238			%\end{itemize}
239
240			\subsection{Thermodynamics}
241			\label{app:thermodynamics}
242
243			In its original formulation the sea ice model \citep{menemenlis05}
244			uses simple thermodynamics following the appendix of
245			\citet{semtner76}. This formulation does not allow storage of heat,
246			that is, the heat capacity of ice is zero. Upward conductive heat flux
247			is parameterized assuming a linear temperature profile and together
248			with a constant ice conductivity. It is expressed as
249			$(K/h)(T_{w}-T_{0})$, where $K$ is the ice conductivity, $h$ the ice
250			thickness, and $T_{w}-T_{0}$ the difference between water and ice
251			surface temperatures. This type of model is often refered to as a
252			``zero-layer'' model. The surface heat flux is computed in a similar
253			way to that of \citet{parkinson79} and \citet{manabe79}.
254
255			The conductive heat flux depends strongly on the ice thickness $h$.
256			However, the ice thickness in the model represents a mean over a
257			potentially very heterogeneous thickness distribution. In order to
258			parameterize a sub-grid scale distribution for heat flux
259			computations, the mean ice thickness $h$ is split into seven thickness
260			categories $H_{n}$ that are equally distributed between $2h$ and a
261			minimum imposed ice thickness of $5\text{\,cm}$ by $H_n=
262			\frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each
263			thickness category is area-averaged to give the total heat flux
264			\citep{hibler84}.
265
266			\ml{[This is Ian Fenty's work and we may want to remove this paragraph
267			from the paper]: %
268			The atmospheric heat flux is balanced by an oceanic heat flux from
269			below. The oceanic flux is proportional to
270			$\rho\,c_{p}\left(T_{w}-T_{fr}\right)$ where $\rho$ and $c_{p}$ are
271			the density and heat capacity of sea water and $T_{fr}$ is the local
272			freezing point temperature that is a function of salinity. Contrary to
273			\citet{menemenlis05}, this flux is not assumed to instantaneously melt
274			or create ice, but a time scale of three days is used to relax $T_{w}$
275			to the freezing point.}
276			%
277			The parameterization of lateral and vertical growth of sea ice follows
278			that of \citet{hibler79, hibler80}.
279
280			On top of the ice there is a layer of snow that modifies the heat flux
281			and the albedo \citep{zhang98}. Snow modifies the effective
282			conductivity according to
283			\[\frac{K}{h} \rightarrow \frac{1}{\frac{h_{s}}{K_{s}}+\frac{h}{K}},\]
284			where $K_s$ is the conductivity of snow and $h_s$ the snow thickness.
285			If enough snow accumulates so that its weight submerges the ice and
286			the snow is flooded, a simple mass conserving parameterization of
287			snowice formation (a flood-freeze algorithm following Archimedes'
288			principle) turns snow into ice until the ice surface is back at $z=0$
289			\citep{leppaeranta83}.
290
291			Effective ice thickness (ice volume per unit area,
292			$c\cdot{h}$), concentration $c$ and effective snow thickness
293			($c\cdot{h}_{s}$) are advected by ice velocities:
294			\begin{equation}
295			\label{eq:advection}
296			\frac{\partial{X}}{\partial{t}} = - \nabla\cdot\left(\vek{u}\,X\right) +
297			\Gamma_{X} + D_{X}
298			\end{equation}
299			where $\Gamma_X$ are the thermodynamic source terms and $D_{X}$ the
300			diffusive terms for quantities $X=(c\cdot{h}), c, (c\cdot{h}_{s})$.
301			%
302			From the various advection scheme that are available in the MITgcm
303			\citep{mitgcm02}, we choose flux-limited schemes
304			\citep[multidimensional 2nd and 3rd-order advection scheme with flux
305			limiter][]{roe85, hundsdorfer94} to preserve sharp gradients and edges
306			that are typical of sea ice distributions and to rule out unphysical
307			over- and undershoots (negative thickness or concentration). These
308			scheme conserve volume and horizontal area and are unconditionally
309			stable, so that we can set $D_{X}=0$. \ml{[do we need to proove that?
310			can we proove that? citation?]}
311
312			There is considerable doubt about the reliability of such a simple
313			thermodynamic model---\citet{semtner84} found significant errors in
314			phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in
315			such models---, so that today many sea ice models employ more complex
316			thermodynamics. A popular thermodynamics model of \citet{winton00} is
317			based on the 3-layer model of \citet{semtner76} and treats brine
318			content by means of enthalphy conservation. This model requires in
319			addition to ice-thickness and compactness (fractional area) additional
320			state variables to be advected by ice velocities, namely enthalphy of
321			the two ice layers and the thickness of the overlying snow layer.
322			\ml{[Jean-Michel, your turf: ]Care must be taken in advecting these
323			quantities in order to ensure conservation of enthalphy. Currently
324			this can only be accomplished with a 2nd-order advection scheme with
325			flux limiter \citep{roe85}.}
326
327			Further documentation and model code can be found at
328	dimitri	1.11	\url{http://mitgcm.org}.
329	dimitri	1.1
330	dimitri	1.12
331	mlosch	1.9	%\subsection{C-grid}
332			%\begin{itemize}
333			%\item no-slip vs. free-slip for both lsr and evp;
334			% "diagnostics" to look at and use for comparison
335			% \begin{itemize}
336			% \item ice transport through Fram Strait/Denmark Strait/Davis
337			% Strait/Bering strait (these are general)
338			% \item ice transport through narrow passages, e.g.\ Nares-Strait
339			% \end{itemize}
340			%\item compare different advection schemes (if lsr turns out to be more
341			% effective, then with lsr otherwise I prefer evp), eg.
342			% \begin{itemize}
343			% \item default 2nd-order with diff1=0.002
344			% \item 1st-order upwind with diff1=0.
345			% \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)
346			% \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)
347			% \end{itemize}
348			% That should be enough. Here, total ice mass and location of ice edge
349			% is interesting. However, this comparison can be done in an idealized
350			% domain, may not require full Arctic Domain?
351			%\item
352			%Do a little study on the parameters of LSR and EVP
353			%\begin{enumerate}
354			%\item convergence of LSR, how many iterations do you need to get a
355			% true elliptic yield curve
356			%\item vary deltaTevp and the relaxation parameter for EVP and see when
357			% the EVP solution breaks down (relative to the forcing time scale).
358			% For this, it is essential that the evp solver gives use "stripeless"
359			% solutions, that is your dtevp = 1sec solutions/or 10sec solutions
360			% with SEAICE\_evpDampC = 615.
361			%\end{enumerate}
362	dimitri	1.2
363	mlosch	1.9	%\end{itemize}
364	mlosch	1.3
365			%%% Local Variables:
366			%%% mode: latex
367			%%% TeX-master: "ceaice"
368			%%% End: