--- MITgcm_contrib/articles/ceaice/ceaice_model.tex 2008/02/26 19:27:26 1.1 +++ MITgcm_contrib/articles/ceaice/ceaice_model.tex 2008/04/29 14:03:31 1.8 @@ -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 parameterizaion 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 @@ -205,20 +226,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} @@ -227,13 +242,64 @@ uses simple thermodynamics following the appendix of \citet{semtner76}. This formulation does not allow storage of heat (heat capacity of ice is zero, and this type of model is often refered -to as a ``zero-layer'' model). Upward heat flux is parameterized +to as a ``zero-layer'' model). Upward conductive heat flux is parameterized assuming a linear temperature profile and together with a constant ice conductivity. It is expressed as $(K/h)(T_{w}-T_{0})$, where $K$ is the ice conductivity, $h$ the ice thickness, and $T_{w}-T_{0}$ the difference between water and ice surface temperatures. The surface -heat budget is computed in a similar way to that of -\citet{parkinson79} and \citet{manabe79}. +heat flux is computed in a similar way to that of \citet{parkinson79} +and \citet{manabe79}. + +The conductive heat flux depends strongly on the ice thickness $h$. +However, the ice thickness in the model represents a mean over a +potentially very heterogeneous thickness distribution. In order to +parameterize 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 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 +$\rho\,c_{p}\left(T_{w}-T_{fr}\right)$ where $\rho$ and $c_{p}$ are +the density and heat capacity of sea water and $T_{fr}$ is the local +freezing point temperature that is a function of salinity. Contrary to +\citet{menemenlis05}, this flux is not assumed to instantaneously melt +or create ice, but a time scale of three days is used to relax $T_{w}$ +to the freezing point. + +The parameterization of lateral and vertical growth of sea ice follows +that of \citet{hibler79, hibler80}. + +On top of the ice there is a layer of snow that modifies the heat flux +and the albedo \citep{zhang98}. If enough snow accumulates so that its +weight submerges the ice and the snow is flooded, a simple mass +conserving parameterization of snowice formation (a flood-freeze +algorithm following Archimedes' principle) turns snow into ice until +the ice surface is back at $z=0$ \citep{leppaeranta83}. + +Effective ice thickness (ice volume per unit area, +$c\cdot{h}$), concentration $c$ and effective snow thickness +($c\cdot{h}_{s}$) are advected by ice velocities: +\begin{equation} + \label{eq:advection} + \frac{\partial{X}}{\partial{t}} = - \nabla\cdot\left(\vek{u}\,X\right) + + \Gamma_{X} + D_{X} +\end{equation} +where $\Gamma_X$ are the thermodynamic source terms and $D_{X}$ the +diffusive terms for quantities $X=(c\cdot{h}), c, (c\cdot{h}_{s})$. +% +From the various advection scheme that are available in the MITgcm +\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 @@ -282,4 +348,10 @@ solutions, that is your dtevp = 1sec solutions/or 10sec solutions with SEAICE\_evpDampC = 615. \end{enumerate} + \end{itemize} + +%%% Local Variables: +%%% mode: latex +%%% TeX-master: "ceaice" +%%% End: