/[MITgcm]/MITgcm_contrib/articles/ceaice/ceaice_model.tex

Diff of /MITgcm_contrib/articles/ceaice/ceaice_model.tex

Parent Directory | Revision Log | View Revision Graph Revision Graph | View Patch Patch

-revision 1.3 by mlosch,
Thu Feb 28 16:34:56 2008 UTC
+revision 1.12 by dimitri,
Fri Aug 15 14:16:35 2008 UTC
 Line 1
- \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 state estimation, 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}, and EVP
+   \citep{hunke97};
+ \item ice-ocean stress can be formulated as in \citet{hibler87} or as in \citet{cam08};
+ \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 more stable automatic differentiation of the code.
+ \end{itemize}
+ The sea ice model is tightly coupled to the ocean compontent of the MITgcm
+ \citep{mar97a}.  Heat, fresh water fluxes and surface stresses are computed
+ from the atmospheric state and modified by the ice model at every time step.
  \subsection{Dynamics}
- \label{sec:dynamics}
+ \label{app:dynamics}
  The momentum equation of the sea-ice model is
  \begin{equation}
-Line 19 
 $\vtau_{air}$ and $\vtau_{ocean}$ are th
+Line 44 
 $\vtau_{air}$ and $\vtau_{ocean}$ are th
  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
+ $\phi(0) = g\eta + p_{a}/\rho_{0} + mg/\rho_{0}$ is the sea surface
- in response to ocean dynamics ($g\eta$) and to atmospheric pressure
+ height potential in response to ocean dynamics ($g\eta$), to
- loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density);
+ atmospheric pressure loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a
- and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress
+ reference density) and a term due to snow and ice loading \citep{campin08};
- tensor $\sigma_{ij}$.
+ and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice
- When using the rescaled vertical coordinate system, z$^\ast$, of
+ stress tensor $\sigma_{ij}$. %
- \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*}
-Line 43 
 wind/current vectors.
+Line 66 
 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, 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}
-Line 77 
 minor axis $e$ equal to $2$; they are gi
+Line 100 
 minor axis $e$ equal to $2$; they are gi
      \left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right)
      (1+e^{-2}) +  4e^{-2}\dot{\epsilon}_{12}^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
-Line 126 
 timestep. According to the recommendatio
+Line 149 
 timestep. According to the recommendatio
  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
-Line 134 
 $\sigma_{12}$. Introducing the divergenc
+Line 158 
 $\sigma_{12}$. Introducing the divergenc
  \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 =
-\dot{\epsilon}_{12}$, respectively, and using the above abbreviations,
+\dot{\epsilon}_{12}$, respectively, and using the above
- the equations can be written as:
+ abbreviations, the equations\refeq{evpequation} can be written as:
  \begin{align}
    \label{eq:evpstresstensor1}
    \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} +
-Line 153 
 $T=E_{0}\Delta{t}$ with the tunable para
+Line 177 
 $T=E_{0}\Delta{t}$ with the tunable para
  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
+ Our discretization differs from \citet{zhang98, zhang03} in the
- \citet{zhang98, zhang03}. Our discretization differs only (but
+ underlying grid, namely the Arakawa C-grid, but is otherwise
- importantly) in the underlying grid, namely the Arakawa C-grid, but is
+ straightforward. The EVP model, in particular, is discretized
- 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
-Line 205 
 advect tracers, is really an average ove
+Line 228 
 advect tracers, is really an average ove
  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.
