/[MITgcm]/manual/s_phys_pkgs/text/seaice.tex

Diff of /manual/s_phys_pkgs/text/seaice.tex

Parent Directory | Revision Log | View Revision Graph Revision Graph | View Patch Patch

-revision 1.6 by heimbach,
Tue Jan 15 23:58:53 2008 UTC
+revision 1.22 by mlosch,
Wed Jan 21 18:18:22 2015 UTC
 Line 16 
 Patrick Heimbach, Chris Hill and Jinlun
  %----------------------------------------------------------------------
  \subsubsection{Introduction
- \label{sec:pkg:exf:intro}}
+ \label{sec:pkg:seaice:intro}}
  Package ``seaice'' provides a dynamic and thermodynamic interactive
 Line 24 
 sea-ice model.
  CPP options enable or disable different aspects of the package
  (Section \ref{sec:pkg:seaice:config}).
- Runtime options, flags, filenames and field-related dates/times are
+ Run-Time options, flags, filenames and field-related dates/times are
- set in \texttt{data.seaice}
+ set in \code{data.seaice}
  (Section \ref{sec:pkg:seaice:runtime}).
  A description of key subroutines is given in Section
  \ref{sec:pkg:seaice:subroutines}.
  Input fields, units and sign conventions are summarized in
  Section \ref{sec:pkg:seaice:fields_units}, and available diagnostics
- output is listed in Section \ref{sec:pkg:seaice:fields_diagnostics}.
+ output is listed in Section \ref{sec:pkg:seaice:diagnostics}.
  %----------------------------------------------------------------------
 Line 46 
 As with all MITgcm packages, SEAICE can
  \begin{itemize}
  %
  \item
- using the \texttt{packages.conf} file by adding \texttt{seaice} to it,
+ using the \code{packages.conf} file by adding \code{seaice} to it,
  %
  \item
- or using \texttt{genmake2} adding
+ or using \code{genmake2} adding
- \texttt{-enable=seaice} or \texttt{-disable=seaice} switches
+ \code{-enable=seaice} or \code{-disable=seaice} switches
  %
  \item
  \textit{required packages and CPP options}: \\
- SEAICE requires the external forcing package \texttt{exf} to be enabled;
+ SEAICE requires the external forcing package \code{exf} to be enabled;
  no additional CPP options are required.
  %
  \end{itemize}
- (see Section \ref{sect:buildingCode}).
+ (see Section \ref{sec:buildingCode}).
  Parts of the SEAICE code can be enabled or disabled at compile time
