articles/ceaice/ceaice_forward.tex

\section{Forward sensitivity experiments}
\label{sec:forward}

This section presents results from global and regional coupled ocean and sea
ice simulations that exercise various capabilities of the MITgcm sea ice
model.  The first set of results is from a global, eddy-permitting, ocean and
sea ice configuration.  The second set of results is from a regional Arctic
configuration, which is used to compare the B-grid and C-grid dynamic solvers
and various other capabilities of the MITgcm sea ice model.

\subsection{Global Ocean and Sea Ice Simulation}
\label{sec:global}

The global ocean and sea ice results presented below were carried out as part
of the Estimating the Circulation and Climate of the Ocean, Phase II (ECCO2)
project.  ECCO2 aims to produce increasingly accurate syntheses of all
available global-scale ocean and sea-ice data at resolutions that start to
resolve ocean eddies and other narrow current systems, which transport heat,
carbon, and other properties within the ocean \citep{menemenlis05}.  The
particular ECCO2 simulation discussed next is a baseline 28-year (1979-2006)
integration, labeled cube76, which has not yet been constrained by oceanic and
by sea ice data.  A cube-sphere grid projection is employed, which permits
relatively even grid spacing throughout the domain and which avoids polar
singularities \citep{adcroft04:_cubed_sphere}. Each face of the cube comprises
510 by 510 grid cells for a mean horizontal grid spacing of 18 km. 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. The model employs the
partial-cell formulation of
\citet{adcroft97:_shaved_cells}, which permits accurate representation of the
bathymetry. Bathymetry is from the S2004 (Smith, unpublished) blend of the
\citet{smi97} and the General Bathymetric Charts of the Oceans (GEBCO) one
arc-minute bathymetric grid (see Fig.~\ref{fig:CubeBathymetry}).
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{jac95} is
used.

\begin{figure}[h]
  \centering
  \includegraphics[width=\textwidth]{\fpath/CubeBathymetry}
  \caption{Bathymetry of the global cubed sphere model configuration.  The
    solid lines indicate domain boundaries for the regional Arctic
    configuration discussed in Section~\ref{sec:arctic}.} 
  \label{fig:CubeBathymetry}
\end{figure}

The ocean model is coupled to the sea-ice model discussed in
\refsec{model} using the following specific options.  The
zero-heat-capacity thermodynamics formulation of \citet{hibler80} is
used to compute sea ice thickness and concentration.  Snow cover and
sea ice salinity are prognostic.  Open water, dry ice, wet ice, dry
snow, and wet snow albedo are, respectively, 0.15, 0.88, 0.79, 0.97,
and 0.83. Ice mechanics follow the viscous plastic rheology of
\citet{hibler79} and the ice momentum equation is solved numerically
using the C-grid implementation of the \citet{zhang97}'s LSOR dynamics
model discussed hereinabove.  The ice is coupled to the ocean using
the rescaled vertical coordinate system, z$^\ast$, of \citet{cam08},
that is, sea ice does not float above the ocean model but rather
deforms the ocean's model surface level.

This particular ECCO2 simulation is initialized from temperature and salinity
fields derived from the Polar science center Hydrographic Climatology (PHC)
3.0 \citep{ste01a}. Surface boundary conditions for the period January 1979 to
July 2002 are derived from the European Centre for Medium-Range Weather
Forecasts (ECMWF) 40 year re-analysis (ERA-40) \citep{upp05}.  Surface
boundary conditions after September 2002 are derived from the ECMWF
operational analysis.  There is a one month transition period, August 2002,
during which the ERA-40 contribution decreases linearly from 1 to 0 and the
ECMWF analysis contribution increases linearly from 0 to 1.  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 \citet{pau77}.  Low frequency
precipitation has been adjusted using the pentad (5-day) data from the Global
Precipitation Climatology Project \citep[GPCP][]{huf01}.  The time-mean river
run-off from \citet{lar01} is applied globally, except in the Arctic Ocean
where monthly mean river runoff based on the Arctic Runoff Data Base (ARDB)
and prepared by P. Winsor (personnal communication, 2007) is specificied.
Additionally, there is a relaxation to the monthly-mean climatological sea
surface salinity values from PHC 3.0, a relaxation time scale of 101 days.

