--- MITgcm_contrib/articles/ceaice/ceaice_model.tex	2008/02/28 16:34:56	1.3
+++ MITgcm_contrib/articles/ceaice/ceaice_model.tex	2008/02/28 20:09:23	1.4
@@ -1,6 +1,27 @@
-\section{Model}
+\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 simulations, 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{zha97}, EVP
+  \citep{hunke97};
+\item ice-ocean stress can be formulated as in \citet{hibler87};
+\item ice variables are advected by sophisticated advection schemes;
+\item growth and melt parameterizaion have been refined and extended
+  in order to allow for automatic differentiation of the code.
+\end{itemize}
+The model equations and their numerical realization are summarized
+below.
+
 \subsection{Dynamics}
 \label{sec:dynamics}
 
@@ -43,7 +64,7 @@
 
 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, zhang98}:
+viscous-plastic (VP) constitutive law \citep{hibler79, zhang97}:
 \begin{equation}
   \label{eq:vpequation}
   \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij} 
@@ -205,20 +226,14 @@
 velocity and the ice velocity leading to an inconsistency as the ice
 temperature and salinity are different from the oceanic variables.
 
-Sea ice distributions are characterized by sharp gradients and edges.
-For this reason, we employ positive, multidimensional 2nd and 3rd-order
-advection scheme with flux limiter \citep{roe85, hundsdorfer94} for the
-thermodynamic variables discussed in the next section.
-
-\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}
+%\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{sec:thermodynamics}
@@ -242,8 +257,9 @@
 computations, the mean ice thickness $h$ is split into seven thickness
 categories $H_{n}$ that are equally distributed between $2h$ and
 minimum imposed ice thickness of $5\text{\,cm}$ by $H_n=
-\frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat flux for all thickness
-categories is averaged to give the total heat flux.
+\frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each
+thickness category area averaged to give the total heat flux. \ml{[I
+  don't have citation for that, anyone?]}
 
 The atmospheric heat flux is balanced by an oceanic heat flux from
 below.  The oceanic flux is proportional to
@@ -266,13 +282,27 @@
 
 Effective ich thickness (ice volume per unit area,
 $c\cdot{h}$), concentration $c$ and effective snow thickness
-($c\cdot{h}_{snow}$) are advected by ice velocities as described in
-\refsec{dynamics}. From the various advection scheme that are
-available in the MITgcm \citep{mitgcm02}, we choose flux-limited
-schemes to preserve sharp gradients and edges and to rule out
-unphysical over- and undershoots (negative thickness or
-concentration). These scheme conserve volume and horizontal area.
-\ml{[do we need to proove that? can we proove that? citation?]}
+($c\cdot{h}_{s}$) are advected by ice velocities:
+\begin{align}
+  \frac{\partial(c\,{h})}{\partial{t}} &= - \nabla\left(\vek{u}\,c\,{h}\right) +
+  \Gamma_{h} + D_{h} \\
+  \frac{\partial{c}}{\partial{t}} &= - \nabla\left(\vek{u}\,c\right) +
+  \Gamma_{c} + D_{c} \\
+  \frac{\partial(c\,{h}_{s})}{\partial{t}} &= - \nabla\left(\vek{u}\,c\,{h}_{s}\right) +
+  \Gamma_{h_{s}} + D_{h_{s}} 
+\end{align}
+where $\Gamma_X$ are the thermodynamic source terms and $D_{X}$ the
+diffusive terms for quantity $X=h, c, 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