- via CPP preprocessor flags. These options are set in either
+ via CPP preprocessor flags. These options are set in
- \texttt{SEAICE\_OPTIONS.h} or in \texttt{ECCO\_CPPOPTIONS.h}.
+ \code{SEAICE\_OPTIONS.h}.
- Table \ref{tab:pkg:seaice:cpp} summarizes these options.
+ Table \ref{tab:pkg:seaice:cpp} summarizes the most important ones.
- \begin{table}[h!]
+ \begin{table}[!ht]
  \centering
    \label{tab:pkg:seaice:cpp}
    {\footnotesize
-     \begin{tabular}{|l|l|}
+     \begin{tabular}{|l|p{10cm}|}
        \hline
        \textbf{CPP option}  &  \textbf{Description}  \\
        \hline \hline
-         \texttt{SEAICE\_DEBUG} &
+         \code{SEAICE\_DEBUG} &
            Enhance STDOUT for debugging \\
-         \texttt{SEAICE\_ALLOW\_DYNAMICS} &
+         \code{SEAICE\_ALLOW\_DYNAMICS} &
            sea-ice dynamics code \\
-         \texttt{SEAICE\_CGRID} &
+         \code{SEAICE\_CGRID} &
-           LSR solver on C-grid (rather than original B-grid \\
+           LSR solver on C-grid (rather than original B-grid) \\
-         \texttt{SEAICE\_ALLOW\_EVP} &
+         \code{SEAICE\_ALLOW\_EVP} &
-           use EVP rather than LSR rheology solver \\
+           enable use of EVP rheology solver \\
-         \texttt{SEAICE\_EXTERNAL\_FLUXES} &
+         \code{SEAICE\_ALLOW\_JFNK} &
+           enable use of JFNK rheology solver \\
+         \code{SEAICE\_EXTERNAL\_FLUXES} &
            use EXF-computed fluxes as starting point \\
-         \texttt{SEAICE\_MULTICATEGORY} &
+         \code{SEAICE\_ZETA\_SMOOTHREG} &
-           enable 8-category thermodynamics \\
+           use differentialable regularization for viscosities \\
-         \texttt{SEAICE\_VARIABLE\_FREEZING\_POINT} &
+         \code{SEAICE\_VARIABLE\_FREEZING\_POINT} &
-           enable linear dependence of the freezing point on salinity \\
+           enable linear dependence of the freezing point on salinity
-         \texttt{ALLOW\_SEAICE\_FLOODING} &
+           (by default undefined)\\
+         \code{ALLOW\_SEAICE\_FLOODING} &
            enable snow to ice conversion for submerged sea-ice \\
-         \texttt{SEAICE\_SALINITY} &
+         \code{SEAICE\_VARIABLE\_SALINITY} &
-           enable "salty" sea-ice \\
+           enable sea-ice with variable salinity (by default undefined) \\
-         \texttt{SEAICE\_CAP\_HEFF} &
+         \code{SEAICE\_SITRACER} &
-           enable capping of sea-ice thickness to MAX\_HEFF \\
+           enable sea-ice tracer package (by default undefined) \\
+         \code{SEAICE\_BICE\_STRESS} &
+           B-grid only for backward compatiblity: turn on ice-stress on
+           ocean\\
+         \code{EXPLICIT\_SSH\_SLOPE} &
+           B-grid only for backward compatiblity: use ETAN for tilt
+           computations rather than geostrophic velocities \\
        \hline
      \end{tabular}
    }
-   \caption{~}
+   \caption{Some of the most relevant CPP preprocessor flags in the
+     \code{seaice}-package.}
  \end{table}
  %----------------------------------------------------------------------
-Line 104 
 Table \ref{tab:pkg:seaice:cpp} summarize
+Line 114 
 Table \ref{tab:pkg:seaice:cpp} summarize
  \subsubsection{Run-time parameters
  \label{sec:pkg:seaice:runtime}}
- Run-time parameters are set in files
+ Run-time parameters (see Table~\ref{tab:pkg:seaice:runtimeparms}) are set
- \texttt{data.pkg} (read in \texttt{packages\_readparms.F}),
+ in files \code{data.pkg} (read in \code{packages\_readparms.F}), and
- and \texttt{data.seaice} (read in \texttt{seaice\_readparms.F}).
+ \code{data.seaice} (read in \code{seaice\_readparms.F}).
  \paragraph{Enabling the package}
  ~ \\
  %
- A package is switched on/off at runtime by setting
+ A package is switched on/off at run-time by setting
- (e.g. for SEAICE) \texttt{useSEAICE = .TRUE.} in \texttt{data.pkg}.
+ (e.g. for SEAICE) \code{useSEAICE = .TRUE.} in \code{data.pkg}.
  \paragraph{General flags and parameters}
  ~ \\
  %
- \input{part6/seaice-parms.tex}
+ Table~\ref{tab:pkg:seaice:runtimeparms} lists most run-time parameters.
+ \input{s_phys_pkgs/text/seaice-parms.tex}
+ \paragraph{Input fields and units\label{sec:pkg:seaice:fields_units}}
+ \begin{description}
+ \item[\code{HeffFile}:] Initial sea ice thickness averaged over grid cell
+   in meters; initializes variable \code{HEFF};
+ \item[\code{AreaFile}:] Initial fractional sea ice cover, range $[0,1]$;
+   initializes variable \code{AREA};
+ \item[\code{HsnowFile}:] Initial snow thickness on sea ice averaged
+   over grid cell in meters; initializes variable \code{HSNOW};
+ \item[\code{HsaltFile}:] Initial salinity of sea ice averaged over grid
+   cell in g/m$^2$; initializes variable \code{HSALT};
+ \end{description}
  %----------------------------------------------------------------------
  \subsubsection{Description
-Line 127 
 A package is switched on/off at runtime
+Line 148 
 A package is switched on/off at runtime
  [TO BE CONTINUED/MODIFIED]
- Sea-ice model thermodynamics are based on Hibler
+ % Sea-ice model thermodynamics are based on Hibler
- \cite{hib80}, that is, a 2-category model that simulates ice thickness
+ % \cite{hib80}, that is, a 2-category model that simulates ice thickness
- and concentration.  Snow is simulated as per Zhang et al.
+ % and concentration.  Snow is simulated as per Zhang et al.
- \cite{zha98a}.  Although recent years have seen an increased use of
+ % \cite{zha98a}.  Although recent years have seen an increased use of
- multi-category thickness distribution sea-ice models for climate
+ % multi-category thickness distribution sea-ice models for climate
- studies, the Hibler 2-category ice model is still the most widely used
+ % studies, the Hibler 2-category ice model is still the most widely used
- model and has resulted in realistic simulation of sea-ice variability
+ % model and has resulted in realistic simulation of sea-ice variability
- on regional and global scales.  Being less complicated, compared to
+ % on regional and global scales.  Being less complicated, compared to
- multi-category models, the 2-category model permits easier application
+ % multi-category models, the 2-category model permits easier application
- of adjoint model optimization methods.
+ % of adjoint model optimization methods.
- Note, however, that the Hibler 2-category model and its variants use a
+ % Note, however, that the Hibler 2-category model and its variants use a
- so-called zero-layer thermodynamic model to estimate ice growth and
+ % so-called zero-layer thermodynamic model to estimate ice growth and
- decay.  The zero-layer thermodynamic model assumes that ice does not
+ % decay.  The zero-layer thermodynamic model assumes that ice does not
- store heat and, therefore, tends to exaggerate the seasonal
+ % store heat and, therefore, tends to exaggerate the seasonal
- variability in ice thickness.  This exaggeration can be significantly
+ % variability in ice thickness.  This exaggeration can be significantly
- reduced by using Semtner's \cite{sem76} three-layer thermodynamic
+ % reduced by using Semtner's \cite{sem76} three-layer thermodynamic
- model that permits heat storage in ice.  Recently, the three-layer
+ % model that permits heat storage in ice.  Recently, the three-layer
- thermodynamic model has been reformulated by Winton \cite{win00}.  The
+ % thermodynamic model has been reformulated by Winton \cite{win00}.  The
- reformulation improves model physics by representing the brine content
+ % reformulation improves model physics by representing the brine content
- of the upper ice with a variable heat capacity.  It also improves
+ % of the upper ice with a variable heat capacity.  It also improves
- model numerics and consumes less computer time and memory.  The Winton
+ % model numerics and consumes less computer time and memory.  The Winton
- sea-ice thermodynamics have been ported to the MIT GCM; they currently
+ % sea-ice thermodynamics have been ported to the MIT GCM; they currently
- reside under pkg/thsice.  At present pkg/thsice is not fully
+ % reside under pkg/thsice. The package pkg/thsice is fully
- compatible with pkg/seaice and with pkg/exf.  But the eventual
+ % compatible with pkg/seaice and with pkg/exf. When turned on togeter
- objective is to have fully compatible and interchangeable
+ % with pkg/seaice, the zero-layer thermodynamics are replaced by the by
- thermodynamic packages for sea-ice, so that it becomes possible to use
+ % Winton thermodynamics
- Hibler dynamics with Winton thermodyanmics.
+ The MITgcm sea ice model (MITgcm/sim) is based on a variant of the
+ viscous-plastic (VP) dynamic-thermodynamic sea ice model \citep{zhang97}
+ first introduced by \citet{hib79, hib80}. In order to adapt this model
+ to the requirements of coupled ice-ocean state estimation, many
+ important aspects of the original code have been modified and
+ improved \citep{losch10:_mitsim}:
+ \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 three different solution methods for solving the nonlinear
+   momentum equations have been adopted: LSOR \citep{zhang97}, EVP
+   \citep{hun97}, JFNK \citep{lemieux10,losch14:_jfnk};
+ \item ice-ocean stress can be formulated as in \citet{hibler87} or as in
+   \citet{cam08};
+ \item ice variables are advected by sophisticated, conservative
+   advection schemes 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 sea ice model is tightly coupled to the ocean compontent of the
+ MITgcm.  Heat, fresh water fluxes and surface stresses are computed
+ from the atmospheric state and -- by default -- modified by the ice
+ model at every time step.
  The ice dynamics models that are most widely used for large-scale
- climate studies are the viscous-plastic (VP) model \cite{hib79}, the
+ climate studies are the viscous-plastic (VP) model \citep{hib79}, the
- cavitating fluid (CF) model \cite{fla92}, and the
+ cavitating fluid (CF) model \citep{fla92}, and the
- elastic-viscous-plastic (EVP) model \cite{hun97}.  Compared to the VP
+ elastic-viscous-plastic (EVP) model \citep{hun97}.  Compared to the VP
  model, the CF model does not allow ice shear in calculating ice
  motion, stress, and deformation.  EVP models approximate VP by adding
  an elastic term to the equations for easier adaptation to parallel
  computers.  Because of its higher accuracy in plastic solution and
  relatively simpler formulation, compared to the EVP model, we decided
- to use the VP model as the dynamic component of our ice model.  To do
+ to use the VP model as the default dynamic component of our ice
- this we extended the alternating-direction-implicit (ADI) method of
+ model. To do this we extended the line successive over relaxation
- Zhang and Rothrock \cite{zha00} for use in a parallel configuration.
+ (LSOR) method of \citet{zhang97} for use in a parallel
+ configuration. An EVP model and a free-drift implemtation can be
+ selected with runtime flags.
+ \paragraph{Compatibility with ice-thermodynamics package \code{thsice}\label{sec:pkg:seaice:thsice}}~\\
+ %
+ Note, that by default the \code{seaice}-package includes the orginial
+ so-called zero-layer thermodynamics following \citet{hib80} with a
+ snow cover as in \citet{zha98a}. The zero-layer thermodynamic model
+ assumes that ice does not store heat and, therefore, tends to
+ exaggerate the seasonal variability in ice thickness.  This
+ exaggeration can be significantly reduced by using
+ \citeauthor{sem76}'s~[\citeyear{sem76}] three-layer thermodynamic
+ model that permits heat storage in ice. Recently, the three-layer thermodynamic model has been reformulated by
+ \citet{win00}.  The reformulation improves model physics by
+ representing the brine content of the upper ice with a variable heat
+ capacity.  It also improves model numerics and consumes less computer
+ time and memory.
+ The Winton sea-ice thermodynamics have been ported to the MIT GCM;
+ they currently reside under \code{pkg/thsice}.  The package
+ \code{thsice} is described in section~\ref{sec:pkg:thsice}; it is
+ fully compatible with the packages \code{seaice} and \code{exf}. When
+ turned on together with \code{seaice}, the zero-layer thermodynamics
+ are replaced by the Winton thermodynamics. In order to use the
+ \code{seaice}-package with the thermodynamics of \code{thsice},
+ compile both packages and turn both package on in \code{data.pkg}; see
+ an example in \code{global\_ocean.cs32x15/input.icedyn}. Note, that
+ once \code{thsice} is turned on, the variables and diagnostics
+ associated to the default thermodynamics are meaningless, and the
+ diagnostics of \code{thsice} have to be used instead.
+ \paragraph{Surface forcing\label{sec:pkg:seaice:surfaceforcing}}~\\
+ %
  The sea ice model requires the following input fields: 10-m winds, 2-m
  air temperature and specific humidity, downward longwave and shortwave
  radiations, precipitation, evaporation, and river and glacier runoff.
  The sea ice model also requires surface temperature from the ocean
- model and third level horizontal velocity which is used as a proxy for
+ model and the top level horizontal velocity.  Output fields are
- surface geostrophic velocity.  Output fields are surface wind stress,
+ surface wind stress, evaporation minus precipitation minus runoff, net
- evaporation minus precipitation minus runoff, net surface heat flux,
+ surface heat flux, and net shortwave flux.  The sea-ice model is
- and net shortwave flux.  The sea-ice model is global: in ice-free
+ global: in ice-free regions bulk formulae are used to estimate oceanic
- regions bulk formulae are used to estimate oceanic forcing from the
+ forcing from the atmospheric fields.
- atmospheric fields.
+ \paragraph{Dynamics\label{sec:pkg:seaice:dynamics}}~\\
+ %
+ \newcommand{\vek}[1]{\ensuremath{\vec{\mathbf{#1}}}}
+ \newcommand{\vtau}{\vek{\mathbf{\tau}}}
+ 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{cam08};
+ 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.
+ \paragraph{Viscous-Plastic (VP) Rheology\label{sec:pkg:seaice:VPrheology}}~\\
+ %
+ 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{hib79, 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\,\exp\{-C^{*}\cdot(1-c)\},
+ \label{eq:icestrength}
+ \end{equation}
+ with the constants $P^{*}$ (run-time parameter \code{SEAICE\_strength}) and
+ $C^{*}=20$. 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 +
+\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}$ (for numerical reasons, run-time parameter
+ \code{SEAICE\_EPS} with a default value of
+ $10^{-10}\text{\,s}^{-1}$) and a maximum $\zeta_{\max} =
+ P_{\max}/\Delta^*$, where
+ $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. (There is also
+ the option of bounding $\zeta$ from below by setting run-time
+ parameter \code{SEAICE\_zetaMin} $>0$, but this is generally not
+ recommended). 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.
+ Defining the CPP-flag \code{SEAICE\_ZETA\_SMOOTHREG} in
+ \code{SEAICE\_OPTIONS.h} before compiling replaces the method for
+ bounding $\zeta$ by a smooth (differentiable) expression:
+ \begin{equation}
+   \label{eq:zetaregsmooth}
+   \begin{split}
+   \zeta &= \zeta_{\max}\tanh\left(\frac{P}{2\,\min(\Delta,\Delta_{\min})
+       \,\zeta_{\max}}\right)\\
+   &= \frac{P}{2\Delta^*}
+   \tanh\left(\frac{\Delta^*}{\min(\Delta,\Delta_{\min})}\right)
+   \end{split}
+ \end{equation}
+ where $\Delta_{\min}=10^{-20}\text{\,s}^{-1}$ is chosen to avoid divisions
+ by zero.
+ \paragraph{LSR and  JFNK solver \label{sec:pkg:seaice:LSRJFNK}}~\\
+ %
+ % By default, the VP-model is integrated by a Pickwith the
+ % semi-implicit line successive over relaxation (LSOR)-solver of
+ % \citet{zhang97}, 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.
+ \newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}}
+ %
+ In the matrix notation, the discretized momentum equations can be
+ written as
+ \begin{equation}
+   \label{eq:matrixmom}
+   \mat{A}(\vek{x})\,\vek{x} = \vek{b}(\vek{x}).
+ \end{equation}
+ The solution vector $\vek{x}$ consists of the two velocity components
+ $u$ and $v$ that contain the velocity variables at all grid points and
+ at one time level. The standard (and default) method for solving
+ Eq.\,(\ref{eq:matrixmom}) in the sea ice component of the
+ \mbox{MITgcm}, as in many sea ice models, is an iterative Picard
+ solver: in the $k$-th iteration a linearized form
+ $\mat{A}(\vek{x}^{k-1})\,\vek{x}^{k} = \vek{b}(\vek{x}^{k-1})$ is
+ solved (in the case of the MITgcm it is a Line Successive (over)
+ Relaxation (LSR) algorithm \citep{zhang97}).  Picard solvers converge
+ slowly, but generally the iteration is terminated after only a few
+ non-linear steps \citep{zhang97, lemieux09} and the calculation
+ continues with the next time level. This method is the default method
+ in the MITgcm. The number of non-linear iteration steps or pseudo-time
+ steps can be controlled by the runtime parameter
+ \code{NPSEUDOTIMESTEPS} (default is 2).
+ In order to overcome the poor convergence of the Picard-solver,
+ \citet{lemieux10} introduced a Jacobian-free Newton-Krylov solver for
+ the sea ice momentum equations. This solver is also implemented in the
+ MITgcm \citep{losch14:_jfnk}. The Newton method transforms minimizing
+ the residual $\vek{F}(\vek{x}) = \mat{A}(\vek{x})\,\vek{x} -
+ \vek{b}(\vek{x})$ to finding the roots of a multivariate Taylor
+ expansion of the residual \vek{F} around the previous ($k-1$) estimate
+ $\vek{x}^{k-1}$:
+ \begin{equation}
+   \label{eq:jfnktaylor}
+   \vek{F}(\vek{x}^{k-1}+\delta\vek{x}^{k}) =
+   \vek{F}(\vek{x}^{k-1}) + \vek{F}'(\vek{x}^{k-1})\,\delta\vek{x}^{k}
+ \end{equation}
+ with the Jacobian $\mat{J}\equiv\vek{F}'$. The root
+ $\vek{F}(\vek{x}^{k-1}+\delta\vek{x}^{k})=0$ is found by solving
+ \begin{equation}
+   \label{eq:jfnklin}
+   \mat{J}(\vek{x}^{k-1})\,\delta\vek{x}^{k} = -\vek{F}(\vek{x}^{k-1})
+ \end{equation}
+ for $\delta\vek{x}^{k}$. The next ($k$-th) estimate is given by
+ $\vek{x}^{k}=\vek{x}^{k-1}+a\,\delta\vek{x}^{k}$. In order to avoid
+ overshoots the factor $a$ is iteratively reduced in a line search
+ ($a=1, \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \ldots$) until
+ $\|\vek{F}(\vek{x}^k)\| < \|\vek{F}(\vek{x}^{k-1})\|$, where
+ $\|\cdot\|=\int\cdot\,dx^2$ is the $L_2$-norm. In practice, the line
+ search is stopped at $a=\frac{1}{8}$. The line search starts after
+ $\code{SEAICE\_JFNK\_lsIter}$ non-linear Newton iterations (off by
+ default).
+ Forming the Jacobian $\mat{J}$ explicitly is often avoided as ``too
+ error prone and time consuming'' \citep{knoll04:_jfnk}. Instead,
+ Krylov methods only require the action of \mat{J} on an arbitrary
+ vector \vek{w} and hence allow a matrix free algorithm for solving
+ Eq.\,(\ref{eq:jfnklin}) \citep{knoll04:_jfnk}. The action of \mat{J}
+ can be approximated by a first-order Taylor series expansion:
+ \begin{equation}
+   \label{eq:jfnkjacvecfd}
+   \mat{J}(\vek{x}^{k-1})\,\vek{w} \approx
+   \frac{\vek{F}(\vek{x}^{k-1}+\epsilon\vek{w}) - \vek{F}(\vek{x}^{k-1})}
+   {\epsilon}
+ \end{equation}
+ or computed exactly with the help of automatic differentiation (AD)
+ tools. \code{SEAICE\_JFNKepsilon} sets the step size
+ $\epsilon$.
+ We use the Flexible Generalized Minimum RESidual method
+ \citep[FGMRES,][]{saad93:_fgmres} with right-hand side preconditioning
+ to solve Eq.\,(\ref{eq:jfnklin}) iteratively starting from a first
+ guess of $\delta\vek{x}^{k}_{0} = 0$. For the preconditioning matrix
+ \mat{P} we choose a simplified form of the system matrix
+ $\mat{A}(\vek{x}^{k-1})$ \citep{lemieux10} where $\vek{x}^{k-1}$ is
+ the estimate of the previous Newton step $k-1$. The transformed
+ equation\,(\ref{eq:jfnklin}) becomes
+ \begin{equation}
+   \label{eq:jfnklinpc}
+   \mat{J}(\vek{x}^{k-1})\,\mat{P}^{-1}\delta\vek{z} =
+   -\vek{F}(\vek{x}^{k-1}),
+   \quad\text{with}\quad \delta\vek{z}=\mat{P}\delta\vek{x}^{k}.
+ \end{equation}
+ The Krylov method iteratively improves the approximate solution
+ to~(\ref{eq:jfnklinpc}) in subspace ($\vek{r}_0$,
+ $\mat{J}\mat{P}^{-1}\vek{r}_0$, $(\mat{J}\mat{P}^{-1})^2\vek{r}_0$,
+ \ldots, $(\mat{J}\mat{P}^{-1})^m\vek{r}_0$) with increasing $m$;
+ $\vek{r}_0 = -\vek{F}(\vek{x}^{k-1})
+ -\mat{J}(\vek{x}^{k-1})\,\delta\vek{x}^{k}_{0}$
+ %-\vek{F}(\vek{x}^{k-1})
+ %-\mat{J}(\vek{x}^{k-1})\,\mat{P}^{-1}\delta\vek{z}$
+ is the initial residual of
+ (\ref{eq:jfnklin}); $\vek{r}_0=-\vek{F}(\vek{x}^{k-1})$ with the first
+ guess $\delta\vek{x}^{k}_{0}=0$. We allow a Krylov-subspace of
+ dimension~$m=50$ and we do not use restarts. The preconditioning operation
+ involves applying $\mat{P}^{-1}$ to the basis vectors $\vek{v}_0,
+ \vek{v}_1, \vek{v}_2, \ldots, \vek{v}_m$ of the Krylov subspace. This
+ operation is approximated by solving the linear system
+ $\mat{P}\,\vek{w}=\vek{v}_i$. Because $\mat{P} \approx
+ \mat{A}(\vek{x}^{k-1})$, we can use the LSR-algorithm \citep{zhang97}
+ already implemented in the Picard solver. Each preconditioning
+ operation uses a fixed number of 10~LSR-iterations avoiding any
+ termination criterion. More details and results can be found in
+ \citet{lemieux10, losch14:_jfnk}.
+ To use the JFNK-solver set \code{SEAICEuseJFNK = .TRUE.} in the
+ namelist file \code{data.seaice}; \code{SEAICE\_ALLOW\_JFNK} needs to
+ be defined in \code{SEAICE\_OPTIONS.h} and we recommend using a smooth
+ regularization of $\zeta$ by defining \code{SEAICE\_ZETA\_SMOOTHREG}
+ (see above) for better convergence.  The non-linear Newton iteration
+ is terminated when the $L_2$-norm of the residual is reduced by
+ $\gamma_{\mathrm{nl}}$ (runtime parameter \code{JFNKgamma\_nonlin =
+.e-4} will already lead to expensive simulations) with respect to
+ the initial norm: $\|\vek{F}(\vek{x}^k)\| <
+ \gamma_{\mathrm{nl}}\|\vek{F}(\vek{x}^0)\|$.  Within a non-linear
+ iteration, the linear FGMRES solver is terminated when the residual is
+ smaller than $\gamma_k\|\vek{F}(\vek{x}^{k-1})\|$ where $\gamma_k$ is
+ determined by
+ \begin{equation}
+   \label{eq:jfnkgammalin}
+   \gamma_k =
+   \begin{cases}
+     \gamma_0 &\text{for $\|\vek{F}(\vek{x}^{k-1})\| \geq r$},  \\
+     \max\left(\gamma_{\min},
+     \frac{\|\vek{F}(\vek{x}^{k-1})\|}{\|\vek{F}(\vek{x}^{k-2})\|}\right)
+ %    \phi\left(\frac{\|\vek{F}(\vek{x}^{k-1})\|}{\|\vek{F}(\vek{x}^{k-2})\|}\right)^\alpha\right)
+     &\text{for $\|\vek{F}(\vek{x}^{k-1})\| < r$,}
+   \end{cases}
+ \end{equation}
+ so that the linear tolerance parameter $\gamma_k$ decreases with the
+ non-linear Newton step as the non-linear solution is approached. This
+ inexact Newton method is generally more robust and computationally
+ more efficient than exact methods \citep[e.g.,][]{knoll04:_jfnk}.
+ % \footnote{The general idea behind
+ %   inexact Newton methods is this: The Krylov solver is ``only''
+ %   used to find an intermediate solution of the linear
+ %   equation\,(\ref{eq:jfnklin}) that is used to improve the approximation of
+ %   the actual equation\,(\ref{eq:matrixmom}). With the choice of a
+ %   relatively weak lower limit for FGMRES convergence
+ %   $\gamma_{\min}$ we make sure that the time spent in the FGMRES
+ %   solver is reduced at the cost of more Newton iterations. Newton
+ %   iterations are cheaper than Krylov iterations so that this choice
+ %   improves the overall efficiency.}
+ Typical parameter choices are
+ $\gamma_0=\code{JFNKgamma\_lin\_max}=0.99$,
+ $\gamma_{\min}=\code{JFNKgamma\_lin\_min}=0.1$, and $r =
+ \code{JFNKres\_tFac}\times\|\vek{F}(\vek{x}^{0})\|$ with
+ $\code{JFNKres\_tFac} = \frac{1}{2}$. We recommend a maximum number of
+ non-linear iterations $\code{SEAICEnewtonIterMax} = 100$ and a maximum
+ number of Krylov iterations $\code{SEAICEkrylovIterMax} = 50$, because
+ the Krylov subspace has a fixed dimension of 50.
+ \paragraph{Elastic-Viscous-Plastic (EVP) Dynamics\label{sec:pkg:seaice:EVPdynamics}}~\\
+ %
+ \citet{hun97}'s introduced an elastic contribution to the strain
+ rate in order to regularize Eq.~\ref{eq: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{hun97}, the
+ EVP-model should be stepped forward in time 120 times
+ ($\code{SEAICE\_deltaTevp} = \code{SEAICIE\_deltaTdyn}/120$) 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{hun97}. 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 =
+\dot{\epsilon}_{12}$, respectively, and using the above
+ abbreviations, the equations~\ref{eq: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}$.  $E_{0} = \frac{1}{3}$ is the default
+ value in the code and close to what \citet{hun97} and
+ \citet{hun01} recommend.
+ To use the EVP solver, make sure that both \code{SEAICE\_CGRID} and
+ \code{SEAICE\_ALLOW\_EVP} are defined in \code{SEAICE\_OPTIONS.h}
+ (default). The solver is turned on by setting the sub-cycling time
+ step \code{SEAICE\_deltaTevp} to a value larger than zero. The
+ choice of this time step is under debate. \citet{hun97} recommend
+ order(120) time steps for the EVP solver within one model time step
+ $\Delta{t}$ (\code{deltaTmom}). One can also choose order(120) time
+ steps within the forcing time scale, but then we recommend adjusting
+ the damping time scale $T$ accordingly, by setting either
+ \code{SEAICE\_elasticParm} ($E_{0}$), so that
+ $E_{0}\Delta{t}=\mbox{forcing time scale}$, or directly
+ \code{SEAICE\_evpTauRelax} ($T$) to the forcing time scale.
+ \paragraph{More stable variant of Elastic-Viscous-Plastic Dynamics:  EVP*\label{sec:pkg:seaice:EVPstar}}~\\
+ %
+ The genuine EVP schemes appears to give noisy solutions \citep{hun01,
+   lemieux12, bouillon13}. This has lead to a modified EVP or EVP*
+ \citep{lemieux12, bouillon13, kimmritz15}; here, refer to these
+ variants by EVP*. The main idea is to modify the ``natural''
+ time-discretization of the momentum equations:
+ \begin{equation}
+   \label{eq:evpstar}
+   m\frac{D\vec{u}}{Dt} \approx m\frac{u^{p+1}-u^{n}}{\Delta{t}}
+   + \beta^{*}\frac{u^{p+1}-u^{p}}{\Delta{t}_{\mathrm{EVP}}}
+ \end{equation}
+ where $n$ is the previous time step index, and $p$ is the previous
+ sub-cycling index. The term allows the definition of a residual
+ $|u^{p+1}-u^{p}|$ that, as $u^{p+1} \rightarrow u^{n+1}$, converges to
+ $0$ and a re-interpretation of EVP as a pure iterative solver where
+ the sub-cycling has lost all time-relation \citep{bouillon13,
+   kimmritz15}. Using the terminology of \citet{kimmritz15}, the
+ evolution equations of stress $\sigma_{ij}$ and momentum $\vec{u}$ can
+ be written as:
+ \begin{align}
+   \label{eq:evpstarsigma}
+   \sigma_{ij}^{p+1}&=\sigma_{ij}^p+\frac{1}{\alpha}
+   \Big(\sigma_{ij}(\vec{u}^p)-\sigma_{ij}^p\Big),
+   \phantom{\int}\\
+   \label{eq:evpstarmom}
+   \vec{u}^{p+1}&=\vec{u}^p+\frac{1}{\beta}
+   \Big(\frac{\Delta t}{m}\nabla \cdot{\bf \sigma}^{p+1}+
+   \frac{\Delta t}{m}\vec{R}^{p+1/2}+\vec{u}_n-\vec{u}^p\Big).
+ \end{align}
+ $\vec{R}$ contains all terms in the momentum equations except for the
+ rheology terms and the time derivative, $\alpha$ and $\beta$ are free
+ parameters (\code{SEAICE\_evpAlpha}, \code{SEAICE\_evpBeta}) that
+ replace the time stepping parameters \code{SEAICE\_deltaTevp}
+ ($\Delta{T}_{\mathrm{EVP}}$), \code{SEAICE\_elasticParm} ($E_{0}$), or
+ \code{SEAICE\_evpTauRelax} ($T$). $\alpha$ and $\beta$ determine the
+ speed of convergence and the stability. Usually, it makes sense to use
+ $\alpha = \beta$, and \code{SEAICEnEVPstarSteps} $>> \alpha = \beta$
+ \citep{kimmritz15}.
+ In order to use EVP* in the MITgcm, set \code{SEAICEuseEVPstar =
+   .TRUE.,} in \code{data.seaice}. \code{SEAICEuseEVPrev =.TRUE.,} uses
+ the actual form of equations (\ref{eq:evpstarsigma}) and
+ (\ref{eq:evpstarmom}) with fewer implicit terms and the factor of
+ $e^{2}$ dropped in the stress equations (\ref{eq:evpstresstensor2})
+ and (\ref{eq:evpstresstensor12}).  This turns out to improve
+ convergence \citep{bouillon13}.
+ Note, that for historical reasons, \code{SEAICE\_deltaTevp} needs to
+ be set to some value in order to use also EVP*. Also note, that
+ probably because of the C-grid staggering of velocities and stresses,
+ EVP* does not converge as successfully as in \citet{kimmritz15}.
+ \paragraph{Truncated ellipse method (TEM) for yield curve \label{sec:pkg:seaice:TEM}}~\\
+ %
+ 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}
+ To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in
+ \code{SEAICE\_OPTIONS.h} and turn it on with
+ \code{SEAICEuseTEM} in \code{data.seaice}.
+ \paragraph{Ice-Ocean stress \label{sec:pkg:seaice:iceoceanstress}}~\\
+ %
+ Moving sea ice exerts a stress on the ocean which is the opposite of
+ the stress $\vtau_{ocean}$ in Eq.~\ref{eq: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.
+ To turn on the stress formulation of \citet{hibler87}, set
+ \code{useHB87StressCoupling=.TRUE.} in \code{data.seaice}.
+ % Our discretization differs from \citet{zhang97, 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.
+ \paragraph{Finite-volume discretization of the stress tensor
+   divergence\label{sec:pkg:seaice:discretization}}~\\
+ %
+ On an Arakawa C~grid, ice thickness and concentration and thus ice
+ strength $P$ and bulk and shear viscosities $\zeta$ and $\eta$ are
+ naturally defined a C-points in the center of the grid
+ cell. Discretization requires only averaging of $\zeta$ and $\eta$ to
+ vorticity or Z-points (or $\zeta$-points, but here we use Z in order
+ avoid confusion with the bulk viscosity) at the bottom left corner of
+ the cell to give $\overline{\zeta}^{Z}$ and $\overline{\eta}^{Z}$. In
+ the following, the superscripts indicate location at Z or C points,
+ distance across the cell (F), along the cell edge (G), between
+ $u$-points (U), $v$-points (V), and C-points (C). The control volumes
+ of the $u$- and $v$-equations in the grid cell at indices $(i,j)$ are
+ $A_{i,j}^{w}$ and $A_{i,j}^{s}$, respectively. With these definitions
+ (which follow the model code documentation except that $\zeta$-points
+ have been renamed to Z-points), the strain rates are discretized as:
+ \begin{align}
+   \dot{\epsilon}_{11} &= \partial_{1}{u}_{1} + k_{2}u_{2} \\ \notag
+   => (\epsilon_{11})_{i,j}^C &= \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}}
+    + k_{2,i,j}^{C}\frac{v_{i,j+1}+v_{i,j}}{2} \\
+   \dot{\epsilon}_{22} &= \partial_{2}{u}_{2} + k_{1}u_{1} \\\notag
+   => (\epsilon_{22})_{i,j}^C &= \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}}
+    + k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\
+    \dot{\epsilon}_{12} = \dot{\epsilon}_{21} &= \frac{1}{2}\biggl(
+    \partial_{1}{u}_{2} + \partial_{2}{u}_{1} - k_{1}u_{2} - k_{2}u_{1}
+    \biggr) \\ \notag
+   => (\epsilon_{12})_{i,j}^Z &= \frac{1}{2}
+   \biggl( \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^V}
+    + \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^U} \\\notag
+   &\phantom{=\frac{1}{2}\biggl(}
+    - k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2}
+    - k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2}
+    \biggr),
+ \end{align}
+ so that the diagonal terms of the strain rate tensor are naturally
+ defined at C-points and the symmetric off-diagonal term at
+ Z-points. No-slip boundary conditions ($u_{i,j-1}+u_{i,j}=0$ and
+ $v_{i-1,j}+v_{i,j}=0$ across boundaries) are implemented via
+ ``ghost-points''; for free slip boundary conditions
+ $(\epsilon_{12})^Z=0$ on boundaries.
+ For a spherical polar grid, the coefficients of the metric terms are
+ $k_{1}=0$ and $k_{2}=-\tan\phi/a$, with the spherical radius $a$ and
+ the latitude $\phi$; $\Delta{x}_1 = \Delta{x} = a\cos\phi
+ \Delta\lambda$, and $\Delta{x}_2 = \Delta{y}=a\Delta\phi$. For a
+ general orthogonal curvilinear grid, $k_{1}$ and
+ $k_{2}$ can be approximated by finite differences of the cell widths:
+ \begin{align}
+   k_{1,i,j}^{C} &= \frac{1}{\Delta{y}_{i,j}^{F}}
+   \frac{\Delta{y}_{i+1,j}^{G}-\Delta{y}_{i,j}^{G}}{\Delta{x}_{i,j}^{F}} \\
+   k_{2,i,j}^{C} &= \frac{1}{\Delta{x}_{i,j}^{F}}
+   \frac{\Delta{x}_{i,j+1}^{G}-\Delta{x}_{i,j}^{G}}{\Delta{y}_{i,j}^{F}} \\
+   k_{1,i,j}^{Z} &= \frac{1}{\Delta{y}_{i,j}^{U}}
+   \frac{\Delta{y}_{i,j}^{C}-\Delta{y}_{i-1,j}^{C}}{\Delta{x}_{i,j}^{V}} \\
+   k_{2,i,j}^{Z} &= \frac{1}{\Delta{x}_{i,j}^{V}}
+   \frac{\Delta{x}_{i,j}^{C}-\Delta{x}_{i,j-1}^{C}}{\Delta{y}_{i,j}^{U}}
+ \end{align}
+ The stress tensor is given by the constitutive viscous-plastic
+ relation $\sigma_{\alpha\beta} = 2\eta\dot{\epsilon}_{\alpha\beta} +
+ [(\zeta-\eta)\dot{\epsilon}_{\gamma\gamma} - P/2
+ ]\delta_{\alpha\beta}$ \citep{hib79}. The stress tensor divergence
+ $(\nabla\sigma)_{\alpha} = \partial_\beta\sigma_{\beta\alpha}$, is
+ discretized in finite volumes \citep[see
+ also][]{losch10:_mitsim}. This conveniently avoids dealing with
+ further metric terms, as these are ``hidden'' in the differential cell
+ widths. For the $u$-equation ($\alpha=1$) we have:
+ \begin{align}
+   (\nabla\sigma)_{1}: \phantom{=}&
+   \frac{1}{A_{i,j}^w}
+   \int_{\mathrm{cell}}(\partial_1\sigma_{11}+\partial_2\sigma_{21})\,dx_1\,dx_2
+   \\\notag
+   =& \frac{1}{A_{i,j}^w} \biggl\{
+   \int_{x_2}^{x_2+\Delta{x}_2}\sigma_{11}dx_2\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
+   + \int_{x_1}^{x_1+\Delta{x}_1}\sigma_{21}dx_1\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
+   \biggr\} \\ \notag
+   \approx& \frac{1}{A_{i,j}^w} \biggl\{
+   \Delta{x}_2\sigma_{11}\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
+   + \Delta{x}_1\sigma_{21}\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
+   \biggr\} \\ \notag
+   =& \frac{1}{A_{i,j}^w} \biggl\{
+   (\Delta{x}_2\sigma_{11})_{i,j}^C -
+   (\Delta{x}_2\sigma_{11})_{i-1,j}^C
+   \\\notag
+   \phantom{=}& \phantom{\frac{1}{A_{i,j}^w} \biggl\{}
+   + (\Delta{x}_1\sigma_{21})_{i,j+1}^Z - (\Delta{x}_1\sigma_{21})_{i,j}^Z
+   \biggr\}
+ \end{align}
+ with
+ \begin{align}
+   (\Delta{x}_2\sigma_{11})_{i,j}^C =& \phantom{+}
+   \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
+   \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag
+   &+ \Delta{y}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
+   k_{2,i,j}^C \frac{v_{i,j+1}+v_{i,j}}{2} \\ \notag
+   \phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
+   \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} \\ \notag
+   \phantom{=}& + \Delta{y}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
+   k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag
+   \phantom{=}& - \Delta{y}_{i,j}^{F} \frac{P}{2} \\
+   (\Delta{x}_1\sigma_{21})_{i,j}^Z =& \phantom{+}
+   \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
+   \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}} \\ \notag
+   & + \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
+   \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\ \notag
+   & - \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
+   k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2} \\ \notag
+   & - \Delta{x}_{i,j}^{V}\overline{\eta}^{Z}_{i,j}
+   k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2}
+ \end{align}
+ Similarly, we have for the $v$-equation ($\alpha=2$):
+ \begin{align}
+   (\nabla\sigma)_{2}: \phantom{=}&
+   \frac{1}{A_{i,j}^s}
+   \int_{\mathrm{cell}}(\partial_1\sigma_{12}+\partial_2\sigma_{22})\,dx_1\,dx_2
+   \\\notag
+   =& \frac{1}{A_{i,j}^s} \biggl\{
+   \int_{x_2}^{x_2+\Delta{x}_2}\sigma_{12}dx_2\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
+   + \int_{x_1}^{x_1+\Delta{x}_1}\sigma_{22}dx_1\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
+   \biggr\} \\ \notag
+   \approx& \frac{1}{A_{i,j}^s} \biggl\{
+   \Delta{x}_2\sigma_{12}\biggl|_{x_{1}}^{x_{1}+\Delta{x}_{1}}
+   + \Delta{x}_1\sigma_{22}\biggl|_{x_{2}}^{x_{2}+\Delta{x}_{2}}
+   \biggr\} \\ \notag
+   =& \frac{1}{A_{i,j}^s} \biggl\{
+   (\Delta{x}_2\sigma_{12})_{i+1,j}^Z - (\Delta{x}_2\sigma_{12})_{i,j}^Z
+   \\ \notag
+   \phantom{=}& \phantom{\frac{1}{A_{i,j}^s} \biggl\{}
+   + (\Delta{x}_1\sigma_{22})_{i,j}^C - (\Delta{x}_1\sigma_{22})_{i,j-1}^C
+   \biggr\}
+ \end{align}
+ with
+ \begin{align}
+   (\Delta{x}_1\sigma_{12})_{i,j}^Z =& \phantom{+}
+   \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
+   \frac{u_{i,j}-u_{i,j-1}}{\Delta{y}_{i,j}^{U}}
+   \\\notag &
+   + \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
+   \frac{v_{i,j}-v_{i-1,j}}{\Delta{x}_{i,j}^{V}} \\\notag
+   &- \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
+   k_{2,i,j}^{Z}\frac{u_{i,j}+u_{i,j-1}}{2}
+   \\\notag &
+   - \Delta{y}_{i,j}^{U}\overline{\eta}^{Z}_{i,j}
+   k_{1,i,j}^{Z}\frac{v_{i,j}+v_{i-1,j}}{2} \\ \notag
+   (\Delta{x}_2\sigma_{22})_{i,j}^C =& \phantom{+}
+   \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
+   \frac{u_{i+1,j}-u_{i,j}}{\Delta{x}_{i,j}^{F}} \\ \notag
+   &+ \Delta{x}_{i,j}^{F}(\zeta - \eta)^{C}_{i,j}
+   k_{2,i,j}^{C} \frac{v_{i,j+1}+v_{i,j}}{2} \\ \notag
+   & + \Delta{x}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
+   \frac{v_{i,j+1}-v_{i,j}}{\Delta{y}_{i,j}^{F}} \\ \notag
+   & + \Delta{x}_{i,j}^{F}(\zeta + \eta)^{C}_{i,j}
+   k_{1,i,j}^{C}\frac{u_{i+1,j}+u_{i,j}}{2} \\ \notag
+   & -\Delta{x}_{i,j}^{F} \frac{P}{2}
+ \end{align}
+ Again, no slip boundary conditions are realized via ghost points and
+ $u_{i,j-1}+u_{i,j}=0$ and $v_{i-1,j}+v_{i,j}=0$ across boundaries. For
+ free slip boundary conditions the lateral stress is set to zeros. In
+ analogy to $(\epsilon_{12})^Z=0$ on boundaries, we set
+ $\sigma_{21}^{Z}=0$, or equivalently $\eta_{i,j}^{Z}=0$, on boundaries.
+ \paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}}~\\
+ %
+ In its original formulation the sea ice model \citep{menemenlis05}
+ uses simple thermodynamics following the appendix of
+ \citet{sem76}. 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}. To use this thickness category parameterization set
+ \code{\#define SEAICE\_MULTICATEGORY}; note that this requires
+ different restart files and switching this flag on in the middle of an
+ integration is not possible.
+ 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. This flux
+ is not assumed to instantaneously melt or create ice, but a time scale
+ of three days (run-time parameter \code{SEAICE\_gamma\_t}) is used
+ to relax $T_{w}$ to the freezing point.
+ %
+ The parameterization of lateral and vertical growth of sea ice follows
+ that of \citet{hib79, hib80}; the so-called lead closing parameter
+ $h_{0}$ (run-time parameter \code{HO}) has a default value of
+.5~meters.
+ On top of the ice there is a layer of snow that modifies the heat flux
+ and the albedo \citep{zha98a}. 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}. The flood-freeze algorithm is enabled with the CPP-flag
+ \code{SEAICE\_ALLOW\_FLOODING} and turned on with run-time parameter
+ \code{SEAICEuseFlooding=.true.}.
+ \paragraph{Advection of thermodynamic variables\label{sec:pkg:seaice:advection}}~\\
+ %
+ 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, we
+ recommend flux-limited schemes \citep[multidimensional 2nd and
+rd-order advection scheme with flux limiter][]{roe:85, 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 schemes conserve volume
+ and horizontal area and are unconditionally stable, so that we can set
+ $D_{X}=0$. Run-timeflags: \code{SEAICEadvScheme} (default=2, is the
+ historic 2nd-order, centered difference scheme), \code{DIFF1} =
+ $D_{X}/\Delta{x}$
+ (default=0.004).
+ The MITgcm sea ice model provides the option to use
+ the thermodynamics model of \citet{win00}, which in turn is based on
+ the 3-layer model of \citet{sem76} and which treats brine content by
+ means of enthalpy conservation; the corresponding package
+ \code{thsice} is described in section~\ref{sec:pkg:thsice}. This
+ scheme requires additional state variables, namely the enthalpy of the
+ two ice layers (instead of effective ice salinity), to be advected by
+ ice velocities.
+ %
+ The internal sea ice temperature is inferred from ice enthalpy.  To
+ avoid unphysical (negative) values for ice thickness and
+ concentration, a positive 2nd-order advection scheme with a SuperBee
+ flux limiter \citep{roe:85} should be used to advect all
+ sea-ice-related quantities of the \citet{win00} thermodynamic model
+ (runtime flag \code{thSIceAdvScheme=77} and
+ \code{thSIce\_diffK}=$D_{X}$=0 in \code{data.ice}, defaults are 0).  Because of the
+ non-linearity of the advection scheme, care must be taken in advecting
+ these quantities: when simply using ice velocity to advect enthalpy,
+ the total energy (i.e., the volume integral of enthalpy) is not
+ conserved. Alternatively, one can advect the energy content (i.e.,
+ product of ice-volume and enthalpy) but then false enthalpy extrema
+ can occur, which then leads to unrealistic ice temperature.  In the
+ currently implemented solution, the sea-ice mass flux is used to
+ advect the enthalpy in order to ensure conservation of enthalpy and to
+ prevent false enthalpy extrema. %
  %----------------------------------------------------------------------
  \subsubsection{Key subroutines
  \label{sec:pkg:seaice:subroutines}}
- Top-level routine: \texttt{exf\_getforcing.F}
+ Top-level routine: \code{seaice\_model.F}
  {\footnotesize
  \begin{verbatim}
-Line 197 
 c  seaice_model (TOP LEVEL ROUTINE)
+Line 977 
 c  seaice_model (TOP LEVEL ROUTINE)
  c  |
  c  |-- #ifdef SEAICE_CGRID
  c  |     SEAICE_DYNSOLVER
+ c  |     |
+ c  |     |-- < compute proxy for geostrophic velocity >
+ c  |     |
+ c  |     |-- < set up mass per unit area and Coriolis terms >
+ c  |     |
+ c  |     |-- < dynamic masking of areas with no ice >
+ c  |     |
+ c  |     |
  c  |   #ELSE
  c  |     DYNSOLVER
  c  |   #ENDIF
  c  |
- c  ...
+ c  |-- if ( useOBCS )
+ c  |     OBCS_APPLY_UVICE
+ c  |
+ c  |-- if ( SEAICEadvHeff .OR. SEAICEadvArea .OR. SEAICEadvSnow .OR. SEAICEadvSalt )
+ c  |     SEAICE_ADVDIFF
+ c  |
+ c  |-- if ( usePW79thermodynamics )
+ c  |     SEAICE_GROWTH
+ c  |
+ c  |-- if ( useOBCS )
+ c  |     if ( SEAICEadvHeff ) OBCS_APPLY_HEFF
+ c  |     if ( SEAICEadvArea ) OBCS_APPLY_AREA
+ c  |     if ( SEAICEadvSALT ) OBCS_APPLY_HSALT
+ c  |     if ( SEAICEadvSNOW ) OBCS_APPLY_HSNOW
+ c  |
+ c  |-- < do various exchanges >
+ c  |
+ c  |-- < do additional diagnostics >
+ c  |
+ c  o
  \end{verbatim}
  }
-Line 209 
 c  ...
+Line 1017 
 c  ...
  %----------------------------------------------------------------------
- \subsubsection{EXF diagnostics
+ \subsubsection{SEAICE diagnostics
  \label{sec:pkg:seaice:diagnostics}}
  Diagnostics output is available via the diagnostics package
-Line 217 
 Diagnostics output is available via the
+Line 1025 
 Diagnostics output is available via the
  Available output fields are summarized in
  Table \ref{tab:pkg:seaice:diagnostics}.
- \begin{table}[h!]
+ \input{s_phys_pkgs/text/seaice_diags.tex}
- \centering
- \label{tab:pkg:seaice:diagnostics}
- {\footnotesize
- \begin{verbatim}
- ---------+----+----+----------------+-----------------
-  <-Name->|Levs|grid|<--  Units   -->|<- Tile (max=80c)
- ---------+----+----+----------------+-----------------
-  SIarea  |  1 |SM  |m^2/m^2         |SEAICE fractional ice-covered area [0 to 1]
-  SIheff  |  1 |SM  |m               |SEAICE effective ice thickness
-  SIuice  |  1 |UU  |m/s             |SEAICE zonal ice velocity, >0 from West to East
-  SIvice  |  1 |VV  |m/s             |SEAICE merid. ice velocity, >0 from South to North
-  SIhsnow |  1 |SM  |m               |SEAICE snow thickness
-  SIhsalt |  1 |SM  |g/m^2           |SEAICE effective salinity
-  SIatmFW |  1 |SM  |m/s             |Net freshwater flux from the atmosphere (+=down)
-  SIuwind |  1 |SM  |m/s             |SEAICE zonal 10-m wind speed, >0 increases uVel
-  SIvwind |  1 |SM  |m/s             |SEAICE meridional 10-m wind speed, >0 increases uVel
-  SIfu    |  1 |UU  |N/m^2           |SEAICE zonal surface wind stress, >0 increases uVel
-  SIfv    |  1 |VV  |N/m^2           |SEAICE merid. surface wind stress, >0 increases vVel
-  SIempmr |  1 |SM  |m/s             |SEAICE upward freshwater flux, > 0 increases salt
-  SIqnet  |  1 |SM  |W/m^2           |SEAICE upward heatflux, turb+rad, >0 decreases theta
-  SIqsw   |  1 |SM  |W/m^2           |SEAICE upward shortwave radiat., >0 decreases theta
-  SIpress |  1 |SM  |m^2/s^2         |SEAICE strength (with upper and lower limit)
-  SIzeta  |  1 |SM  |m^2/s           |SEAICE nonlinear bulk viscosity
-  SIeta   |  1 |SM  |m^2/s           |SEAICE nonlinear shear viscosity
-  SIsigI  |  1 |SM  |no units        |SEAICE normalized principle stress, component one
-  SIsigII |  1 |SM  |no units        |SEAICE normalized principle stress, component two
-  SIthdgrh|  1 |SM  |m/s             |SEAICE thermodynamic growth rate of effective ice thickness
-  SIsnwice|  1 |SM  |m/s             |SEAICE ice formation rate due to flooding
-  SIuheff |  1 |UU  |m^2/s           |Zonal Transport of effective ice thickness
-  SIvheff |  1 |VV  |m^2/s           |Meridional Transport of effective ice thickness
-  ADVxHEFF|  1 |UU  |m.m^2/s         |Zonal      Advective Flux of eff ice thickn
-  ADVyHEFF|  1 |VV  |m.m^2/s         |Meridional Advective Flux of eff ice thickn
-  DFxEHEFF|  1 |UU  |m.m^2/s         |Zonal      Diffusive Flux of eff ice thickn
-  DFyEHEFF|  1 |VV  |m.m^2/s         |Meridional Diffusive Flux of eff ice thickn
-  ADVxAREA|  1 |UU  |m^2/m^2.m^2/s   |Zonal      Advective Flux of fract area
-  ADVyAREA|  1 |VV  |m^2/m^2.m^2/s   |Meridional Advective Flux of fract area
-  DFxEAREA|  1 |UU  |m^2/m^2.m^2/s   |Zonal      Diffusive Flux of fract area
-  DFyEAREA|  1 |VV  |m^2/m^2.m^2/s   |Meridional Diffusive Flux of fract area
-  ADVxSNOW|  1 |UU  |m.m^2/s         |Zonal      Advective Flux of eff snow thickn
-  ADVySNOW|  1 |VV  |m.m^2/s         |Meridional Advective Flux of eff snow thickn
-  DFxESNOW|  1 |UU  |m.m^2/s         |Zonal      Diffusive Flux of eff snow thickn
-  DFyESNOW|  1 |VV  |m.m^2/s         |Meridional Diffusive Flux of eff snow thickn
-  ADVxSSLT|  1 |UU  |psu.m^2/s       |Zonal      Advective Flux of seaice salinity
-  ADVySSLT|  1 |VV  |psu.m^2/s       |Meridional Advective Flux of seaice salinity
-  DFxESSLT|  1 |UU  |psu.m^2/s       |Zonal      Diffusive Flux of seaice salinity
-  DFyESSLT|  1 |VV  |psu.m^2/s       |Meridional Diffusive Flux of seaice salinity
- \end{verbatim}
- }
- \caption{~}
- \end{table}
  %\subsubsection{Package Reference}
-Line 276 
 Table \ref{tab:pkg:seaice:diagnostics}.
+Line 1033 
 Table \ref{tab:pkg:seaice:diagnostics}.
  \label{sec:pkg:seaice:experiments}
  \begin{itemize}
- \item{Labrador Sea experiment in lab\_sea verification directory. }
+ \item{Labrador Sea experiment in \code{lab\_sea} verification directory. }
+ \item \code{seaice\_obcs}, based on \code{lab\_sea}
+ \item \code{offline\_exf\_seaice/input.seaicetd}, based on \code{lab\_sea}
+ \item \code{global\_ocean.cs32x15/input.icedyn} and
+   \code{global\_ocean.cs32x15/input.seaice}, global
+   cubed-sphere-experiment with combinations of \code{seaice} and
+   \code{thsice}
  \end{itemize}
+ %%% Local Variables:
+ %%% mode: latex
+ %%% TeX-master: "../../manual"
+ %%% End:

 Legend:



Removed from v.1.6
 


changed lines


 
Added in v.1.22
 Legend:



Removed from v.1.6
 


changed lines


 
Added in v.1.22
-Removed from v.1.6
+Added in v.1.22

	ViewVC Help
Powered by ViewVC 1.1.22