/[MITgcm]/MITgcm_contrib/articles/ceaice/ceaice.tex
ViewVC logotype

Diff of /MITgcm_contrib/articles/ceaice/ceaice.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.9 by mlosch, Mon Jan 21 08:06:00 2008 UTC revision 1.16 by mlosch, Tue Feb 26 19:14:36 2008 UTC
# Line 52  Line 52 
52  \maketitle  \maketitle
53    
54  \begin{abstract}  \begin{abstract}
55    Some blabla  As part of ongoing efforts to obtain a best possible synthesis of most
56    available, global-scale, ocean and sea ice data, a dynamic and thermodynamic
57    sea-ice model has been coupled to the Massachusetts Institute of Technology
58    general circulation model (MITgcm).  Ice mechanics follow a viscous plastic
59    rheology and the ice momentum equations are solved numerically using either
60    line successive relaxation (LSR) or elastic-viscous-plastic (EVP) dynamic
61    models.  Ice thermodynamics are represented using either a zero-heat-capacity
62    formulation or a two-layer formulation that conserves enthalpy.  The model
63    includes prognostic variables for snow and for sea-ice salinity.  The above
64    sea ice model components were borrowed from current-generation climate models
65    but they were reformulated on an Arakawa C-grid in order to match the MITgcm
66    oceanic grid and they were modified in many ways to permit efficient and
67    accurate automatic differentiation.  This paper describes the MITgcm sea ice
68    model; it presents example Arctic and Antarctic results from a realistic,
69    eddy-permitting, global ocean and sea-ice configuration; it compares B-grid
70    and C-grid dynamic solvers in a regional Arctic configuration; and it presents
71    example results from coupled ocean and sea-ice adjoint-model integrations.
72    
73  \end{abstract}  \end{abstract}
74    
75  \section{Introduction}  \section{Introduction}
76  \label{sec:intro}  \label{sec:intro}
77    
78  more blabla  The availability of an adjoint model as a powerful research tool
79    complementary to an ocean model was a major design requirement early
80  \section{Model}  on in the development of the MIT general circulation model (MITgcm)
81  \label{sec:model}  [Marshall et al. 1997a, Marotzke et al. 1999, Adcroft et al. 2002]. It
82    was recognized that the adjoint model permitted computing the
83    gradients of various scalar-valued model diagnostics, norms or,
84    generally, objective functions with respect to external or independent
85    parameters very efficiently. The information associtated with these
86    gradients is useful in at least two major contexts. First, for state
87    estimation problems, the objective function is the sum of squared
88    differences between observations and model results weighted by the
89    inverse error covariances. The gradient of such an objective function
90    can be used to reduce this measure of model-data misfit to find the
91    optimal model solution in a least-squares sense.  Second, the
92    objective function can be a key oceanographic quantity such as
93    meridional heat or volume transport, ocean heat content or mean
94    surface temperature index. In this case the gradient provides a
95    complete set of sensitivities of this quantity to all independent
96    variables simultaneously. These sensitivities can be used to address
97    the cause of, say, changing net transports accurately.
98    
99    References to existing sea-ice adjoint models, explaining that they are either
100    for simplified configurations, for ice-only studies, or for short-duration
101    studies to motivate the present work.
102    
103  Traditionally, probably for historical reasons and the ease of  Traditionally, probably for historical reasons and the ease of
104  treating the Coriolis term, most standard sea-ice models are  treating the Coriolis term, most standard sea-ice models are
# Line 71  sea ice-model to a C-grid ocean model, t Line 108  sea ice-model to a C-grid ocean model, t
108  and fresh-water pose no difficulty for a B-grid sea-ice model  and fresh-water pose no difficulty for a B-grid sea-ice model
109  \citep[e.g.,][]{timmermann02a}. However, surface stress is defined at  \citep[e.g.,][]{timmermann02a}. However, surface stress is defined at
110  velocities points and thus needs to be interpolated between a B-grid  velocities points and thus needs to be interpolated between a B-grid
111  sea-ice model and a C-grid ocean model. While the smoothing implicitly  sea-ice model and a C-grid ocean model. Smoothing implicitly
112  associated with this interpolation may mask grid scale noise, it may  associated with this interpolation may mask grid scale noise and may
113  in two-way coupling lead to a computational mode as will be shown. By  contribute to stabilizing the solution. On the other hand, by
114  choosing a C-grid for the sea-ice model, we circumvene this difficulty  smoothing the stress signals are damped which could lead to reduced
115  altogether and render the stress coupling as consistent as the  variability of the system. By choosing a C-grid for the sea-ice model,
116  buoyancy coupling.  we circumvent this difficulty altogether and render the stress
117    coupling as consistent as the buoyancy coupling.
118    
119  A further advantage of the C-grid formulation is apparent in narrow  A further advantage of the C-grid formulation is apparent in narrow
120  straits. In the limit of only one grid cell between coasts there is no  straits. In the limit of only one grid cell between coasts there is no
121  flux allowed for a B-grid (with no-slip lateral boundary counditions),  flux allowed for a B-grid (with no-slip lateral boundary counditions),
122  whereas the C-grid formulation allows a flux of sea-ice through this  and models have used topographies artificially widened straits to
123  passage for all types of lateral boundary conditions. We (will)  avoid this problem \citep{holloway07}. The C-grid formulation on the
124  demonstrate this effect in the Candian archipelago.  other hand allows a flux of sea-ice through narrow passages if
125    free-slip along the boundaries is allowed. We demonstrate this effect
126    in the Candian archipelago.
127    
128    Talk about problems that make the sea-ice-ocean code very sensitive and
129    changes in the code that reduce these sensitivities.
130    
131    This paper describes the MITgcm sea ice
132    model; it presents example Arctic and Antarctic results from a realistic,
133    eddy-permitting, global ocean and sea-ice configuration; it compares B-grid
134    and C-grid dynamic solvers in a regional Arctic configuration; and it presents
135    example results from coupled ocean and sea-ice adjoint-model integrations.
136    
137    \section{Model}
138    \label{sec:model}
139    
140  \subsection{Dynamics}  \subsection{Dynamics}
141  \label{sec:dynamics}  \label{sec:dynamics}
142    
143  The momentum equations of the sea-ice model are standard with  The momentum equation of the sea-ice model is
144  \begin{equation}  \begin{equation}
145    \label{eq:momseaice}    \label{eq:momseaice}
146    m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} +    m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} +
147    \vtau_{ocean} - m \nabla{\phi(0)} + \vek{F},    \vtau_{ocean} - m \nabla{\phi(0)} + \vek{F},
148  \end{equation}  \end{equation}
149  where $\vek{u} = u\vek{i}+v\vek{j}$ is the ice velocity vectory, $m$  where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area;
150  the ice mass per unit area, $f$ the Coriolis parameter, $g$ is the  $\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector;
151  gravity accelation, $\nabla\phi$ is the gradient (tilt) of the sea  $\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$
152  surface height potential beneath the ice. $\phi$ is the sum of  directions, respectively;
153  atmpheric pressure $p_{a}$ and loading due to ice and snow  $f$ is the Coriolis parameter;
154  $(m_{i}+m_{s})g$. $\vtau_{air}$ and $\vtau_{ocean}$ are the wind and  $\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses,
155  ice-ocean stresses, respectively.  $\vek{F}$ is the interaction force  respectively;
156  and $\vek{i}$, $\vek{j}$, and $\vek{k}$ are the unit vectors in the  $g$ is the gravity accelation;
157  $x$, $y$, and $z$ directions.  Advection of sea-ice momentum is  $\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;
158  neglected. The wind and ice-ocean stress terms are given by  $\phi(0) = g\eta + p_{a}/\rho_{0}$ is the sea surface height potential
159    in response to ocean dynamics ($g\eta$) and to atmospheric pressure
160    loading ($p_{a}/\rho_{0}$, where $\rho_{0}$ is a reference density);
161    and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress
162    tensor $\sigma_{ij}$.
163    When using the rescaled vertical coordinate system, z$^\ast$, of
164    \citet{cam08}, $\phi(0)$ also includes a term due to snow and ice
165    loading, $mg/\rho_{0}$.
166    Advection of sea-ice momentum is neglected. The wind and ice-ocean stress
167    terms are given by
168  \begin{align*}  \begin{align*}
169    \vtau_{air} =& \rho_{air} |\vek{U}_{air}|R_{air}(\vek{U}_{air}) \\    \vtau_{air}   = & \rho_{air}  C_{air}   |\vek{U}_{air}  -\vek{u}|
170    \vtau_{ocean} =& \rho_{ocean} |\vek{U}_{ocean}-\vek{u}|                     R_{air}  (\vek{U}_{air}  -\vek{u}), \\
171      \vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}|
172                     R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\                     R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\
173  \end{align*}  \end{align*}
174  where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere  where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere
175  and surface currents of the ocean, respectively. $C_{air/ocean}$ are  and surface currents of the ocean, respectively; $C_{air/ocean}$ are
176  air and ocean drag coefficients, $\rho_{air/ocean}$ reference  air and ocean drag coefficients; $\rho_{air/ocean}$ are reference
177  densities, and $R_{air/ocean}$ rotation matrices that act on the  densities; and $R_{air/ocean}$ are rotation matrices that act on the
178  wind/current vectors. $\vek{F} = \nabla\cdot\sigma$ is the divergence  wind/current vectors.
179  of the interal stress tensor $\sigma_{ij}$.  
180    For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can
181  For an isotropic system this stress tensor can be related to the ice  be related to the ice strain rate and strength by a nonlinear
182  strain rate and strength by a nonlinear viscous-plastic (VP)  viscous-plastic (VP) constitutive law \citep{hibler79, zhang98}:
 constitutive law \citep{hibler79, zhang98}:  