Vertical mixing follows \citet{lar94} but with meridionally and vertically
varying background vertical diffusivity; at the surface, vertical diffusivity
is $4.4\times 10^{-6}$~m$^2$~s$^{-1}$ at the Equator, $3.6\times
10^{-6}$~m$^2$~s$^{-1}$ north of 70$^\circ$N, and $1.9\times
10^{-5}$~m$^2$~s$^{-1}$ south of 30$^\circ$S and between 30$^\circ$N and
60$^\circ$N , with sinusoidally varying values in between these latitudes;
vertically, diffusivity increases to $1.1\times 10^{-4}$~m$^2$~s$^{-1}$ at a a
depth of 6150 m as per \citet{bry79}.  A high order monotonicity-preserving
advection scheme \citep{dar04} is employed and there is no explicit horizontal
diffusivity.  Horizontal viscosity follows \citet{lei96} but modified to sense
the divergent flow as per \citet{kem08}.

\ml{[Dimitris, here you need to either provide figures, so that I can
  write text, or you can provide both figures and text. I guess, one
  figure, showing the northern and southern hemisphere in summer and
  winter is fine (four panels), as we are showing so many figures in
  the next section.]}


\subsection{Arctic Domain with Open Boundaries}
\label{sec:arctic}

A series of forward sensitivity experiments have been carried out on
an Arctic Ocean domain with open boundaries.  The objective is to
compare the old B-grid LSR dynamic solver with the new C-grid LSR and
EVP solvers.  Additional experiments are carried out to illustrate
the differences between different ice advection schemes, ocean-ice
stress formulations and the two main options for sea ice
thermodynamics in the MITgcm.

The Arctic domain of integration is illustrated in
\reffig{arctic_topog}.  It is carved out from, and obtains open
boundary conditions from, the global cubed-sphere configuration
described above.  The horizontal domain size is 420 by 384 grid boxes.
\begin{figure*}
%\includegraphics*[width=0.44\linewidth,viewport=139 210 496 606,clip]{\fpath/topography}
%\includegraphics*[width=0.44\linewidth,viewport=0 0 496 606,clip]{\fpath/topography}
%\includegraphics*[width=0.44\linewidth]{\fpath/topography}
%\includegraphics*[width=0.46\linewidth]{\fpath/archipelago}
\includegraphics*[width=\linewidth]{\fpath/topography}
\caption{Left: Bathymetry and domain boudaries of Arctic
  Domain; the dashed line marks the boundaries of the inset on the
  right hand side. The letters in the inset label sections in the
  Canadian Archipelago, where ice transport is evaluated: 
  A: Nares Strait; %
  B: \ml{Meighen Island}; %
  C: Prince Gustaf Adolf Sea; %
  D: \ml{Brock Island}; %
  E: M'Clure Strait; %
  F: Amundsen Gulf; %
  G: Lancaster Sound; %
  H: Barrow Strait \ml{W.}; %
  I: Barrow Strait \ml{E.}; %
  J: Barrow Strait \ml{N.}; %
  K: Fram Strait. %
  The sections A through F comprise the total inflow into the Canadian
  Archipelago. \ml{[May still need to check the geography.]}
  \label{fig:arctic_topog}}
\end{figure*}

The main dynamic difference from cube sphere is that the Arctic domain
configuration does not use rescaled vertical coordinates (z$^\ast$)
and the surface boundary conditions for freshwater input are
different, because those features are not supported by the open
boundary code.
%
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.

