/[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.11 by mlosch, Mon Feb 25 19:30:56 2008 UTC revision 1.14 by dimitri, Tue Feb 26 00:13:20 2008 UTC
# Line 52  Line 52 
52  \maketitle  \maketitle
53    
54  \begin{abstract}  \begin{abstract}
   
55  As part of ongoing efforts to obtain a best possible synthesis of most  As part of ongoing efforts to obtain a best possible synthesis of most
56  available, global-scale, ocean and sea ice data, dynamic and thermodynamic  available, global-scale, ocean and sea ice data, a dynamic and thermodynamic
57  sea-ice model components have been incorporated in the Massachusetts Institute  sea-ice model has been coupled to the Massachusetts Institute of Technology
58  of Technology general circulation model (MITgcm).  Sea-ice dynamics use either  general circulation model (MITgcm).  Ice mechanics follow a viscous plastic
59  a visco-plastic rheology solved with a line successive relaxation (LSR)  rheology and the ice momentum equations are solved numerically using either
60  technique, reformulated on an Arakawa C-grid in order to match the oceanic and  line successive relaxation (LSR) or elastic-viscous-plastic (EVP) dynamic
61  atmospheric grids of the MITgcm, and modified to permit efficient and accurate  models.  Ice thermodynamics are represented using either a zero-heat-capacity
62  automatic differentiation of the coupled ocean and sea-ice model  formulation or a two-layer formulation that conserves enthalpy.  The model
63  configurations.  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
79    tool complementary to an ocean model was a major design
80  \section{Model}  requirement early on in the development of the MIT general
81  \label{sec:model}  circulation model (MITgcm) [Marshall et al. 1997a,
82    Marotzke et al. 1999, Adcroft et al. 2002]. It was recognized
83    that the adjoint permitted very efficient computation
84    of gradients of various scalar-valued model diagnostics,
85    norms or, generally, objective functions with respect
86    to external or independent parameters. Such gradients
87    arise in at least two major contexts. If the objective function
88    is the sum of squared model vs. obervation differences
89    weighted by e.g. the inverse error covariances, the gradient
90    of the objective function can be used to optimize this measure
91    of model vs. data misfit in a least-squares sense. One
92    is then solving a problem of statistical state estimation.
93    If the objective function is a key oceanographic quantity
94    such as meridional heat or volume transport, ocean heat
95    content or mean surface temperature index, the gradient
96    provides a complete set of sensitivities of this quantity
97    with respect to all independent variables simultaneously.
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
105  discretized on Arakawa-B-grids \citep[e.g.,][]{hibler79, harder99,  discretized on Arakawa-B-grids \citep[e.g.,][]{hibler79, harder99,
106    kreyscher00, zhang98, hunke97}. From the perspective of coupling a  kreyscher00, zhang98, hunke97}. From the perspective of coupling a
107  sea ice-model to a C-grid ocean model, the exchange of fluxes of heat  sea ice-model to a C-grid ocean model, the exchange of fluxes of heat
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
# Line 84  velocities points and thus needs to be i Line 111  velocities points and thus needs to be i
111  sea-ice model and a C-grid ocean model. While the smoothing implicitly  sea-ice model and a C-grid ocean model. While the smoothing implicitly
112  associated with this interpolation may mask grid scale noise, it may  associated with this interpolation may mask grid scale noise, it may
113  in two-way coupling lead to a computational mode as will be shown. By  in two-way coupling lead to a computational mode as will be shown. By
114  choosing a C-grid for the sea-ice model, we circumvene this difficulty  choosing a C-grid for the sea-ice model, we circumvent this difficulty
115  altogether and render the stress coupling as consistent as the  altogether and render the stress coupling as consistent as the
116  buoyancy coupling.  buoyancy coupling.
117    
# Line 92  A further advantage of the C-grid formul Line 119  A further advantage of the C-grid formul
119  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
120  flux allowed for a B-grid (with no-slip lateral boundary counditions),  flux allowed for a B-grid (with no-slip lateral boundary counditions),
121  whereas the C-grid formulation allows a flux of sea-ice through this  whereas the C-grid formulation allows a flux of sea-ice through this
122  passage for all types of lateral boundary conditions. We (will)  passage for all types of lateral boundary conditions. We
123  demonstrate this effect in the Candian archipelago.  demonstrate this effect in the Candian archipelago.
124    
125    Talk about problems that make the sea-ice-ocean code very sensitive and
126    changes in the code that reduce these sensitivities.
127    
128    This paper describes the MITgcm sea ice
129    model; it presents example Arctic and Antarctic results from a realistic,
130    eddy-permitting, global ocean and sea-ice configuration; it compares B-grid
131    and C-grid dynamic solvers in a regional Arctic configuration; and it presents
132    example results from coupled ocean and sea-ice adjoint-model integrations.
133    
134    \section{Model}
135    \label{sec:model}
136    
137  \subsection{Dynamics}  \subsection{Dynamics}
138  \label{sec:dynamics}  \label{sec:dynamics}
139    
140  The momentum equations of the sea-ice model are standard with  The momentum equation of the sea-ice model is
141  \begin{equation}  \begin{equation}
142    \label{eq:momseaice}    \label{eq:momseaice}
143    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} +
144    \vtau_{ocean} - m \nabla{\phi(0)} + \vek{F},    \vtau_{ocean} - mg \nabla{\phi(0)} + \vek{F},
145  \end{equation}  \end{equation}
146  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;
147  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;
148  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$
149  surface height potential beneath the ice. $\phi$ is the sum of  directions, respectively;
150  atmpheric pressure $p_{a}$ and loading due to ice and snow  $f$ is the Coriolis parameter;
151  $(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,
152  ice-ocean stresses, respectively.  $\vek{F}$ is the interaction force  respectively;
153  and $\vek{i}$, $\vek{j}$, and $\vek{k}$ are the unit vectors in the  $g$ is the gravity accelation;
154  $x$, $y$, and $z$ directions.  Advection of sea-ice momentum is  $\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height;
155  neglected. The wind and ice-ocean stress terms are given by  $\phi(0)$ is the sea surface height potential in response to ocean dynamics
156    and to atmospheric pressure loading;
157    and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress
158    tensor $\sigma_{ij}$.
159    When using the rescaled vertical coordinate system, z$^\ast$, of
160    \citet{cam08}, $\phi(0)$ also includes a term due to snow and ice loading, $mg$.
161    Advection of sea-ice momentum is neglected. The wind and ice-ocean stress
162    terms are given by
163  \begin{align*}  \begin{align*}
164    \vtau_{air} =& \rho_{air} |\vek{U}_{air}|R_{air}(\vek{U}_{air}) \\    \vtau_{air}   = & \rho_{air}  C_{air}   |\vek{U}_{air}  -\vek{u}|
165    \vtau_{ocean} =& \rho_{ocean} |\vek{U}_{ocean}-\vek{u}|                     R_{air}  (\vek{U}_{air}  -\vek{u}), \\
166      \vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}|
167                     R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\                     R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\
168  \end{align*}  \end{align*}
169  where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere  where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere
170  and surface currents of the ocean, respectively. $C_{air/ocean}$ are  and surface currents of the ocean, respectively; $C_{air/ocean}$ are
171  air and ocean drag coefficients, $\rho_{air/ocean}$ reference  air and ocean drag coefficients; $\rho_{air/ocean}$ are reference
172  densities, and $R_{air/ocean}$ rotation matrices that act on the  densities; and $R_{air/ocean}$ are rotation matrices that act on the
173  wind/current vectors. $\vek{F} = \nabla\cdot\sigma$ is the divergence  wind/current vectors.
 of the interal stress tensor $\sigma_{ij}$.  
