/[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.15 by mlosch,
Mon Feb 28 16:27:56 2011 UTC
+revision 1.25 by mlosch,
Tue Mar 29 14:50:54 2016 UTC
 Line 61 
 no additional CPP options are required.
  (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
- \code{SEAICE\_OPTIONS.h} or in \code{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. For
+ more options see the default \code{pkg/seaice/SEAICE\_OPTIONS.h}.
  \begin{table}[!ht]
  \centering
-Line 80 
 Table \ref{tab:pkg:seaice:cpp} summarize
+Line 81 
 Table \ref{tab:pkg:seaice:cpp} summarize
          \code{SEAICE\_CGRID} &
            LSR solver on C-grid (rather than original B-grid) \\
          \code{SEAICE\_ALLOW\_EVP} &
-           use EVP rather than LSR rheology solver \\
+           enable use of EVP rheology solver \\
+         \code{SEAICE\_ALLOW\_JFNK} &
+           enable use of JFNK rheology solver \\
          \code{SEAICE\_EXTERNAL\_FLUXES} &
            use EXF-computed fluxes as starting point \\
-         \code{SEAICE\_MULTICATEGORY} &
+         \code{SEAICE\_ZETA\_SMOOTHREG} &
-           enable 8-category thermodynamics (by default undefined)\\
+           use differentialable regularization for viscosities \\
          \code{SEAICE\_VARIABLE\_FREEZING\_POINT} &
            enable linear dependence of the freezing point on salinity
            (by default undefined)\\
          \code{ALLOW\_SEAICE\_FLOODING} &
            enable snow to ice conversion for submerged sea-ice \\
-         \code{SEAICE\_SALINITY} &
+         \code{SEAICE\_VARIABLE\_SALINITY} &
-           enable "salty" sea-ice (by default undefined) \\
+           enable sea-ice with variable salinity (by default undefined) \\
-         \code{SEAICE\_AGE} &
+         \code{SEAICE\_SITRACER} &
-           enable "age tracer" sea-ice (by default undefined) \\
+           enable sea-ice tracer package (by default undefined) \\
-         \code{SEAICE\_CAP\_HEFF} &
-           enable capping of sea-ice thickness to MAX\_HEFF \\ \hline
          \code{SEAICE\_BICE\_STRESS} &
            B-grid only for backward compatiblity: turn on ice-stress on
            ocean\\
-Line 105 
 Table \ref{tab:pkg:seaice:cpp} summarize
+Line 106 
 Table \ref{tab:pkg:seaice:cpp} summarize
        \hline
      \end{tabular}
    }
-   \caption{~}
+   \caption{Some of the most relevant CPP preprocessor flags in the
+     \code{seaice}-package.}
  \end{table}
  %----------------------------------------------------------------------
-Line 113 
 Table \ref{tab:pkg:seaice:cpp} summarize
+Line 115 
 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
- \code{data.pkg} (read in \code{packages\_readparms.F}),
+ in files \code{data.pkg} (read in \code{packages\_readparms.F}), and
- and \code{data.seaice} (read in \code{seaice\_readparms.F}).
+ \code{data.seaice} (read in \code{seaice\_readparms.F}).
  \paragraph{Enabling the package}
  ~ \\
-Line 139 
 Table~\ref{tab:pkg:seaice:runtimeparms}
+Line 141 
 Table~\ref{tab:pkg:seaice:runtimeparms}
    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};
- \item[\code{IceAgeFile}:] Initial ice age of sea ice averaged over grid
-   cell in seconds; initializes variable \code{ICEAGE};
  \end{description}
  %----------------------------------------------------------------------
 Line 182 
 viscous-plastic (VP) dynamic-thermodynam
  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:
+ 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 two different solution methods for solving the nonlinear
+ \item three different solution methods for solving the nonlinear
-   momentum equations have been adopted: LSOR \citep{zhang97}, and EVP
+   momentum equations have been adopted: LSOR \citep{zhang97}, EVP
-   \citep{hun97};
+   \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
 Line 214 
 relatively simpler formulation, compared
  to use the VP model as the default dynamic component of our ice
  model. To do this we extended the line successive over relaxation
  (LSOR) method of \citet{zhang97} for use in a parallel
- configuration.
+ configuration. An EVP model and a free-drift implemtation can be
+ selected with runtime flags.
- Note, that by default the seaice-package includes the orginial
+ \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
+ \citeauthor{sem76}'s~[\citeyear{sem76}] three-layer thermodynamic
- 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
- thermodynamic model has been reformulated by \citet{win00}.  The
+ \citet{win00}.  The reformulation improves model physics by
- reformulation improves model physics by representing the brine content
+ representing the brine content of the upper ice with a variable heat
- of the upper ice with a variable heat capacity.  It also improves
+ capacity.  It also improves model numerics and consumes less computer
- model numerics and consumes less computer time and memory.  The Winton
+ time and memory.
- sea-ice thermodynamics have been ported to the MIT GCM; they currently
- reside under pkg/thsice. The package pkg/thsice is fully compatible
+ The Winton sea-ice thermodynamics have been ported to the MIT GCM;
- with pkg/seaice and with pkg/exf. When turned on together with
+ they currently reside under \code{pkg/thsice}.  The package
- pkg/seaice, the zero-layer thermodynamics are replaced by the Winton
+ \code{thsice} is described in section~\ref{sec:pkg:thsice}; it is
- thermodynamics.
+ 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.
-Line 244 
 surface heat flux, and net shortwave flu
+Line 257 
 surface heat flux, and net shortwave flu
  global: in ice-free regions bulk formulae are used to estimate oceanic
  forcing from the atmospheric fields.
- \paragraph{Dynamics\label{sec:pkg:seaice:dynamics}}
+ \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
-Line 283 
 air and ocean drag coefficients; $\rho_{
+Line 296 
 air and ocean drag coefficients; $\rho_{
  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}:
-Line 302 
 The ice strain rate is given by
+Line 317 
 The ice strain rate is given by
  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)]},