183  \begin{equation}  \begin{equation}
184    \label{eq:vpequation}    \label{eq:vpequation}
185    \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij}    \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij}
# Line 158  The bulk viscosities are bounded above b Line 219  The bulk viscosities are bounded above b
219  $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a  $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a
220  maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where  maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where
221  $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress  $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress
222  tensor compuation the replacement pressure $P = 2\,\Delta\zeta$  tensor computation the replacement pressure $P = 2\,\Delta\zeta$
223  \citep{hibler95} is used so that the stress state always lies on the  \citep{hibler95} is used so that the stress state always lies on the
224  elliptic yield curve by definition.  elliptic yield curve by definition.
225    
# Line 166  In the so-called truncated ellipse metho Line 227  In the so-called truncated ellipse metho
227  is capped to suppress any tensile stress \citep{hibler97, geiger98}:  is capped to suppress any tensile stress \citep{hibler97, geiger98}:
228  \begin{equation}  \begin{equation}
229    \label{eq:etatem}    \label{eq:etatem}
230    \eta = \min(\frac{\zeta}{e^2}    \eta = \min\left(\frac{\zeta}{e^2},
231    \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}    \frac{\frac{P}{2}-\zeta(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})}
232    {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2    {\sqrt{(\dot{\epsilon}_{11}+\dot{\epsilon}_{22})^2
233        +4\dot{\epsilon}_{12}^2}}        +4\dot{\epsilon}_{12}^2}}\right).
234  \end{equation}  \end{equation}
235    
236  In the current implementation, the VP-model is integrated with the  In the current implementation, the VP-model is integrated with the
237  semi-implicit line successive over relaxation (LSOR)-solver of  semi-implicit line successive over relaxation (LSOR)-solver of
238  \citet{zhang98}, which allows for long time steps that, in our case,  \citet{zhang98}, which allows for long time steps that, in our case,
239  is limited by the explicit treatment of the Coriolis term. The  are limited by the explicit treatment of the Coriolis term. The
240  explicit treatment of the Coriolis term does not represent a severe  explicit treatment of the Coriolis term does not represent a severe
241  limitation because it restricts the time step to approximately the  limitation because it restricts the time step to approximately the
242  same length as in the ocean model where the Coriolis term is also  same length as in the ocean model where the Coriolis term is also
243  treated explicitly.  treated explicitly.
244    
245  \citet{hunke97}'s introduced an elastic contribution to the strain  \citet{hunke97}'s introduced an elastic contribution to the strain
246  rate elatic-viscous-plastic in order to regularize  rate in order to regularize Eq.\refeq{vpequation} in such a way that
247  Eq.\refeq{vpequation} in such a way that the resulting  the resulting elastic-viscous-plastic (EVP) and VP models are
248  elatic-viscous-plastic (EVP) and VP models are identical at steady  identical at steady state,
 state,  
