s_phys_pkgs/text/seaice.tex

% $Header: /u/gcmpack/manual/part6/seaice.tex,v 1.7 2008/01/17 22:32:38 heimbach Exp $
% $Name:  $

%%EH3  Copied from "MITgcm/pkg/seaice/seaice_description.tex"
%%EH3  which was written by Dimitris M.


\subsection{SEAICE Package}
\label{sec:pkg:seaice}
\begin{rawhtml}
<!-- CMIREDIR:package_seaice: -->
\end{rawhtml}

Authors: Martin Losch, Dimitris Menemenlis, An Nguyen, Jean-Michel Campin,
Patrick Heimbach, Chris Hill and Jinlun Zhang

%----------------------------------------------------------------------
\subsubsection{Introduction
\label{sec:pkg:exf:intro}}


Package ``seaice'' provides a dynamic and thermodynamic interactive
sea-ice model. 

CPP options enable or disable different aspects of the package
(Section \ref{sec:pkg:seaice:config}).
Run-Time options, flags, filenames and field-related dates/times are
set in \texttt{data.seaice}
(Section \ref{sec:pkg:seaice:runtime}).
A description of key subroutines is given in Section
\ref{sec:pkg:seaice:subroutines}.
Input fields, units and sign conventions are summarized in
Section \ref{sec:pkg:seaice:fields_units}, and available diagnostics
output is listed in Section \ref{sec:pkg:seaice:fields_diagnostics}.

%----------------------------------------------------------------------

\subsubsection{SEAICE configuration, compiling \& running}

\paragraph{Compile-time options
\label{sec:pkg:seaice:config}}
~

As with all MITgcm packages, SEAICE can be turned on or off at compile time
%
\begin{itemize}
%
\item
using the \texttt{packages.conf} file by adding \texttt{seaice} to it,
%
\item
or using \texttt{genmake2} adding
\texttt{-enable=seaice} or \texttt{-disable=seaice} switches
%
\item
\textit{required packages and CPP options}: \\
SEAICE requires the external forcing package \texttt{exf} to be enabled;
no additional CPP options are required.
%
\end{itemize}
(see Section \ref{sect:buildingCode}).

Parts of the SEAICE code can be enabled or disabled at compile time
via CPP preprocessor flags. These options are set in either
\texttt{SEAICE\_OPTIONS.h} or in \texttt{ECCO\_CPPOPTIONS.h}.
Table \ref{tab:pkg:seaice:cpp} summarizes these options.

\begin{table}[h!]
\centering
  \label{tab:pkg:seaice:cpp}
  {\footnotesize
    \begin{tabular}{|l|p{10cm}|}
      \hline 
      \textbf{CPP option}  &  \textbf{Description}  \\
      \hline \hline
        \texttt{SEAICE\_DEBUG} & 
          Enhance STDOUT for debugging \\
        \texttt{SEAICE\_ALLOW\_DYNAMICS} & 
          sea-ice dynamics code \\
        \texttt{SEAICE\_CGRID} & 
          LSR solver on C-grid (rather than original B-grid) \\
        \texttt{SEAICE\_ALLOW\_EVP} & 
          use EVP rather than LSR rheology solver \\
        \texttt{SEAICE\_EXTERNAL\_FLUXES} & 
          use EXF-computed fluxes as starting point \\
        \texttt{SEAICE\_MULTICATEGORY} & 
          enable 8-category thermodynamics (by default undefined)\\
        \texttt{SEAICE\_VARIABLE\_FREEZING\_POINT} & 
          enable linear dependence of the freezing point on salinity
          (by default undefined)\\
        \texttt{ALLOW\_SEAICE\_FLOODING} & 
          enable snow to ice conversion for submerged sea-ice \\
        \texttt{SEAICE\_SALINITY} & 
          enable "salty" sea-ice (by default undefined) \\
        \texttt{SEAICE\_AGE} & 
          enable "age tracer" sea-ice (by default undefined) \\
        \texttt{SEAICE\_CAP\_HEFF} & 
          enable capping of sea-ice thickness to MAX\_HEFF \\ \hline
        \texttt{SEAICE\_BICE\_STRESS} &
          B-grid only for backward compatiblity: turn on ice-stress on
          ocean\\
        \texttt{EXPLICIT\_SSH\_SLOPE} &
          B-grid only for backward compatiblity: use ETAN for tilt
          computations rather than geostrophic velocities \\
      \hline
    \end{tabular}
  }
  \caption{~}
\end{table}

%----------------------------------------------------------------------

\subsubsection{Run-time parameters
\label{sec:pkg:seaice:runtime}}

