--- MITgcm_contrib/articles/ceaice/ceaice.tex	2007/11/07 14:35:09	1.1
+++ MITgcm_contrib/articles/ceaice/ceaice.tex	2008/02/26 00:13:20	1.14
@@ -1,10 +1,12 @@
+% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.14 2008/02/26 00:13:20 dimitri Exp $
+% $Name:  $
 \documentclass[12pt]{article}
-\usepackage{epsfig}
-\usepackage{graphics}
+
+\usepackage[]{graphicx}
 \usepackage{subfigure}
 
 \usepackage[round,comma]{natbib}
-\bibliographystyle{agu04}
+\bibliographystyle{bib/agu04}
 
 \usepackage{amsmath,amssymb}
 \newcommand\bmmax{10} \newcommand\hmmax{10}
@@ -35,7 +37,10 @@
 \newlength{\mediumfigwidth}\setlength{\mediumfigwidth}{39pc}
 %\newlength{\widefigwidth}\setlength{\widefigwidth}{39pc}
 \newlength{\widefigwidth}\setlength{\widefigwidth}{\textwidth}
-\newcommand{\fpath}{.}
+\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}
@@ -47,21 +52,58 @@
 \maketitle
 
 \begin{abstract}
-  Some blabla
+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}
 
-more blabla
-
-\section{Model}
-\label{sec:model}
+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 permitted very efficient computation
+of gradients of various scalar-valued model diagnostics,
+norms or, generally, objective functions with respect
+to external or independent parameters. Such gradients
+arise in at least two major contexts. If the objective function
+is the sum of squared model vs. obervation differences
+weighted by e.g. the inverse error covariances, the gradient
+of the objective function can be used to optimize this measure
+of model vs. data misfit in a least-squares sense. One
+is then solving a problem of statistical state estimation.
+If the objective function is a key oceanographic quantity
+such as meridional heat or volume transport, ocean heat
+content or mean surface temperature index, the gradient
+provides a complete set of sensitivities of this quantity
+with respect to all independent variables simultaneously.
+
+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
+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
@@ -69,7 +111,7 @@
 sea-ice model and a C-grid ocean model. While the smoothing implicitly
 associated with this interpolation may mask grid scale noise, it may
 in two-way coupling lead to a computational mode as will be shown. By
-choosing a C-grid for the sea-ice model, we circumvene this difficulty
+choosing a C-grid for the sea-ice model, we circumvent this difficulty
 altogether and render the stress coupling as consistent as the
 buoyancy coupling.
 
@@ -77,39 +119,58 @@
 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)
+passage for all types of lateral boundary conditions. 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 equations of the sea-ice model are standard with
+The momentum equation of the sea-ice model is
 \begin{equation}
   \label{eq:momseaice}
   m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} +
-  \vtau_{ocean} - m \nabla{\phi(0)} + \vek{F},
+  \vtau_{ocean} - mg \nabla{\phi(0)} + \vek{F},
 \end{equation}
-where $\vek{u} = u\vek{i}+v\vek{j}$ is the ice velocity vectory, $m$
-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
+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} |\vek{U}_{air}|R_{air}(\vek{U}_{air}) \\ 
-  \vtau_{ocean} =& \rho_{ocean} |\vek{U}_{ocean}-\vek{u}| 
+  \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}$ 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}$. 
+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)
@@ -127,17 +188,21 @@
     \frac{\partial{u_{i}}}{\partial{x_{j}}} +
     \frac{\partial{u_{j}}}{\partial{x_{i}}}\right).
 \end{equation*}
-The pressure $P$, a measure of ice strength, depends on both thickness
-$h$ and compactness (concentration) $c$: \[P =
-P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]},\] 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:
+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 =& \frac{P}{2\Delta} \\
-  \eta =& \frac{P}{2\Delta{e}^2} \\
+  \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)
@@ -145,6 +210,24 @@
     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,
@@ -155,9 +238,9 @@
 treated explicitly.
 
 \citet{hunke97}'s introduced an elastic contribution to the strain
-rate elatic-viscous-plastic in order to regularize
+rate elastic-viscous-plastic in order to regularize
 Eq.\refeq{vpequation} in such a way that the resulting
-elatic-viscous-plastic (EVP) and VP models are identical at steady
+elastic-viscous-plastic (EVP) and VP models are identical at steady
 state,
 \begin{equation}
   \label{eq:evpequation}
@@ -183,7 +266,7 @@
 \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,
+2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations,
 the equations can be written as:
 \begin{align}
   \label{eq:evpstresstensor1}
@@ -213,8 +296,8 @@
 $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
+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
@@ -223,7 +306,7 @@
 cartesian.  However, for last-glacial-maximum or snowball-earth-like
 simulations the question of metric terms needs to be reconsidered.
 Either, one includes these terms as in \citet{zhang03}, or one finds a
