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

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

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

-revision 1.10 by dimitri,
Mon Feb 25 16:50:56 2008 UTC
+revision 1.13 by dimitri,
Mon Feb 25 23:45:46 2008 UTC
 Line 54
  \begin{abstract}
  As part of ongoing efforts to obtain a best possible synthesis of most
- available, global-scale, ocean and sea ice data, dynamic and thermodynamic
+ available, global-scale, ocean and sea ice data, a dynamic and thermodynamic
- sea-ice model components have been incorporated in the Massachusetts Institute
+ sea-ice model has been coupled to the Massachusetts Institute of Technology
- of Technology general circulation model (MITgcm).  Sea-ice dynamics use either
+ general circulation model (MITgcm).  Ice mechanics follow a viscous plastic
- a visco-plastic rheology solved with a line successive relaxation (LSR)
+ rheology and the ice momentum equations are solved numerically using either
- technique, reformulated on an Arakawa C-grid in order to match the oceanic and
+ line successive relaxation (LSR) or elastic-viscous-plastic (EVP) dynamic
- atmospheric grids of the MITgcm, and modified to permit efficient and accurate
+ models.  Ice thermodynamics are represented using either a zero-heat-capacity
- automatic differentiation of the coupled ocean and sea-ice model
+ formulation or a two-layer formulation that conserves enthalpy.  The model
- configurations.
+ includes prognostic variables for snow and for sea-ice salinity.  The above
+ sea ice model components were borrowed from current-generation climate models
+ but they were reformulated on an Arakawa C-grid in order to match the MITgcm
+ oceanic grid and they were modified in many ways to permit efficient and
+ accurate automatic differentiation.  This paper describes the MITgcm sea ice
+ model; it presents example Arctic and Antarctic results from a realistic,
+ eddy-permitting, global ocean and sea-ice configuration; it compares B-grid
+ and C-grid dynamic solvers in a regional Arctic configuration; and it presents
+ example results from coupled ocean and sea-ice adjoint-model integrations.
  \end{abstract}
  \section{Introduction}
  \label{sec:intro}
- more blabla
- \section{Model}
- \label{sec:model}
  Traditionally, probably for historical reasons and the ease of
  treating the Coriolis term, most standard sea-ice models are
  discretized on Arakawa-B-grids \citep[e.g.,][]{hibler79, harder99,
-Line 84 
 velocities points and thus needs to be i
+Line 87 
 velocities points and thus needs to be i
  sea-ice model and a C-grid ocean model. While the smoothing implicitly
  associated with this interpolation may mask grid scale noise, it may
  in two-way coupling lead to a computational mode as will be shown. By
- choosing a C-grid for the sea-ice model, we circumvene this difficulty
+ choosing a C-grid for the sea-ice model, we circumvent this difficulty
  altogether and render the stress coupling as consistent as the
  buoyancy coupling.
-Line 95 
 whereas the C-grid formulation allows a
+Line 98 
 whereas the C-grid formulation allows a
  passage for all types of lateral boundary conditions. We (will)
  demonstrate this effect in the Candian archipelago.
+ \section{Model}
+ \label{sec:model}
  \subsection{Dynamics}
  \label{sec:dynamics}
- The momentum equations of the sea-ice model are standard with
+ 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},
+   \vtau_{ocean} - mg \nabla{\phi(0)} + \vek{F},
  \end{equation}
- where $\vek{u} = u\vek{i}+v\vek{j}$ is the ice velocity vectory, $m$
+ where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area;
- the ice mass per unit area, $f$ the Coriolis parameter, $g$ is the
+ $\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector;
- gravity accelation, $\nabla\phi$ is the gradient (tilt) of the sea
+ $\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$
- surface height potential beneath the ice. $\phi$ is the sum of
+ directions, respectively;
- atmpheric pressure $p_{a}$ and loading due to ice and snow
+ $f$ is the Coriolis parameter;
- $(m_{i}+m_{s})g$. $\vtau_{air}$ and $\vtau_{ocean}$ are the wind and
+ $\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses,
- ice-ocean stresses, respectively.  $\vek{F}$ is the interaction force
+ respectively;
- and $\vek{i}$, $\vek{j}$, and $\vek{k}$ are the unit vectors in the
+ $g$ is the gravity accelation;
- $x$, $y$, and $z$ directions.  Advection of sea-ice momentum is
+ $\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;
- neglected. The wind and ice-ocean stress terms are given by
+ $\phi(0)$ is the sea surface height potential in response to ocean dynamics
+ and to atmospheric pressure loading;
+ 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$.
+ Advection of sea-ice momentum is neglected. The wind and ice-ocean stress
+ terms are given by
  \begin{align*}