+   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
-Line 334 
 recommended). For stress tensor computat
+Line 349 
 recommended). For stress tensor computat
  = 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$
+ Defining the CPP-flag \code{SEAICE\_ZETA\_SMOOTHREG} in
- is capped to suppress any tensile stress \citep{hibler97, geiger98}:
+ \code{SEAICE\_OPTIONS.h} before compiling replaces the method for
+ bounding $\zeta$ by a smooth (differentiable) expression:
  \begin{equation}
-   \label{eq:etatem}
+   \label{eq:zetaregsmooth}
-   \eta = \min\left(\frac{\zeta}{e^2},
+   \begin{split}
-   \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
+   \zeta &= \zeta_{\max}\tanh\left(\frac{P}{2\,\min(\Delta,\Delta_{\min})
-   {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2
+       \,\zeta_{\max}}\right)\\
-       +4\dot{\epsilon}_{12}^2}}\right).
+   &= \frac{P}{2\Delta^*}
+   \tanh\left(\frac{\Delta^*}{\min(\Delta,\Delta_{\min})}\right)
+   \end{split}
  \end{equation}
- To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in
+ where $\Delta_{\min}=10^{-20}\text{\,s}^{-1}$ is chosen to avoid divisions
- \code{SEAICE\_OPTIONS.h} and turn it on with
+ by zero.
- \code{SEAICEuseTEM=.TRUE.} in \code{data.seaice}.
- In the current implementation, the VP-model is integrated with the
+ \paragraph{LSR and  JFNK solver \label{sec:pkg:seaice:LSRJFNK}}~\\
- semi-implicit line successive over relaxation (LSOR)-solver of
+ %
- \citet{zhang97}, which allows for long time steps that, in our case,
+ % By default, the VP-model is integrated by a Pickwith the
- are limited by the explicit treatment of the Coriolis term. The
+ % semi-implicit line successive over relaxation (LSOR)-solver of
- explicit treatment of the Coriolis term does not represent a severe
+ % \citet{zhang97}, which allows for long time steps that, in our case,
- limitation because it restricts the time step to approximately the
+ % are limited by the explicit treatment of the Coriolis term. The
- same length as in the ocean model where the Coriolis term is also
+ % explicit treatment of the Coriolis term does not represent a severe
- treated explicitly.
+ % 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.
+ Setting \code{SEAICEuseStrImpCpl = .TRUE.,} turns on ``strength
+ implicit coupling'' \citep{hutchings04} in the LSR-solver and in the
+ LSR-preconditioner for the JFNK-solver. In this mode, the different
+ contributions of the stress divergence terms are re-ordered in order
+ to increase the diagonal dominance of the system
+ matrix. Unfortunately, the convergence rate of the LSR solver is
+ increased only slightly, while the JFNK-convergence appears to be
+ unaffected.
+ \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
-Line 373 
 identical at steady state,
