--- MITgcm_contrib/articles/ceaice/ceaice.tex 2008/01/21 08:06:00 1.9 +++ 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.9 2008/01/21 08:06:00 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,21 +52,58 @@ \maketitle \begin{abstract} - Some blabla +As part of ongoing efforts to obtain a best possible synthesis of most +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 @@ -74,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. @@ -82,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) @@ -158,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. @@ -182,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} @@ -210,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} @@ -240,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 @@ -250,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}. @@ -311,142 +367,6 @@ state variables to be advected by ice velocities, namely enthalphy of the two ice layers and the thickness of the overlying snow layer. -\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?]} \subsection{C-grid} \begin{itemize} @@ -493,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