+ %\subparagraph{boundary conditions: no-slip, free-slip, half-slip}
- For this reason, we employ positive, multidimensional 2nd and 3rd-order
+ %\begin{itemize}
- advection scheme with flux limiter \citep{roe85, hundsdorfer94} for the
+ %\item transition from existing B-Grid to C-Grid
- thermodynamic variables discussed in the next section.
+ %\item boundary conditions: no-slip, free-slip, half-slip
+ %\item fancy (multi dimensional) advection schemes
- \subparagraph{boundary conditions: no-slip, free-slip, half-slip}
+ %\item VP vs.\ EVP \citep{hunke97}
+ %\item ocean stress formulation \citep{hibler87}
- \begin{itemize}
+ %\end{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}
+ \label{app:thermodynamics}
- In the original formulation the sea ice model \citep{menemenlis05}
+ In its 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
+ \citet{semtner76}. This formulation does not allow storage of heat,
- (heat capacity of ice is zero, and this type of model is often refered
+ that is, the heat capacity of ice is zero. Upward conductive heat flux
- to as a ``zero-layer'' model). Upward conductive heat flux is parameterized
+ is parameterized assuming a linear temperature profile and together
- assuming a linear temperature profile and together with a constant ice
+ with a constant ice conductivity. It is expressed as
- conductivity. It is expressed as $(K/h)(T_{w}-T_{0})$, where $K$ is
+ $(K/h)(T_{w}-T_{0})$, where $K$ is the ice conductivity, $h$ the ice
- the ice conductivity, $h$ the ice thickness, and $T_{w}-T_{0}$ the
+ thickness, and $T_{w}-T_{0}$ the difference between water and ice
- difference between water and ice surface temperatures. The surface
+ surface temperatures. This type of model is often refered to as a
- heat flux is computed in a similar way to that of \citet{parkinson79}
+ ``zero-layer'' model. The surface heat flux is computed in a similar
- and \citet{manabe79}.
+ 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
+ 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
+ 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 flux for all thickness
+ \frac{2n-1}{7}\,h$ for $n\in[1,7]$. The heat fluxes computed for each
- categories is averaged to give the total heat flux.
+ thickness category is area-averaged to give the total heat flux
+ \citep{hibler84}.
+ \ml{[This is Ian Fenty's work and we may want to remove this paragraph
+   from the paper]: %
  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
-Line 252 
 the density and heat capacity of sea wat
+Line 272 
 the density and heat capacity of sea wat
  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.
+ 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
+ and the albedo \citep{zhang98}. Snow modifies the effective
- weight submerges the ice and the snow is flooded, a simple mass
+ conductivity according to
- conserving parameterization of snowice formation (a flood-freeze
+ \[\frac{K}{h} \rightarrow \frac{1}{\frac{h_{s}}{K_{s}}+\frac{h}{K}},\]
- algorithm following Archimedes' principle) turns snow into ice until
+ where $K_s$ is the conductivity of snow and $h_s$ the snow thickness.
- the ice surface is back at $z=0$ \citep{leppaeranta83}.
+ 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 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
+ ($c\cdot{h}_{s}$) are advected by ice velocities:
- \refsec{dynamics}. From the various advection scheme that are
+ \begin{equation}
- available in the MITgcm \citep{mitgcm02}, we choose flux-limited
+   \label{eq:advection}
- schemes to preserve sharp gradients and edges and to rule out
+   \frac{\partial{X}}{\partial{t}} = - \nabla\cdot\left(\vek{u}\,X\right) +
- unphysical over- and undershoots (negative thickness or
+   \Gamma_{X} + D_{X}
- concentration). These scheme conserve volume and horizontal area.
+ \end{equation}
- \ml{[do we need to proove that? can we proove that? citation?]}
+ 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
-Line 289 
 the two ice layers and the thickness of
+Line 324 
 the two ice layers and the thickness of
    this can only be accomplished with a 2nd-order advection scheme with
    flux limiter \citep{roe85}.}
+ Further documentation and model code can be found at
+ \url{http://mitgcm.org}.
- \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}
+ %\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}
  %%% Local Variables:
  %%% mode: latex

 Legend:



Removed from v.1.3
 


changed lines


 
Added in v.1.12
 Legend:



Removed from v.1.3
 


changed lines


 
Added in v.1.12
-Removed from v.1.3
+Added in v.1.12

	ViewVC Help
Powered by ViewVC 1.1.22