--- MITgcm_contrib/articles/ceaice/ceaice_model.tex	2008/06/04 13:33:45	1.9
+++ MITgcm_contrib/articles/ceaice/ceaice_model.tex	2008/08/15 14:16:35	1.12
@@ -5,25 +5,29 @@
 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
+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}, EVP
+  momentum equations have been adopted: LSOR \citep{zhang97}, and 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 parameterization have been refined and extended
-  in order to allow for automatic differentiation of the code.
+\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 model equations and their numerical realization are summarized
-below.
+
+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{sec:dynamics}
+\label{app:dynamics}
 
 The momentum equation of the sea-ice model is
 \begin{equation}
@@ -40,14 +44,12 @@
 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}$ is the sea surface height potential
-in response to ocean dynamics ($g\eta$) and to atmospheric pressure
-loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density);
-and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress
-tensor $\sigma_{ij}$.
-When using the rescaled vertical coordinate system, z$^\ast$, of
-\citet{cam08}, $\phi(0)$ also includes a term due to snow and ice
-loading, $mg/\rho_{0}$. 
+$\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*}
@@ -98,7 +100,7 @@
     \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}}
+  \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
@@ -175,10 +177,9 @@
 the external (long) timestep $\Delta{t}$. \citet{hunke97} recommend
 $E_{0} = \frac{1}{3}$.
 
-For details of the spatial discretization, the reader is referred to
-\citet{zhang98, zhang03}. Our discretization differs only (but
-importantly) in the underlying grid, namely the Arakawa C-grid, but is
-otherwise straightforward. The EVP model, in particular, is discretized
+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
@@ -237,31 +238,33 @@
 %\end{itemize}
 
 \subsection{Thermodynamics}
-\label{sec:thermodynamics}
+\label{app:thermodynamics}
 
-In the original formulation the sea ice model \citep{menemenlis05}
+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
-(heat capacity of ice is zero, and this type of model is often refered
-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 flux is computed in a similar way to that of \citet{parkinson79}
-and \citet{manabe79}.
+\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 this sub-grid scale distribution for heat flux
+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=
+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 area averaged to give the total heat flux. \ml{[I
-  don't have citation for that, anyone?]}
+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
@@ -269,17 +272,21 @@
 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.
-
+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}.
+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
@@ -317,6 +324,9 @@
   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}