--- MITgcm_contrib/articles/ceaice/ceaice.tex	2008/01/21 08:06:00	1.9
+++ MITgcm_contrib/articles/ceaice/ceaice.tex	2008/02/26 19:14:36	1.16
@@ -1,4 +1,4 @@
-% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.9 2008/01/21 08:06:00 mlosch Exp $
+% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.16 2008/02/26 19:14:36 mlosch Exp $
 % $Name:  $
 \documentclass[12pt]{article}
 
@@ -52,16 +52,53 @@
 \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 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
@@ -71,54 +108,78 @@
 and fresh-water pose no difficulty for a B-grid sea-ice model
 \citep[e.g.,][]{timmermann02a}. However, surface stress is defined at
 velocities points and thus needs to be interpolated between a B-grid
-sea-ice model and a C-grid ocean model. While the smoothing implicitly
-associated with this interpolation may mask grid scale noise, it may
-in two-way coupling lead to a computational mode as will be shown. By
-choosing a C-grid for the sea-ice model, we circumvene this difficulty
-altogether and render the stress coupling as consistent as the
-buoyancy coupling.
+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),
-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.
+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 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},
 \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) = 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} |\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}$. 
-
-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}:
+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} 
@@ -158,7 +219,7 @@
 $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a
 maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where
 $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress
-tensor compuation the replacement pressure $P = 2\,\Delta\zeta$
+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.
 
@@ -166,26 +227,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 elatic-viscous-plastic in order to regularize
-Eq.\refeq{vpequation} in such a way that the resulting
-elatic-viscous-plastic (EVP) and VP models are identical at steady
-state,
+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}} +
@@ -210,7 +270,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}
@@ -232,7 +292,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
@@ -240,19 +300,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 a curvilinear grid
-on the sphere. However, for now we can neglect these metric terms
-because they are very small on the cubed sphere grids used in this
-paper; in particular, only near the edges of the cubed sphere grid, we
-expect them to be non-zero, but these edges are at approximately
-35\degS\ or 35\degN\ which are never covered by sea-ice in our
-simulations.  Everywhere else the coordinate system is locally nearly
-cartesian.  However, for last-glacial-maximum or snowball-earth-like
-simulations the question of metric terms needs to be reconsidered.
-Either, one includes these terms as in \citet{zhang03}, or one finds a
-vector-invariant formulation fo the sea-ice internal stress term that
-does not require any metric terms, as it is done in the ocean dynamics
-of the MITgcm \citep{adcroft04:_cubed_sphere}.
+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
@@ -271,9 +342,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}
 
@@ -310,143 +381,11 @@
 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}.}
 