The model is integrated from Jan~01, 1992 to Mar~31, 2000,
with three different dynamical solvers, two different boundary
conditions, different stress coupling, rheology, and advection
schemes. \reftab{experiments} gives an overview over the experiments
discussed in this section.
\begin{table}[t]
  \caption{Overview over model simulations in \refsec{arctic}.
    \label{tab:experiments}}
  \begin{tabular}{p{.3\linewidth}p{.65\linewidth}}
    experiment name & description \\ \hline
    B-LSR-ns       &  the original LSOR solver of \citet{zhang97} on an
  Arakawa B-grid, implying no-slip lateral boundary conditions
  ($\vek{u}=0$ exactly) \\
    C-LSR-ns       &  the LSOR solver discretized on a C-grid with no-slip lateral
  boundary conditions (implemented via ghost-points) \\
    C-LSR-fs       &  the LSOR solver on a C-grid with free-slip lateral boundary
  conditions \\
    C-EVP-ns       &  the EVP solver of \citet{hunke01} on a C-grid with
  no-slip lateral boundary conditions and $\Delta{t}_\mathrm{evp} =
  150\text{\,s}$ \\
    C-EVP-ns10     &  the EVP solver of \citet{hunke01} on a C-grid with
  no-slip lateral boundary conditions and $\Delta{t}_\mathrm{evp} =
  10\text{\,s}$ \\
    C-LSR-ns HB87  &  C-LSR-ns with ocean-ice stress coupling according
  to \citet{hibler87}\\
    C-LSR-ns TEM   &  C-LSR-ns with a truncated ellispe method (TEM)
  rheology \citep{hibler97} \\
    C-LSR-ns WTD   &   C-LSR-ns with 3-layer thermodynamics following
  \citet{winton00} \\
    C-LSR-ns DST3FL& C-LSR-ns with a third-order flux limited
  direct-space-time advection scheme for thermodynamic variables
  \citep{hundsdorfer94}
  \end{tabular}
\end{table}
%\begin{description}
%\item[B-LSR-ns:] the original LSOR solver of \citet{zhang97} on an
%  Arakawa B-grid, implying no-slip lateral boundary conditions
%  ($\vek{u}=0$ exactly);
%\item[C-LSR-ns:] the LSOR solver discretized on a C-grid with no-slip lateral
%  boundary conditions (implemented via ghost-points);
%\item[C-LSR-fs:] the LSOR solver on a C-grid with free-slip lateral boundary
%  conditions;
%\item[C-EVP-ns:] the EVP solver of \citet{hunke01} on a C-grid with
%  no-slip lateral boundary conditions and $\Delta{t}_\mathrm{evp} =
%  150\text{\,s}$; 
%\item[C-EVP-fs:] the EVP solver on a C-grid with free-slip lateral
%  boundary conditions  and $\Delta{t}_\mathrm{evp} = 150\text{\,s}$;
%\item[C-LSR-ns DST3FL:] C-LSR-ns with a third-order flux limited
%  direct-space-time advection scheme \citep{hundsdorfer94};
%\item[C-LSR-ns TEM:] C-LSR-ns with a truncated ellispe method (TEM)
%  rheology \citep{hibler97};
%\item[C-LSR-ns HB87:] C-LSR-ns with ocean-ice stress coupling according
%  to \citet{hibler87};
%\item[C-LSR-ns WTD:] C-LSR-ns with 3-layer thermodynamics following
%  \citet{winton00};
%%\item[C-EVP-ns damp:] C-EVP-ns with additional damping to reduce small
%%  scale noise \citep{hunke01};
%\item[C-EVP-ns10:] the EVP solver of \citet{hunke01} on a C-grid with
%  no-slip lateral boundary conditions and $\Delta{t}_\mathrm{evp} =
%  10\text{\,s}$.
%\end{description}
Both LSOR and EVP solvers solve the same viscous-plastic rheology, so
that differences between runs B-LSR-ns, C-LSR-ns, and C-EVP-ns can be
interpreted as pure model error. Lateral boundary conditions on a
coarse grid (coarse compared to the roughness of the true coast line) are
unclear, so that comparing the no-slip solutions to the free-slip
solutions gives another measure of uncertainty in sea ice modeling.
The remaining experiments explore further sensitivities of the system
to different physics (change in rheology, advection and diffusion
properties, stress coupling, and thermodynamics) and different time
steps for the EVP solutions: \citet{hunke01} uses 120 subcycling steps
for the EVP solution. We use two interpretations of this choice where
the EVP model is subcycled 120 times within a (short) model timestep
of 1200\,s resulting in a very long and expensive integration
($\Delta{t}_\mathrm{evp}=10\text{\,s}$) and 120 times within the
forcing timescale of 6\,h ($\Delta{t}_\mathrm{evp}=150\text{\,s}$).