Run-time parameters are set in files 
\texttt{data.pkg} (read in \texttt{packages\_readparms.F}),
and \texttt{data.seaice} (read in \texttt{seaice\_readparms.F}).

\paragraph{Enabling the package}
~ \\
%
A package is switched on/off at run-time by setting
(e.g. for SEAICE) \texttt{useSEAICE = .TRUE.} in \texttt{data.pkg}.

\paragraph{General flags and parameters}
~ \\
%
Table~\ref{tab:pkg:seaice:runtimeparms} lists most run-time parameters.
\input{part6/seaice-parms.tex}


%----------------------------------------------------------------------
\subsubsection{Description
\label{sec:pkg:seaice:descr}}

[TO BE CONTINUED/MODIFIED]

% Sea-ice model thermodynamics are based on Hibler
% \cite{hib80}, that is, a 2-category model that simulates ice thickness
% and concentration.  Snow is simulated as per Zhang et al.
% \cite{zha98a}.  Although recent years have seen an increased use of
% multi-category thickness distribution sea-ice models for climate
% studies, the Hibler 2-category ice model is still the most widely used
% model and has resulted in realistic simulation of sea-ice variability
% on regional and global scales.  Being less complicated, compared to
% multi-category models, the 2-category model permits easier application
% of adjoint model optimization methods.

% Note, however, that the Hibler 2-category model and its variants use a
% so-called zero-layer thermodynamic model to estimate ice growth and
% decay.  The zero-layer thermodynamic model assumes that ice does not
% store heat and, therefore, tends to exaggerate the seasonal
% variability in ice thickness.  This exaggeration can be significantly
% reduced by using Semtner's \cite{sem76} three-layer thermodynamic
% model that permits heat storage in ice.  Recently, the three-layer
% thermodynamic model has been reformulated by Winton \cite{win00}.  The
% reformulation improves model physics by representing the brine content
% of the upper ice with a variable heat capacity.  It also improves
% model numerics and consumes less computer time and memory.  The Winton
% sea-ice thermodynamics have been ported to the MIT GCM; they currently
% reside under pkg/thsice. The package pkg/thsice is fully
% compatible with pkg/seaice and with pkg/exf. When turned on togeter
% with pkg/seaice, the zero-layer thermodynamics are replaced by the by
% Winton thermodynamics

The MITgcm sea ice model (MITgcm/sim) is based on a variant of the
viscous-plastic (VP) dynamic-thermodynamic sea ice model \citep{zhang97}
first introduced by \citet{hib79, hib80}. 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{hun97};
\item ice-ocean stress can be formulated as in \citet{hibler87} or as in
  \citet{cam08}; 
\item ice variables are advected by sophisticated, conservative
  advection schemes 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.  Heat, fresh water fluxes and surface stresses are computed
from the atmospheric state and -- by default -- modified by the ice
model at every time step.

The ice dynamics models that are most widely used for large-scale
climate studies are the viscous-plastic (VP) model \citep{hib79}, the
cavitating fluid (CF) model \citep{fla92}, and the
elastic-viscous-plastic (EVP) model \citep{hun97}.  Compared to the VP
model, the CF model does not allow ice shear in calculating ice
motion, stress, and deformation.  EVP models approximate VP by adding
an elastic term to the equations for easier adaptation to parallel
computers.  Because of its higher accuracy in plastic solution and
relatively simpler formulation, compared to the EVP model, we decided
to use the VP model as the default dynamic component of our ice
model. To do this we extended the line successive over relaxation
(LSOR) method of \citet{zhang97} for use in a parallel
configuration.

Note, that by default the seaice-package includes the orginial
so-called zero-layer thermodynamics following \citet{hib80} with a
snow cover as in \citet{zha98a}. The zero-layer thermodynamic model
assumes that ice does not store heat and, therefore, tends to
exaggerate the seasonal variability in ice thickness.  This
exaggeration can be significantly reduced by using
\citeauthor{sem76}'s~[\citeyear{sem76}] three-layer thermodynamic model
that permits heat storage in ice.  Recently, the three-layer
thermodynamic model has been reformulated by \citet{win00}.  The
reformulation improves model physics by representing the brine content
of the upper ice with a variable heat capacity.  It also improves
model numerics and consumes less computer time and memory.  The Winton
sea-ice thermodynamics have been ported to the MIT GCM; they currently
reside under pkg/thsice. The package pkg/thsice is fully compatible
with pkg/seaice and with pkg/exf. When turned on together with
pkg/seaice, the zero-layer thermodynamics are replaced by the Winton
thermodynamics.

