--- MITgcm_contrib/articles/ceaice/ceaice.tex 2008/02/25 19:30:56 1.11 +++ MITgcm_contrib/articles/ceaice/ceaice.tex 2008/02/26 00:13:20 1.14 @@ -1,4 +1,4 @@ -% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.11 2008/02/25 19:30:56 mlosch Exp $ +% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.14 2008/02/26 00:13:20 dimitri Exp $ % $Name: $ \documentclass[12pt]{article} @@ -52,31 +52,58 @@ \maketitle \begin{abstract} - As part of ongoing efforts to obtain a best possible synthesis of most -available, global-scale, ocean and sea ice data, dynamic and thermodynamic -sea-ice model components have been incorporated in the Massachusetts Institute -of Technology general circulation model (MITgcm). Sea-ice dynamics use either -a visco-plastic rheology solved with a line successive relaxation (LSR) -technique, reformulated on an Arakawa C-grid in order to match the oceanic and -atmospheric grids of the MITgcm, and modified to permit efficient and accurate -automatic differentiation of the coupled ocean and sea-ice model -configurations. +available, global-scale, ocean and sea ice data, a dynamic and thermodynamic +sea-ice model has been coupled to the Massachusetts Institute of Technology +general circulation model (MITgcm). Ice mechanics follow a viscous plastic +rheology and the ice momentum equations are solved numerically using either +line successive relaxation (LSR) or elastic-viscous-plastic (EVP) dynamic +models. Ice thermodynamics are represented using either a zero-heat-capacity +formulation or a two-layer formulation that conserves enthalpy. The model +includes prognostic variables for snow and for sea-ice salinity. The above +sea ice model components were borrowed from current-generation climate models +but they were reformulated on an Arakawa C-grid in order to match the MITgcm +oceanic grid and they were modified in many ways to permit efficient and +accurate automatic differentiation. This paper describes the MITgcm sea ice +model; it presents example Arctic and Antarctic results from a realistic, +eddy-permitting, global ocean and sea-ice configuration; it compares B-grid +and C-grid dynamic solvers in a regional Arctic configuration; and it presents +example results from coupled ocean and sea-ice adjoint-model integrations. \end{abstract} \section{Introduction} \label{sec:intro} -more blabla - -\section{Model} -\label{sec:model} +The availability of an adjoint model as a powerful research +tool complementary to an ocean model was a major design +requirement early on in the development of the MIT general +circulation model (MITgcm) [Marshall et al. 1997a, +Marotzke et al. 1999, Adcroft et al. 2002]. It was recognized +that the adjoint permitted very efficient computation +of gradients of various scalar-valued model diagnostics, +norms or, generally, objective functions with respect +to external or independent parameters. Such gradients +arise in at least two major contexts. If the objective function +is the sum of squared model vs. obervation differences +weighted by e.g. the inverse error covariances, the gradient +of the objective function can be used to optimize this measure +of model vs. data misfit in a least-squares sense. One +is then solving a problem of statistical state estimation. +If the objective function is a key oceanographic quantity +such as meridional heat or volume transport, ocean heat +content or mean surface temperature index, the gradient +provides a complete set of sensitivities of this quantity +with respect to all independent variables simultaneously. + +References to existing sea-ice adjoint models, explaining that they are either +for simplified configurations, for ice-only studies, or for short-duration +studies to motivate the present work. Traditionally, probably for historical reasons and the ease of treating the Coriolis term, most standard sea-ice models are discretized on Arakawa-B-grids \citep[e.g.,][]{hibler79, harder99, - kreyscher00, zhang98, hunke97}. From the perspective of coupling a +kreyscher00, zhang98, hunke97}. From the perspective of coupling a sea ice-model to a C-grid ocean model, the exchange of fluxes of heat and fresh-water pose no difficulty for a B-grid sea-ice model \citep[e.g.,][]{timmermann02a}. However, surface stress is defined at @@ -84,7 +111,7 @@ sea-ice model and a C-grid ocean model. While the smoothing implicitly associated with this interpolation may mask grid scale noise, it may in two-way coupling lead to a computational mode as will be shown. By -choosing a C-grid for the sea-ice model, we circumvene this difficulty +choosing a C-grid for the sea-ice model, we circumvent this difficulty altogether and render the stress coupling as consistent as the buoyancy coupling. @@ -92,39 +119,58 @@ straits. In the limit of only one grid cell between coasts there is no flux allowed for a B-grid (with no-slip lateral boundary counditions), whereas the C-grid formulation allows a flux of sea-ice through this -passage for all types of lateral boundary conditions. We (will) +passage for all types of lateral boundary conditions. We demonstrate this effect in the Candian archipelago. +Talk about problems that make the sea-ice-ocean code very sensitive and +changes in the code that reduce these sensitivities. + +This paper describes the MITgcm sea ice +model; it presents example Arctic and Antarctic results from a realistic, +eddy-permitting, global ocean and sea-ice configuration; it compares B-grid +and C-grid dynamic solvers in a regional Arctic configuration; and it presents +example results from coupled ocean and sea-ice adjoint-model integrations. + +\section{Model} +\label{sec:model} + \subsection{Dynamics} \label{sec:dynamics} -The momentum equations of the sea-ice model are standard with +The momentum equation of the sea-ice model is \begin{equation} \label{eq:momseaice} m \frac{D\vek{u}}{Dt} = -mf\vek{k}\times\vek{u} + \vtau_{air} + - \vtau_{ocean} - m \nabla{\phi(0)} + \vek{F}, + \vtau_{ocean} - mg \nabla{\phi(0)} + \vek{F}, \end{equation} -where $\vek{u} = u\vek{i}+v\vek{j}$ is the ice velocity vectory, $m$ -the ice mass per unit area, $f$ the Coriolis parameter, $g$ is the -gravity accelation, $\nabla\phi$ is the gradient (tilt) of the sea -surface height potential beneath the ice. $\phi$ is the sum of -atmpheric pressure $p_{a}$ and loading due to ice and snow -$(m_{i}+m_{s})g$. $\vtau_{air}$ and $\vtau_{ocean}$ are the wind and -ice-ocean stresses, respectively. $\vek{F}$ is the interaction force -and $\vek{i}$, $\vek{j}$, and $\vek{k}$ are the unit vectors in the -$x$, $y$, and $z$ directions. Advection of sea-ice momentum is -neglected. The wind and ice-ocean stress terms are given by +where $m=m_{i}+m_{s}$ is the ice and snow mass per unit area; +$\vek{u}=u\vek{i}+v\vek{j}$ is the ice velocity vector; +$\vek{i}$, $\vek{j}$, and $\vek{k}$ are unit vectors in the $x$, $y$, and $z$ +directions, respectively; +$f$ is the Coriolis parameter; +$\vtau_{air}$ and $\vtau_{ocean}$ are the wind-ice and ocean-ice stresses, +respectively; +$g$ is the gravity accelation; +$\nabla\phi(0)$ is the gradient (or tilt) of the sea surface height; +$\phi(0)$ is the sea surface height potential in response to ocean dynamics +and to atmospheric pressure loading; +and $\vek{F}=\nabla\cdot\sigma$ is the divergence of the internal ice stress +tensor $\sigma_{ij}$. +When using the rescaled vertical coordinate system, z$^\ast$, of +\citet{cam08}, $\phi(0)$ also includes a term due to snow and ice loading, $mg$. +Advection of sea-ice momentum is neglected. The wind and ice-ocean stress +terms are given by \begin{align*} - \vtau_{air} =& \rho_{air} |\vek{U}_{air}|R_{air}(\vek{U}_{air}) \\ - \vtau_{ocean} =& \rho_{ocean} |\vek{U}_{ocean}-\vek{u}| + \vtau_{air} = & \rho_{air} C_{air} |\vek{U}_{air} -\vek{u}| + R_{air} (\vek{U}_{air} -\vek{u}), \\ + \vtau_{ocean} = & \rho_{ocean}C_{ocean} |\vek{U}_{ocean}-\vek{u}| R_{ocean}(\vek{U}_{ocean}-\vek{u}), \\ \end{align*} where $\vek{U}_{air/ocean}$ are the surface winds of the atmosphere -and surface currents of the ocean, respectively. $C_{air/ocean}$ are -air and ocean drag coefficients, $\rho_{air/ocean}$ reference -densities, and $R_{air/ocean}$ rotation matrices that act on the -wind/current vectors. $\vek{F} = \nabla\cdot\sigma$ is the divergence -of the interal stress tensor $\sigma_{ij}$. +and surface currents of the ocean, respectively; $C_{air/ocean}$ are +air and ocean drag coefficients; $\rho_{air/ocean}$ are reference +densities; and $R_{air/ocean}$ are rotation matrices that act on the +wind/current vectors. For an isotropic system this stress tensor can be related to the ice strain rate and strength by a nonlinear viscous-plastic (VP) @@ -168,7 +214,7 @@ $\Delta_{\min}=10^{-11}\text{\,s}^{-1}$ (for numerical reasons) and a maximum $\zeta_{\max} = P_{\max}/\Delta^*$, where $\Delta^*=(5\times10^{12}/2\times10^4)\text{\,s}^{-1}$. For stress -tensor compuation the replacement pressure $P = 2\,\Delta\zeta$ +tensor computation the replacement pressure $P = 2\,\Delta\zeta$ \citep{hibler95} is used so that the stress state always lies on the elliptic yield curve by definition. @@ -192,9 +238,9 @@ treated explicitly. \citet{hunke97}'s introduced an elastic contribution to the strain -rate elatic-viscous-plastic in order to regularize +rate elastic-viscous-plastic in order to regularize Eq.\refeq{vpequation} in such a way that the resulting -elatic-viscous-plastic (EVP) and VP models are identical at steady +elastic-viscous-plastic (EVP) and VP models are identical at steady state, \begin{equation} \label{eq:evpequation} @@ -220,7 +266,7 @@ \dot{\epsilon}_{11}+\dot{\epsilon}_{22}$, and the horizontal tension and shearing strain rates, $D_T = \dot{\epsilon}_{11}-\dot{\epsilon}_{22}$ and $D_S = -2\dot{\epsilon}_{12}$, respectively and using the above abbreviations, +2\dot{\epsilon}_{12}$, respectively, and using the above abbreviations, the equations can be written as: \begin{align} \label{eq:evpstresstensor1} @@ -250,8 +296,8 @@ $P$ at vorticity points. For a general curvilinear grid, one needs in principle to take metric -terms into account that arise in the transformation a curvilinear grid -on the sphere. However, for now we can neglect these metric terms +terms into account that arise in the transformation of a curvilinear grid +on the sphere. For now, however, we can neglect these metric terms because they are very small on the cubed sphere grids used in this paper; in particular, only near the edges of the cubed sphere grid, we expect them to be non-zero, but these edges are at approximately @@ -260,7 +306,7 @@ cartesian. However, for last-glacial-maximum or snowball-earth-like simulations the question of metric terms needs to be reconsidered. Either, one includes these terms as in \citet{zhang03}, or one finds a -vector-invariant formulation fo the sea-ice internal stress term that +vector-invariant formulation for the sea-ice internal stress term that does not require any metric terms, as it is done in the ocean dynamics of the MITgcm \citep{adcroft04:_cubed_sphere}. @@ -367,13 +413,18 @@ \subsection{Arctic Domain with Open Boundaries} \label{sec:arctic} -The Arctic domain of integration is illustrated in Fig.~\ref{???}. It +The Arctic domain of integration is illustrated in Fig.~\ref{fig:arctic1}. It is carved out from, and obtains open boundary conditions from, the global cubed-sphere configuration of the Estimating the Circulation and Climate of the Ocean, Phase II (ECCO2) project \citet{menemenlis05}. The domain size is 420 by 384 grid boxes horizontally with mean horizontal grid spacing of 18 km. +\begin{figure} +%\centerline{{\includegraphics*[width=0.44\linewidth]{\fpath/arctic1.eps}}} +\caption{Bathymetry of Arctic Domain.\label{fig:arctic1}} +\end{figure} + There are 50 vertical levels ranging in thickness from 10 m near the surface to approximately 450 m at a maximum model depth of 6150 m. Bathymetry is from the National Geophysical Data Center (NGDC) 2-minute gridded global relief