A principle difficulty in comparing the solutions obtained with
different realizations of the model dynamics lies in the non-linear
feedback of the ice dynamics and thermodynamics. Already after a few
months the solutions have diverged so far from each other that
comparing velocities only makes sense within the first 3~months of the
integration while the ice distribution is still close to the initial
conditions. At the end of the integration, the differences between the
model solutions can be interpreted as cumulated model uncertainties. 

\reffig{iceveloc} shows ice velocities averaged over Janunary,
February, and March (JFM) of 1992 for the C-LSR-ns solution; also
shown are the differences between B-grid and C-grid, LSR and EVP, and
no-slip and free-slip solution. The velocity field of the C-LSR-ns
solution (\reffig{iceveloc}a) roughly resembles the drift velocities
of some of the AOMIP (Arctic Ocean Model Intercomparison Project)
models in a cyclonic circulation regime (CCR) \citep[their
Figure\,6]{martin07} with a Beaufort Gyre and a transpolar drift
shifted eastwards towards Alaska.

\newcommand{\subplotwidth}{0.44\textwidth}
%\newcommand{\subplotwidth}{0.3\textwidth}
\begin{figure}[tp]
  \centering
  \subfigure[{\footnotesize C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_C-LSR-ns}}
  \subfigure[{\footnotesize B-LSR-ns $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_B-LSR-ns-C-LSR-ns}}
  \\
  \subfigure[{\footnotesize C-LSR-fs $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_C-LSR-fs-C-LSR-ns}}
  \subfigure[{\footnotesize C-EVP-ns $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_C-EVP-ns150-C-LSR-ns}}
%  \\
%  \subfigure[{\footnotesize C-LSR-ns TEM $-$ C-LSR-ns}]
%  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_TEM-C-LSR-ns}}
%  \subfigure[{\footnotesize C-LSR-ns HB87 $-$ C-LSR-ns}]
%  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_HB87-C-LSR-ns}}
%  \\
%  \subfigure[{\footnotesize  C-LSR-ns WTD $-$ C-LSR-ns}]
%  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_ThSIce-C-LSR-ns}}
%  \subfigure[{\footnotesize C-LSR-ns DST3FL $-$ C-LSR-ns}]
%  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_adv33-C-LSR-ns}}
  \caption{(a) Ice drift velocity of the C-LSR-ns solution averaged
    over the first 3 months of integration [cm/s]; (b)-(h) difference
    between solutions with B-grid, free lateral slip, EVP-solver,
    truncated ellipse method (TEM), different ice-ocean stress
    formulation (HB87), different thermodynamics (WTD), different
    advection for thermodynamic variables (DST3FL) and the C-LSR-ns
    reference solution [cm/s]; color indicates speed (or differences
    of speed), vectors indicate direction only.}
  \label{fig:iceveloc}
\end{figure}
\addtocounter{figure}{-1}
\setcounter{subfigure}{4}
\begin{figure}[t]
  \subfigure[{\footnotesize C-LSR-ns TEM $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_TEM-C-LSR-ns}}
  \subfigure[{\footnotesize C-LSR-ns HB87 $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_HB87-C-LSR-ns}}
  \\
  \subfigure[{\footnotesize  C-LSR-ns WTD $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_ThSIce-C-LSR-ns}}
  \subfigure[{\footnotesize C-LSR-ns DST3FL $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMuv_adv33-C-LSR-ns}}
  \caption{continued}
\end{figure}
The difference beween runs C-LSR-ns and B-LSR-ns (\reffig{iceveloc}b)
is most pronounced along the coastlines, where the discretization
differs most between B and C-grids: On a B-grid the tangential
velocity lies on the boundary (and is thus zero through the no-slip
boundary conditions), whereas on the C-grid it is half a cell width
away from the boundary, thus allowing more flow. The B-LSR-ns solution
has less ice drift through the Fram Strait and especially the along
Greenland's east coast; also, the flow through Baffin Bay and Davis
Strait into the Labrador Sea is reduced with respect the C-LSR-ns
solution.  \ml{[Do we expect this? Say something about that]}
%
Compared to the differences between B and C-grid solutions, the
C-LSR-fs ice drift field differs much less from the C-LSR-ns solution
(\reffig{iceveloc}c).  As expected the differences are largest along
coastlines: because of the free-slip boundary conditions, flow is
faster in the C-LSR-fs solution, for example, along the east coast
of Greenland, the north coast of Alaska, and the east Coast of Baffin
Island.

