--- MITgcm_contrib/articles/ceaice/ceaice_model.tex 2008/02/28 16:34:56 1.3 +++ MITgcm_contrib/articles/ceaice/ceaice_model.tex 2008/06/04 13:33:45 1.9 @@ -1,6 +1,27 @@ -\section{Model} +\section{Model Formulation} \label{sec:model} +The MITgcm sea ice model (MITsim) is based on a variant of the +viscous-plastic (VP) dynamic-thermodynamic sea ice model +\citep{zhang97} first introduced by \citet{hibler79, hibler80}. In +order to adapt this model to the requirements of coupled +ice-ocean simulations, many important aspects of the original code have +been modified and improved: +\begin{itemize} +\item the code has been rewritten for an Arakawa C-grid, both B- and + C-grid variants are available; the C-grid code allows for no-slip + and free-slip lateral boundary conditions; +\item two different solution methods for solving the nonlinear + momentum equations have been adopted: LSOR \citep{zhang97}, EVP + \citep{hunke97}; +\item ice-ocean stress can be formulated as in \citet{hibler87}; +\item ice variables are advected by sophisticated advection schemes; +\item growth and melt parameterization have been refined and extended + in order to allow for automatic differentiation of the code. +\end{itemize} +The model equations and their numerical realization are summarized +below. + \subsection{Dynamics} \label{sec:dynamics} @@ -43,7 +64,7 @@ For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can be related to the ice strain rate and strength by a nonlinear -viscous-plastic (VP) constitutive law \citep{hibler79, zhang98}: +viscous-plastic (VP) constitutive law \citep{hibler79, zhang97}: \begin{equation} \label{eq:vpequation} \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij} @@ -77,7 +98,7 @@ \left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right) (1+e^{-2}) + 4e^{-2}\dot{\epsilon}_{12}^2 + 2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2}) - \right]^{-\frac{1}{2}} + \right]^{\frac{1}{2}} \end{align*} The bulk viscosities are bounded above by imposing both a minimum $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a @@ -126,7 +147,8 @@ EVP-model is stepped forward in time 120 times within the physical ocean model time step (although this parameter is under debate), to allow for elastic waves to disappear. Because the scheme does not -require a matrix inversion it is fast in spite of the small timestep +require a matrix inversion it is fast in spite of the small internal +timestep and simple to implement on parallel computers \citep{hunke97}. For completeness, we repeat the equations for the components of the stress tensor $\sigma_{1} = \sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and @@ -134,8 +156,8 @@ \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension and shearing strain rates, $D_T = \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S = -2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations, -the equations can be written as: +2\dot{\epsilon}_{12}$, respectively, and using the above +abbreviations, the equations\refeq{evpequation} can be written as: \begin{align} \label{eq:evpstresstensor1} \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} + @@ -205,20 +227,14 @@ velocity and the ice velocity leading to an inconsistency as the ice temperature and salinity are different from the oceanic variables. -Sea ice distributions are characterized by sharp gradients and edges. -For this reason, we employ positive, multidimensional 2nd and 3rd-order -advection scheme with flux limiter \citep{roe85, hundsdorfer94} for the -thermodynamic variables discussed in the next section. - -\subparagraph{boundary conditions: no-slip, free-slip, half-slip} - -\begin{itemize} -\item transition from existing B-Grid to C-Grid -\item boundary conditions: no-slip, free-slip, half-slip -\item fancy (multi dimensional) advection schemes -\item VP vs.\ EVP \citep{hunke97} -\item ocean stress formulation \citep{hibler87} -\end{itemize} +%\subparagraph{boundary conditions: no-slip, free-slip, half-slip} +%\begin{itemize} +%\item transition from existing B-Grid to C-Grid +%\item boundary conditions: no-slip, free-slip, half-slip +%\item fancy (multi dimensional) advection schemes +%\item VP vs.\ EVP \citep{hunke97} +%\item ocean stress formulation \citep{hibler87} +%\end{itemize} \subsection{Thermodynamics} \label{sec:thermodynamics} @@ -241,9 +257,10 @@ parameterize this sub-grid scale distribution for heat flux computations, the mean ice thickness $h$ is split into seven thickness categories $H_{n}$ that are equally distributed between $2h$ and -minimum imposed ice thickness of $5\text{\,cm}$ by $H_n= -\frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat flux for all thickness -categories is averaged to give the total heat flux. +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 area averaged to give the total heat flux. \ml{[I + don't have citation for that, anyone?]} The atmospheric heat flux is balanced by an oceanic heat flux from below. The oceanic flux is proportional to @@ -264,15 +281,26 @@ algorithm following Archimedes' principle) turns snow into ice until the ice surface is back at $z=0$ \citep{leppaeranta83}. -Effective ich thickness (ice volume per unit area, +Effective ice thickness (ice volume per unit area, $c\cdot{h}$), concentration $c$ and effective snow thickness -($c\cdot{h}_{snow}$) are advected by ice velocities as described in -\refsec{dynamics}. From the various advection scheme that are -available in the MITgcm \citep{mitgcm02}, we choose flux-limited -schemes to preserve sharp gradients and edges and to rule out -unphysical over- and undershoots (negative thickness or -concentration). These scheme conserve volume and horizontal area. -\ml{[do we need to proove that? can we proove that? citation?]} +($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 +\citep{mitgcm02}, we choose flux-limited schemes +\citep[multidimensional 2nd and 3rd-order advection scheme with flux +limiter][]{roe85, 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 and horizontal area and are unconditionally +stable, so that we can set $D_{X}=0$. \ml{[do we need to proove that? + can we proove that? citation?]} There is considerable doubt about the reliability of such a simple thermodynamic model---\citet{semtner84} found significant errors in @@ -290,39 +318,39 @@ flux limiter \citep{roe85}.} -\subsection{C-grid} -\begin{itemize} -\item no-slip vs. free-slip for both lsr and evp; - "diagnostics" to look at and use for comparison - \begin{itemize} - \item ice transport through Fram Strait/Denmark Strait/Davis - Strait/Bering strait (these are general) - \item ice transport through narrow passages, e.g.\ Nares-Strait - \end{itemize} -\item compare different advection schemes (if lsr turns out to be more - effective, then with lsr otherwise I prefer evp), eg. - \begin{itemize} - \item default 2nd-order with diff1=0.002 - \item 1st-order upwind with diff1=0. - \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me) - \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.) - \end{itemize} - That should be enough. Here, total ice mass and location of ice edge - is interesting. However, this comparison can be done in an idealized - domain, may not require full Arctic Domain? -\item -Do a little study on the parameters of LSR and EVP -\begin{enumerate} -\item convergence of LSR, how many iterations do you need to get a - true elliptic yield curve -\item vary deltaTevp and the relaxation parameter for EVP and see when - the EVP solution breaks down (relative to the forcing time scale). - For this, it is essential that the evp solver gives use "stripeless" - solutions, that is your dtevp = 1sec solutions/or 10sec solutions - with SEAICE\_evpDampC = 615. -\end{enumerate} +%\subsection{C-grid} +%\begin{itemize} +%\item no-slip vs. free-slip for both lsr and evp; +% "diagnostics" to look at and use for comparison +% \begin{itemize} +% \item ice transport through Fram Strait/Denmark Strait/Davis +% Strait/Bering strait (these are general) +% \item ice transport through narrow passages, e.g.\ Nares-Strait +% \end{itemize} +%\item compare different advection schemes (if lsr turns out to be more +% effective, then with lsr otherwise I prefer evp), eg. +% \begin{itemize} +% \item default 2nd-order with diff1=0.002 +% \item 1st-order upwind with diff1=0. +% \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me) +% \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.) +% \end{itemize} +% That should be enough. Here, total ice mass and location of ice edge +% is interesting. However, this comparison can be done in an idealized +% domain, may not require full Arctic Domain? +%\item +%Do a little study on the parameters of LSR and EVP +%\begin{enumerate} +%\item convergence of LSR, how many iterations do you need to get a +% true elliptic yield curve +%\item vary deltaTevp and the relaxation parameter for EVP and see when +% the EVP solution breaks down (relative to the forcing time scale). +% For this, it is essential that the evp solver gives use "stripeless" +% solutions, that is your dtevp = 1sec solutions/or 10sec solutions +% with SEAICE\_evpDampC = 615. +%\end{enumerate} -\end{itemize} +%\end{itemize} %%% Local Variables: %%% mode: latex