+Line 556 
 identical at steady state,
  %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 is stepped forward in time 120 times within the physical
+ EVP-model should be stepped forward in time 120 times
- ocean model time step (although this parameter is under debate), to
+ ($\code{SEAICE\_deltaTevp} = \code{SEAICIE\_deltaTdyn}/120$) within
- allow for elastic waves to disappear.  Because the scheme does not
+ the physical ocean model time step (although this parameter is under
- require a matrix inversion it is fast in spite of the small internal
+ debate), to allow for elastic waves to disappear.  Because the scheme
- timestep and simple to implement on parallel computers
+ 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
-Line 398 
 abbreviations, the equations~\ref{eq:evp
+Line 582 
 abbreviations, the equations~\ref{eq:evp
    \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
+ Here, the elastic parameter $E$ is redefined in terms of a damping
- $T$ for elastic waves \[E=\frac{\zeta}{T}.\]
+ timescale $T$ for elastic waves \[E=\frac{\zeta}{T}.\]
- $T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and
+ $T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and the external
- the external (long) timestep $\Delta{t}$. \citet{hun97} recommend
+ (long) timestep $\Delta{t}$.  $E_{0} = \frac{1}{3}$ is the default
- $E_{0} = \frac{1}{3}$ (which is the default value in the code).
+ 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}
-Line 417 
 the damping time scale $T$ accordingly,
+Line 602 
 the damping time scale $T$ accordingly,
  $E_{0}\Delta{t}=\mbox{forcing time scale}$, or directly
  \code{SEAICE\_evpTauRelax} ($T$) to the forcing time scale.
+ \paragraph{More stable variants of Elastic-Viscous-Plastic Dynamics:
+   EVP* , mEVP, and aEVP \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, we refer to these
+ variants by modified EVP (mEVP) and adaptive EVP (aEVP)
+ \citep{kimmritz16}. 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 extra ``intertial'' term
+ $m\,(u^{p+1}-u^{n})/\Delta{t})$ allows the definition of a residual
+ $|u^{p+1}-u^{p}|$ that, as $u^{p+1} \rightarrow u^{n+1}$, converges to
+ $0$. In this way EVP can be re-interpreted as a pure iterative solver
+ where the sub-cycling has no association with time-relation (through
+ $\Delta{t}_{\mathrm{EVP}}$) \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}+\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} $\gg
+ (\alpha,\,\beta)$ \citep{kimmritz15}. Currently, there is no
+ termination criterion and the number of mEVP iterations is fixed to
+ \code{SEAICEnEVPstarSteps}.
+ In order to use mEVP in the MITgcm, set \code{SEAICEuseEVPstar =
+   .TRUE.,} in \code{data.seaice}. If \code{SEAICEuseEVPrev =.TRUE.,}
+ the actual form of equations (\ref{eq:evpstarsigma}) and
+ (\ref{eq:evpstarmom}) is used with fewer implicit terms and the factor
+ of $e^{2}$ dropped in the stress equations (\ref{eq:evpstresstensor2})
+ and (\ref{eq:evpstresstensor12}). Although this modifies the original
+ EVP-equations, it turns out to improve convergence \citep{bouillon13}.
+ Another variant is the aEVP scheme \citep{kimmritz16}, where the value
+ of $\alpha$ is set dynamically based on the stability criterion
+ \begin{equation}
+   \label{eq:aevpalpha}
+   \alpha = \beta = \max\left( \tilde{c}\pi\sqrt{c \frac{\zeta}{A_{c}}
+     \frac{\Delta{t}}{\max(m,10^{-4}\text{\,kg})}},\alpha_{\min} \right)
+ \end{equation}
+ with the grid cell area $A_c$ and the ice and snow mass $m$.  This
+ choice sacrifices speed of convergence for stability with the result
+ that aEVP converges quickly to VP where $\alpha$ can be small and more
+ slowly in areas where the equations are stiff. In practice, aEVP leads
+ to an overall better convergence than mEVP \citep{kimmritz16}.
+ %
+ To use aEVP in the MITgcm set \code{SEAICEaEVPcoeff} $= \tilde{c}$;
+ this also sets the default values of \code{SEAICEaEVPcStar} ($c=4$)
+ and \code{SEAICEaEVPalphaMin} ($\alpha_{\min}=5$). Good convergence
+ has been obtained with setting these values \citep{kimmritz16}:
+ \code{SEAICEaEVPcoeff = 0.5, SEAICEnEVPstarSteps = 500,
+   SEAICEuseEVPstar = .TRUE., SEAICEuseEVPrev = .TRUE.}
+ Note, that probably because of the C-grid staggering of velocities and
+ stresses, mEVP may not converge as successfully as in
+ \citet{kimmritz15}, and that convergence at very high resolution
+ (order 5\,km) has not been studied yet.
+ \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{\max(\Delta_{\min}^{2},(\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
-Line 446 
 To turn on the stress formulation of \ci
+Line 726 
 To turn on the stress formulation of \ci
  % $P$ at vorticity points.
  \paragraph{Finite-volume discretization of the stress tensor
-   divergence\label{sec:pkg:seaice:discretization}}
+   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
-Line 508 
 relation $\sigma_{\alpha\beta} = 2\eta\d
+Line 789 
 relation $\sigma_{\alpha\beta} = 2\eta\d
  [(\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. This conveniently avoids dealing with
+ 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}
-Line 607 
 free slip boundary conditions the latera
+Line 889 
 free slip boundary conditions the latera
  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}}
+ \paragraph{Thermodynamics\label{sec:pkg:seaice:thermodynamics}}~\\
+ %
+ \noindent\textbf{NOTE: THIS SECTION IS TERRIBLY OUT OF DATE}\\
  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,
-Line 624 
 way to that of \citet{parkinson79} and \
+Line 907 
 way to that of \citet{parkinson79} and \
  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
+ parameterize a sub-grid scale distribution for heat flux computations,
- computations, the mean ice thickness $h$ is split into seven thickness
+ the mean ice thickness $h$ is split into $N$ thickness categories
- categories $H_{n}$ that are equally distributed between $2h$ and a
+ $H_{n}$ that are equally distributed between $2h$ and a minimum
- minimum imposed ice thickness of $5\text{\,cm}$ by $H_n=
+ imposed ice thickness of $5\text{\,cm}$ by $H_n= \frac{2n-1}{7}\,h$
- \frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each
+ for $n\in[1,N]$. The heat fluxes computed for each thickness category
- thickness category is area-averaged to give the total heat flux
+ is area-averaged to give the total heat flux \citep{hibler84}. To use
- \citep{hibler84}. To use this thickness category parameterization set
+ this thickness category parameterization set \code{SEAICE\_multDim} to
- \code{\#define SEAICE\_MULTICATEGORY}; note that this requires
+ the number of desired categories (7 is a good guess, for anything
- different restart files and switching this flag on in the middle of an
+ larger than 7 modify \code{SEAICE\_SIZE.h}) in
- integration is not possible.
+ \code{data.seaice}; note that this requires different restart files
+ and switching this flag on in the middle of an integration is not
+ advised. In order to include the same distribution for snow, set
+ \code{SEAICE\_useMultDimSnow = .TRUE.}; only then, the
+ parameterization of always having a fraction of thin ice is efficient
+ and generally thicker ice is produced \citep{castro-morales14}.
  The atmospheric heat flux is balanced by an oceanic heat flux from
  below.  The oceanic flux is proportional to
-Line 662 
 principle) turns snow into ice until the
+Line 951 
 principle) turns snow into ice until the
  \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:
-Line 681 
 distributions and to rule out unphysical
+Line 972 
 distributions and to rule out unphysical
  (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}
+ historic 2nd-order, centered difference scheme), \code{DIFF1} =
+ $D_{X}/\Delta{x}$
  (default=0.004).
