--- MITgcm_contrib/articles/ceaice/ceaice_model.tex	2008/07/03 18:10:31	1.10
+++ MITgcm_contrib/articles/ceaice/ceaice_model.tex	2008/08/15 14:16:35	1.12
@@ -5,29 +5,328 @@
 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-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 automatic differentiation of the code.
+  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{marshall97:_finit_volum_incom_navier_stokes, mitgcm02}.
-Heat, fresh water fluxes and surface stresses are computed from the
-atmospheric state and modified by the ice model at every time step.
-The model equations and details of their numerical realization are summarized
-in the appendix. Further documentation and model code can be found at
-\url{http://mitgcm.org}. 
+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}