Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch
--- MITgcm_contrib/articles/ceaice/ceaice_model.tex 2008/06/04 13:33:45 1.9
+++ MITgcm_contrib/articles/ceaice/ceaice_model.tex 2008/07/03 18:10:31 1.10
@@ -15,308 +15,19 @@
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
+\item ice variables \ml{can be} advected by sophisticated, \ml{conservative}
+ advection schemes \ml{with flux limiting};
+\item growth and melt parameterizations 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}
-
-The momentum equation of the sea-ice model is
-\begin{equation}
- \label{eq:momseaice}
- m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} +
- \vtau_{ocean} - m \nabla{\phi(0)} + \vek{F},
-\end{equation}
-where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area;
-$\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector;
-$\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$
-directions, respectively;
-$f$ is the Coriolis parameter;
-$\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses,
-respectively;
-$g$ is the gravity accelation;
-$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;
-$\phi(0) = g\eta + p_{a}/\rho_{0}$ is the sea surface height potential
-in response to ocean dynamics ($g\eta$) and to atmospheric pressure
-loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density);
-and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress
-tensor $\sigma_{ij}$.
-When using the rescaled vertical coordinate system, z$^\ast$, of
-\citet{cam08}, $\phi(0)$ also includes a term due to snow and ice
-loading, $mg/\rho_{0}$.
-Advection of sea-ice momentum is neglected. The wind and ice-ocean stress
-terms are given by
-\begin{align*}
- \vtau_{air} = & \rho_{air} C_{air} |\vek{U}_{air} -\vek{u}|
- R_{air} (\vek{U}_{air} -\vek{u}), \\
- \vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}|
- R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\
-\end{align*}
-where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere
-and surface currents of the ocean, respectively; $C_{air/ocean}$ are
-air and ocean drag coefficients; $\rho_{air/ocean}$ are reference
-densities; and $R_{air/ocean}$ are rotation matrices that act on the
-wind/current vectors.
-
-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, zhang97}:
-\begin{equation}
- \label{eq:vpequation}
- \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij}
- + \left[\zeta(\dot{\epsilon}_{ij},P) -
- \eta(\dot{\epsilon}_{ij},P)\right]\dot{\epsilon}_{kk}\delta_{ij}
- - \frac{P}{2}\delta_{ij}.
-\end{equation}
-The ice strain rate is given by
-\begin{equation*}
- \dot{\epsilon}_{ij} = \frac{1}{2}\left(
- \frac{\partial{u_{i}}}{\partial{x_{j}}} +
- \frac{\partial{u_{j}}}{\partial{x_{i}}}\right).
-\end{equation*}
-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)]},
-\label{eq:icestrength}
-\end{equation}
-with the constants $P^{*}$ and $C^{*}$. The nonlinear bulk and shear
-viscosities $\eta$ and $\zeta$ are functions of ice strain rate
-invariants and ice strength such that the principal components of the
-stress lie on an elliptical yield curve with the ratio of major to
-minor axis $e$ equal to $2$; they are given by:
-\begin{align*}
- \zeta =& \min\left(\frac{P_{\max}}{2\max(\Delta,\Delta_{\min})},
- \zeta_{\max}\right) \\
- \eta =& \frac{\zeta}{e^2} \\
- \intertext{with the abbreviation}
- \Delta = & \left[
- \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}}
-\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
-maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where
-$\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress
-tensor computation the replacement pressure $P = 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}:
-\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).
-\end{equation}
-
-In the current implementation, the VP-model is integrated with the
-semi-implicit line successive over relaxation (LSOR)-solver of
-\citet{zhang98}, 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.
-
-\citet{hunke97}'s introduced an elastic contribution to the strain
-rate in order to regularize Eq.\refeq{vpequation} in such a way that
-the resulting elastic-viscous-plastic (EVP) and VP models are
-identical at steady state,
-\begin{equation}
- \label{eq:evpequation}
- \frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} +
- \frac{1}{2\eta}\sigma_{ij}
- + \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij}
- + \frac{P}{4\zeta}\delta_{ij}
- = \dot{\epsilon}_{ij}.
-\end{equation}
-%In the EVP model, equations for the components of the stress tensor
-%$\sigma_{ij}$ are solved explicitly. Both model formulations will be
-%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{hunke97}, 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
-\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
-$\sigma_{12}$. Introducing the divergence $D_D =
-\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\refeq{evpequation} can be written as:
-\begin{align}
- \label{eq:evpstresstensor1}
- \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +
- \frac{P}{2T} &= \frac{P}{2T\Delta} D_D \\
- \label{eq:evpstresstensor2}
- \frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T}
- &= \frac{P}{2T\Delta} D_T \\
- \label{eq:evpstresstensor12}
- \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{hunke97} recommend
-$E_{0} = \frac{1}{3}$.
-
-For details of the spatial discretization, the reader is referred to
-\citet{zhang98, zhang03}. Our discretization differs only (but
-importantly) in the underlying grid, namely the Arakawa C-grid, but is
-otherwise straightforward. The EVP model, in particular, is discretized
-naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the
-center points and $\sigma_{12}$ on the corner (or vorticity) points of
-the grid. With this choice all derivatives are discretized as central
-differences and averaging is only involved in computing $\Delta$ and
-$P$ at vorticity points.
-
-For a general curvilinear grid, one needs in principle to take metric
-terms into account that arise in the transformation of a curvilinear
-grid on the sphere. For now, however, we can neglect these metric
-terms because they are very small on the \ml{[modify following
- section3:] cubed sphere grids used in this paper; in particular,
-only near the edges of the cubed sphere grid, we expect them to be
-non-zero, but these edges are at approximately 35\degS\ or 35\degN\
-which are never covered by sea-ice in our simulations. Everywhere
-else the coordinate system is locally nearly cartesian.} However, for
-last-glacial-maximum or snowball-earth-like simulations the question
-of metric terms needs to be reconsidered. Either, one includes these
-terms as in \citet{zhang03}, or one finds a vector-invariant
-formulation for the sea-ice internal stress term that does not require
-any metric terms, as it is done in the ocean dynamics of the MITgcm
-\citep{adcroft04:_cubed_sphere}.
-
-Lateral boundary conditions are naturally ``no-slip'' for B-grids, as
-the tangential velocities points lie on the boundary. For C-grids, the
-lateral boundary condition for tangential velocities is realized via
-``ghost points'', allowing alternatively no-slip or free-slip
-conditions. In ocean models free-slip boundary conditions in
-conjunction with piecewise-constant (``castellated'') coastlines have
-been shown to reduce in effect to no-slip boundary conditions
-\citep{adcroft98:_slippery_coast}; for sea-ice models the effects of
-lateral boundary conditions have not yet been studied.
-
-Moving sea ice exerts a stress on the ocean which is the opposite of
-the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is
-applied directly to the surface layer of the ocean model. An
-alternative ocean stress formulation is given by \citet{hibler87}.
-Rather than applying $\vtau_{ocean}$ directly, the stress is derived
-from integrating over the ice thickness to the bottom of the oceanic
-surface layer. In the resulting equation for the \emph{combined}
-ocean-ice momentum, the interfacial stress cancels and the total
-stress appears as the sum of windstress and divergence of internal ice
-stresses: $\delta(z) (\vtau_{air} + \vek{F})/\rho_0$, \citep[see also
-Eq.\,2 of][]{hibler87}. The disadvantage of this formulation is that
-now the velocity in the surface layer of the ocean that is used to
-advect tracers, is really an average over the ocean surface
-velocity and the ice velocity leading to an inconsistency as the ice
-temperature and salinity are different from the oceanic variables.
-
-%\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}
-
-In the original formulation the sea ice model \citep{menemenlis05}
-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 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 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
-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
-$\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
-phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in
-such models---, so that today many sea ice models employ more complex
-thermodynamics. A popular thermodynamics model of \citet{winton00} is
-based on the 3-layer model of \citet{semtner76} and treats brine
-content by means of enthalphy conservation. This model requires in
-addition to ice-thickness and compactness (fractional area) additional
-state variables to be advected by ice velocities, namely enthalphy of
-the two ice layers and the thickness of the overlying snow layer.
-\ml{[Jean-Michel, your turf: ]Care must be taken in advecting these
- quantities in order to ensure conservation of enthalphy. Currently
- this can only be accomplished with a 2nd-order advection scheme with
- flux limiter \citep{roe85}.}
+The sea ice model is tightly coupled to the ocean compontent of the
+MITgcm \citep{marshall97:_finit_volum_incom_navier_stokes, mitgcm02}.
+Heat, fresh water fluxes and surface stresses are computed from the
+atmospheric state and modified by the ice model at every time step.
+The model equations and details of their numerical realization are summarized
+in the appendix. Further documentation and model code can be found at
+\url{http://mitgcm.org}.
%\subsection{C-grid}
%\begin{itemize}
| ViewVC Help | |
| Powered by ViewVC 1.1.22 |