The C-EVP-ns solution with $\Delta{t}_\mathrm{evp}=150\text{\,s}$ is
very different from the C-LSR-ns solution (\reffig{iceveloc}d). The
EVP-approximation of the VP-dynamics allows for increased drift by
over 2\,cm/s in the Beaufort Gyre and the transarctic drift.
%\ml{Also the Beaufort Gyre is moved towards Alaska in the C-EVP-ns
%solution. [Really?, No]} 
In general, drift velocities are biased towards higher values in the
EVP solutions.
% as can be seen from a histogram of the differences in
%\reffig{drifthist}.
%\begin{figure}[htbp]
%  \centering
%  \includegraphics[width=\textwidth]{\fpath/drifthist_C-EVP-ns-C-LSR-ns}
%  \caption{Histogram of drift velocity differences for C-LSR-ns and
%    C-EVP-ns solution [cm/s].} 
%  \label{fig:drifthist}
%\end{figure}

Compared to the other parameters, the ice rheology TEM
(\reffig{iceveloc}(e)) has a very small effect on the solution. In
general the ice drift tends to be increased, because there is no
tensile stress and ice can be ``pulled appart'' at no cost.
Consequently, the largest effect on drift velocity can be observed
near the ice edge in the Labrador Sea. In contrast, in the run with
the ice-ocean stress formulation of \citet{hibler87},
\reffig{iceveloc}(f) the drift is stronger almost everywhere in the
computational domain. The increase is mostly aligned with the general
direction of the flow, implying that the different stress formulation
reduces the deceleration of drift by the ocean.

The 3-layer thermodynamics following \citet{winton00} requires
additional information on initial conditions for enthalphy. These
different initial conditions make a comparison of the first months
difficult to interpret. The drift in the Beaufort Gyre is slightly
reduced relative to the reference run C-LSR-ns, but the drift through
the Fram Strait is increased. The drift velocities near the ice edge
are very different, because the ice extend is already larger in
\mbox{C-LSR-ns~WTD}; inward from the ice egde, this run has smaller
drift velocities, because the ice motion is more contrained by a
larger ice extent than in \mbox{C-LSR-ns}, where the ice at the same
place is drifting nearly freely.

A more sophisticated advection scheme (\mbox{C-LSR-ns DST3FL},
\reffig{iceveloc}(h)) has its largest effect along the ice edge, where
the gradients of thickness and concentration are largest. Everywhere
else the effect is very small and can mostly be attributed to smaller
numerical diffusion (and to the absense of explicitly diffusion for
numerical stability).

\reffig{icethick}a shows the effective thickness (volume per unit
area) of the C-LSR-ns solution, averaged over January, February, March
of year 2000. By this time of the integration, the differences in the
ice drift velocities have led to the evolution of very different ice
thickness distributions, which are shown in \reffig{icethick}b--d, and
concentrations (not shown).
\begin{figure}[tp]
  \centering
  \subfigure[{\footnotesize C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMheff2000_C-LSR-ns}}
  \subfigure[{\footnotesize B-LSR-ns $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMheff2000_B-LSR-ns-C-LSR-ns}}
  \\
  \subfigure[{\footnotesize C-LSR-fs $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMheff2000_C-LSR-fs-C-LSR-ns}}
  \subfigure[{\footnotesize C-EVP-ns $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMheff2000_C-EVP-ns150-C-LSR-ns}}
  \caption{(a) Effective thickness (volume per unit area) of the
    C-LSR-ns solution, averaged over the months Janurary through March
    2000 [m]; (b)-(h) difference between solutions with B-grid, free
    lateral slip, EVP-solver, truncated ellipse method (TEM),
    different ice-ocean stress formulation (HB87), different
    thermodynamics (WTD), different advection for thermodynamic
    variables (DST3FL) and the C-LSR-ns reference solution [m].}
  \label{fig:icethick}
