--- MITgcm_contrib/articles/ceaice/ceaice.tex	2008/02/26 19:14:36	1.16
+++ MITgcm_contrib/articles/ceaice/ceaice.tex	2008/02/26 19:27:26	1.17
@@ -1,4 +1,4 @@
-% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.16 2008/02/26 19:14:36 mlosch Exp $
+% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.17 2008/02/26 19:27:26 dimitri Exp $
 % $Name:  $
 \documentclass[12pt]{article}
 
@@ -51,701 +51,17 @@
 
 \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}
-
-The availability of an adjoint model as a powerful research tool
-complementary to an ocean model was a major design requirement early
-on in the development of the MIT general circulation model (MITgcm)
-[Marshall et al. 1997a, Marotzke et al. 1999, Adcroft et al. 2002]. It
-was recognized that the adjoint model permitted computing the
-gradients of various scalar-valued model diagnostics, norms or,
-generally, objective functions with respect to external or independent
-parameters very efficiently. The information associtated with these
-gradients is useful in at least two major contexts. First, for state
-estimation problems, the objective function is the sum of squared
-differences between observations and model results weighted by the
-inverse error covariances. The gradient of such an objective function
-can be used to reduce this measure of model-data misfit to find the
-optimal model solution in a least-squares sense.  Second, the
-objective function can be a key oceanographic quantity such as
-meridional heat or volume transport, ocean heat content or mean
-surface temperature index. In this case the gradient provides a
-complete set of sensitivities of this quantity to all independent
-variables simultaneously. These sensitivities can be used to address
-the cause of, say, changing net transports accurately.
-
-References to existing sea-ice adjoint models, explaining that they are either
-for simplified configurations, for ice-only studies, or for short-duration
-studies to motivate the present work.
-
-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. Smoothing implicitly
-associated with this interpolation may mask grid scale noise and may
-contribute to stabilizing the solution. On the other hand, by
-smoothing the stress signals are damped which could lead to reduced
-variability of the system. 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),
-and models have used topographies artificially widened straits to
-avoid this problem \citep{holloway07}. The C-grid formulation on the
-other hand allows a flux of sea-ice through narrow passages if
-free-slip along the boundaries is allowed. We demonstrate this effect
-in the Candian archipelago.
-
-Talk about problems that make the sea-ice-ocean code very sensitive and
-changes in the code that reduce these sensitivities.
-
-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.
-
-\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} - m \nabla{\phi(0)} + \vek{F},
-\end{equation}
-where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area;
-$\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector;
-$\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$
-directions, respectively;
-$f$ is the Coriolis parameter;
-$\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses,
-respectively;
-$g$ is the gravity accelation;
-$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;
-$\phi(0) = g\eta + p_{a}/\rho_{0}$ is the sea surface height potential
-in response to ocean dynamics ($g\eta$) and to atmospheric pressure
-loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density);
-and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress
-tensor $\sigma_{ij}$.
-When using the rescaled vertical coordinate system, z$^\ast$, of
-\citet{cam08}, $\phi(0)$ also includes a term due to snow and ice
-loading, $mg/\rho_{0}$. 
-Advection of sea-ice momentum is neglected. The wind and ice-ocean stress
-terms are given by
-\begin{align*}
-  \vtau_{air}   = & \rho_{air}  C_{air}   |\vek{U}_{air}  -\vek{u}|
-                   R_{air}  (\vek{U}_{air}  -\vek{u}), \\ 
-  \vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}| 
-                   R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\ 
-\end{align*}
-where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere
-and surface currents of the ocean, respectively; $C_{air/ocean}$ are
-air and ocean drag coefficients; $\rho_{air/ocean}$ are reference
-densities; and $R_{air/ocean}$ are rotation matrices that act on the
-wind/current vectors.
-
-For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can
-be related to the ice strain rate and strength by a nonlinear
-viscous-plastic (VP) constitutive law \citep{hibler79, 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\left(\frac{\zeta}{e^2},
-  \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
-  {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2
-      +4\dot{\epsilon}_{12}^2}}\right).
-\end{equation}
-
-In the current implementation, the VP-model is integrated with the
-semi-implicit line successive over relaxation (LSOR)-solver of
-\citet{zhang98}, which allows for long time steps that, in our case,
-are limited by the explicit treatment of the Coriolis term. The
-explicit treatment of the Coriolis term does not represent a severe
-limitation because it restricts the time step to approximately the
-same length as in the ocean model where the Coriolis term is also
-treated explicitly.
-
-\citet{hunke97}'s introduced an elastic contribution to the strain
-rate in order to regularize Eq.\refeq{vpequation} in such a way that
-the resulting elastic-viscous-plastic (EVP) and VP models are
-identical at steady state,
-\begin{equation}
-  \label{eq:evpequation}
-  \frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} +
-  \frac{1}{2\eta}\sigma_{ij} 
-  + \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij}  
-  + \frac{P}{4\zeta}\delta_{ij}
-  = \dot{\epsilon}_{ij}. 
-\end{equation}
-%In the EVP model, equations for the components of the stress tensor
-%$\sigma_{ij}$ are solved explicitly. Both model formulations will be
-%used and compared the present sea-ice model study.
-The EVP-model uses an explicit time stepping scheme with a short
-timestep. According to the recommendation of \citet{hunke97}, the
-EVP-model is stepped forward in time 120 times within the physical
-ocean model time step (although this parameter is under debate), to
-allow for elastic waves to disappear.  Because the scheme does not
-require a matrix inversion it is fast in spite of the small 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 \ml{[modify following
-  section3:] cubed sphere grids used in this paper; in particular,
-only near the edges of the cubed sphere grid, we expect them to be
-non-zero, but these edges are at approximately 35\degS\ or 35\degN\ 
-which are never covered by sea-ice in our simulations.  Everywhere
-else the coordinate system is locally nearly cartesian.}  However, for
-last-glacial-maximum or snowball-earth-like simulations the question
-of metric terms needs to be reconsidered.  Either, one includes these
-terms as in \citet{zhang03}, or one finds a vector-invariant
-formulation for the sea-ice internal stress term that does not require
-any metric terms, as it is done in the ocean dynamics of the MITgcm
-\citep{adcroft04:_cubed_sphere}.
-
-Lateral boundary conditions are naturally ``no-slip'' for B-grids, as
-the tangential velocities points lie on the boundary. For C-grids, the
-lateral boundary condition for tangential velocities is realized via
-``ghost points'', allowing alternatively no-slip or free-slip
-conditions. In ocean models free-slip boundary conditions in
-conjunction with piecewise-constant (``castellated'') coastlines have
-been shown to reduce in effect to no-slip boundary conditions
-\citep{adcroft98:_slippery_coast}; for sea-ice models the effects of
-lateral boundary conditions have not yet been studied.
-
-Moving sea ice exerts a stress on the ocean which is the opposite of
-the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is
-applied directly to the surface layer of the ocean model. An
-alternative ocean stress formulation is given by \citet{hibler87}.
-Rather than applying $\vtau_{ocean}$ directly, the stress is derived
-from integrating over the ice thickness to the bottom of the oceanic
-surface layer. In the resulting equation for the \emph{combined}
-ocean-ice momentum, the interfacial stress cancels and the total
-stress appears as the sum of windstress and divergence of internal ice
-stresses: $\delta(z) (\vtau_{air} + \vek{F})/\rho_0$, \citep[see also
-Eq.\,2 of][]{hibler87}. The disadvantage of this formulation is that
-now the velocity in the surface layer of the ocean that is used to
-advect tracers, is really an average over the ocean surface
-velocity and the ice velocity leading to an inconsistency as the ice
-temperature and salinity are different from the oceanic variables.
-
-Sea ice distributions are characterized by sharp gradients and edges.
-For this reason, we employ positive, multidimensional 2nd and 3rd-order
-advection scheme with flux limiter \citep{roe85, 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.
-\ml{[Jean-Michel, your turf: ]Care must be taken in advecting these
-  quantities in order to ensure conservation of enthalphy. Currently
-  this can only be accomplished with a 2nd-order advection scheme with
-  flux limiter \citep{roe85}.}
-
-
-\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?
+\input{ceaice_abstract.tex}
+
+\input{ceaice_intro.tex}
+
+\input{ceaice_model.tex}
+
+\input{ceaice_forward.tex}
+
+\input{ceaice_adjoint.tex}
+
+\input{ceaice_concl.tex}
 
 \paragraph{Acknowledgements}
 We thank Jinlun Zhang for providing the original B-grid code and many