174    
175  For an isotropic system this stress tensor can be related to the ice  For an isotropic system this stress tensor can be related to the ice
176  strain rate and strength by a nonlinear viscous-plastic (VP)  strain rate and strength by a nonlinear viscous-plastic (VP)
# Line 168  The bulk viscosities are bounded above b Line 214  The bulk viscosities are bounded above b
214  $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a  $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a
215  maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where  maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where
216  $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress  $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress
217  tensor compuation the replacement pressure $P = 2\,\Delta\zeta$  tensor computation the replacement pressure $P = 2\,\Delta\zeta$
218  \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
219  elliptic yield curve by definition.  elliptic yield curve by definition.
220    
# Line 192  same length as in the ocean model where Line 238  same length as in the ocean model where
238  treated explicitly.  treated explicitly.
239    
240  \citet{hunke97}'s introduced an elastic contribution to the strain  \citet{hunke97}'s introduced an elastic contribution to the strain
241  rate elatic-viscous-plastic in order to regularize  rate elastic-viscous-plastic in order to regularize
242  Eq.\refeq{vpequation} in such a way that the resulting  Eq.\refeq{vpequation} in such a way that the resulting
243  elatic-viscous-plastic (EVP) and VP models are identical at steady  elastic-viscous-plastic (EVP) and VP models are identical at steady
244  state,  state,
245  \begin{equation}  \begin{equation}
246    \label{eq:evpequation}    \label{eq:evpequation}
# Line 220  $\sigma_{12}$. Introducing the divergenc Line 266  $\sigma_{12}$. Introducing the divergenc
266  \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension  \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension
267  and shearing strain rates, $D_T =  and shearing strain rates, $D_T =
268  \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =  \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S =
269  2\dot{\epsilon}_{12}$, respectively and using the above abbreviations,  2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations,
270  the equations can be written as:  the equations can be written as:
271  \begin{align}  \begin{align}
272    \label{eq:evpstresstensor1}    \label{eq:evpstresstensor1}
# Line 250  differences and averaging is only involv Line 296  differences and averaging is only involv
296  $P$ at vorticity points.  $P$ at vorticity points.
297    
298  For a general curvilinear grid, one needs in principle to take metric  For a general curvilinear grid, one needs in principle to take metric
299  terms into account that arise in the transformation a curvilinear grid  terms into account that arise in the transformation of a curvilinear grid
300  on the sphere. However, for now we can neglect these metric terms  on the sphere. For now, however, we can neglect these metric terms
301  because they are very small on the cubed sphere grids used in this  because they are very small on the cubed sphere grids used in this
302  paper; in particular, only near the edges of the cubed sphere grid, we  paper; in particular, only near the edges of the cubed sphere grid, we
303  expect them to be non-zero, but these edges are at approximately  expect them to be non-zero, but these edges are at approximately
# Line 260  simulations.  Everywhere else the coordi Line 306  simulations.  Everywhere else the coordi
306  cartesian.  However, for last-glacial-maximum or snowball-earth-like  cartesian.  However, for last-glacial-maximum or snowball-earth-like
307  simulations the question of metric terms needs to be reconsidered.  simulations the question of metric terms needs to be reconsidered.
308  Either, one includes these terms as in \citet{zhang03}, or one finds a  Either, one includes these terms as in \citet{zhang03}, or one finds a
309  vector-invariant formulation fo the sea-ice internal stress term that  vector-invariant formulation for the sea-ice internal stress term that
310  does not require any metric terms, as it is done in the ocean dynamics  does not require any metric terms, as it is done in the ocean dynamics
311  of the MITgcm \citep{adcroft04:_cubed_sphere}.  of the MITgcm \citep{adcroft04:_cubed_sphere}.
312    
# Line 367  differences between the two main options Line 413  differences between the two main options
413  \subsection{Arctic Domain with Open Boundaries}  \subsection{Arctic Domain with Open Boundaries}
414  \label{sec:arctic}  \label{sec:arctic}
415    
416  The Arctic domain of integration is illustrated in Fig.~\ref{???}.  It  The Arctic domain of integration is illustrated in Fig.~\ref{fig:arctic1}.  It
417  is carved out from, and obtains open boundary conditions from, the  is carved out from, and obtains open boundary conditions from, the
418  global cubed-sphere configuration of the Estimating the Circulation  global cubed-sphere configuration of the Estimating the Circulation
419  and Climate of the Ocean, Phase II (ECCO2) project  and Climate of the Ocean, Phase II (ECCO2) project
420  \citet{menemenlis05}.  The domain size is 420 by 384 grid boxes  \citet{menemenlis05}.  The domain size is 420 by 384 grid boxes
421  horizontally with mean horizontal grid spacing of 18 km.  horizontally with mean horizontal grid spacing of 18 km.
422    
423    \begin{figure}
424    %\centerline{{\includegraphics*[width=0.44\linewidth]{\fpath/arctic1.eps}}}
425    \caption{Bathymetry of Arctic Domain.\label{fig:arctic1}}
426    \end{figure}
427    
428  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
429  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
430  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.11  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22