249  \begin{equation}  \begin{equation}
250    \label{eq:evpequation}    \label{eq:evpequation}
251    \frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} +    \frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} +
# Line 210  $\sigma_{12}$. Introducing the divergenc Line 270  $\sigma_{12}$. Introducing the divergenc
270  \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension  \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension
271  and shearing strain rates, $D_T =  and shearing strain rates, $D_T =
272  \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =  \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =
273  2\dot{\epsilon}_{12}$, respectively and using the above abbreviations,  2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations,
274  the equations can be written as:  the equations can be written as:
275  \begin{align}  \begin{align}
276    \label{eq:evpstresstensor1}    \label{eq:evpstresstensor1}
# Line 232  $E_{0} = \frac{1}{3}$. Line 292  $E_{0} = \frac{1}{3}$.
292  For details of the spatial discretization, the reader is referred to  For details of the spatial discretization, the reader is referred to
293  \citet{zhang98, zhang03}. Our discretization differs only (but  \citet{zhang98, zhang03}. Our discretization differs only (but
294  importantly) in the underlying grid, namely the Arakawa C-grid, but is  importantly) in the underlying grid, namely the Arakawa C-grid, but is
295  otherwise straightforward. The EVP model in particular is discretized  otherwise straightforward. The EVP model, in particular, is discretized
296  naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the  naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the
297  center points and $\sigma_{12}$ on the corner (or vorticity) points of  center points and $\sigma_{12}$ on the corner (or vorticity) points of
298  the grid. With this choice all derivatives are discretized as central  the grid. With this choice all derivatives are discretized as central
# Line 240  differences and averaging is only involv Line 300  differences and averaging is only involv
300  $P$ at vorticity points.  $P$ at vorticity points.
301    
302  For a general curvilinear grid, one needs in principle to take metric  For a general curvilinear grid, one needs in principle to take metric
303  terms into account that arise in the transformation a curvilinear grid  terms into account that arise in the transformation of a curvilinear
304  on the sphere. However, for now we can neglect these metric terms  grid on the sphere. For now, however, we can neglect these metric
305  because they are very small on the cubed sphere grids used in this  terms because they are very small on the \ml{[modify following
306  paper; in particular, only near the edges of the cubed sphere grid, we    section3:] cubed sphere grids used in this paper; in particular,
307  expect them to be non-zero, but these edges are at approximately  only near the edges of the cubed sphere grid, we expect them to be
308  35\degS\ or 35\degN\ which are never covered by sea-ice in our  non-zero, but these edges are at approximately 35\degS\ or 35\degN\
309  simulations.  Everywhere else the coordinate system is locally nearly  which are never covered by sea-ice in our simulations.  Everywhere
310  cartesian.  However, for last-glacial-maximum or snowball-earth-like  else the coordinate system is locally nearly cartesian.}  However, for
311  simulations the question of metric terms needs to be reconsidered.  last-glacial-maximum or snowball-earth-like simulations the question
312  Either, one includes these terms as in \citet{zhang03}, or one finds a  of metric terms needs to be reconsidered.  Either, one includes these
313  vector-invariant formulation fo the sea-ice internal stress term that  terms as in \citet{zhang03}, or one finds a vector-invariant
314  does not require any metric terms, as it is done in the ocean dynamics  formulation for the sea-ice internal stress term that does not require
315  of the MITgcm \citep{adcroft04:_cubed_sphere}.  any metric terms, as it is done in the ocean dynamics of the MITgcm
316    \citep{adcroft04:_cubed_sphere}.
317    
318    Lateral boundary conditions are naturally ``no-slip'' for B-grids, as
319    the tangential velocities points lie on the boundary. For C-grids, the
320    lateral boundary condition for tangential velocities is realized via
321    ``ghost points'', allowing alternatively no-slip or free-slip
322    conditions. In ocean models free-slip boundary conditions in
323    conjunction with piecewise-constant (``castellated'') coastlines have
324    been shown to reduce in effect to no-slip boundary conditions
325    \citep{adcroft98:_slippery_coast}; for sea-ice models the effects of
326    lateral boundary conditions have not yet been studied.
327    
328  Moving sea ice exerts a stress on the ocean which is the opposite of  Moving sea ice exerts a stress on the ocean which is the opposite of
329  the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is  the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is
# Line 271  velocity and the ice velocity leading to Line 342  velocity and the ice velocity leading to
342  temperature and salinity are different from the oceanic variables.  temperature and salinity are different from the oceanic variables.
343    
344  Sea ice distributions are characterized by sharp gradients and edges.  Sea ice distributions are characterized by sharp gradients and edges.
345  For this reason, we employ a positive 3rd-order advection scheme  For this reason, we employ positive, multidimensional 2nd and 3rd-order
346  \citep{hundsdorfer94} for the thermodynamic variables discussed in the  advection scheme with flux limiter \citep{roe85, hundsdorfer94} for the
347  next section.  thermodynamic variables discussed in the next section.
348    
349  \subparagraph{boundary conditions: no-slip, free-slip, half-slip}  \subparagraph{boundary conditions: no-slip, free-slip, half-slip}
350    
# Line 310  content by means of enthalphy conservati Line 381  content by means of enthalphy conservati
381  addition to ice-thickness and compactness (fractional area) additional  addition to ice-thickness and compactness (fractional area) additional
382  state variables to be advected by ice velocities, namely enthalphy of  state variables to be advected by ice velocities, namely enthalphy of
383  the two ice layers and the thickness of the overlying snow layer.  the two ice layers and the thickness of the overlying snow layer.
384    \ml{[Jean-Michel, your turf: ]Care must be taken in advecting these
385      quantities in order to ensure conservation of enthalphy. Currently
386      this can only be accomplished with a 2nd-order advection scheme with
387      flux limiter \citep{roe85}.}
388    
 \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?]}  