-vector-invariant formulation fo the sea-ice internal stress term that
+vector-invariant formulation for the sea-ice internal stress term that
 does not require any metric terms, as it is done in the ocean dynamics
 of the MITgcm \citep{adcroft04:_cubed_sphere}.
 
@@ -284,26 +367,6 @@
 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}
-
-\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
-\end{itemize}
-
-\subsection{B-grid vs.\ C-grid}
-Comparison between:
-\begin{itemize}
-\item B-grid, lsr, no-slip
-\item C-grid, lsr, no-slip 
-\item C-grid, evp, no-slip
-\end{itemize}
-all without ice-ocean stress, because ice-ocean stress does not work
-for B-grid.
 
 \subsection{C-grid}
 \begin{itemize}
@@ -350,45 +413,55 @@
 \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 \cite{men05a}.  The domain size is 420 by
-384 grid boxes horizontally with mean horizontal grid spacing of 18 km.  
+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
-\cite{adc97}, which permits accurate representation of the bathymetry. The
+\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 \cite{jac95}.  The ocean model is
+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 Large and Pond [1981, 1982] 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 Large et al. [1994] with background vertical
-diffusivity of 1.5 × 10-5 m2 s-1 and viscosity of 10-3 m2 s-1. A third order,
-direct-space-time advection scheme with flux limiter is employed and there is
-no explicit horizontal diffusivity. Horizontal viscosity follows Leith [1996]
-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.
+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,
@@ -417,6 +490,7 @@
 \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
@@ -428,18 +502,12 @@
 \begin{itemize}
 \item advection schemes: along the ice-edge and regions with large
   gradients
-\item C-grid: more transport through narrow straits for no slip
-  conditons, less for free slip
+\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 sensitivity experiment}
-\label{sec:adjoint}
-
-Adjoint sensitivity experiment on 1/2-res setup
- Sensitivity of sea ice volume flow through Fram Strait
-
 \section{Adjoint sensiivities of the MITsim}
 
 \subsection{The adjoint of MITsim}
@@ -516,7 +584,7 @@
 checkpointing loop.
 Again, an initial code adjustment is required to support TAFs 
 checkpointing capability.
-The code adjustments are sufficiently simply so as not to cause
+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.
@@ -539,9 +607,10 @@
 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
-(see Menemenlis et al. 2005). It covers the entire Arctic,
+\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 
@@ -552,46 +621,41 @@
 (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 Large and Yeager (2004).
-Derivation of air-sea fluxes in the presence of sea-ice is handled
-by the ice model as described in Section XXX.
+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.
-
-Fig. XXX(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 Fig. XXX(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$.
+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.
-Fig. XXX (a--d) depicts ice export sensitivities to
+\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,
@@ -600,21 +664,21 @@
 \begin{figure}[t!]
 \centerline{
 \subfigure[{\footnotesize -12 months}]
-{\includegraphics*[width=0.44\linewidth]{figs/run_4yr_ADJheff_arc_lev1_tim072_cmax2.0E+02.eps}}
+{\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]{figs/run_4yr_ADJheff_arc_lev1_tim145_cmax2.0E+02.eps}}
+{\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]{figs/run_4yr_ADJheff_arc_lev1_tim218_cmax2.0E+02.eps}}
+{\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]{figs/run_4yr_ADJheff_arc_lev1_tim292_cmax2.0E+02.eps}}
+{\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.
@@ -625,23 +689,23 @@
 \begin{figure}[t!]
 \centerline{
 \subfigure[{\footnotesize -12 months}]
-{\includegraphics*[width=0.44\linewidth]{figs/run_4yr_ADJtheta_arc_lev1_tim072_cmax5.0E+01.eps}}
+{\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]{figs/run_4yr_ADJtheta_arc_lev1_tim145_cmax5.0E+01.eps}}
+{\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]{figs/run_4yr_ADJtheta_arc_lev1_tim218_cmax5.0E+01.eps}}
+{\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]{figs/run_4yr_ADJtheta_arc_lev1_tim292_cmax5.0E+01.eps}}
+{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim292_cmax5.0E+01.eps}}
 }
-\caption{Same as Fig. XXX but for sea surface temperature
+\caption{Same as \reffig{4yradjheff} but for sea surface temperature
 \label{fig:4yradjthetalev1}}
 \end{figure}
 
@@ -666,7 +730,7 @@
 
 \paragraph{Acknowledgements}
 We thank Jinlun Zhang for providing the original B-grid code and many
-helpful discussions.
+helpful discussions. ML thanks Elizabeth Hunke for multiple explanations.
 
 \bibliography{bib/journal_abrvs,bib/seaice,bib/genocean,bib/maths,bib/mitgcmuv,bib/fram}
 
@@ -676,3 +740,54 @@
 %%% 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