- There is considerable doubt about the reliability of a ``zero-layer''
+ The MITgcm sea ice model provides the option to use
- 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. 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
-Line 700 
 ice velocities.
+Line 988 
 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} is used in this study to advect all
+ flux limiter \citep{roe:85} should be used to advect all
- sea-ice-related quantities of the \citet{win00} thermodynamic model.
+ sea-ice-related quantities of the \citet{win00} thermodynamic model
- Because of the non-linearity of the advection scheme, care must be
+ (runtime flag \code{thSIceAdvScheme=77} and
- taken in advecting these quantities: when simply using ice velocity to
+ \code{thSIce\_diffK}=$D_{X}$=0 in \code{data.ice}, defaults are 0).  Because of the
- advect enthalpy, the total energy (i.e., the volume integral of
+ non-linearity of the advection scheme, care must be taken in advecting
- enthalpy) is not conserved. Alternatively, one can advect the energy
+ these quantities: when simply using ice velocity to advect enthalpy,
- content (i.e., product of ice-volume and enthalpy) but then false
+ the total energy (i.e., the volume integral of enthalpy) is not
- enthalpy extrema can occur, which then leads to unrealistic ice
+ conserved. Alternatively, one can advect the energy content (i.e.,
- temperature.  In the currently implemented solution, the sea-ice mass
+ product of ice-volume and enthalpy) but then false enthalpy extrema
- flux is used to advect the enthalpy in order to ensure conservation of
+ can occur, which then leads to unrealistic ice temperature.  In the
- enthalpy and to prevent false enthalpy extrema. %
+ currently implemented solution, the sea-ice mass flux is used to
- In order to use the \code{seaice}-package with the more sophisticated
+ advect the enthalpy in order to ensure conservation of enthalpy and to
- thermodynamics of \code{thsice}, compile both packages and turn both
+ prevent false enthalpy extrema. %
- package on in \code{data.pkg}; see an example in
- \code{global\_ocean.cs32x15/input.icedyn}.
  %----------------------------------------------------------------------
-Line 780 
 Diagnostics output is available via the
+Line 1066 
 Diagnostics output is available via the
  Available output fields are summarized in
  Table \ref{tab:pkg:seaice:diagnostics}.
- \begin{table}[!ht]
+ \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  |kg/m^2/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  |kg/m^2/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{Available diagnostics of the seaice-package}
- \end{table}
  %\subsubsection{Package Reference}
-Line 839 
 Table \ref{tab:pkg:seaice:diagnostics}.
+Line 1074 
 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}

 Legend:



Removed from v.1.15
 


changed lines


 
Added in v.1.25
 Legend:



Removed from v.1.15
 


changed lines


 
Added in v.1.25
-Removed from v.1.15
+Added in v.1.25

	ViewVC Help
Powered by ViewVC 1.1.22