-   \vtau_{air} =& \rho_{air} |\vek{U}_{air}|R_{air}(\vek{U}_{air}) \\
+   \vtau_{air}   = & \rho_{air}  C_{air}   |\vek{U}_{air}  -\vek{u}|
-   \vtau_{ocean} =& \rho_{ocean} |\vek{U}_{ocean}-\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
+ and surface currents of the ocean, respectively; $C_{air/ocean}$ are
- air and ocean drag coefficients, $\rho_{air/ocean}$ reference
+ air and ocean drag coefficients; $\rho_{air/ocean}$ are reference
- densities, and $R_{air/ocean}$ rotation matrices that act on the
+ densities; and $R_{air/ocean}$ are rotation matrices that act on the
- wind/current vectors. $\vek{F} = \nabla\cdot\sigma$ is the divergence
+ wind/current vectors.
- of the interal stress tensor $\sigma_{ij}$.
  For an isotropic system this stress tensor can be related to the ice
  strain rate and strength by a nonlinear viscous-plastic (VP)
-Line 168 
 The bulk viscosities are bounded above b
+Line 181 
 The bulk viscosities are bounded above b
  $\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 compuation the replacement pressure $P = 2\,\Delta\zeta$
+ 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.
-Line 192 
 same length as in the ocean model where
+Line 205 
 same length as in the ocean model where
  treated explicitly.
  \citet{hunke97}'s introduced an elastic contribution to the strain
- rate elatic-viscous-plastic in order to regularize
+ rate elastic-viscous-plastic in order to regularize
  Eq.\refeq{vpequation} in such a way that the resulting
- elatic-viscous-plastic (EVP) and VP models are identical at steady
+ elastic-viscous-plastic (EVP) and VP models are identical at steady
  state,
  \begin{equation}
    \label{eq:evpequation}
-Line 220 
 $\sigma_{12}$. Introducing the divergenc
+Line 233 
 $\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 abbreviations,
  the equations can be written as:
  \begin{align}
    \label{eq:evpstresstensor1}
-Line 250 
 differences and averaging is only involv
+Line 263 
 differences and averaging is only involv
  $P$ at vorticity points.
  For a general curvilinear grid, one needs in principle to take metric
- terms into account that arise in the transformation a curvilinear grid
+ terms into account that arise in the transformation of a curvilinear grid
- on the sphere. However, for now we can neglect these metric terms
+ on the sphere. For now, however, we can neglect these metric terms
  because they are very small on the 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
-Line 260 
 simulations.  Everywhere else the coordi
+Line 273 
 simulations.  Everywhere else the coordi
  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 fo the sea-ice internal stress term that
+ 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}.
-Line 321 
 addition to ice-thickness and compactnes
+Line 334 
 addition to ice-thickness and compactnes
  state variables to be advected by ice velocities, namely enthalphy of
  the two ice layers and the thickness of the overlying snow layer.
