--- MITgcm_contrib/articles/ceaice/ceaice.tex	2008/01/14 15:46:54	1.5
+++ MITgcm_contrib/articles/ceaice/ceaice.tex	2008/01/21 08:06:00	1.9
@@ -1,3 +1,5 @@
+% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.9 2008/01/21 08:06:00 mlosch Exp $
+% $Name:  $
 \documentclass[12pt]{article}
 
 \usepackage[]{graphicx}
@@ -134,7 +136,7 @@
 both thickness $h$ and compactness (concentration) $c$:
 \begin{equation}
   P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]},
-\label{icestrength}
+\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
@@ -160,6 +162,16 @@
 \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,
@@ -376,12 +388,12 @@
 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 has only a small
+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 velocity field
-appears to be a little noisy. This noise has been address by
+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
@@ -390,55 +402,51 @@
 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 variable
-atmospheric wind whose mean directed diagnonally to the north-east
+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 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}.
+(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{Hunke's test case.}
+  \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 confinded by the boundary, it resists any further
+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}$ in the
-unregularized case (middle row), at the cost of reduced viscosities
+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).
 
-
-%\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.
+\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}
@@ -485,45 +493,50 @@
 \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{???}.  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.
 
 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,
@@ -552,6 +565,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
@@ -563,8 +577,8 @@
 \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}
@@ -645,7 +659,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.
@@ -668,9 +682,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 
@@ -681,46 +696,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,
@@ -770,7 +780,7 @@
 -48 months}]
 {\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}
 
@@ -795,10 +805,9 @@
 
 \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}
-\bibliography{journal_abrvs,seaice,genocean,maths,mitgcmuv,bib/fram}
+\bibliography{bib/journal_abrvs,bib/seaice,bib/genocean,bib/maths,bib/mitgcmuv,bib/fram}
 
 \end{document}