--- MITgcm_contrib/articles/ceaice/ceaice.tex	2008/02/25 23:45:46	1.13
+++ MITgcm_contrib/articles/ceaice/ceaice.tex	2008/02/26 17:21:48	1.15
@@ -1,4 +1,4 @@
-% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.13 2008/02/25 23:45:46 dimitri Exp $
+% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.15 2008/02/26 17:21:48 mlosch Exp $
 % $Name:  $
 \documentclass[12pt]{article}
 
@@ -52,7 +52,6 @@
 \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
@@ -76,10 +75,35 @@
 \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 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
@@ -95,9 +119,18 @@
 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}
 
@@ -108,7 +141,7 @@
 \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},
+  \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;
@@ -119,12 +152,14 @@
 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;
+$\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$.
+\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*}
@@ -139,9 +174,9 @@
 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}:
+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} 
@@ -189,26 +224,25 @@
 is capped to suppress any tensile stress \citep{hibler97, geiger98}:
 \begin{equation}
   \label{eq:etatem}
-  \eta = \min(\frac{\zeta}{e^2}
+  \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}}
+      +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,
-is limited by the explicit treatment of the Coriolis term. The
+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 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,
+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}} +
@@ -255,7 +289,7 @@
 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
+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
@@ -263,19 +297,30 @@
 $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}.
+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
@@ -294,9 +339,9 @@
 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.
+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}
 
@@ -333,6 +378,10 @@
 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}