articles/ceaice/ceaice.tex

% $Header: /u/gcmpack/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.12 2008/02/25 22:06:17 dimitri 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, a dynamic and thermodynamic
sea-ice model has been coupled to the Massachusetts Institute of Technology
general circulation model (MITgcm).  Ice mechanics follow a viscous plastic
rheology and the ice momentum equations are solved numerically using either
line successive relaxation (LSR) or elastic-viscous-plastic (EVP) dynamic
models.  Ice thermodynamics are represented using either a zero-heat-capacity
formulation or a two-layer formulation that conserves enthalpy.  The model
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}

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 circumvent 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.

\section{Model}
\label{sec:model}

\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} - mg \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)$ 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}  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 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 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(\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 elastic-viscous-plastic 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 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 of a curvilinear grid
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
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}.

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.


\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{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
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