389    
390  \subsection{C-grid}  \subsection{C-grid}
391  \begin{itemize}  \begin{itemize}
# Line 493  differences between the two main options Line 432  differences between the two main options
432  \subsection{Arctic Domain with Open Boundaries}  \subsection{Arctic Domain with Open Boundaries}
433  \label{sec:arctic}  \label{sec:arctic}
434    
435  The Arctic domain of integration is illustrated in Fig.~\ref{???}.  It  The Arctic domain of integration is illustrated in Fig.~\ref{fig:arctic1}.  It
436  is carved out from, and obtains open boundary conditions from, the  is carved out from, and obtains open boundary conditions from, the
437  global cubed-sphere configuration of the Estimating the Circulation  global cubed-sphere configuration of the Estimating the Circulation
438  and Climate of the Ocean, Phase II (ECCO2) project  and Climate of the Ocean, Phase II (ECCO2) project
439  \citet{menemenlis05}.  The domain size is 420 by 384 grid boxes  \citet{menemenlis05}.  The domain size is 420 by 384 grid boxes
440  horizontally with mean horizontal grid spacing of 18 km.  horizontally with mean horizontal grid spacing of 18 km.
441    
442    \begin{figure}
443    %\centerline{{\includegraphics*[width=0.44\linewidth]{\fpath/arctic1.eps}}}
444    \caption{Bathymetry of Arctic Domain.\label{fig:arctic1}}
445    \end{figure}
446    
447  There are 50 vertical levels ranging in thickness from 10 m near the surface  There are 50 vertical levels ranging in thickness from 10 m near the surface
448  to approximately 450 m at a maximum model depth of 6150 m. Bathymetry is from  to approximately 450 m at a maximum model depth of 6150 m. Bathymetry is from
449  the National Geophysical Data Center (NGDC) 2-minute gridded global relief  the National Geophysical Data Center (NGDC) 2-minute gridded global relief

Legend:
Removed from v.1.9  
changed lines
  Added in v.1.16

  ViewVC Help
Powered by ViewVC 1.1.22