--- MITgcm_contrib/articles/ceaice/ceaice_model.tex	2008/02/28 16:34:56	1.3
+++ MITgcm_contrib/articles/ceaice/ceaice_model.tex	2008/07/03 18:10:31	1.10
@@ -1,328 +1,67 @@
-\section{Model}
+\section{Model Formulation}
 \label{sec:model}
 
-\subsection{Dynamics}
-\label{sec: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}$ 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}$. 
-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, zhang98}:
-\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 timestep
-\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 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}$.
-
-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
-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.
-
-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}
-
+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 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}
+\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
+  \citep{hunke97};
+\item ice-ocean stress can be formulated as in \citet{hibler87};
+\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.
 \end{itemize}
 
-\subsection{Thermodynamics}
-\label{sec:thermodynamics}
+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}. 
+
+%\subsection{C-grid}
+%\begin{itemize}
+%\item no-slip vs. free-slip for both lsr and evp; 
+%  "diagnostics" to look at and use for comparison
+%  \begin{itemize}
+%  \item ice transport through Fram Strait/Denmark Strait/Davis
+%    Strait/Bering strait (these are general) 
+%  \item ice transport through narrow passages, e.g.\ Nares-Strait
+%  \end{itemize}
+%\item compare different advection schemes (if lsr turns out to be more
+%  effective, then with lsr otherwise I prefer evp), eg. 
+%  \begin{itemize}
+%  \item default 2nd-order with diff1=0.002
+%  \item 1st-order upwind with diff1=0.
+%  \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)
+%  \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)
+%  \end{itemize}
+%  That should be enough. Here, total ice mass and location of ice edge
+%  is interesting. However, this comparison can be done in an idealized
+%  domain, may not require full Arctic Domain?
+%\item
+%Do a little study on the parameters of LSR and EVP
+%\begin{enumerate}
+%\item convergence of LSR, how many iterations do you need to get a
+%  true elliptic yield curve 
+%\item vary deltaTevp and the relaxation parameter for EVP and see when
+%  the EVP solution breaks down (relative to the forcing time scale).
+%  For this, it is essential that the evp solver gives use "stripeless"
+%  solutions, that is your dtevp = 1sec solutions/or 10sec solutions
+%  with SEAICE\_evpDampC = 615.
+%\end{enumerate}
 
-In the 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}.
-
-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 flux for all thickness
-categories is averaged to give the total heat flux.
-
-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}_{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?]}
-
-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}.}
-
-
-\subsection{C-grid}
-\begin{itemize}
-\item no-slip vs. free-slip for both lsr and evp; 
-  "diagnostics" to look at and use for comparison
-  \begin{itemize}
-  \item ice transport through Fram Strait/Denmark Strait/Davis
-    Strait/Bering strait (these are general) 
-  \item ice transport through narrow passages, e.g.\ Nares-Strait
-  \end{itemize}
-\item compare different advection schemes (if lsr turns out to be more
-  effective, then with lsr otherwise I prefer evp), eg. 
-  \begin{itemize}
-  \item default 2nd-order with diff1=0.002
-  \item 1st-order upwind with diff1=0.
-  \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me)
-  \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.)
-  \end{itemize}
-  That should be enough. Here, total ice mass and location of ice edge
-  is interesting. However, this comparison can be done in an idealized
-  domain, may not require full Arctic Domain?
-\item
-Do a little study on the parameters of LSR and EVP
-\begin{enumerate}
-\item convergence of LSR, how many iterations do you need to get a
-  true elliptic yield curve 
-\item vary deltaTevp and the relaxation parameter for EVP and see when
-  the EVP solution breaks down (relative to the forcing time scale).
-  For this, it is essential that the evp solver gives use "stripeless"
-  solutions, that is your dtevp = 1sec solutions/or 10sec solutions
-  with SEAICE\_evpDampC = 615.
-\end{enumerate}
-
-\end{itemize}
+%\end{itemize}
 
 %%% Local Variables: 
 %%% mode: latex