--- MITgcm_contrib/articles/ceaice/ceaice_model.tex	2008/02/26 19:27:26	1.1
+++ MITgcm_contrib/articles/ceaice/ceaice_model.tex	2008/02/28 20:45:29	1.5
@@ -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}
@@ -227,13 +242,64 @@
 uses simple thermodynamics following the appendix of
 \citet{semtner76}. This formulation does not allow storage of heat
 (heat capacity of ice is zero, and this type of model is often refered
-to as a ``zero-layer'' model). Upward heat flux is parameterized
+to as a ``zero-layer'' model). 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. The surface
-heat budget is computed in a similar way to that of
-\citet{parkinson79} and \citet{manabe79}.
+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 this 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
+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 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
+$\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}. 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 ich 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
@@ -282,4 +348,10 @@
   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: 