\end{figure}
\addtocounter{figure}{-1}
\setcounter{subfigure}{4}
\begin{figure}[t]
  \subfigure[{\footnotesize C-LSR-ns TEM $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMheff2000_TEM-C-LSR-ns}}
  \subfigure[{\footnotesize C-EVP-ns HB87 $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMheff2000_HB87-C-LSR-ns}}
  \\
  \subfigure[{\footnotesize C-LSR-ns WTD $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMheff2000_ThSIce-C-LSR-ns}}
  \subfigure[{\footnotesize C-EVP-ns DST3FL $-$ C-LSR-ns}]
  {\includegraphics[width=\subplotwidth]{\fpath/JFMheff2000_adv33-C-LSR-ns}}
  \caption{continued}
\end{figure}
%
The generally weaker ice drift velocities in the B-LSR-ns solution,
when compared to the C-LSR-ns solution, in particular through the
narrow passages in the Canadian Archipelago, lead to a larger build-up
of ice north of Greenland and the Archipelago by 2\,m effective
thickness and more in the B-grid solution (\reffig{icethick}b). But
the ice volume in not larger everywhere: further west, there are
patches of smaller ice volume in the B-grid solution, most likely
because the Beaufort Gyre is weaker and hence not as effective in
transporting ice westwards. There are also dipoles of ice volume
differences with more ice on the upstream side of island groups and
less ice in their lee, such as Franz-Josef-Land and
Severnaya Semlya\ml{/or Nordland?},
because ice tends to flow along coasts less easily in the B-LSR-ns
solution.

Imposing a free-slip boundary condition in C-LSR-fs leads to a much
smaller differences to C-LSR-ns in the central Arctic than the
transition from the B-grid to the C-grid (\reffig{icethick}c), except
in the Canadian Archipelago. There it reduces the effective ice
thickness by 2\,m and more where the ice is thick and the straits are
narrow.  Dipoles of ice thickness differences can also be observed
around islands, because the free-slip solution allows more flow around
islands than the no-slip solution. Everywhere else the ice volume is
affected only slightly by the different boundary condition.
%
The C-EVP-ns solution has generally stronger drift velocities than the
C-LSR-ns solution. Consequently, more ice can be moved from the
eastern part of the Arctic, where ice volumes are smaller, to the
western Arctic (\reffig{icethick}d). Within the Canadian Archipelago,
more drift leads to faster ice export and reduced effective ice
thickness. With a shorter time step of
$\Delta{t}_\mathrm{evp}=10\text{\,s}$ the EVP solution seems to
converge to the LSOR solution (not shown). Only in the narrow straits
in the Archipelago the ice thickness is not affected by the shorter
time step and the ice is still thinner by 2\,m and more, as in the EVP
solution with $\Delta{t}_\mathrm{evp}=150\text{\,s}$.

In year 2000, there more ice everywhere in the domain in
\mbox{C-LSR-ns WTD}, \reffig{icethick}(g). This difference, which is
even more pronounced in summer (not shown), can be attributed to
direct effects of the different thermodynamics in this run. The
remaining runs have the largest differences in effective ice thickness
long the north coasts of Greenland and Ellesmere Island. Although the
effects of TEM and \citet{hibler87}'s ice-ocean stress formulation are
so different on the initial ice velocities, both runs have similarly
reduced ice thicknesses in this area. The 3rd-order advection scheme
has an opposite effect of similar magnitude, point towards more
implicit lateral stress with this numerical scheme.