- \section{Funnel Experiments}
- \label{sec:funnel}
- For a first/detailed comparison between the different variants of the
- MIT sea ice model an idealized geometry of a periodic channel,
-\,km long and 500\,m wide on a non-rotating plane, with converging
- walls forming a symmetric funnel and a narrow strait of 40\,km width
- is used. The horizontal resolution is 5\,km throughout the domain
- making the narrow strait 8 grid points wide. The ice model is
- initialized with a complete ice cover of 50\,cm uniform thickness. The
- ice model is driven by a constant along channel eastward ocean current
- of 25\,cm/s that does not see the walls in the domain. All other
- ice-ocean-atmosphere interactions are turned off, in particular there
- is no feedback of ice dynamics on the ocean current. All thermodynamic
- processes are turned off so that ice thickness variations are only
- caused by convergent or divergent ice flow. Ice volume (effective
- thickness) and concentration are advected with a third-order scheme
- with a flux limiter \citep{hundsdorfer94} to avoid undershoots. This
- scheme is unconditionally stable and does not require additional
- diffusion. The time step used here is 1\,h.
- \reffig{funnelf0} compares the dynamic fields ice concentration $c$,
- effective thickness $h_{eff} = h\cdot{c}$, and velocities $(u,v)$ for
- five different cases at steady state (after 10\,years of integration):
- \begin{description}
- \item[B-LSRns:] LSR solver with no-slip boundary conditions on a B-grid,
- \item[C-LSRns:] LSR solver with no-slip boundary conditions on a C-grid,
- \item[C-LSRfs:] LSR solver with free-slip boundary conditions on a C-grid,
- \item[C-EVPns:] EVP solver with no-slip boundary conditions on a C-grid,
- \item[C-EVPfs:] EVP solver with free-slip boundary conditions on a C-grid,
- \end{description}
- \ml{[We have not implemented the EVP solver on a B-grid.]}
- \begin{figure*}[htbp]
-   \includegraphics[width=\widefigwidth]{\fpath/all_086280}
-   \caption{Ice concentration, effective thickness [m], and ice
-     velocities [m/s]
-     for 5 different numerical solutions.}
-   \label{fig:funnelf0}
- \end{figure*}
- At a first glance, the solutions look similar. This is encouraging as
- the details of discretization and numerics should not affect the
- solutions to first order. In all cases the ice-ocean stress pushes the
- ice cover eastwards, where it converges in the funnel. In the narrow
- channel the ice moves quickly (nearly free drift) and leaves the
- channel as narrow band.
- A close look reveals interesting differences between the B- and C-grid
- results. The zonal velocity in the narrow channel is nearly the free
- drift velocity ( = ocean velocity) of 25\,cm/s for the C-grid
- solutions, regardless of the boundary conditions, while it is just
- above 20\,cm/s for the B-grid solution. The ice accelerates to
-\,cm/s after it exits the channel. Concentrating on the solutions
- B-LSRns and C-LSRns, the ice volume (effective thickness) along the
- boundaries in the narrow channel is larger in the B-grid case although
- the ice concentration is reduces in the C-grid case. The combined
- effect leads to a larger actual ice thickness at smaller
- concentrations in the C-grid case. However, since the effective
- thickness determines the ice strength $P$ in Eq\refeq{icestrength},
- the ice strength and thus the bulk and shear viscosities are larger in
- the B-grid case leading to more horizontal friction. This circumstance
- might explain why the no-slip boundary conditions in the B-grid case
- appear to be more effective in reducing the flow within the narrow
- channel, than in the C-grid case. Further, the viscosities are also
- sensitive to details of the velocity gradients. Via $\Delta$, these
- gradients enter the viscosities in the denominator so that large
- gradients tend to reduce the viscosities. This again favors more flow
- along the boundaries in the C-grid case: larger velocities
- (\reffig{funnelf0}) on grid points that are closer to the boundary by
- a factor $\frac{1}{2}$ than in the B-grid case because of the stagger
- nature of the C-grid lead numerically to larger tangential gradients
- across the boundary; these in turn make the viscosities smaller for
- less tangential friction and allow more tangential flow along the
- boundaries.
- The above argument can also be invoked to explain the small
- differences between the free-slip and no-slip solutions on the C-grid.
- Because of the non-linearities in the ice viscosities, in particular
- along the boundaries, the no-slip boundary conditions have only a small
- impact on the solution.
- The difference between LSR and EVP solutions is largest in the
- effective thickness and meridional velocity fields. The EVP velocity
- fields appears to be a little noisy. This noise has been address by
- \citet{hunke01}. It can be dealt with by reducing EVP's internal time
- step (increasing the number of iterations along with the computational
- cost) or by regularizing the bulk and shear viscosities. We revisit
- the latter option by reproducing some of the results of
- \citet{hunke01}, namely the experiment described in her section~4, for
- our C-grid no-slip cases: in a square domain with a few islands the
- ice model is initialized with constant ice thickness and linearly
- increasing ice concentration to the east. The model dynamics are
- forced with a constant anticyclonic ocean gyre and by variable
- atmospheric wind whose mean direction is diagnonal to the north-east
- corner of the domain; ice volume and concentration are held constant
- (no thermodynamics and no advection by ice velocity).
- \reffig{hunke01} shows the ice velocity field, its divergence, and the
- bulk viscosity $\zeta$ for the cases C-LRSns and C-EVPns, and for a
- C-EVPns case, where \citet{hunke01}'s regularization has been
- implemented; compare to Fig.\,4 in \citet{hunke01}. The regularization
- contraint limits ice strength and viscosities as a function of damping
- time scale, resolution and EVP-time step, effectively allowing the
- elastic waves to damp out more quickly \citep{hunke01}.
- \begin{figure*}[htbp]
-   \includegraphics[width=\widefigwidth]{\fpath/hun12days}
-   \caption{Ice flow, divergence and bulk viscosities of three
-     experiments with \citet{hunke01}'s test case: C-LSRns (top),
-     C-EVPns (middle), and C-EVPns with damping described in
-     \citet{hunke01} (bottom).}
-   \label{fig:hunke01}
- \end{figure*}
- In the far right (``east'') side of the domain the ice concentration
- is close to one and the ice should be nearly rigid. The applied wind
- tends to push ice toward the upper right corner. Because the highly
- compact ice is confined by the boundary, it resists any further
- compression and exhibits little motion in the rigid region on the
- right hand side. The C-LSRns solution (top row) allows high
- viscosities in the rigid region suppressing nearly all flow.
- \citet{hunke01}'s regularization for the C-EVPns solution (bottom row)
- clearly suppresses the noise present in $\nabla\cdot\vek{u}$ and
- $\log_{10}\zeta$ in the
- unregularized case (middle row), at the cost of reduced viscosities.
- These reduced viscosities lead to small but finite ice velocities
- which in turn can have a strong effect on solutions in the limit of
- nearly rigid regimes (arching and blocking, not shown).
- \ml{[Say something about performance? This is tricky, as the
-   perfomance depends strongly on the configuration. A run with slowly
-   changing forcing is favorable for LSR, because then only very few
-   iterations are required for convergences while EVP uses its fixed
-   number of internal timesteps. If the forcing in changing fast, LSR
-   needs far more iterations while EVP still uses the fixed number of
-   internal timesteps. I have produces runs where for slow forcing LSR
-   is much faster than EVP and for fast forcing, LSR is much slower
-   than EVP. EVP is certainly more efficient in terms of vectorization
-   and MFLOPS on our SX8, but is that a criterion?]}
  \subsection{C-grid}
  \begin{itemize}