The sea ice model requires the following input fields: 10-m winds, 2-m
air temperature and specific humidity, downward longwave and shortwave
radiations, precipitation, evaporation, and river and glacier runoff.
The sea ice model also requires surface temperature from the ocean
model and the top level horizontal velocity.  Output fields are
surface wind stress, evaporation minus precipitation minus runoff, net
surface heat flux, and net shortwave flux.  The sea-ice model is
global: in ice-free regions bulk formulae are used to estimate oceanic
forcing from the atmospheric fields.

\subsubsection{Dynamics}
\label{sec:pkg:seaice:dynamics}

\newcommand{\vek}[1]{\ensuremath{\vec{\mathbf{#1}}}}
\newcommand{\vtau}{\vek{\mathbf{\tau}}}
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{cam08};
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{hib79, 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^{*}$ (run-time parameter \texttt{SEAICE\_strength}) 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*}
The bulk viscosities are bounded above by imposing both a minimum
$\Delta_{\min}$ (for numerical reasons, run-time parameter
\texttt{SEAICE\_EPS} with a default value of
$10^{-10}\text{\,s}^{-1}$) and a maximum $\zeta_{\max} =
P_{\max}/\Delta^*$, where
$\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. (There is also
the option of bounding $\zeta$ from below by setting run-time
parameter \texttt{SEAICE\_zetaMin} $>0$, but this is generally not
recommended). 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}
To enable this method, set \texttt{\#define SEAICE\_ALLOW\_TEM} in
\texttt{SEAICE\_OPTIONS.h} and turn it on with
\texttt{SEAICEuseTEM=.TRUE.} in \texttt{data.seaice}. 

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{hun97}'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}
%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{hun97}, 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{hun97}. 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}$. \citet{hun97} recommend
$E_{0} = \frac{1}{3}$ (which is the default value in the code).

To use the EVP solver, make sure that both \texttt{SEAICE\_CGRID} and
\texttt{SEAICE\_ALLOW\_EVP} are defined in \texttt{SEAICE\_OPTIONS.h}
(default). The solver is turned on by setting the sub-cycling time
step \texttt{SEAICE\_deltaTevp} to a value larger than zero. The
choice of this time step is under debate. \citet{hun97} recommend
order(120) time steps for the EVP solver within one model time step
$\Delta{t}$ (\texttt{deltaTmom}). One can also choose order(120) time
steps within the forcing time scale, but then we recommend adjusting
the damping time scale $T$ accordingly, by setting either
\texttt{SEAICE\_elasticParm} ($E_{0}$), so that
$E_{0}\Delta{t}=\mbox{forcing time scale}$, or directly
\texttt{SEAICE\_evpTauRelax} ($T$) to the forcing time scale.

Moving sea ice exerts a stress on the ocean which is the opposite of
the stress $\vtau_{ocean}$ in Eq.~\ref{eq: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.
To turn on the stress formulation of \citet{hibler87}, set
\texttt{useHB87StressCoupling=.TRUE.} in \texttt{data.seaice}.


% Our discretization differs from \citet{zhang97, 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.

\subsubsection{Finite-volume discretization of the stress tensor
  divergence}
\label{sec:pkg:seaice: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 a C-points in the center of the grid
cell. Discretization requires only averaging of $\zeta$ and $\eta$ to
vorticity or Z-points (or $\zeta$-points, but here we use Z in order
avoid confusion with the bulk viscosity) 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 except that $\zeta$-points
have been renamed to Z-points), the strain rates are discretized as:
\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}
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, $k_{1}$ and
$k_{2}$ can be approximated by finite differences of the cell widths:
\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}

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{hib79}. 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{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}

Similarly, we have for the $v$-equation ($\alpha=2$):
\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}

Again, no slip boundary conditions are realized via ghost points and
$u_{i,j-1}+u_{i,j}=0$ and $v_{i-1,j}+v_{i,j}=0$ across boundaries. For
free slip boundary conditions the lateral stress is set to zeros. In
analogy to $(\epsilon_{12})^Z=0$ on boundaries, we set
$\sigma_{21}^{Z}=0$, or equivalently $\eta_{i,j}^{Z}=0$, on boundaries.

\subsubsection{Thermodynamics}
\label{sec:pkg:seaice:thermodynamics}

In its original formulation the sea ice model \citep{menemenlis05}
uses simple thermodynamics following the appendix of
\citet{sem76}. 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}. To use this thickness category parameterization set
\texttt{\#define SEAICE\_MULTICATEGORY}; note that this requires
different restart files and switching this flag on in the middle of an
integration is not possible.

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. This flux
is not assumed to instantaneously melt or create ice, but a time scale
of three days (run-time parameter \texttt{SEAICE\_gamma\_t}) is used
to relax $T_{w}$ to the freezing point.
%
The parameterization of lateral and vertical growth of sea ice follows
that of \citet{hib79, hib80}; the so-called lead closing parameter
$h_{0}$ (run-time parameter \texttt{HO}) has a default value of
0.5~meters.

On top of the ice there is a layer of snow that modifies the heat flux
and the albedo \citep{zha98a}. 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}. The flood-freeze algorithm is enabled with the CPP-flag
\texttt{SEAICE\_ALLOW\_FLOODING} and turned on with run-time parameter
\texttt{SEAICEuseFlooding=.true.}.

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, we
choose flux-limited schemes \citep[multidimensional 2nd and 3rd-order
advection scheme with flux limiter][]{roe:85, 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$. Run-timeflags: \texttt{SEAICEadvScheme} (default=2),
\texttt{DIFF1} (default=0.004).

There is considerable doubt about the reliability of a ``zero-layer''
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. The MITgcm sea ice model provides the option to use
the thermodynamics model of \citet{win00}, which in turn is based
on the 3-layer model of \citet{sem76} and which treats brine
content by means of enthalpy conservation. This scheme requires
additional state variables, namely the enthalpy of the two ice layers
(instead of effective ice salinity), to be advected by ice velocities.
%
The internal sea ice temperature is inferred from ice enthalpy.  To
avoid unphysical (negative) values for ice thickness and
concentration, a positive 2nd-order advection scheme with a SuperBee
flux limiter \citep{roe:85} is used in this study to advect all
sea-ice-related quantities of the \citet{win00} thermodynamic
model.  Because of the non-linearity of the advection scheme, care
must be taken in advecting these quantities: when simply using ice
velocity to advect enthalpy, the total energy (i.e., the volume
integral of enthalpy) is not conserved. Alternatively, one can advect
the energy content (i.e., product of ice-volume and enthalpy) but then
false enthalpy extrema can occur, which then leads to unrealistic ice
temperature.  In the currently implemented solution, the sea-ice mass
flux is used to advect the enthalpy in order to ensure conservation of
enthalpy and to prevent false enthalpy extrema.

%----------------------------------------------------------------------

\subsubsection{Key subroutines
\label{sec:pkg:seaice:subroutines}}

Top-level routine: \texttt{seaice\_model.F}

{\footnotesize
\begin{verbatim}

C     !CALLING SEQUENCE:
c ...
c  seaice_model (TOP LEVEL ROUTINE)
c  |
c  |-- #ifdef SEAICE_CGRID
c  |     SEAICE_DYNSOLVER
c  |     |
c  |     |-- < compute proxy for geostrophic velocity >
c  |     |
c  |     |-- < set up mass per unit area and Coriolis terms >
c  |     |
c  |     |-- < dynamic masking of areas with no ice >
c  |     |
c  |     |

c  |   #ELSE
c  |     DYNSOLVER
c  |   #ENDIF
c  |
c  |-- if ( useOBCS ) 
c  |     OBCS_APPLY_UVICE
c  |
c  |-- if ( SEAICEadvHeff .OR. SEAICEadvArea .OR. SEAICEadvSnow .OR. SEAICEadvSalt )
c  |     SEAICE_ADVDIFF
c  |
c  |-- if ( usePW79thermodynamics ) 
c  |     SEAICE_GROWTH
c  |
c  |-- if ( useOBCS ) 
c  |     if ( SEAICEadvHeff ) OBCS_APPLY_HEFF
c  |     if ( SEAICEadvArea ) OBCS_APPLY_AREA
c  |     if ( SEAICEadvSALT ) OBCS_APPLY_HSALT
c  |     if ( SEAICEadvSNOW ) OBCS_APPLY_HSNOW
c  |
c  |-- < do various exchanges >
c  |
c  |-- < do additional diagnostics >
c  |
c  o

\end{verbatim}
}


%----------------------------------------------------------------------

\subsubsection{SEAICE diagnostics
\label{sec:pkg:seaice:diagnostics}}

Diagnostics output is available via the diagnostics package
(see Section \ref{sec:pkg:diagnostics}).
Available output fields are summarized in 
Table \ref{tab:pkg:seaice:diagnostics}.

\begin{table}[h!]
\centering
\label{tab:pkg:seaice:diagnostics}
{\footnotesize
\begin{verbatim}
---------+----+----+----------------+-----------------
 <-Name->|Levs|grid|<--  Units   -->|<- Tile (max=80c)
---------+----+----+----------------+-----------------
 SIarea  |  1 |SM  |m^2/m^2         |SEAICE fractional ice-covered area [0 to 1]
 SIheff  |  1 |SM  |m               |SEAICE effective ice thickness
 SIuice  |  1 |UU  |m/s             |SEAICE zonal ice velocity, >0 from West to East
 SIvice  |  1 |VV  |m/s             |SEAICE merid. ice velocity, >0 from South to North
 SIhsnow |  1 |SM  |m               |SEAICE snow thickness
 SIhsalt |  1 |SM  |g/m^2           |SEAICE effective salinity
 SIatmFW |  1 |SM  |kg/m^2/s        |Net freshwater flux from the atmosphere (+=down)
 SIuwind |  1 |SM  |m/s             |SEAICE zonal 10-m wind speed, >0 increases uVel
 SIvwind |  1 |SM  |m/s             |SEAICE meridional 10-m wind speed, >0 increases uVel
 SIfu    |  1 |UU  |N/m^2           |SEAICE zonal surface wind stress, >0 increases uVel
 SIfv    |  1 |VV  |N/m^2           |SEAICE merid. surface wind stress, >0 increases vVel
 SIempmr |  1 |SM  |kg/m^2/s        |SEAICE upward freshwater flux, > 0 increases salt
 SIqnet  |  1 |SM  |W/m^2           |SEAICE upward heatflux, turb+rad, >0 decreases theta
 SIqsw   |  1 |SM  |W/m^2           |SEAICE upward shortwave radiat., >0 decreases theta
 SIpress |  1 |SM  |m^2/s^2         |SEAICE strength (with upper and lower limit)
 SIzeta  |  1 |SM  |m^2/s           |SEAICE nonlinear bulk viscosity
 SIeta   |  1 |SM  |m^2/s           |SEAICE nonlinear shear viscosity
 SIsigI  |  1 |SM  |no units        |SEAICE normalized principle stress, component one
 SIsigII |  1 |SM  |no units        |SEAICE normalized principle stress, component two
 SIthdgrh|  1 |SM  |m/s             |SEAICE thermodynamic growth rate of effective ice thickness
 SIsnwice|  1 |SM  |m/s             |SEAICE ice formation rate due to flooding
 SIuheff |  1 |UU  |m^2/s           |Zonal Transport of effective ice thickness
 SIvheff |  1 |VV  |m^2/s           |Meridional Transport of effective ice thickness
 ADVxHEFF|  1 |UU  |m.m^2/s         |Zonal      Advective Flux of eff ice thickn
 ADVyHEFF|  1 |VV  |m.m^2/s         |Meridional Advective Flux of eff ice thickn
 DFxEHEFF|  1 |UU  |m.m^2/s         |Zonal      Diffusive Flux of eff ice thickn
 DFyEHEFF|  1 |VV  |m.m^2/s         |Meridional Diffusive Flux of eff ice thickn
 ADVxAREA|  1 |UU  |m^2/m^2.m^2/s   |Zonal      Advective Flux of fract area
 ADVyAREA|  1 |VV  |m^2/m^2.m^2/s   |Meridional Advective Flux of fract area
 DFxEAREA|  1 |UU  |m^2/m^2.m^2/s   |Zonal      Diffusive Flux of fract area
 DFyEAREA|  1 |VV  |m^2/m^2.m^2/s   |Meridional Diffusive Flux of fract area
 ADVxSNOW|  1 |UU  |m.m^2/s         |Zonal      Advective Flux of eff snow thickn
 ADVySNOW|  1 |VV  |m.m^2/s         |Meridional Advective Flux of eff snow thickn
 DFxESNOW|  1 |UU  |m.m^2/s         |Zonal      Diffusive Flux of eff snow thickn
 DFyESNOW|  1 |VV  |m.m^2/s         |Meridional Diffusive Flux of eff snow thickn
 ADVxSSLT|  1 |UU  |psu.m^2/s       |Zonal      Advective Flux of seaice salinity
 ADVySSLT|  1 |VV  |psu.m^2/s       |Meridional Advective Flux of seaice salinity
 DFxESSLT|  1 |UU  |psu.m^2/s       |Zonal      Diffusive Flux of seaice salinity
 DFyESSLT|  1 |VV  |psu.m^2/s       |Meridional Diffusive Flux of seaice salinity
\end{verbatim}
}
\caption{Available diagnostics of the seaice-package}
\end{table}


%\subsubsection{Package Reference}

\subsubsection{Experiments and tutorials that use seaice}
\label{sec:pkg:seaice:experiments}

\begin{itemize}
\item{Labrador Sea experiment in lab\_sea verification directory. }
\end{itemize}


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "../manual"
%%% End: 