The observed difference of order 2\,m and less are smaller than the
differences that were observed between different hindcast models and climate
models in \citet{gerdes07}. There the range of sea ice volume of
different sea ice-ocean models (which shared very similar forcing
fields) was on the order of $10,000\text{km$^{3}$}$; this range was
even larger for coupled climate models. Here, the range (and the
averaging period) is smaller than $4,000\text{km$^{3}$}$ except for
the run \mbox{C-LSR-ns~WTD} where the more complicated thermodynamics
leads to generally thicker ice (\reffig{icethick} and
\reftab{icevolume}).
\begin{table}[htbp]
  \begin{tabular}{lr@{\hspace{5ex}}r@{$\pm$}rr@{$\pm$}rr@{$\pm$}r}
    model run & ice volume 
    & \multicolumn{6}{c}{ice transport [$\text{flux$\pm$ std.,
        km$^{3}$\,y$^{-1}$}$]}\\ 
    & [$\text{km$^{3}$}$] 
    & \multicolumn{2}{c}{FS}
    & \multicolumn{2}{c}{NI}
    & \multicolumn{2}{c}{LS} \\ \hline
    B-LSR-ns       & 23,824 & 2126 & 1278 &   34 &  122 &   43 &   76 \\
    C-LSR-ns       & 24,769 & 2196 & 1253 &   70 &  224 &   77 &  110 \\
    C-LSR-fs       & 23,286 & 2236 & 1289 &   80 &  276 &   91 &   85 \\
    C-EVP-ns       & 27,056 & 3050 & 1652 &  352 &  735 &  256 &  151 \\
    C-EVP-ns10     & 22,633 & 2174 & 1260 &  186 &  496 &  133 &  128 \\
    C-LSR-ns HB87  & 23,060 & 2256 & 1327 &   64 &  230 &   77 &  114 \\
    C-LSR-ns TEM   & 23,529 & 2222 & 1258 &   60 &  242 &   87 &  112 \\
    C-LSR-ns WTD   & 31,634 & 2761 & 1563 &   23 &  140 &   94 &   63 \\
    C-LSR-ns DST3FL& 24,023 & 2191 & 1261 &   88 &  251 &   84 &  129 
  \end{tabular}
  \caption{Arctic ice volume averaged over Jan--Mar 2000, in
    $\text{km$^{3}$}$. Mean ice transport and standard deviation for the
    period Jan 1992 -- Dec 1999 through the Fram Strait (FS), the
    total northern inflow into the Canadian Archipelago (NI), and the
    export through Lancaster Sound (LS), in $\text{km$^{3}$\,y$^{-1}$}$.} 
  \label{tab:icevolume}
\end{table}

The difference in ice volume and ice drift velocities between the
different experiments has consequences for the ice transport out of
the Arctic. Although by far the most exported ice drifts through the
Fram Strait (approximately $2300\pm610\text{\,km$^3$\,y$^{-1}$}$), a
considerable amount (order $160\text{\,km$^3$\,y$^{-1}$}$) ice is
exported through the Canadian Archipelago \citep[and references
therein]{serreze06}. Note, that ice transport estimates are associated
with large uncertainties; also note that tuning an Arctic sea
ice-ocean model to reproduce observations is not our goal, but we use
the published numbers as an orientation.