-Line 503 
 differences between the two main options
+Line 380 
 differences between the two main options
  \subsection{Arctic Domain with Open Boundaries}
  \label{sec:arctic}
- The Arctic domain of integration is illustrated in Fig.~\ref{???}.  It
+ The Arctic domain of integration is illustrated in Fig.~\ref{fig:arctic1}.  It
  is carved out from, and obtains open boundary conditions from, the
  global cubed-sphere configuration of the Estimating the Circulation
  and Climate of the Ocean, Phase II (ECCO2) project
  \citet{menemenlis05}.  The domain size is 420 by 384 grid boxes
  horizontally with mean horizontal grid spacing of 18 km.
+ \begin{figure}
+ %\centerline{{\includegraphics*[width=0.44\linewidth]{\fpath/arctic1.eps}}}
+ \caption{Bathymetry of Arctic Domain.\label{fig:arctic1}}
+ \end{figure}
  There are 50 vertical levels ranging in thickness from 10 m near the surface
  to approximately 450 m at a maximum model depth of 6150 m. Bathymetry is from
  the National Geophysical Data Center (NGDC) 2-minute gridded global relief

 Legend:



Removed from v.1.10
 


changed lines


 
Added in v.1.13
 Legend:



Removed from v.1.10
 


changed lines


 
Added in v.1.13
-Removed from v.1.10
+Added in v.1.13

	ViewVC Help
Powered by ViewVC 1.1.22