--- MITgcm_contrib/articles/ceaice/ceaice.tex 2008/02/25 16:50:56 1.10 +++ MITgcm_contrib/articles/ceaice/ceaice.tex 2008/02/28 16:34:56 1.19 @@ -1,4 +1,4 @@ -% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.10 2008/02/25 16:50:56 dimitri Exp $ +% $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/articles/ceaice/ceaice.tex,v 1.19 2008/02/28 16:34:56 mlosch Exp $ % $Name: $ \documentclass[12pt]{article} @@ -12,6 +12,8 @@ \newcommand\bmmax{10} \newcommand\hmmax{10} \usepackage{bm} +\usepackage{url} + % math abbreviations \newcommand{\vek}[1]{\ensuremath{\mathbf{#1}}} \newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}} @@ -51,772 +53,29 @@ \maketitle -\begin{abstract} +\input{ceaice_abstract.tex} + +\input{ceaice_intro.tex} + +\input{ceaice_model.tex} -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. - -\end{abstract} - -\section{Introduction} -\label{sec:intro} - -more blabla - -\section{Model} -\label{sec:model} - -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 -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 -velocities points and thus needs to be interpolated between a B-grid -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 -altogether and render the stress coupling as consistent as the -buoyancy coupling. - -A further advantage of the C-grid formulation is apparent in narrow -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) -demonstrate this effect in the Candian archipelago. - -\subsection{Dynamics} -\label{sec:dynamics} - -The momentum equations of the sea-ice model are standard with -\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}, -\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 -\begin{align*} - \vtau_{air} =& \rho_{air} |\vek{U}_{air}|R_{air}(\vek{U}_{air}) \\ - \vtau_{ocean} =& \rho_{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}$. - -For an isotropic system this stress tensor can be related to the ice -strain rate and strength by a nonlinear viscous-plastic (VP) -constitutive law \citep{hibler79, zhang98}: -\begin{equation} - \label{eq:vpequation} - \sigma_{ij}=2\eta(\dot{\epsilon}_{ij},P)\dot{\epsilon}_{ij} - + \left[\zeta(\dot{\epsilon}_{ij},P) - - \eta(\dot{\epsilon}_{ij},P)\right]\dot{\epsilon}_{kk}\delta_{ij} - - \frac{P}{2}\delta_{ij}. -\end{equation} -The ice strain rate is given by -\begin{equation*} - \dot{\epsilon}_{ij} = \frac{1}{2}\left( - \frac{\partial{u_{i}}}{\partial{x_{j}}} + - \frac{\partial{u_{j}}}{\partial{x_{i}}}\right). -\end{equation*} -The maximum ice pressure $P_{\max}$, a measure of ice strength, depends on -both thickness $h$ and compactness (concentration) $c$: -\begin{equation} - P_{\max} = P^{*}c\,h\,e^{[C^{*}\cdot(1-c)]}, -\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 -invariants and ice strength such that the principal components of the -stress lie on an elliptical yield curve with the ratio of major to -minor axis $e$ equal to $2$; they are given by: -\begin{align*} - \zeta =& \min\left(\frac{P_{\max}}{2\max(\Delta,\Delta_{\min})}, - \zeta_{\max}\right) \\ - \eta =& \frac{\zeta}{e^2} \\ - \intertext{with the abbreviation} - \Delta = & \left[ - \left(\dot{\epsilon}_{11}^2+\dot{\epsilon}_{22}^2\right) - (1+e^{-2}) + 4e^{-2}\dot{\epsilon}_{12}^2 + - 2\dot{\epsilon}_{11}\dot{\epsilon}_{22} (1-e^{-2}) - \right]^{-\frac{1}{2}} -\end{align*} -The bulk viscosities are bounded above by imposing both a minimum -$\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$ -\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, -is limited by the explicit treatment of the Coriolis term. The -explicit treatment of the Coriolis term does not represent a severe -limitation because it restricts the time step to approximately the -same length as in the ocean model where the Coriolis term is also -treated explicitly. - -\citet{hunke97}'s introduced an elastic contribution to the strain -rate elatic-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 -state, -\begin{equation} - \label{eq:evpequation} - \frac{1}{E}\frac{\partial\sigma_{ij}}{\partial{t}} + - \frac{1}{2\eta}\sigma_{ij} - + \frac{\eta - \zeta}{4\zeta\eta}\sigma_{kk}\delta_{ij} - + \frac{P}{4\zeta}\delta_{ij} - = \dot{\epsilon}_{ij}. -\end{equation} -%In the EVP model, equations for the components of the stress tensor -%$\sigma_{ij}$ are solved explicitly. Both model formulations will be -%used and compared the present sea-ice model study. -The EVP-model uses an explicit time stepping scheme with a short -timestep. According to the recommendation of \citet{hunke97}, the -EVP-model is stepped forward in time 120 times within the physical -ocean model time step (although this parameter is under debate), to -allow for elastic waves to disappear. Because the scheme does not -require a matrix inversion it is fast in spite of the small timestep -\citep{hunke97}. For completeness, we repeat the equations for the -components of the stress tensor $\sigma_{1} = -\sigma_{11}+\sigma_{22}$, $\sigma_{2}= \sigma_{11}-\sigma_{22}$, and -$\sigma_{12}$. Introducing the divergence $D_D = -\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, -the equations can be written as: -\begin{align} - \label{eq:evpstresstensor1} - \frac{\partial\sigma_{1}}{\partial{t}} + \frac{\sigma_{1}}{2T} + - \frac{P}{2T} &= \frac{P}{2T\Delta} D_D \\ - \label{eq:evpstresstensor2} - \frac{\partial\sigma_{2}}{\partial{t}} + \frac{\sigma_{2} e^{2}}{2T} - &= \frac{P}{2T\Delta} D_T \\ - \label{eq:evpstresstensor12} - \frac{\partial\sigma_{12}}{\partial{t}} + \frac{\sigma_{12} e^{2}}{2T} - &= \frac{P}{4T\Delta} D_S -\end{align} -Here, the elastic parameter $E$ is redefined in terms of a damping timescale -$T$ for elastic waves \[E=\frac{\zeta}{T}.\] -$T=E_{0}\Delta{t}$ with the tunable parameter $E_0<1$ and -the external (long) timestep $\Delta{t}$. \citet{hunke97} recommend -$E_{0} = \frac{1}{3}$. - -For details of the spatial discretization, the reader is referred to -\citet{zhang98, zhang03}. Our discretization differs only (but -importantly) in the underlying grid, namely the Arakawa C-grid, but is -otherwise straightforward. The EVP model in particular is discretized -naturally on the C-grid with $\sigma_{1}$ and $\sigma_{2}$ on the -center points and $\sigma_{12}$ on the corner (or vorticity) points of -the grid. With this choice all derivatives are discretized as central -differences and averaging is only involved in computing $\Delta$ and -$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 -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 -35\degS\ or 35\degN\ which are never covered by sea-ice in our -simulations. Everywhere else the coordinate system is locally nearly -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 -does not require any metric terms, as it is done in the ocean dynamics -of the MITgcm \citep{adcroft04:_cubed_sphere}. - -Moving sea ice exerts a stress on the ocean which is the opposite of -the stress $\vtau_{ocean}$ in Eq.\refeq{momseaice}. This stess is -applied directly to the surface layer of the ocean model. An -alternative ocean stress formulation is given by \citet{hibler87}. -Rather than applying $\vtau_{ocean}$ directly, the stress is derived -from integrating over the ice thickness to the bottom of the oceanic -surface layer. In the resulting equation for the \emph{combined} -ocean-ice momentum, the interfacial stress cancels and the total -stress appears as the sum of windstress and divergence of internal ice -stresses: $\delta(z) (\vtau_{air} + \vek{F})/\rho_0$, \citep[see also -Eq.\,2 of][]{hibler87}. The disadvantage of this formulation is that -now the velocity in the surface layer of the ocean that is used to -advect tracers, is really an average over the ocean surface -velocity and the ice velocity leading to an inconsistency as the ice -temperature and salinity are different from the oceanic variables. - -Sea ice distributions are characterized by sharp gradients and edges. -For this reason, we employ a positive 3rd-order advection scheme -\citep{hundsdorfer94} for the thermodynamic variables discussed in the -next section. - -\subparagraph{boundary conditions: no-slip, free-slip, half-slip} - -\begin{itemize} -\item transition from existing B-Grid to C-Grid -\item boundary conditions: no-slip, free-slip, half-slip -\item fancy (multi dimensional) advection schemes -\item VP vs.\ EVP \citep{hunke97} -\item ocean stress formulation \citep{hibler87} -\end{itemize} - -\subsection{Thermodynamics} -\label{sec:thermodynamics} - -In the original formulation the sea ice model \citep{menemenlis05} -uses simple thermodynamics following the appendix of -\citet{semtner76}. This formulation does not allow storage of heat -(heat capacity of ice is zero, and this type of model is often refered -to as a ``zero-layer'' model). Upward heat flux is parameterized -assuming a linear temperature profile and together with a constant ice -conductivity. It is expressed as $(K/h)(T_{w}-T_{0})$, where $K$ is -the ice conductivity, $h$ the ice thickness, and $T_{w}-T_{0}$ the -difference between water and ice surface temperatures. The surface -heat budget is computed in a similar way to that of -\citet{parkinson79} and \citet{manabe79}. - -There is considerable doubt about the reliability of such a simple -thermodynamic model---\citet{semtner84} found significant errors in -phase (one month lead) and amplitude ($\approx$50\%\,overestimate) in -such models---, so that today many sea ice models employ more complex -thermodynamics. A popular thermodynamics model of \citet{winton00} is -based on the 3-layer model of \citet{semtner76} and treats brine -content by means of enthalphy conservation. This model requires in -addition to ice-thickness and compactness (fractional area) additional -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} -\item no-slip vs. free-slip for both lsr and evp; - "diagnostics" to look at and use for comparison - \begin{itemize} - \item ice transport through Fram Strait/Denmark Strait/Davis - Strait/Bering strait (these are general) - \item ice transport through narrow passages, e.g.\ Nares-Strait - \end{itemize} -\item compare different advection schemes (if lsr turns out to be more - effective, then with lsr otherwise I prefer evp), eg. - \begin{itemize} - \item default 2nd-order with diff1=0.002 - \item 1st-order upwind with diff1=0. - \item DST3FL (SEAICEadvScheme=33 with diff1=0., should work, works for me) - \item 2nd-order wit flux limiter (SEAICEadvScheme=77 with diff1=0.) - \end{itemize} - That should be enough. Here, total ice mass and location of ice edge - is interesting. However, this comparison can be done in an idealized - domain, may not require full Arctic Domain? -\item -Do a little study on the parameters of LSR and EVP -\begin{enumerate} -\item convergence of LSR, how many iterations do you need to get a - true elliptic yield curve -\item vary deltaTevp and the relaxation parameter for EVP and see when - the EVP solution breaks down (relative to the forcing time scale). - For this, it is essential that the evp solver gives use "stripeless" - solutions, that is your dtevp = 1sec solutions/or 10sec solutions - with SEAICE\_evpDampC = 615. -\end{enumerate} -\end{itemize} - -\section{Forward sensitivity experiments} -\label{sec:forward} - -A second series of forward sensitivity experiments have been carried out on an -Arctic Ocean domain with open boundaries. Once again the objective is to -compare the old B-grid LSR dynamic solver with the new C-grid LSR and EVP -solvers. One additional experiment is carried out to illustrate the -differences between the two main options for sea ice thermodynamics in the MITgcm. - -\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 -\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 -\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 \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 \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, -0.76, 0.94, and 0.8. - -\begin{itemize} -\item Configuration -\item OBCS from cube -\item forcing -\item 1/2 and full resolution -\item with a few JFM figs from C-grid LSR no slip - ice transport through Canadian Archipelago - thickness distribution - ice velocity and transport -\end{itemize} - -\begin{itemize} -\item Arctic configuration -\item ice transport through straits and near boundaries -\item focus on narrow straits in the Canadian Archipelago -\end{itemize} - -\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 -\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 - thickness distribution differences - ice velocity and transport differences -\end{itemize} - -We anticipate small differences between the different models due to: -\begin{itemize} -\item advection schemes: along the ice-edge and regions with large - gradients -\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} - -\section{Adjoint sensiivities of the MITsim} - -\subsection{The adjoint of MITsim} - -The ability to generate tangent linear and adjoint model components -of the MITsim has been a main design task. -For the ocean the adjoint capability has proven to be an -invaluable tool for sensitivity analysis as well as state estimation. -In short, the adjoint enables very efficient computation of the gradient -of scalar-valued model diagnostics (called cost function or objective function) -with respect to many model "variables". -These variables can be two- or three-dimensional fields of initial -conditions, model parameters such as mixing coefficients, or -time-varying surface or lateral (open) boundary conditions. -When combined, these variables span a potentially high-dimensional -(e.g. O(10$^8$)) so-called control space. Performing parameter perturbations -to assess model sensitivities quickly becomes prohibitive at these scales. -Alternatively, (time-varying) sensitivities of the objective function -to any element of the control space can be computed very efficiently in -one single adjoint -model integration, provided an efficient adjoint model is available. - -[REFERENCES] - - -The adjoint operator (ADM) is the transpose of the tangent linear operator (TLM) -of the full (in general nonlinear) forward model, i.e. the MITsim. -The TLM maps perturbations of elements of the control space -(e.g. initial ice thickness distribution) -via the model Jacobian -to a perturbation in the objective function -(e.g. sea-ice export at the end of the integration interval). -\textit{Tangent} linearity ensures that the derivatives are evaluated -with respect to the underlying model trajectory at each point in time. -This is crucial for nonlinear trajectories and the presence of different -regimes (e.g. effect of the seaice growth term at or away from the -freezing point of the ocean surface). -Ensuring tangent linearity can be easily achieved by integrating -the full model in sync with the TLM to provide the underlying model state. -Ensuring \textit{tangent} adjoints is equally crucial, but much more -difficult to achieve because of the reverse nature of the integration: -the adjoint accumulates sensitivities backward in time, -starting from a unit perturbation of the objective function. -The adjoint model requires the model state in reverse order. -This presents one of the major complications in deriving an -exact, i.e. \textit{tangent} adjoint model. - -Following closely the development and maintenance of TLM and ADM -components of the MITgcm we have relied heavily on the -autmomatic differentiation (AD) tool -"Transformation of Algorithms in Fortran" (TAF) -developed by Fastopt (Giering and Kaminski, 1998) -to derive TLM and ADM code of the MITsim. -Briefly, the nonlinear parent model is fed to the AD tool which produces -derivative code for the specified control space and objective function. -Following this approach has (apart from its evident success) -several advantages: -(1) the adjoint model is the exact adjoint operator of the parent model, -(2) the adjoint model can be kept up to date with respect to ongoing -development of the parent model, and adjustments to the parent model -to extend the automatically generated adjoint are incremental changes -only, rather than extensive re-developments, -(3) the parallel structure of the parent model is preserved -by the adjoint model, ensuring efficient use in high performance -computing environments. - -Some initial code adjustments are required to support dependency analysis -of the flow reversal and certain language limitations which may lead -to irreducible flow graphs (e.g. GOTO statements). -The problem of providing the required model state in reverse order -at the time of evaluating nonlinear or conditional -derivatives is solved via balancing -storing vs. recomputation of the model state in a multi-level -checkpointing loop. -Again, an initial code adjustment is required to support TAFs -checkpointing capability. -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. - -[HIGHLIGHT COUPLED NATURE OF THE ADJOINT!] - -\subsection{Special considerations} - -* growth term(?) - -* small active denominators - -* dynamic solver (implicit function theorem) - -* approximate adjoints - - -\subsection{An example: sensitivities of sea-ice export through Fram Strait} - -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 -\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 -North Atlantic current to sea-ice variability and export. -The horizontal resolution varies between XX and YY km -with 50 unevenly spaced vertical levels. -The adjoint models run efficiently on 80 processors -(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 -\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. - -\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. -\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, -the circulation around Svalbard, and ... - -\begin{figure}[t!] -\centerline{ -\subfigure[{\footnotesize -12 months}] -{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim072_cmax2.0E+02.eps}} -%\includegraphics*[width=.3\textwidth]{H_c.bin_res_100_lev1.pdf} -% -\subfigure[{\footnotesize -24 months}] -{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim145_cmax2.0E+02.eps}} -} - -\centerline{ -\subfigure[{\footnotesize --36 months}] -{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim218_cmax2.0E+02.eps}} -% -\subfigure[{\footnotesize --48 months}] -{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJheff_arc_lev1_tim292_cmax2.0E+02.eps}} -} -\caption{Sensitivity of sea-ice export through Fram Strait in December 2005 to -sea-ice thickness at various prior times. -\label{fig:4yradjheff}} -\end{figure} - - -\begin{figure}[t!] -\centerline{ -\subfigure[{\footnotesize -12 months}] -{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim072_cmax5.0E+01.eps}} -%\includegraphics*[width=.3\textwidth]{H_c.bin_res_100_lev1.pdf} -% -\subfigure[{\footnotesize -24 months}] -{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim145_cmax5.0E+01.eps}} -} - -\centerline{ -\subfigure[{\footnotesize --36 months}] -{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim218_cmax5.0E+01.eps}} -% -\subfigure[{\footnotesize --48 months}] -{\includegraphics*[width=0.44\linewidth]{\fpath/run_4yr_ADJtheta_arc_lev1_tim292_cmax5.0E+01.eps}} -} -\caption{Same as \reffig{4yradjheff} but for sea surface temperature -\label{fig:4yradjthetalev1}} -\end{figure} - - - -\section{Discussion and conclusion} -\label{sec:concl} - -The story of the paper could be: -B-grid ice model + C-grid ocean model does not work properly for us, -therefore C-grid ice model with advantages: -\begin{enumerate} -\item stress coupling simpler (no interpolation required) -\item different boundary conditions -\item advection schemes carry over trivially from main code -\end{enumerate} -LSR/EVP solutions are similar with appropriate bcs, evp parameters as -a function of forcing time scale (when does VP solution break -down). Same for LSR solver, provided that it works (o: -Which scheme is more efficient for the resolution/time stepping -parameters that we use here. What about other resolutions? +\input{ceaice_forward.tex} + +\input{ceaice_adjoint.tex} + +\input{ceaice_concl.tex} \paragraph{Acknowledgements} We thank Jinlun Zhang for providing the original B-grid code and many helpful discussions. ML thanks Elizabeth Hunke for multiple explanations. +This work is a contribution to Estimating the Circulation and Climate of the +Ocean, Phase II (ECCO2). The ECCO2 project (http://ecco2.org/) is sponsored +by the NASA Modeling Analysis and Prediction (MAP) program. D. Menemenlis +carried out this work at the Jet Propulsion Laboratory, California Institute +of Technology under contract with the National Aeronautics and Space +Administration. + \bibliography{bib/journal_abrvs,bib/seaice,bib/genocean,bib/maths,bib/mitgcmuv,bib/fram} \end{document}