-\section{Funnel Experiments}
-\label{sec:funnel}
-
-For a first/detailed comparison between the different variants of the
-MIT sea ice model an idealized geometry of a periodic channel,
-1000\,km long and 500\,m wide on a non-rotating plane, with converging
-walls forming a symmetric funnel and a narrow strait of 40\,km width
-is used. The horizontal resolution is 5\,km throughout the domain
-making the narrow strait 8 grid points wide. The ice model is
-initialized with a complete ice cover of 50\,cm uniform thickness. The
-ice model is driven by a constant along channel eastward ocean current
-of 25\,cm/s that does not see the walls in the domain. All other
-ice-ocean-atmosphere interactions are turned off, in particular there
-is no feedback of ice dynamics on the ocean current. All thermodynamic
-processes are turned off so that ice thickness variations are only
-caused by convergent or divergent ice flow. Ice volume (effective
-thickness) and concentration are advected with a third-order scheme
-with a flux limiter \citep{hundsdorfer94} to avoid undershoots. This
-scheme is unconditionally stable and does not require additional
-diffusion. The time step used here is 1\,h.
-
-\reffig{funnelf0} compares the dynamic fields ice concentration $c$,
-effective thickness $h_{eff} = h\cdot{c}$, and velocities $(u,v)$ for
-five different cases at steady state (after 10\,years of integration): 
-\begin{description}
-\item[B-LSRns:] LSR solver with no-slip boundary conditions on a B-grid, 
-\item[C-LSRns:] LSR solver with no-slip boundary conditions on a C-grid,
-\item[C-LSRfs:] LSR solver with free-slip boundary conditions on a C-grid,
-\item[C-EVPns:] EVP solver with no-slip boundary conditions on a C-grid,
-\item[C-EVPfs:] EVP solver with free-slip boundary conditions on a C-grid,
-\end{description}
-\ml{[We have not implemented the EVP solver on a B-grid.]}
-\begin{figure*}[htbp]
-  \includegraphics[width=\widefigwidth]{\fpath/all_086280}
-  \caption{Ice concentration, effective thickness [m], and ice
-    velocities [m/s]
-    for 5 different numerical solutions.} 
-  \label{fig:funnelf0}
-\end{figure*}
-At a first glance, the solutions look similar. This is encouraging as
-the details of discretization and numerics should not affect the
-solutions to first order. In all cases the ice-ocean stress pushes the
-ice cover eastwards, where it converges in the funnel. In the narrow
-channel the ice moves quickly (nearly free drift) and leaves the
-channel as narrow band.
-
-A close look reveals interesting differences between the B- and C-grid
-results. The zonal velocity in the narrow channel is nearly the free
-drift velocity ( = ocean velocity) of 25\,cm/s for the C-grid
-solutions, regardless of the boundary conditions, while it is just
-above 20\,cm/s for the B-grid solution. The ice accelerates to
-25\,cm/s after it exits the channel. Concentrating on the solutions
-B-LSRns and C-LSRns, the ice volume (effective thickness) along the
-boundaries in the narrow channel is larger in the B-grid case although
-the ice concentration is reduces in the C-grid case. The combined
-effect leads to a larger actual ice thickness at smaller
-concentrations in the C-grid case. However, since the effective
-thickness determines the ice strength $P$ in Eq\refeq{icestrength},
-the ice strength and thus the bulk and shear viscosities are larger in
-the B-grid case leading to more horizontal friction. This circumstance
-might explain why the no-slip boundary conditions in the B-grid case
-appear to be more effective in reducing the flow within the narrow
-channel, than in the C-grid case. Further, the viscosities are also
-sensitive to details of the velocity gradients. Via $\Delta$, these
-gradients enter the viscosities in the denominator so that large
-gradients tend to reduce the viscosities. This again favors more flow
-along the boundaries in the C-grid case: larger velocities
-(\reffig{funnelf0}) on grid points that are closer to the boundary by
-a factor $\frac{1}{2}$ than in the B-grid case because of the stagger
-nature of the C-grid lead numerically to larger tangential gradients
-across the boundary; these in turn make the viscosities smaller for
-less tangential friction and allow more tangential flow along the
-boundaries.
-
-The above argument can also be invoked to explain the small
-differences between the free-slip and no-slip solutions on the C-grid.
-Because of the non-linearities in the ice viscosities, in particular
-along the boundaries, the no-slip boundary conditions have only a small
-impact on the solution.
-
-The difference between LSR and EVP solutions is largest in the
-effective thickness and meridional velocity fields. The EVP velocity
-fields appears to be a little noisy. This noise has been address by
-\citet{hunke01}. It can be dealt with by reducing EVP's internal time
-step (increasing the number of iterations along with the computational
-cost) or by regularizing the bulk and shear viscosities. We revisit
-the latter option by reproducing some of the results of
-\citet{hunke01}, namely the experiment described in her section~4, for
-our C-grid no-slip cases: in a square domain with a few islands the
-ice model is initialized with constant ice thickness and linearly
-increasing ice concentration to the east. The model dynamics are
-forced with a constant anticyclonic ocean gyre and by variable
-atmospheric wind whose mean direction is diagnonal to the north-east
-corner of the domain; ice volume and concentration are held constant
-(no thermodynamics and no advection by ice velocity).
-\reffig{hunke01} shows the ice velocity field, its divergence, and the
-bulk viscosity $\zeta$ for the cases C-LRSns and C-EVPns, and for a
-C-EVPns case, where \citet{hunke01}'s regularization has been
-implemented; compare to Fig.\,4 in \citet{hunke01}. The regularization
-contraint limits ice strength and viscosities as a function of damping
-time scale, resolution and EVP-time step, effectively allowing the
-elastic waves to damp out more quickly \citep{hunke01}.
-\begin{figure*}[htbp]
-  \includegraphics[width=\widefigwidth]{\fpath/hun12days}
-  \caption{Ice flow, divergence and bulk viscosities of three
-    experiments with \citet{hunke01}'s test case: C-LSRns (top),
-    C-EVPns (middle), and C-EVPns with damping described in
-    \citet{hunke01} (bottom).}
-  \label{fig:hunke01}
-\end{figure*}
-
-In the far right (``east'') side of the domain the ice concentration
-is close to one and the ice should be nearly rigid. The applied wind
-tends to push ice toward the upper right corner. Because the highly
-compact ice is confined by the boundary, it resists any further
-compression and exhibits little motion in the rigid region on the
-right hand side. The C-LSRns solution (top row) allows high
-viscosities in the rigid region suppressing nearly all flow.
-\citet{hunke01}'s regularization for the C-EVPns solution (bottom row)
-clearly suppresses the noise present in $\nabla\cdot\vek{u}$ and
-$\log_{10}\zeta$ in the
-unregularized case (middle row), at the cost of reduced viscosities.
-These reduced viscosities lead to small but finite ice velocities
-which in turn can have a strong effect on solutions in the limit of
-nearly rigid regimes (arching and blocking, not shown).
-
-\ml{[Say something about performance? This is tricky, as the
-  perfomance depends strongly on the configuration. A run with slowly
-  changing forcing is favorable for LSR, because then only very few
-  iterations are required for convergences while EVP uses its fixed
-  number of internal timesteps. If the forcing in changing fast, LSR
-  needs far more iterations while EVP still uses the fixed number of
-  internal timesteps. I have produces runs where for slow forcing LSR
-  is much faster than EVP and for fast forcing, LSR is much slower
-  than EVP. EVP is certainly more efficient in terms of vectorization
-  and MFLOPS on our SX8, but is that a criterion?]}
 
 \subsection{C-grid}
 \begin{itemize}
@@ -493,13 +432,18 @@
 \subsection{Arctic Domain with Open Boundaries}
 \label{sec:arctic}
 
-The Arctic domain of integration is illustrated in Fig.~\ref{???}.  It
+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