\reffig{archipelago} shows a time series of daily averaged, smoothed
with monthly running means, ice transports through various straits in
the Canadian Archipelago and the Fram Strait for the different model
solutions and \reftab{icevolume} summarizes the time series. The
export through Fram Strait agrees with the observations in all model
solutions (annual averages range from $2110$ to
$2300\text{\,km$^3$\,y$^{-1}$}$, except for \mbox{C-LSR-ns~WTD} with
$2760\text{\,km$^3$\,y$^{-1}$}$ and the EVP solution with the long
time step of 150\,s with nearly $3000\text{\,km$^3$\,y$^{-1}$}$),
while the export through the Candian Archipelago is smaller than
generally thought. For example, the ice transport through Lancaster
Sound is lower (annual averages are $43$ to
$256\text{\,km$^3$\,y$^{-1}$}$) than in \citet{dey81} who estimates an
inflow into Baffin Bay of $370$ to $537\text{\,km$^3$\,y$^{-1}$}$, but
a flow of only $102$ to $137\text{\,km$^3$\,y$^{-1}$}$ further
upstream in Barrow Strait in the 1970ies from satellite images.
Generally, the EVP solutions have the highest maximum (export out of
the Artic) and lowest minimum (import into the Artic) fluxes as the
drift velocities are largest in these solutions.  In the extreme of
the Nares Strait, which is only a few grid points wide in our
configuration, both B- and C-grid LSOR solvers lead to practically no
ice transport, while the C-EVP solutions allow up to
$600\text{\,km$^3$\,y$^{-1}$}$ in summer; \citet{tang04} report $300$
to $350\text{\,km$^3$\,y$^{-1}$}$.  As as consequence, the import into
the Candian Archipelago is larger in all EVP solutions
%(range: $539$ to $773\text{\,km$^3$\,y$^{-1}$}$)
than in the LSOR solutions.
%get the order of magnitude right (range: $132$ to
%$165\text{\,km$^3$\,y$^{-1}$}$); 
The B-LSR-ns solution is even smaller by another factor of two than the
C-LSR solutions (an exception is the WTD solution, where larger ice thickness
tends to block the transport).
%underestimates the ice transport with $34\text{\,km$^3$\,y$^{-1}$}$.
\begin{figure}
%\centerline{{\includegraphics*[width=0.6\linewidth]{\fpath/Jan1992xport}}}
%\centerline{{\includegraphics*[width=0.6\linewidth]{\fpath/ice_export}}}
\centerline{{\includegraphics[width=\linewidth]{\fpath/ice_export}}}
\caption{Transport through Canadian Archipelago for different solver
  flavors. The letters refer to the labels of the sections in
  \reffig{arctic_topog}; positive values are flux out of the Arctic;
  legend abbreviations are explained in \reftab{experiments}.
\label{fig:archipelago}}
\end{figure}

%\ml{[Transport to narrow straits, area?, more runs, TEM, advection
%  schemes, Winton TD, discussion about differences in terms of model
%  error? that's tricky as it means refering to Tremblay, thus our ice
%  models are all erroneous!]}

In summary, we find that different dynamical solvers can yield very
different solutions. In constrast to that, the differences between
free-slip and no-slip solutions \emph{with the same solver} are
considerably smaller (the difference for the EVP solver is not shown,
but similar to that for the LSOR solver). Albeit smaller, the
differences between free and no-slip solutions in ice drift can lead
to equally large differences in ice volume, especially in the Canadian
Archipelago over the integration time. At first, this observation
seems counterintuitive, as we expect that the solution
\emph{technique} should not affect the \emph{solution} to a higher
degree than actually modifying the equations. A more detailed study on
these differences is beyond the scope of this paper, but at this point
we may speculate, that the large difference between B-grid, C-grid,
LSOR, and EVP solutions stem from incomplete convergence of the
solvers due to linearization and due to different methods of
linearization \citep[and Bruno Tremblay, personal
communication]{hunke01}: if the convergence of the non-linear momentum
equations is not complete for all linearized solvers, then one can
imagine that each solver stops at a different point in velocity-space
thus leading to different solutions for the ice drift velocities. If
this were true, this tantalizing circumstance would have a dramatic
impact on sea-ice modeling in general, and we would need to improve
the solution techniques for dynamic sea ice models, most likely at a very
high compuational cost (Bruno Tremblay, personal communication). Further,
we observe that the EVP solutions tends to produce effectively
``weaker'' ice that yields more easily to stress. The fast response to
changing wind was also observed by \citet{hunke99}, their Fig.\,10--12,
where the EVP model adjusts quickly to a cyclonic wind pattern, while
the LSOR solution does not. This property of the EVP solutions allows
larger ice transports through narrow straits, where the implicit
solver LSOR forms rigid ice. The underlying reasons for this striking
difference need further exploration.

% THIS is now almost all in the text:
%\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: B-LSR-ns
%\item C-grid LSR no-slip: C-LSR-ns
%\item C-grid LSR slip:    C-LSR-fs
%\item C-grid EVP no-slip: C-EVP-ns
%\item C-grid EVP slip:    C-EVP-fs
%\item C-grid LSR + TEM (truncated ellipse method, no tensile stress,
%  new flag): C-LSR-ns+TEM
%\item C-grid LSR with different advection scheme: 33 vs 77, vs. default?
%\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}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "ceaice"
%%% End: 