--- manual/s_phys_pkgs/text/seaice.tex 2015/01/21 18:18:22 1.22 +++ 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.22 2015/01/21 18:18:22 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" @@ -63,7 +63,8 @@ Parts of the SEAICE code can be enabled or disabled at compile time via CPP preprocessor flags. These options are set in \code{SEAICE\_OPTIONS.h}. -Table \ref{tab:pkg:seaice:cpp} summarizes the most important ones. +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 @@ -527,6 +528,15 @@ 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 @@ -592,12 +602,14 @@ $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}}~\\ +\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, refer to these -variants by EVP*. The main idea is to modify the ``natural'' +\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} @@ -605,13 +617,14 @@ + \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 +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$ 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: +$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} @@ -620,31 +633,52 @@ \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). + \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 +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}. +$\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 EVP* in the MITgcm, set \code{SEAICEuseEVPstar = - .TRUE.,} in \code{data.seaice}. \code{SEAICEuseEVPrev =.TRUE.,} uses +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}) 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}. - +(\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}}~\\ % @@ -654,8 +688,8 @@ \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). + {\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 @@ -857,6 +891,7 @@ \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, @@ -872,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