--- manual/s_phys_pkgs/text/seaice.tex 2009/05/14 15:35:17 1.9 +++ manual/s_phys_pkgs/text/seaice.tex 2016/03/29 14:50:54 1.25 @@ -1,4 +1,4 @@ -% $Header: /home/ubuntu/mnt/e9_copy/manual/s_phys_pkgs/text/seaice.tex,v 1.9 2009/05/14 15:35:17 mlosch Exp $ +% $Header: /home/ubuntu/mnt/e9_copy/manual/s_phys_pkgs/text/seaice.tex,v 1.25 2016/03/29 14:50:54 mlosch Exp $ % $Name: $ %%EH3 Copied from "MITgcm/pkg/seaice/seaice_description.tex" @@ -16,7 +16,7 @@ %---------------------------------------------------------------------- \subsubsection{Introduction -\label{sec:pkg:exf:intro}} +\label{sec:pkg:seaice:intro}} Package ``seaice'' provides a dynamic and thermodynamic interactive @@ -31,7 +31,7 @@ \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}. %---------------------------------------------------------------------- @@ -58,14 +58,15 @@ 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 -\code{SEAICE\_OPTIONS.h} or in \code{ECCO\_CPPOPTIONS.h}. -Table \ref{tab:pkg:seaice:cpp} summarizes these options. +via CPP preprocessor flags. These options are set in +\code{SEAICE\_OPTIONS.h}. +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}[h!] +\begin{table}[!ht] \centering \label{tab:pkg:seaice:cpp} {\footnotesize @@ -80,22 +81,22 @@ \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} & - enable 8-category thermodynamics (by default undefined)\\ + \code{SEAICE\_ZETA\_SMOOTHREG} & + 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} & - enable "salty" sea-ice (by default undefined) \\ - \code{SEAICE\_AGE} & - enable "age tracer" sea-ice (by default undefined) \\ - \code{SEAICE\_CAP\_HEFF} & - enable capping of sea-ice thickness to MAX\_HEFF \\ \hline + \code{SEAICE\_VARIABLE\_SALINITY} & + enable sea-ice with variable salinity (by default undefined) \\ + \code{SEAICE\_SITRACER} & + enable sea-ice tracer package (by default undefined) \\ \code{SEAICE\_BICE\_STRESS} & B-grid only for backward compatiblity: turn on ice-stress on ocean\\ @@ -105,7 +106,8 @@ \hline \end{tabular} } - \caption{~} + \caption{Some of the most relevant CPP preprocessor flags in the + \code{seaice}-package.} \end{table} %---------------------------------------------------------------------- @@ -113,9 +115,9 @@ \subsubsection{Run-time parameters \label{sec:pkg:seaice:runtime}} -Run-time parameters are set in files -\code{data.pkg} (read in \code{packages\_readparms.F}), -and \code{data.seaice} (read in \code{seaice\_readparms.F}). +Run-time parameters (see Table~\ref{tab:pkg:seaice:runtimeparms}) are set +in files \code{data.pkg} (read in \code{packages\_readparms.F}), and +\code{data.seaice} (read in \code{seaice\_readparms.F}). \paragraph{Enabling the package} ~ \\ @@ -127,9 +129,19 @@ ~ \\ % Table~\ref{tab:pkg:seaice:runtimeparms} lists most run-time parameters. -\input{part6/seaice-parms.tex} - +\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 @@ -170,14 +182,14 @@ 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 - momentum equations have been adopted: LSOR \citep{zhang97}, and EVP - \citep{hun97}; +\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 @@ -202,26 +214,39 @@ 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 -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 pkg/thsice. The package pkg/thsice is fully compatible -with pkg/seaice and with pkg/exf. When turned on together with -pkg/seaice, the zero-layer thermodynamics are replaced by the Winton -thermodynamics. +\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. @@ -232,8 +257,8 @@ 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 @@ -271,6 +296,8 @@ 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}: @@ -290,7 +317,7 @@ 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 @@ -322,28 +349,196 @@ = 2\,\Delta\zeta$ \citep{hibler95} is used so that the stress state always lies on the elliptic yield curve by definition. -In the so-called truncated ellipse method the shear viscosity $\eta$ -is capped to suppress any tensile stress \citep{hibler97, geiger98}: +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: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). + \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} -To enable this method, set \code{\#define SEAICE\_ALLOW\_TEM} in -\code{SEAICE\_OPTIONS.h} and turn it on with -\code{SEAICEuseTEM=.TRUE.} in \code{data.seaice}. +where $\Delta_{\min}=10^{-20}\text{\,s}^{-1}$ is chosen to avoid divisions +by zero. -In the current implementation, the VP-model is integrated with 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. +\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 = + 1.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 @@ -361,11 +556,12 @@ %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 -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 +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 @@ -386,11 +582,12 @@ \frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T} &= \frac{P}{4T\Delta} D_S \end{align} -Here, the elastic parameter $E$ is redefined in terms of a damping timescale -$T$ for elastic waves \[E=\frac{\zeta}{T}.\] -$T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and -the external (long) timestep $\Delta{t}$. \citet{hun97} recommend -$E_{0} = \frac{1}{3}$ (which is the default value in the code). +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} @@ -405,6 +602,101 @@ $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 @@ -434,7 +726,8 @@ % $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 @@ -496,7 +789,8 @@ [(\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} @@ -519,7 +813,9 @@ \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\} - \intertext{with} +\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 @@ -561,7 +857,9 @@ \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\} - \intertext{with} +\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}} @@ -591,8 +889,9 @@ 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, @@ -608,16 +907,22 @@ 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. +parameterize a sub-grid scale distribution for heat flux computations, +the mean ice thickness $h$ is split into $N$ 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,N]$. 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{SEAICE\_multDim} to +the number of desired categories (7 is a good guess, for anything +larger than 7 modify \code{SEAICE\_SIZE.h}) in +\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 @@ -646,6 +951,8 @@ \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: @@ -658,40 +965,42 @@ 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 -choose flux-limited schemes \citep[multidimensional 2nd and 3rd-order -advection scheme with flux limiter][]{roe:85, hundsdorfer94} to -preserve sharp gradients and edges that are typical of sea ice +recommend flux-limited schemes \citep[multidimensional 2nd and +3rd-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 scheme conserve volume +(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), -\code{DIFF1} (default=0.004). - -There is considerable doubt about the reliability of a ``zero-layer'' -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. 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. +$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} is used in this study to advect all -sea-ice-related quantities of the \citet{win00} thermodynamic -model. 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. +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. % %---------------------------------------------------------------------- @@ -757,58 +1066,7 @@ Available output fields are summarized in Table \ref{tab:pkg:seaice:diagnostics}. -\begin{table}[h!] -\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} - +\input{s_phys_pkgs/text/seaice_diags.tex} %\subsubsection{Package Reference} @@ -816,11 +1074,17 @@ \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" +%%% TeX-master: "../../manual" %%% End: