articles/ceaice/ceaice.tex

% $Header: /u/gcmpack/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.9 2008/01/21 08:06:00 mlosch Exp $
% $Name:  $
\documentclass[12pt]{article}

\usepackage[]{graphicx}
\usepackage{subfigure}

\usepackage[round,comma]{natbib}
\bibliographystyle{bib/agu04}

\usepackage{amsmath,amssymb}
\newcommand\bmmax{10} \newcommand\hmmax{10}
\usepackage{bm}

% math abbreviations
\newcommand{\vek}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\vtau}{\bm{{\tau}}}

\newcommand{\degree}{\ensuremath{^\circ}}
\newcommand{\degC}{\,\ensuremath{\degree}C}
\newcommand{\degE}{\ensuremath{\degree}\,E}
\newcommand{\degS}{\ensuremath{\degree}\,S}
\newcommand{\degN}{\ensuremath{\degree}\,N}
\newcommand{\degW}{\ensuremath{\degree}\,W}

% cross reference scheme
\newcommand{\reffig}[1]{Figure~\ref{fig:#1}}
\newcommand{\reftab}[1]{Table~\ref{tab:#1}}
\newcommand{\refapp}[1]{Appendix~\ref{app:#1}}
\newcommand{\refsec}[1]{Section~\ref{sec:#1}}
\newcommand{\refeq}[1]{\,(\ref{eq:#1})}
\newcommand{\refeqs}[2]{\,(\ref{eq:#1})--(\ref{eq:#2})}

\newlength{\stdfigwidth}\setlength{\stdfigwidth}{20pc}
%\newlength{\stdfigwidth}\setlength{\stdfigwidth}{\columnwidth}
\newlength{\mediumfigwidth}\setlength{\mediumfigwidth}{39pc}
%\newlength{\widefigwidth}\setlength{\widefigwidth}{39pc}
\newlength{\widefigwidth}\setlength{\widefigwidth}{\textwidth}
\newcommand{\fpath}{figs}

% commenting scheme
\newcommand{\ml}[1]{\textsf{\slshape #1}}

\title{A Dynamic-Thermodynamic Sea ice Model for Ocean Climate
  Estimation on an Arakawa C-Grid}

\author{Martin Losch, Dimitris Menemenlis, Patrick Heimbach, \\
        Jean-Michel Campin, and Chris Hill}
\begin{document}

\maketitle

\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
sea-ice model components have been incorporated in the Massachusetts Institute
of Technology general circulation model (MITgcm).  Sea-ice dynamics use either
a visco-plastic rheology solved with a line successive relaxation (LSR)
technique, reformulated on an Arakawa C-grid in order to match the oceanic and
atmospheric grids of the MITgcm, and modified to permit efficient and accurate
automatic differentiation of the coupled ocean and sea-ice model
configurations.

\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,
  kreyscher00, zhang98, hunke97}. From the perspective of coupling a
sea ice-model to a C-grid ocean model, the exchange of fluxes of heat
and fresh-water pose no difficulty for a B-grid sea-ice model
\citep[e.g.,][]{timmermann02a}. However, surface stress is defined at
velocities points and thus needs to be interpolated between a B-grid
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
altogether and render the stress coupling as consistent as the
buoyancy coupling.

A further advantage of the C-grid formulation is apparent in narrow
straits. In the limit of only one grid cell between coasts there is no
flux allowed for a B-grid (with no-slip lateral boundary counditions),
whereas the C-grid formulation allows a flux of sea-ice through this
passage for all types of lateral boundary conditions. We (will)
demonstrate this effect in the Candian archipelago.

\subsection{Dynamics}
\label{sec:dynamics}

The momentum equations of the sea-ice model are standard with
\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 $\vek{u} = u\vek{i}+v\vek{j}$ is the ice velocity vectory, $m$
the ice mass per unit area, $f$ the Coriolis parameter, $g$ is the
gravity accelation, $\nabla\phi$ is the gradient (tilt) of the sea
surface height potential beneath the ice. $\phi$ is the sum of
atmpheric pressure $p_{a}$ and loading due to ice and snow
$(m_{i}+m_{s})g$. $\vtau_{air}$ and $\vtau_{ocean}$ are the wind and
ice-ocean stresses, respectively.  $\vek{F}$ is the interaction force
and $\vek{i}$, $\vek{j}$, and $\vek{k}$ are the unit vectors in the
$x$, $y$, and $z$ directions.  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_{ocean} =& \rho_{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}$ reference
densities, and $R_{air/ocean}$ rotation matrices that act on the
wind/current vectors. $\vek{F} = \nabla\cdot\sigma$ is the divergence
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)
constitutive law \citep{hibler79, zhang98}:
\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 compuation 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(\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}}
\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,
is 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 elatic-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
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 timestep
\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 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 a curvilinear grid
on the sphere. However, for now 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
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 fo 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}.

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.

Sea ice distributions are characterized by sharp gradients and edges.
For this reason, we employ a positive 3rd-order advection scheme
\citep{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}

\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 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}.

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.

\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,
1000\,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
25\,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}
\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}

\section{Forward sensitivity experiments}
\label{sec:forward}

A second series of forward sensitivity experiments have been carried out on an
Arctic Ocean domain with open boundaries.  Once again the objective is to
compare the old B-grid LSR dynamic solver with the new C-grid LSR and EVP
solvers.  One additional experiment is carried out to illustrate the
differences between the two main options for sea ice thermodynamics in the MITgcm.

\subsection{Arctic Domain with Open Boundaries}
\label{sec:arctic}

The Arctic domain of integration is illustrated in Fig.~\ref{???}.  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.

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
data (ETOPO2) and the model employs the partial-cell formulation of
\citet{adcroft97:_shaved_cells}, which permits accurate representation of the bathymetry. The
model is integrated in a volume-conserving configuration using a finite volume
discretization with C-grid staggering of the prognostic variables. In the
ocean, the non-linear equation of state of \citet{jackett95}.  The ocean model is
coupled to a sea-ice model described hereinabove.  

This particular ECCO2 simulation is initialized from rest using the
January temperature and salinity distribution from the World Ocean
Atlas 2001 (WOA01) [Conkright et al., 2002] and it is integrated for
32 years prior to the 1996--2001 period discussed in the study. Surface
boundary conditions are from the National Centers for Environmental
Prediction and the National Center for Atmospheric Research
(NCEP/NCAR) atmospheric reanalysis [Kistler et al., 2001]. Six-hourly
surface winds, temperature, humidity, downward short- and long-wave
radiations, and precipitation are converted to heat, freshwater, and
wind stress fluxes using the \citet{large81, large82} bulk formulae.
Shortwave radiation decays exponentially as per Paulson and Simpson
[1977]. Additionally the time-mean river run-off from Large and Nurser
[2001] is applied and there is a relaxation to the monthly-mean
climatological sea surface salinity values from WOA01 with a
relaxation time scale of 3 months. Vertical mixing follows
\citet{large94} with background vertical diffusivity of
$1.5\times10^{-5}\text{\,m$^{2}$\,s$^{-1}$}$ and viscosity of
$10^{-3}\text{\,m$^{2}$\,s$^{-1}$}$. A third order, direct-space-time
advection scheme with flux limiter is employed \citep{hundsdorfer94}
and there is no explicit horizontal diffusivity. Horizontal viscosity
follows \citet{lei96} but
modified to sense the divergent flow as per Fox-Kemper and Menemenlis
[in press].  Shortwave radiation decays exponentially as per Paulson
and Simpson [1977].  Additionally, the time-mean runoff of Large and
Nurser [2001] is applied near the coastline and, where there is open
water, there is a relaxation to monthly-mean WOA01 sea surface
salinity with a time constant of 45 days.

Open water, dry
ice, wet ice, dry snow, and wet snow albedo are, respectively, 0.15, 0.85,
0.76, 0.94, and 0.8.

\begin{itemize}
\item Configuration
\item OBCS from cube
\item forcing
\item 1/2 and full resolution
\item with a few JFM figs from C-grid LSR no slip
  ice transport through Canadian Archipelago
  thickness distribution
  ice velocity and transport
\end{itemize}

\begin{itemize}
\item Arctic configuration
\item ice transport through straits and near boundaries
\item focus on narrow straits in the Canadian Archipelago
\end{itemize}

\begin{itemize}
\item B-grid LSR no-slip
\item C-grid LSR no-slip
\item C-grid LSR slip
\item C-grid EVP no-slip
\item C-grid EVP slip
\item C-grid LSR + TEM (truncated ellipse method, no tensile stress, new flag)
\item C-grid LSR no-slip + Winton
\item  speed-performance-accuracy (small)
  ice transport through Canadian Archipelago differences
  thickness distribution differences
  ice velocity and transport differences
\end{itemize}

We anticipate small differences between the different models due to:
\begin{itemize}
\item advection schemes: along the ice-edge and regions with large
  gradients
\item C-grid: less transport through narrow straits for no slip
  conditons, more for free slip
\item VP vs.\ EVP: speed performance, accuracy?
\item ocean stress: different water mass properties beneath the ice
\end{itemize}

\section{Adjoint sensiivities of the MITsim}

\subsection{The adjoint of MITsim}

The ability to generate tangent linear and adjoint model components
of the MITsim has been a main design task.
For the ocean the adjoint capability has proven to be an
invaluable tool for sensitivity analysis as well as state estimation.
In short, the adjoint enables very efficient computation of the gradient
of scalar-valued model diagnostics (called cost function or objective function)
with respect to many model "variables".
These variables can be two- or three-dimensional fields of initial 
conditions, model parameters such as mixing coefficients, or
time-varying surface or lateral (open) boundary conditions.
When combined, these variables span a potentially high-dimensional
(e.g. O(10$^8$)) so-called control space. Performing parameter perturbations
to assess model sensitivities quickly becomes prohibitive at these scales.
Alternatively, (time-varying) sensitivities of the objective function 
to any element of the  control space can be computed very efficiently in 
one single adjoint 
model integration, provided an efficient adjoint model is available.

[REFERENCES]


The adjoint operator (ADM) is the transpose of the tangent linear operator (TLM)
of the full (in general nonlinear) forward model, i.e. the MITsim.
The TLM maps perturbations of elements of the control space
(e.g. initial ice thickness distribution)
via the model Jacobian
to a perturbation in the objective function 
(e.g. sea-ice export at the end of the integration interval).
\textit{Tangent} linearity ensures that the derivatives are evaluated
with respect to the underlying model trajectory at each point in time.
This is crucial for nonlinear trajectories and the presence of different
regimes (e.g. effect of the seaice growth term at or away from the
freezing point of the ocean surface).
Ensuring tangent linearity can be easily achieved by integrating
the full model in sync with the TLM to provide the underlying model state.
Ensuring \textit{tangent} adjoints is equally crucial, but much more
difficult to achieve because of the reverse nature of the integration:
the adjoint accumulates sensitivities backward in time,
starting from a unit perturbation of the objective function.
The adjoint model requires the model state in reverse order.
This presents one of the major complications in deriving an
exact, i.e. \textit{tangent} adjoint model.

Following closely the development and maintenance of TLM and ADM
components of the MITgcm we have relied heavily on the
autmomatic differentiation (AD) tool
"Transformation of Algorithms in Fortran" (TAF)
developed by Fastopt (Giering and Kaminski, 1998)
to derive TLM and ADM code of the MITsim.
Briefly, the nonlinear parent model is fed to the AD tool which produces 
derivative code for the specified control space and objective function.
Following this approach has (apart from its evident success)
several advantages:
(1) the adjoint model is the exact adjoint operator of the parent model,
(2) the adjoint model can be kept up to date with respect to ongoing
development of the parent model, and adjustments to the parent model
to extend the automatically generated adjoint are incremental changes
only, rather than extensive re-developments,
(3) the parallel structure of the parent model is preserved 
by the adjoint model, ensuring efficient use in high performance
computing environments.

Some initial code adjustments are required to support dependency analysis
of the flow reversal and certain language limitations which may lead
to irreducible flow graphs (e.g. GOTO statements).
The problem of providing the required model state in reverse order
at the time of evaluating nonlinear or conditional
derivatives is solved via balancing
storing vs. recomputation of the model state in a multi-level
checkpointing loop.
Again, an initial code adjustment is required to support TAFs 
checkpointing capability.
The code adjustments are sufficiently simple so as not to cause
major limitations to the full nonlinear parent model.
Once in place, an adjoint model of a new model configuration
may be derived in about 10 minutes.

[HIGHLIGHT COUPLED NATURE OF THE ADJOINT!]

\subsection{Special considerations}

* growth term(?)

* small active denominators

* dynamic solver (implicit function theorem)

* approximate adjoints


\subsection{An example: sensitivities of sea-ice export through Fram Strait}

We demonstrate the power of the adjoint method
in the context of investigating sea-ice export sensitivities through Fram Strait
(for details of this study see Heimbach et al., 2007).
%\citep[for details of this study see][]{heimbach07}. %Heimbach et al., 2007).
The domain chosen is a coarsened version of the Arctic face of the
high-resolution cubed-sphere configuration of the ECCO2 project
\citep[see][]{menemenlis05}. It covers the entire Arctic,
extends into the North Pacific such as to cover the entire
ice-covered regions, and comprises parts of the North Atlantic
down to XXN to enable analysis of remote influences of the 
North Atlantic current to sea-ice variability and export.
The horizontal resolution varies between XX and YY km
with 50 unevenly spaced vertical levels.
The adjoint models run efficiently on 80 processors
(benchmarks have been performed both on an SGI Altix as well as an 
IBM SP5 at NASA/ARC).

Following a 1-year spinup, the model has been integrated for four
years between 1992 and 1995. It is forced using realistic 6-hourly
NCEP/NCAR atmospheric state variables. Over the open ocean these are
converted into air-sea fluxes via the bulk formulae of
\citet{large04}.  Derivation of air-sea fluxes in the presence of
sea-ice is handled by the ice model as described in \refsec{model}.
The objective function chosen is sea-ice export through Fram Strait
computed for December 1995.  The adjoint model computes sensitivities
to sea-ice export back in time from 1995 to 1992 along this
trajectory.  In principle all adjoint model variable (i.e., Lagrange
multipliers) of the coupled ocean/sea-ice model are available to
analyze the transient sensitivity behaviour of the ocean and sea-ice
state.  Over the open ocean, the adjoint of the bulk formula scheme
computes sensitivities to the time-varying atmospheric state.  Over
ice-covered parts, the sea-ice adjoint converts surface ocean
sensitivities to atmospheric sensitivities.

\reffig{4yradjheff}(a--d) depict sensitivities of sea-ice export
through Fram Strait in December 1995 to changes in sea-ice thickness
12, 24, 36, 48 months back in time. Corresponding sensitivities to
ocean surface temperature are depicted in
\reffig{4yradjthetalev1}(a--d).  The main characteristics is
consistency with expected advection of sea-ice over the relevant time
scales considered.  The general positive pattern means that an
increase in sea-ice thickness at location $(x,y)$ and time $t$ will
increase sea-ice export through Fram Strait at time $T_e$.  Largest
distances from Fram Strait indicate fastest sea-ice advection over the
time span considered.  The ice thickness sensitivities are in close
correspondence to ocean surface sentivitites, but of opposite sign.
An increase in temperature will incur ice melting, decrease in ice
thickness, and therefore decrease in sea-ice export at time $T_e$.

The picture is fundamentally different and much more complex
for sensitivities to ocean temperatures away from the surface.
\reffig{4yradjthetalev10??}(a--d) depicts ice export sensitivities to
temperatures at roughly 400 m depth.
Primary features are the effect of the heat transport of the North
Atlantic current which feeds into the West Spitsbergen current,
the circulation around Svalbard, and ...

\begin{figure}[t!]
\centerline{
\subfigure[{\footnotesize -12 months}]
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim072_cmax2.0E+02.eps}}
%\includegraphics*[width=.3\textwidth]{H_c.bin_res_100_lev1.pdf}
%
\subfigure[{\footnotesize -24 months}]
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim145_cmax2.0E+02.eps}}
}

\centerline{
\subfigure[{\footnotesize
-36 months}]
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim218_cmax2.0E+02.eps}}
%
\subfigure[{\footnotesize
-48 months}]
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim292_cmax2.0E+02.eps}}
}
\caption{Sensitivity of sea-ice export through Fram Strait in December 2005 to
sea-ice thickness at various prior times.
\label{fig:4yradjheff}}
\end{figure}


\begin{figure}[t!]
\centerline{
\subfigure[{\footnotesize -12 months}]
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim072_cmax5.0E+01.eps}}
%\includegraphics*[width=.3\textwidth]{H_c.bin_res_100_lev1.pdf}
%
\subfigure[{\footnotesize -24 months}]
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim145_cmax5.0E+01.eps}}
}

\centerline{
\subfigure[{\footnotesize
-36 months}] 
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim218_cmax5.0E+01.eps}}
%
\subfigure[{\footnotesize
-48 months}]
{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim292_cmax5.0E+01.eps}}
}
\caption{Same as \reffig{4yradjheff} but for sea surface temperature
\label{fig:4yradjthetalev1}}
\end{figure}


\section{Discussion and conclusion}
\label{sec:concl}

The story of the paper could be:
B-grid ice model + C-grid ocean model does not work properly for us,
therefore C-grid ice model  with advantages: 
\begin{enumerate}
\item stress coupling simpler (no interpolation required)
\item different boundary conditions
\item advection schemes carry over trivially from main code
\end{enumerate}
LSR/EVP solutions are similar with appropriate bcs, evp parameters as
a function of forcing time scale (when does VP solution break
down). Same for LSR solver, provided that it works (o: 
Which scheme is more efficient for the resolution/time stepping
parameters that we use here. What about other resolutions?

\paragraph{Acknowledgements}
We thank Jinlun Zhang for providing the original B-grid code and many
helpful discussions. ML thanks Elizabeth Hunke for multiple explanations.

\bibliography{bib/journal_abrvs,bib/seaice,bib/genocean,bib/maths,bib/mitgcmuv,bib/fram}

\end{document}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 


A Dynamic-Thermodynamic Sea ice Model for Ocean Climate
  Estimation on an Arakawa C-Grid

Introduction

Ice Model:
 Dynamics formulation.
  B-C, LSR, EVP, no-slip, slip
  parallellization
 Thermodynamics formulation.
  0-layer Hibler salinity + snow
  3-layer Winton

Idealized tests
 Funnel Experiments
 Downstream Island tests
  B-grid LSR no-slip
  C-grid LSR no-slip
  C-grid LSR slip
  C-grid EVP no-slip
  C-grid EVP slip

Arctic Setup
 Configuration
 OBCS from cube
 forcing
 1/2 and full resolution
 with a few JFM figs from C-grid LSR no slip
  ice transport through Canadian Archipelago
  thickness distribution
  ice velocity and transport

Arctic forward sensitivity experiments
 B-grid LSR no-slip
 C-grid LSR no-slip
 C-grid LSR slip
 C-grid EVP no-slip
 C-grid EVP slip
 C-grid LSR no-slip + Winton
  speed-performance-accuracy (small)
  ice transport through Canadian Archipelago differences
  thickness distribution differences
  ice velocity and transport differences

Adjoint sensitivity experiment on 1/2-res setup
 Sensitivity of sea ice volume flow through Fram Strait
*** Sensitivity of sea ice volume flow through Canadian Archipelago

Summary and conluding remarks