Index: /issm/trunk/bugs/LarouretalJGR11.tex =================================================================== --- /issm/trunk/bugs/LarouretalJGR11.tex (revision 7954) +++ /issm/trunk/bugs/LarouretalJGR11.tex (revision 7954) @@ -0,0 +1,920 @@ +%Preamble{{{1 +\documentclass[jgrga]{agutex} + +%pictures +\usepackage[dvips]{graphicx} +\setkeys{Gin}{draft=false} + +%Additional packages +\usepackage{amsmath} %for the multline environment +\usepackage{bm} % maths letters, real = IR + +%new commands +\newcommand*{\MAT}[1]{\ensuremath{\mathbf{#1}} } %Matrice name (Bold) +\newcommand*{\VEC}[1]{\ensuremath{\bm{#1}} } %Vector name (Bold Italic) +\newcommand*{\SET}[1]{\ensuremath{\mathcal{#1}} } %Set, Space (italic) +\newcommand*{\R}[0]{{\mathbb R} } %real set IR +\newcommand*{\abs}[1]{\mathopen| #1 \mathclose| } %absolute value |x| + +%}}} +%Authors and affiliations {{{1 +\authorrunninghead{Larour et al.} \titlerunninghead{Ice Sheet System Model \copyright2011. All rights reserved} + +%\authoraddr{E. Larour} Mechanical Division, Thermal and Cryogenics Section, 354, Jet Propulsion +%Laboratory, MS 157-316, 4800 Oak Grove Drive, Pasadena, CA 91109, USA.% (eric.larour@jpl.nasa.gov) +%\authoraddr{M. Morlighem Communication Tracking and Radar Division, Radar Science and Engineering +%Section, 334, Jet Propulsion Laboratory, MS 300-319, 4800 Oak Grove Drive, Pasadena, CA 91109, USA. +%(mathieu.morlighem@jpl.nasa.gov)} +%\authoraddr{E. Rignot} University of California, Irvine, Department of Earth System Science, Croul +%Hall, Irvine, CA 92697-3100, USA.% (erignot@uci.edu) +%\authoraddr{H. Seroussi} Communication Tracking and Radar Division, Radar Science and Engineering +%Section, 334, Jet Propulsion Laboratory, MS 300-319, 4800 Oak Grove Drive, Pasadena, CA 91109, +%USA.% (helene.seroussi@jpl.nasa.gov) + +\begin{document} + +\title{Ice Sheet System Model: a new 2-D/3-D Higer-Order/Full Stokes Ice Flow Model} + +\authors{E. Larour\altaffilmark{1}, H. Seroussi\altaffilmark{1,3}, M. Morlighem\altaffilmark{1,3}, E. Rignot \altaffilmark{2,1}} + +\altaffiltext{1}{Jet Propulsion Laboratory - California Institute of technology, 4800 Oak Grove +Drive MS 300-227, Pasadena, CA 91109-8099, USA.} +\altaffiltext{2}{University of California, Irvine, Department of Earth System Science, Croul Hall, +Irvine, CA 92697-3100, USA.} +\altaffiltext{3}{Laboratoire MSSMat - CNRS U.M.R 8579, \'Ecole Centrale Paris, Grande Voie des +Vignes, 92295 Ch\^atenay-Malabry Cedex, FRANCE.} +%}}} +\begin{abstract}%{{{1 + Most ice sheet numerical models up to present have been mainly relying on the shallow ice + approximation, which is not adequate to represent ice streams, their fluctuations and ongoing + changes in mass balance of ice sheets. Here we present a new finite element, thermomechanical + numerical model of Antarctica named ISSM (Ice Sheet System Model) that includes high order stress + terms, high mesh resolution at the coast and that captures ice stream dynamics and relies on a + heavily parallelized architecture. One of the main model parameters, basal drag, is constrained + along the ice sheet periphery by observations of surface velocity from satellite radar + interferometry using inverse control methods. Our results demonstrates that reliable modeling of + an entire ice continent at high resolution is now possible and provides a pathway for improving + projections of ice sheet evolution in a warming climate. +\end{abstract} %}}} + +\begin{article} +\section{Introduction} %{{{1 +Modeling of ice sheets at a continental scale (Antarctic and Greenland ice sheets) is a major +challenge that must be addressed in order to estimate future sea level rise with a reasonable degree +of confidence. Such modeling will have to account for different sets of physics (multi-model +approach), different geometrical scales (multi-resolution approach) and be computationaly scalable. + +Most ice sheet models were initially based on simplified formulations for ice flow and temperature. +The Shallow Ice Approximation, or SIA,... +%}}} +\section{Ice Flow thermodynamics in ISSM} %{{{1 +Ice thermodynamics is based on the classical conservation laws: momentum balance, mass balance and +energy balance, along with a material constitutive law, and boundary conditions. This section +describes each of these components in detail. + +\subsection{Mechanical model}%{{{2 +ISSM implements several types of ice flow models, which are based on a gradual simplification of the +Full-Stokes (FS) stress-equilibrium equations, according to a set of physical assumptions. The +FS equations are: +\begin{equation} + \nabla \cdot {\bm \sigma} + \rho \;{\bf g} = {\bf 0} + \label{eq:momentumstokes} +\end{equation} + +\begin{equation} + \mbox{Tr} \left( \dot{\bm \varepsilon}\right) + = 0 + \label{eq:incompressibility} +\end{equation} +where $\nabla \cdot {\bm \sigma}$ is the divergence vector of the stress tensor ${\bm \sigma}$, +$\rho$ is the ice density ${\bf g}$ the acceleration due to gravity, $ \dot{\bm \varepsilon}$ is the +strain rate tensor and $\mbox{Tr}$ is the trace operator. + +Eq. \ref{eq:momentumstokes} expresses the balance of stresses, and Eq. \ref{eq:incompressibility} +the incompressibility of ice flow. In Eq. \ref{eq:momentumstokes}, acceleration is neglected, +following a scale analysis by \cite{Reist2005} that shows ice flow acceleration is negligible even +in the most extreme surges or stream flow that may occur in a glacier. This same scale analysis +shows that coriolis effect is also negligible. + +The material consitutive law describes the deformation of ice under stress. For incompressible +fluids, the constitutive law has the following form: +\begin{equation} + {\bm \sigma}'=2\mu\dot{\bm \varepsilon} + \label{eq:Viscosity} +\end{equation} + +where ${\bm \sigma'}={\bm \sigma} + P {\bf I} $ is the deviatoric stress tensor, ${\bf I}$ is the identity +matrix and $\mu$ is the viscosity. The ice viscosity is assumed to be non-linear and follows a +Norton-Hoff law \citep{Glen1955}. \begin{equation} + \mu=B \sigma_e^{\frac{1-n}{n}} + \label{eq:Glen} +\end{equation} +$B$ is the ice hardness, $n$ Glen's law coefficient and $\sigma_e$ the effective stress. $B$ is +highly temperature dependent. In ISSM it can either follow a classical Arrhenius temperature dependence: +\begin{equation} + B\left(T\right) = \left(A_0\left(\exp{\frac{-Q}{RT^\ast}}\right)\right)^{-1/n} +\end{equation} +where $A_0$ is the flow factor, $Q$ is the activation energy for ice creep, $R$ the universal gas +constant and $T^\ast=T - \beta \left(s-z\right)$ is the absolute temperature corrected for the +dependence of melting point on pressure. The values of these constants follow \cite{Payne2000}. +Alternatively, the ice hardness, $B$, can follow the temperature-dependent relation for B given by a +table in \cite{Paterson1994}. + +Using Eq. \ref{eq:Viscosity}, Eq. \ref{eq:momentumstokes} can be rewritten in terms of strain +rate tensor as follows: +\begin{equation} + \nabla \cdot \left( 2 \mu \dot{\bm\varepsilon} \right) - \nabla P + \rho \;{\bf g} = {\bf 0} + \label{eq:momentumstokes2} +\end{equation} + +Using Glens's flow law (Eq. \ref{eq:Glen}), Eq. \ref{eq:momentumstokes2} and Eq. \ref{eq:incompressibility} can then +be written in terms of velocity components: +\begin{multline} + \frac{\partial}{\partial x} \left( 2 \mu \frac{\partial u}{\partial x} \right) + +\frac{\partial}{\partial y} \left( \mu \frac{\partial u}{\partial y} +\mu \frac{\partial v}{\partial x} \right) + \\ + +\frac{\partial}{\partial z} \left( \mu \frac{\partial u}{\partial z} +\mu \frac{\partial w}{\partial x} \right) + -\frac{\partial P}{\partial x}=0 + \label{eq:Stokesx} +\end{multline} +\begin{multline} + \frac{\partial}{\partial x} \left( \mu \frac{\partial u}{\partial y} +\mu \frac{\partial v}{\partial x} \right) + +\frac{\partial}{\partial y} \left( 2 \mu \frac{\partial v}{\partial y} \right) + \\ + +\frac{\partial}{\partial z} \left( \mu \frac{\partial v}{\partial z} +\mu \frac{\partial w}{\partial y} \right) + -\frac{\partial P}{\partial y}=0 + \label{eq:Stokesy} +\end{multline} +\begin{multline} + \frac{\partial}{\partial x} \left( \mu \frac{\partial u}{\partial z} +\mu \frac{\partial w}{\partial x} \right) + +\frac{\partial}{\partial y} \left( \mu \frac{\partial v}{\partial z} +\mu \frac{\partial w}{\partial y} \right) + \\ + +\frac{\partial}{\partial z} \left( 2 \mu \frac{\partial w}{\partial z} \right) + -\frac{\partial P}{\partial z} - \rho g=0 + \label{eq:Stokesz} +\end{multline} +\begin{equation} + \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} + = 0 + \label{eq:Stokesp} +\end{equation} + +where ($u$,$v$,$w$) are the $x$, $y$ and $z$ components of the velocity vector $\textbf{v}$, in the +($x$,$y$,$z$) cartesian coordinate system, with $z$ in the vertical direction. Using this +formulation, pressure $P$ can be seen as a Lagrange multiplier that ensures the validity of the +compressibility/continuity equation Eq. \ref{eq:incompressibility}. Eq. +\ref{eq:Stokesx},\ref{eq:Stokesy},\ref{eq:Stokesz} and \ref{eq:Stokesp} are the basis of the +FS formulation in ISSM, which we solve for using the Finite Element Method. This model +solves for four degrees of freedom ($u$,$v$,$w$,$P$) on each vertex of the finite element mesh, and +is therefore computationally very demanding. + +%Pattyn +A simplified three-dimensional (3D) model from \cite{Blatter1995} and \cite{Pattyn2003} (BP) is derived +from FS by making two assumptions: 1) the horizontal gradients of vertical velocities are negligible +compared to the vertical gradients of horizontal velocities: +\begin{equation} + \dot{\varepsilon}_{xz} = \frac{1}{2} \frac{\partial u}{\partial z}; + \hspace{4em} + \dot{\varepsilon}_{yz} = \frac{1}{2} \frac{\partial v}{\partial z} +\end{equation} +and 2) the bridging effects \citep{Veen1989} are negligible, which reduces the third equation of the +momentum balance (Eq. \ref{eq:momentumstokes2}) to: +\begin{equation} + \frac{\partial}{\partial z} \left( 2 \mu \frac{\partial w}{\partial z} \right) - + \frac{\partial P}{\partial z} - \rho g =0 + \label{eq:vertBP} +\end{equation} +This reduces the FS equations to: +\begin{multline} + \frac{\partial}{\partial x} \left( 2 \mu \frac{\partial u}{\partial x} \right) + +\frac{\partial}{\partial y} \left( \mu \frac{\partial u}{\partial y} +\mu \frac{\partial v}{\partial x} \right) + \\ + +\frac{\partial}{\partial z} \left( \mu \frac{\partial u}{\partial z} \right) + = \rho g \frac{\partial s}{\partial x} +\end{multline} +\begin{multline} + \frac{\partial}{\partial x} \left( \mu \frac{\partial u}{\partial y} +\mu \frac{\partial v}{\partial x} \right) + +\frac{\partial}{\partial y} \left( 2 \mu \frac{\partial v}{\partial y} \right) + \\ + +\frac{\partial}{\partial z} \left( \mu \frac{\partial v}{\partial z} \right) + = \rho g \frac{\partial s}{\partial y} +\end{multline} +The equations of BP only involve two degrees of freedom ($u$ and $v$), as the vertical velocity, +$w$, has been decoupled from the initial system equations. + +%MacAyeal +The third model, Shelfy-stream approximation or SSA \citep{MacAyeal1989}, assumes in addition that +the vertical shear is negligible: +\begin{equation} \dot{\varepsilon}_{xz} = 0 ; + \hspace{6em} + \dot{\varepsilon}_{yz} = 0 +\end{equation} +This assumption reduces the equations to a 2D model, as $u$ and $v$ do not depend on depth $z$: +\begin{multline} + \frac{\partial}{\partial x} \left( + 4 H \bar{\mu} \frac{\partial u}{\partial x} + 2 H \bar{\mu} \frac{\partial v}{\partial y} \right) + \\ + + \frac{\partial}{\partial y} \left( + H \bar{\mu} \frac{\partial u}{\partial y} + H \bar{\mu} \frac{\partial v}{\partial x} \right) + = \rho g H \frac{\partial s}{\partial x} +\end{multline} +\begin{multline} + \frac{\partial}{\partial y} \left( + 4 H \bar{\mu} \frac{\partial v}{\partial y} + 2 H \bar{\mu} \frac{\partial u}{\partial x} \right) + \\ + + \frac{\partial}{\partial x} \left( + H \bar{\mu} \frac{\partial u}{\partial y} + H \bar{\mu} \frac{\partial v}{\partial x} \right) + = \rho g H \frac{\partial s}{\partial y} +\end{multline} +The vertical velocity, $w$, is deduced from the horizontal velocities, $u$ and $v$, using Eq. +\ref{eq:incompressibility} in BP and SSA. + +%SIA +The last model implemented in ISSM is the Shallow Ice Approximation (SIA) \citep{Hutter1983}. In +this 3D model, only the deviatoric stress components $\sigma_{xz}$ and $\sigma_{yz}$ are included. +The horizontal gradients of vertical velocity are also neglected compared to the vertical gradients +of horizontal velocities, so the equations are reduced to: +\begin{equation} + \begin{array}{c} + \displaystyle + \frac{\partial}{\partial z} \left( \mu \frac{\partial u}{\partial z}\right)= \rho g \frac{\partial s}{\partial x} \\ + \\ + \displaystyle + \frac{\partial}{\partial z} \left( \mu \frac{\partial v}{\partial z}\right)= \rho g \frac{\partial s}{\partial y} + \end{array} +\end{equation} + +The upper boundary condition of the ice flow model is a +stress-free surface. A friction law is applied at the ice-bedrock interface. The basal drag is +modeled following \cite{Paterson1994} written in a Coulomb-like law of friction: +\begin{equation} + {\bf \tau_b} = -k^2 N {\bf v_b} + \label{eq:frictionlaw} +\end{equation} +where ${\bf v_b}$ is the basal velocity vector tangential to the glacier base plane, $N$ is the +effective pressure on the glacier base, here equal to $N=\rho g h$, where $h$ is the height of +the ice sheet surface above buoyancy, ${\bf \tau_b}$ is the tangential component of the external +force, ${\bm \sigma} \cdot {\bf n}$, ${\bf n}$ is the outward pointing normal vector and $k^2$ is a +positive constant (i.e. stress opposes the motion). Water pressure is imposed on the ice-sea water +interface. The observed surface velocity is imposed on the remaining boundaries. +%}}} +\subsection{Thermal model}%{{{2 +The ice hardness $B$ of (\ref{eq:Glen}) is highly temperature dependent. A thermal model is hence +required to compute its value. The thermal equation is derived from the energy balance equation and +includes conduction-advection in the three directions: +\begin{equation} + \frac{\partial T}{\partial t} + = + -\textbf{v} \cdot \nabla T + +\frac{k_{th}}{ \rho c } \Delta T + +\frac{\Phi}{\rho c} + \label{eq:thermal} +\end{equation} +where $k_{th}$ is the ice thermal conductivity, $c$ the ice heat capacity, $\Phi$ is the heat +production term (deformational heating) and $\Delta$ is the Laplace operator. The temperature is +kept below the pressure melting point using an iterative penalty-based scheme described in the +following section. + +The following boundary conditions are applied: the surface temperature is the mean annual air +temperature from \cite{Giovinetto1990}. On grounded ice, we imposed a geothermal heat flux +\citep{Maule2005} and a frictional heat flux equal to ${\bf \tau_b} \cdot {\bf v_b}$. On the ice +shelf, basal drag is zero, thermal modeling is unresolved due to the complexity of ice-ocean +interaction and the ice hardness $B$ is inferred using an independent control method. Surface +topography is from a digital elevation model of Antarctica from \cite{Bamber2009}, a firn depth +correction from \cite{Broeke2008} and ice thickness is from \cite{Vaughan2006}. +%}}} +\subsection{Prognostic model}%{{{2 +\subsubsection{Thickness Evolution}%{{{3 +The mass conservation equation relates the ice thickness and the ice flux divergence to the surface and +basal mass balance. This equation is used to compute the ice thickness rate of change: +\begin{equation} + \frac{\partial H}{\partial t} + \nabla \cdot {H {\bf \overline{v}} } = \dot{M}_s - \dot{M}_b + \label{eq:massconservation} +\end{equation} +where $H$ is the ice thickness, ${\bf \overline{v}}=\left(\overline{u},\overline{v}\right)$ the +depth-averaged velocity, $\dot{M}_s$ the surface mass balance (m/yr in ice equivalent, positive for +accumulation, negative for ablation) and $\dot{M}_b$ the basal melting rate (m/yr in ice equivalent, +positive when melting, negative when freezing). The thickness is constrained on the inflow boundary. +The outflow boundary is left free (free-flux boundary condition) + +As this equation is hyperbolic, the standard Galerkin finite element method is generally unstable +and one must use either an artificial diffusion method or a discontinuous Galerkin finite element +method \citep{Johnson1984,Brezzi2004}. +%}}} +\subsubsection{Calving Front Dynamics}%{{{3 +The present-day ice front position is imposed as a mesh boundary and remains fixed with time. We +apply a free-flux boundary conditions: the melting-rate and the calving rate are assumed to be +exactly equal to the velocity at the front termini. We assume that any ice trying to flow beyond +the continental shelf break is removed from the system, no calving process is described. +%}}} +\subsubsection{Grounding Line Dynamics}%{{{3 + +Grounding line migration is a key control of ice flow dynamics +\citep{Schoof2007a,Schoof2007b,Nowicki2007,Nowicki2008,Durand2009,Durand2009a} which must be captured in order for +transient ice flow models to be realistic. In ISSM, we implemented a 3D grounding line migration criterion based +on the hydrostatic equilibrium line. At each time step of the transient ice flow solution, we check the following for every +vertex of our mesh: +\begin{itemize} +\item $b<=ba$ where $b$ is the new computed bed and $ba$ the bathymetry. If this condition is realized for a floating +vertex (i.e. on an ice shelf), we ground the vertex and force $b=ba$. +\item $b>b_{hydro}$ where $b_{hydro}$ is the bed position if the ice column at that location was in hydrostatic equilibrium. +$b_{hydro}=-di*H$ where $di$ is the ice density, and $H$ the new column thickness computed. If this condition is realized for +a grounded vertex (i.e. on the ice sheet), we unground the vertex and force $b=b_{hydro}$. +\end{itemize} +These two criterions are applied with the additional constraint that water cannot penetrate under an ice shelf cavity if there +is no open channel for the water to circulate. Practically, this means that we do not allow for vertices that are not connected +to the grounding line to unground. This condition rests on the assumption that warm water cannot sip under ice sheets, even when +the ice/bed interface is wet, and that the ice sheet represents an impermeable medium. + +This grounding line migration treatment is a first order approximation, and further work is currently being carried out +to implement a physically more realistic migration, similar to that found in \cite{Durand2009}. + + +%}}} +%}}} +\subsection{Control methods}%{{{2 +The basal drag coefficient $k$ in Eq. \ref{eq:frictionlaw} cannot be measured directly and is +inferred using a control method. We use a partial differential equations constrained optimization +algorithm similar to \cite{Vieli2003}, which consists in a gradient minimization of a cost function +that measures the misfit between observed ($u_{obs},v_{obs}$) and modeled ($u,v$) horizontal surface +velocities. The algorithm relies on the adjoint method, which calculates the gradient of the cost +function with respect to the unknown parameters. This cost function is usually taken as: +\begin{equation} + J=\int \int_{\Omega} \frac{1}{2} \left(u - u_{obs} \right)^2 + \frac{1}{2} \left(v - v_{obs} \right)^2 d\Omega + \label{eq:costfunction1} +\end{equation} +This cost function works better in areas of high-velocity than in slow moving regions because the +adjoint state (Lagrange multipliers vector) is larger where the velocity misfit $\left| u-u_{obs} +\right|$ is high, which occurs in regions of high speed. To minimize this effect, we introduce a new +cost function that measures the logarithm of the misfit: +\begin{equation} + J=\int\int_{\Omega} + \left( \log\left(\frac{\sqrt{u^2+v^2}+\varepsilon}{\sqrt{{u_{obs}}^2+{v_{obs} }^2}+\varepsilon} \right) \right)^2 d\Omega +\end{equation} +where $\varepsilon$ is a minimum velocity used to avoid zero velocities, and $\log$ is the natural +logarithm. This cost function enables a robust estimation of the basal drag coefficient after only a +few iterations over the entire model domain. +A Tikhonov regularization term, which penalizes the oscillations of the basal drag coefficient, $k$, is +added to this misfit to stabilize the inversion \citep{Vogel2002}. + +The finite element stiffness matrix is assumed to be independent of the velocity in order to have a +self-adjoint problem. This assumption is not correct as the viscosity, $\mu$, depends on the strain +rate, but it allows an easier calculation of the adjoint state for the three ice flow models and is +widely employed \citep{MacAyeal1993a}. + +This data assimilation technique has successfully been extended to BP and FS \citep{Morlighem2010}. +The major difference with SSA is that only the surface horizontal velocities are taken into account +in the cost function evaluation, while its gradient with respect to $k$ is computed at the base +only. + +At each iteration of the optimization procedure, we recalculate a thermo-mechanical equilibrium +solution and accordingly update the ice hardness $B$ on grounded ice to ensure consistency between +the ice flow and the viscosity $\mu$. +%}}} +%}}} +\section{Numerics} %{{{1 +\subsection{Finite Element Discretization} %{{{2 +ISSM uses the Continuous Galerkin Finite Element Method. Lagrange $P1$ elements (piecewise linear) +are used except fot the full-Stokes equation, for which \emph{Mini-elements} \citep{Gresho2000a} are +used in the finite element implementation of this model to fulfill the compatibility +Ladyzhenskaya-Babu\v{s}ka-Brezzi (LBB) condition. +The thermal model is solved using the Streamline Upwind Petrov-Galerkin (SUPG) \citep{Gresho2000} +formulation of the finite element method to prevent potential numerial oscillations due to dominant +advection terms. The temperature $T$ is forced to remain below the pressure melting point using an +iterative penalty-based scheme as in a contact problem \citep{Courant1943}. +\subsection{Mesh Refinement} %{{{2 + +When applying ice flow models at the continental scale, discretized systems of equations can +involve a large number of dofs, which can prove computationally challenging, even using state of +the art parallel technologies and clusters. In order to minimize the amount of dofs in the system, +different strategies can be implemented, which include among others mesh refinement. Mesh refinement +can be user-driven by relying on space-dependent mesh densities. However, this user-dependent +approach is highly subjective, and rarely used in the ISSM framework. Instead, we rely on metrics to +carry out static anisotropic adaptive meshing. + +Static anisotropic adaptive mesh refinement is based on the use of a metric to distort meshes to +better capture discretization errors. This approach limits the number of elements while maximizing +spatial resolution. The mesh refinement is done once during model setup (static), as opposed to +during the solution run (dynamic). We assume that the mesh remains constant during the model run. +This is certainly the case for short term runs, where major features of the ice flow do not advect +across the mesh (such as shear margin locations, or calving front and grounding line positions), but +for longer transient runs, this approach can be problematic, as major features of the ice flow will +not be properly captured after having advected across the mesh significantly. Further work is needed +to include dynamic anisotropic adaptive mesh refinement for these types of runs, but we believe +our static capability still improves significantly on past models which relied on structured meshes +that remain constant through time. + +Static anisotropic adaptive mesh refinement is based on the fact that interpolation-based a-priori +error estimates of a finite element $P1$ solution (piecewise linear) depend only on its Hessian +\citep{Habashi2000}. This is, provided that the solution is regular enough. For each element $E$ of +a mesh, the error of a $P1$ interpolated field, $u^h$, and the exact field, $u^{exact}$ is bounded +as follows: +\begin{equation} + \abs{u^{exact}\left(\VEC{x}\right)-u_E^h\left(\VEC{x}\right)} + \leq c_d h_E^2 + \sup_{\VEC{x}\in E}\abs{\MAT{H_u}\left(\VEC{x}\right)} +\end{equation} +where +$c_d$ is a constant that depends only on the space dimension, +$h_E$ is the characteristic length of the element, +$\MAT{H_u}\left(\VEC{x}\right)$ is the Hessian +matrix of $u\left(\VEC{x}\right)$, and $\abs{\MAT{H_u} \left(\VEC{x}\right)}$ its spectral norm. If we call $\VEC{e}$ +each edge of each element $E$ of the mesh ($\VEC{e}\in E_K$, where $E_K$ is the set of edges for element $E$), the previous equation becomes: +\begin{equation} + \max_{x \in E}\abs{u^{exact}(x)-u_E^h(x)} + \leq c_d \; + \max_{x \in E} \max_{\VEC{e} \in E_K} + \left(\VEC{e}\cdot\abs{\MAT{H_u}\left(\VEC{x}\right)}\VEC{e}\right) +\end{equation} +In practice the right hand side is delicate to estimate, since the maximum of the exact Hessian matrix +spectrum $\abs{\MAT{H_u}\left(\VEC{x}\right)}$ is unknown. To overcome this difficulty, we build +a metric tensor $\MAT{M}\left(E\right)$ for each element that verifies: +\begin{equation} + \forall \VEC{e} \in E \hspace{2em} + \max_{x \in E} + \left(\VEC{e}\cdot\abs{\MAT{H_u}\left(\VEC{x}\right)}\VEC{e}\right) + \leq + \VEC{e}\cdot\MAT{M}\left(E\right)\VEC{e} +\end{equation} +We also define an error estimate, $\varepsilon_E$, for every element of the mesh as follows: +\begin{equation} + \varepsilon_E= c_d + \max_{\VEC{e} \in E_K} + \left(\VEC{e}\cdot\MAT{M}\left(E\right)\VEC{e}\right) + \label{anisoerror} +\end{equation} +This means that for any element $E$ of the mesh, the interpolation error of $P1$ elements is +strictly proportional to the square of its longest edge in the metric $\MAT{M}\left(E\right)$. Therefore, +controlling of the element edges allows controllin of the interpolation error. + +We implemented an edge-based anisotropic mesh optimization methodology inspired by \cite{Frey2001} +and \cite{Hecht2006} to equi-distribute the interpolation error in each direction over each element +to control the global interpolation error. Let $\varepsilon$ be the maximum authorized error for all +elements of the mesh: +\begin{equation} + \forall \VEC{e} \in E_K \hspace{2em} + \varepsilon = c_d + \left(\VEC{e}\cdot\MAT{M}\left(E\right)\VEC{e}\right) +\end{equation} +all edges must verify: +\begin{equation} + \forall \VEC{e} \in E_K \hspace{2em} + \VEC{e}\cdot\overline{\MAT{M}}\left(E\right)\VEC{e}=1 +\end{equation} +where $\overline{\MAT{M}}\left(E\right)=\MAT{M}\left(E\right) c_d/\varepsilon$ is the final metric +tensor. This equation shows that the interpolation error is $\varepsilon$ if all edges of the mesh +have a length equal to $1$ in the metric $\overline{\MAT{M}}\left(E\right)$. + +In ISSM, we use the hessian matrix of surface velocities as a standard metric. This ensures that the +velocity are efficiently captured. Another metric that can be used is the thickness hessian matrix, +which tends to distort the mesh heavily when transitions in the surface or bed features occur, for +example along deep troughs in the bedrock (as observed on Jakobhsvan Gletscher), or near the +grounding line of ice shelves whose tributaries are steep glaciers (as observed on many glaciers of +the Antarctic Peninsula). Whenever possible, ISSM relies on observations to carry out adaptive mesh +refinement. Combined with the use of control methods, this approach ensures that model runs conform +to observations. In the absence of observations, adaptive mesh refinement has to be based on model +results. This involves iterating on the mesh refinement until a stable solution is reached. + +An example of static adaptive anisotropic mesh refinement is shown on Fig. \ref{fig:pigbamg} for +the Pine Island Glacier, Antarctica. Frame a) shows a 2700 elements mesh refined using InSAR +surface velocities from \cite{Rignot2008}. Shear margins are very well captured by the algorithm, +while the interior of the ice sheet, where ice flow deformation is low is meshed using coarser +sizes. In comparison, we show on frame b) a regular mesh (with elements of approximately equal area) +comprising the same number of elements, on which we also overlay InSAR surface velocities. +Comparison between the two demonstrate how, at equal mesh size, static anisotropic mesh refinement +is able to capture the physics of ice flow more efficiently, resulting in tremendous computational +gains. +%}}} +\subsection{Solvers} %{{{2 +%}}} +%}}} +\section{Verification and Validation} %{{{1 +Nighlty runs are performed every day, and several times a day in development phases. They consists +of almost a hundred test cases with a variety of configurations that ranges from square ice shelves +and ice sheets to simplified versions of Pine Island Glacier and Nioghalvfjerdsfjorden Glacier. 2000 +fields are tested against archived files and an HTML report is sent via email to all developers. +This ensures that no error is introduced in the software when adding new capabilities. + +We also performed several benckmarks, the ISMIP-HOM (Ice Sheet Model Intercomparison Project +for Higher-Order Models) \citep{Pattyn2008} that adresses two-dimensionnal and three dimensional ice +flows, and the EISMINT (European Ice Sheet Modelling INiTiative) \citep{Huybrechts1996,Payne2000} to +test the accuracy of our solutions. The results of these tests show a good agreement with the other +ice flow models. + +We present here the results of several tests. Fig. ? presents the results of the test +A of the ISMIP-HOM benchmark for the higher-order and full-Stokes models. This test involves an ice +slab frozen at the surface and flowing over a bumpy-bed. To create the periodic conditions necessary +for this test, we used the penalty method. The calculated surface velocity from ISSM is comparable +to the velocity found with other software for the benchmark (fig. ?). + +Test C of the ISMIP-HOM benchmark is a similar test except that the bedrock is flat but the basal +drag varies. Periodic boundary conditions are also applied here. Results found with ISSM for both +higher-order and full-Stokes models agree well with results found with other software for the +benchmark experiment (fig. ?) which gives confidence in the implementation of the +sliding law. + +The last test of the ISMIP-HOM benchmark presented here is the test F. This test is a prognostic +experiment: a slab of ice flows over a sloping bed. The initial bedrock and surface are parallel and +the ice thickness is 1000 m, but a gaussian bump is introduced in the center of the basal +topography. The free surface and velocity evolve until a steady-state solution is found. This +experiment is run with and without basal sliding and allows to validate the prognostic components of +the model. Results are presented in fig. ? for both higher-order and full-Stokes +models. + +All the tests are performed using triangular elements on a regular mesh with 50 layers on both +horizontal directions and 20 layers on the horizontal direction. +%}}} +\section{Performance} %{{{1 +Solver for MacAyeal and Pattyn, Stokes is still under development. +Residual Solver (GMRES) \citep{Saad1986} and the external direct solver MUMPS \citep{MUMPS:1,MUMPS:2}. Both solvers are + +Scalability: +Greenland, 2D, 125,000 dofs: 1mn 6s 50 control iterations +Greenland, 3D BP, 125,000 df. 2,500,000 dofs: 1h52 50 control iterations +Greenland, 3D FS, 125,000 df: 5,000,000 dofs 17h54 50 control iterations + + +%}}} +\section{Software Architecture and Management}%{{{1 + +ISSM relies mainly on the C language\citep{Kernighan1988} for the numerical implementation of finite +elements. For management of all objects, the C++ language \citep{Stroustrup1997} is used , which +gives strong emphasis to the use of polymorphic capabilities. This ensures ISSM is flexible (new +finite elements can be added rapidly) and scalable (C is a fast language, for which many compilers +and math kernels are well optimized, such as the Intel compiler and the Math Kernel Library). + +To improve ease of use, ISSM is hosted in MATLAB \citep{matlab}, a common scientific platform. The +C/C++ core is interfaced to the MATLAB environment using the MATLAB External API. This way, pre and +post-processing can be carried out using MATLAB modules (also called MEX functions). Those modules +behave as standard MATLAB modules, but they encapsulate intrinsic ISSM capabilities. The result is +a seamless integration of ISSM's ice flow model within the MATLAB platform. + +ISSM is at its core a parallel architecture. It can be run in serial mode, within MATLAB, but its +focus is mainly on massively parallel computations, run on large clusters such as the NASA Pleiades +cluster. When run on these types of clusters, ISSM relies on its C/C++ core compiled as a standalone +executable. Parallelism is achieved by using the Message Passing Interface +\citep{Gropp:1996:HPI,mpich-user}. This library is flexible enough to allow runs on distributed as +well as shared memory clusters. + +The numerics are implemented with the Portable, Extensible Toolkit for Scientific Computation +package (PETSc) \citep{petsc-efficient,petsc-user-ref,petsc-website}. This library defines objects +such as vectors, matrices and solvers, which are used directly in ISSM. These objects are abstracted +to hide serial vs parallel implementations, so that PETSC can be used in serial or parallel mode the +same way. PETSc also provides access to a wide array of direct and iterative solvers, propietary as +well as external, along with corresponding preconditioners (see section above). + +Apart from MATLAB and PETSC, ISSM relies on a wide array of external packages, of which we only list +here the most important ones: +\begin{itemize} + \item[--] TRIANGLE: Two-Dimensional Quality Mesh Generator and Delaunay Triangulator + \citep{Shewchuk1996,Shewchuk2002}. This package is used for creating unstructured + isotropic two-dimensional triangular meshes. + \item[--] VALGRIND: instrumentation framework for building dynamic analysis tools + \citep{valgrind-website,Nethercote2007,Nethercote2007a}. This package is used to debug memory + leaks in the code, as well as for optimizing memory management. + \item[--] METIS: Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and + Computing Fill-Reducing Orderings of Sparse Matrices \citep{Karypis1998}. This package is + used to partition objects such as elements and vertices across a cluster. The partitioning + scheme used results in partitions with equal numbers of elements for each cluster node. This + kind of partitioning ensures the best load balance, as it results in well-partitioned + stiffness matrices that optimize the most computationally intensive phase of the simulation, + which is the solver phase. +\end{itemize} + +ISSM uses a versioning system for hosting the code \citep{Subversion}, which allows multiple +developpers to synchronize development of new features, and regular maintenance of the code. This +system also allows nightly runs to be carried out, which validate the code numerics, based on the +EISMINT and ISMIP-HOM benchmarks (see section above), as well as approximately 250 verification and +validation tests. These nightly runs are carried out both in serial and parallel mode, which ensures +consistency throughout the code development process. + +Finally, ISSM comes with a full-fledged documentation, which can be compiled at installation, and is +synchronized with the ISSM website \citep{issm-website}. This documentation offers a user guide, +targeted mainly to active users of the code, and a theory guide, targeted to active developpers or +scientists interested in the numerics behind the finite elements. +%}}} +\section{Application to the Greenland ice sheet}%{{{1 + +Transient ice flow modeling at the continental scale is a challenge in terms of computational size +as well as the amount of physics that needs to be captured. In order to spin-up such a large scale ice flow +model, one must either reconstruct with a paleo-run the entire history of ice flow deformation, from a +know past configuration up to present time, or use inverse control methods to infer material properties +and/or basal boundary conditions at present time. Here, we show an application of ISSM to infer the basal +drag of Greenland, using three implemented ice flow formulations (FS, BP and SSA). Differences are shown +between the three formulations, and implications for spinning up the Greenland Ice Sheet are drawn. + +\subsection{Model Setup and Data}%{{{2 + +%}}} +\subsection{Results and discussion}%{{{2 +%}}} +%}}} +\section{Conclusions}%{{{1 + +%}}} +\begin{acknowledgments} %{{{1 + This work was performed at the Jet Propulsion Laboratory, California Institute of Technology, + at the Department of Earth System Science, University of California Irvine, + and at Laboratoire MSSMat, \'Ecole Centrale Paris, + under a contract with the National Aeronautics and Space Administration, Cryospheric Sciences Program and Modeling + Analysis and Prediction Program, and a contract with the Jet Propulsion Laboratory Research Technology and Development Program. + + Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced + Supercomputing (NAS) Division at Ames Research Center. + + The authors would like to thank Frank Pattyn and James Fishbaugh for providing the results and + the figures scripts of the ISMIP-HOM benchmark. + +\end{acknowledgments} +%}}} +%References {{{1 +\begin{thebibliography}{54} +\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi + +\bibitem[{val(2011)}]{valgrind-website} +Valgrind website, 2011, http://valgrind.org. + +\bibitem[{{\it Amestoy et~al.\/}(2001){\it Amestoy, Duff, Koster, and + L'Excellent\/}}]{MUMPS:1} +Amestoy, P.~R., I.~S. Duff, J.~Koster, and J.-Y. L'Excellent, A fully + asynchronous multifrontal solver using distributed dynamic scheduling, {\it + SIAM J. Matrix Anal. Appl.\/}, {\it 23\/}, 15--41, 2001. + +\bibitem[{{\it Amestoy et~al.\/}(2006){\it Amestoy, Guermouche, L'Excellent, + and Pralet\/}}]{MUMPS:2} +Amestoy, P.~R., A.~Guermouche, J.-Y. L'Excellent, and S.~Pralet, Hybrid + scheduling for the parallel solution of linear systems, {\it Parallel + Computing\/}, {\it 32\/}, 136--156, 2006. + +\bibitem[{{\it Balay et~al.\/}(1997){\it Balay, Gropp, McInnes, and + Smith\/}}]{petsc-efficient} +Balay, S., W.~D. Gropp, L.~C. McInnes, and B.~F. Smith, Efficient management of + parallelism in object oriented numerical software libraries, in {\it Modern + Software Tools in Scientific Computing\/}, edited by E.~Arge, A.~M. Bruaset, + and H.~P. Langtangen, pp. 163--202, Birkh{\"{a}}user Press, 1997. + +\bibitem[{{\it Balay et~al.\/}(2008){\it Balay, Buschelman, Eijkhout, Gropp, + Kaushik, Knepley, McInnes, Smith, and Zhang\/}}]{petsc-user-ref} +Balay, S., K.~Buschelman, V.~Eijkhout, W.~D. Gropp, D.~Kaushik, M.~G. Knepley, + L.~C. McInnes, B.~F. Smith, and H.~Zhang, Petsc users manual, {\it Tech. Rep. + ANL-95/11 - Revision 3.0.0\/}, Argonne National Laboratory, 2008. + +\bibitem[{{\it Balay et~al.\/}(2009){\it Balay, Buschelman, Gropp, Kaushik, + Knepley, McInnes, Smith, and Zhang\/}}]{petsc-website} +Balay, S., K.~Buschelman, W.~D. Gropp, D.~Kaushik, M.~G. Knepley, L.~C. + McInnes, B.~F. Smith, and H.~Zhang, Petsc website, 2009, + http://www.mcs.anl.gov/petsc. + +\bibitem[{{\it Bamber et~al.\/}(2009){\it Bamber, Gomez-Dans, and + Griggs\/}}]{Bamber2009} +Bamber, J.~L., J.~L. Gomez-Dans, and J.~A. Griggs, A new 1 km digital elevation + model of the {A}ntarctic derived from combined satellite radar and laser data + - {P}art 1: {D}ata and methods, {\it The Cryosphere\/}, {\it 3\/}, 101--111, + 2009. + +\bibitem[{{\it Blatter\/}(1995)}]{Blatter1995} +Blatter, H., Velocity and stress-fields in grounded glaciers: A simple + algorithm for including deviatoric stress gradients, {\it J. Glaciol.\/}, + {\it 41\/}, 333--344, 1995. + +\bibitem[{{\it Brezzi et~al.\/}(2004){\it Brezzi, Marini, and + Suli\/}}]{Brezzi2004} +Brezzi, F., L.~Marini, and E.~Suli, Discontinuous {G}alerkin methods for + first-order hyperbolic problems, {\it Math. Models Methods Appl. Sci.\/}, + {\it 14\/}, 1893--1903, 2004. + +\bibitem[{{\it Courant\/}(1943)}]{Courant1943} +Courant, R., Variational methods for the solution of problems of equilibrium + and vibrations, {\it Bull. Amer. Math. Soc.\/}, {\it 49\/}, 1--23, 1943. + +\bibitem[{{\it Durand et~al.\/}(2009{\natexlab{a}}){\it Durand, Gagliardini, + de~Fleurian, Zwinger, and Le~Meur\/}}]{Durand2009a} +Durand, G., O.~Gagliardini, B.~de~Fleurian, T.~Zwinger, and E.~Le~Meur, Marine + ice sheet dynamics: Hysteresis and neutral equilibrium, {\it J. Geophys. + Res.\/}, {\it 114\/}, 1--10, 2009{\natexlab{a}}. + +\bibitem[{{\it Durand et~al.\/}(2009{\natexlab{b}}){\it Durand, Gagliardini, + Zwinger, Le~Meur, and Hindmarsh\/}}]{Durand2009} +Durand, G., O.~Gagliardini, T.~Zwinger, E.~Le~Meur, and R.~Hindmarsh, Full + {S}tokes modeling of marine ice sheets: influence of the grid size, {\it Ann. + Glaciol.\/}, {\it 50(52)\/}, 109--114, 2009{\natexlab{b}}. + +\bibitem[{{\it Frey\/}(2001)}]{Frey2001} +Frey, P.~J., Yams, a fully automatic adaptive isotropic surface remeshing + procedure, {\it Tech. Rep. RT-0252\/}, INRIA, Rocquencourt, 2001. + +\bibitem[{{\it Giovinetto et~al.\/}(1990){\it Giovinetto, Waters, and + Bentley\/}}]{Giovinetto1990} +Giovinetto, M., N.~Waters, and C.~Bentley, Dependence of {A}ntarctic surface + mass balance on temperature, elevation, and distance to open ocean, {\it J. + Geophys. Res.\/}, {\it 95\/}, 3517--3531, 1990. + +\bibitem[{{\it Glen\/}(1955)}]{Glen1955} +Glen, J., The creep of polycrystalline ice, {\it Proc. R. Soc. A\/}, {\it + 228\/}, 519--538, 1955. + +\bibitem[{{\it Gresho and Sani\/}(2000{\natexlab{a}})}]{Gresho2000} +Gresho, P.~M., and R.~L. Sani, {\it Incompressible Flow and the Finite Element + Method, Volume1 : Advection-Diffusion\/}, vol.~1, John Wiley and Sons, + Hoboken, New Jersey, 2000{\natexlab{a}}. + +\bibitem[{{\it Gresho and Sani\/}(2000{\natexlab{b}})}]{Gresho2000a} +Gresho, P.~M., and R.~L. Sani, {\it Incompressible Flow and the Finite Element + Method, Volume2 : Isothermal Laminar Flow\/}, vol.~2, John Wiley and Sons, + Hoboken, New Jersey, 2000{\natexlab{b}}. + +\bibitem[{{\it Gropp et~al.\/}(1996){\it Gropp, Lusk, Doss, and + Skjellum\/}}]{Gropp:1996:HPI} +Gropp, W., E.~Lusk, N.~Doss, and A.~Skjellum, A high-performance, portable + implementation of the {MPI} message passing interface standard, {\it Parallel + Computing\/}, {\it 22\/}, 789--828, 1996. + +\bibitem[{{\it Gropp and Lusk\/}(1996)}]{mpich-user} +Gropp, W.~D., and E.~Lusk, {\it User's Guide for {\tt mpich}, a Portable + Implementation of {MPI}\/}, Mathematics and Computer Science Division, + Argonne National Laboratory, 1996, aNL-96/6. + +\bibitem[{{\it Habashi et~al.\/}(2000){\it Habashi, Dompierre, Bourgault, + Ait-Ali-Yahia, Fortin, and Vallet\/}}]{Habashi2000} +Habashi, W., J.~Dompierre, Y.~Bourgault, D.~Ait-Ali-Yahia, M.~Fortin, and + M.~Vallet, Anisotropic mesh adaptation: towards user-independent, + mesh-independent and solver-independent cfd. part i: general principles, {\it + Int. J. Numer. Meth. Fluids\/}, {\it 32\/}, 725--744, 2000. + +\bibitem[{{\it Hecht\/}(2006)}]{Hecht2006} +Hecht, F., {BAMG}: Bi-dimensional anisotropic mesh generator, {\it Tech. + rep.\/}, FreeFem++, 2006. + +\bibitem[{{\it Hutter\/}(1983)}]{Hutter1983} +Hutter, K., {\it Theoretical glaciology: material science of ice and the + mechanics of glaciers and ice sheets\/}, D. Reidel Publishing Co, Dordrecht, + The Netherlands, 1983. + +\bibitem[{{\it Huybrechts et~al.\/}(1996){\it Huybrechts, Payne, and the + {EISMINT} Intercomparison~Group\/}}]{Huybrechts1996} +Huybrechts, P., A.~Payne, and the {EISMINT} Intercomparison~Group, The + {EISMINT} benchmarks for testing ice-sheet models, {\it Ann. Glaciol.\/}, + {\it 23\/}, 1--14, 1996. + +\bibitem[{{\it Johnson et~al.\/}(1984){\it Johnson, Navert, and + Pitkaranta\/}}]{Johnson1984} +Johnson, C., U.~Navert, and J.~Pitkaranta, Finite-element {M}ethods for linear + hyperbolic problems, {\it Comput. Methods Appl. Mech. Engrg.\/}, {\it 45\/}, + 285--312, 1984. + +\bibitem[{{\it Karypis and Kumar\/}(1998)}]{Karypis1998} +Karypis, G., and V.~Kumar, {\it A {S}oftware {P}ackage for {P}artitioning + {U}nstructured {G}raphs, {P}artitioning {M}eshes, and {C}omputing + {F}ill-{R}educing {O}rderings of {S}parse {M}atrices\/}, University of + {M}innesota, 1998. + +\bibitem[{{\it Kernighan\/}(1988)}]{Kernighan1988} +Kernighan, B.~W., {\it The {C} {P}rogramming {L}anguage\/}, 2nd ed., Prentice + {H}all {P}rofessional {T}echnical {R}eference, 1988. + +\bibitem[{{\it Larour et~al.\/}(2009){\it Larour, Rignot, Seroussi, Morlighem, + and Schiermeier\/}}]{issm-website} +Larour, E., E.~Rignot, H.~Seroussi, M.~Morlighem, and J.~Schiermeier, {ISSM} + website, 2009, http://issm.jpl.nasa.gov. + +\bibitem[{{\it MacAyeal\/}(1989)}]{MacAyeal1989} +MacAyeal, D., Large-scale ice flow over a viscous basal sediment - {T}heory and + application to ice stream-{B}, {A}ntarctica, {\it J. Geophys. Res.\/}, {\it + 94\/}, 4071--4087, 1989. + +\bibitem[{{\it MacAyeal\/}(1993)}]{MacAyeal1993a} +MacAyeal, D., A tutorial on the use of control methods in ice-sheet modeling, + {\it J. Glaciol.\/}, {\it 39\/}, 91--98, 1993. + +\bibitem[{{\it Mathworks\/}(2008)}]{matlab} +Mathworks, T., Matlab (2008a, the mathworks, natick, ma, united states), 2008. + +\bibitem[{{\it Maule et~al.\/}(2005){\it Maule, Purucker, Olsen, and + Mosegaard\/}}]{Maule2005} +Maule, C.~F., M.~E. Purucker, N.~Olsen, and K.~Mosegaard, Heat flux anomalies + in {A}ntarctica revealed by satellite magnetic data, {\it Science\/}, {\it + 309\/}, 464--467, 2005. + +\bibitem[{{\it Morlighem et~al.\/}(2010){\it Morlighem, Rignot, Seroussi, + Larour, Ben~Dhia, and Aubry\/}}]{Morlighem2010} +Morlighem, M., E.~Rignot, H.~Seroussi, E.~Larour, H.~Ben~Dhia, and D.~Aubry, + Spatial patterns of basal drag inferred using control methods from a + full-{S}tokes and simpler models for {P}ine {I}sland {G}lacier, {W}est + {A}ntarctica, {\it Geophys. Res. Lett.\/}, {\it 37, L14502\/}, 1--6, 2010. + +\bibitem[{{\it Nethercote and Seward\/}(2007{\natexlab{a}})}]{Nethercote2007} +Nethercote, N., and J.~Seward, Valgrind: a framework for heavyweight dynamic + binary instrumentation, {\it SIGPLAN Not.\/}, {\it 42\/}, 89--100, + 2007{\natexlab{a}}. + +\bibitem[{{\it Nethercote and Seward\/}(2007{\natexlab{b}})}]{Nethercote2007a} +Nethercote, N., and J.~Seward, Valgrind: a framework for heavyweight dynamic + binary instrumentation, in {\it Proceedings of the 2007 ACM SIGPLAN + conference on Programming language design and implementation\/}, PLDI '07, + pp. 89--100, ACM, New York, NY, USA, 2007{\natexlab{b}}. + +\bibitem[{{\it Nowicki\/}(2007)}]{Nowicki2007} +Nowicki, S. M.~J., Modelling the transition zone of marine ice sheets, Ph.D. + thesis, University College London, 2007. + +\bibitem[{{\it Nowicki and Wingham\/}(2008)}]{Nowicki2008} +Nowicki, S. M.~J., and D.~J. Wingham, Conditions for a steady ice sheet-ice + shelf junction, {\it Earth Planet. Sci. Lett.\/}, {\it 265\/}, 246--255, + 2008. + +\bibitem[{{\it Paterson\/}(1994)}]{Paterson1994} +Paterson, W., {\it The Physics of Glaciers\/}, 3rd ed., Pergamon Press, Oxford, + London, New York, 1994. + +\bibitem[{{\it Pattyn\/}(2003)}]{Pattyn2003} +Pattyn, F., A new three-dimensional higher-order thermomechanical ice sheet + model: Basic sensitivity, ice stream development, and ice flow across + subglacial lakes, {\it J. Geophys. Res.\/}, {\it 108\/}, 1--15, 2003. + +\bibitem[{{\it Pattyn et~al.\/}(2008)}]{Pattyn2008} +Pattyn, F., et~al., Benchmark experiments for higher-order and full-stokes ice + sheet models (ismip–hom), {\it The Cryosphere\/}, {\it 2\/}, 95--108, 2008. + +\bibitem[{{\it Payne et~al.\/}(2000)}]{Payne2000} +Payne, A., et~al., Results from the eismint model intercomparison: the effects + of thermomechanical coupling, {\it J. Glaciol.\/}, {\it 46\/}, 227--238, + 2000. + +\bibitem[{{\it Pilato et~al.\/}(2008){\it Pilato, Collins-Sussman, and + Fitzpatrick\/}}]{Subversion} +Pilato, C.~M., B.~Collins-Sussman, and B.~W. Fitzpatrick, {\it {Version Control + with Subversion}\/}, 2 ed., O'Reilly Media, 2008. + +\bibitem[{{\it Reist\/}(2005)}]{Reist2005} +Reist, A., Mathematical analysis and numerical simulation of the motion of a + glacier, Ph.D. thesis, Ecole Polytechnique F\'ed\'erale de Lausanne, 2005. + +\bibitem[{{\it Rignot\/}(2008)}]{Rignot2008} +Rignot, E., {PALSAR} studies of ice sheet motion in {A}ntarctica, in {\it + {ALOS} {PI} Symposium, Nov. 3-7 2008\/}, 2008. + +\bibitem[{{\it Saad and Schultz\/}(1986)}]{Saad1986} +Saad, Y., and M.~H. Schultz, {GMRES: A generalized minimal residual method for + solving nonsymmetric linear systems}, {\it SIAM J. Sci. Stat. Comput.\/}, + {\it 7\/}, 856--869, 1986. + +\bibitem[{{\it Schoof\/}(2007{\natexlab{a}})}]{Schoof2007a} +Schoof, C., Marine ice-sheet dynamics. part 1. the case of rapid sliding, {\it + J. Fluid Mech.\/}, {\it 573\/}, 27--55, 2007{\natexlab{a}}. + +\bibitem[{{\it Schoof\/}(2007{\natexlab{b}})}]{Schoof2007b} +Schoof, C., Ice sheet grounding line dynamics: Steady states, stability, and + hysteresis, {\it J. Geophys. Res.\/}, {\it 112\/}, 1--19, 2007{\natexlab{b}}. + +\bibitem[{{\it Shewchuk\/}(1996)}]{Shewchuk1996} +Shewchuk, J.~R., Triangle: {E}ngineering a {2D} {Q}uality {M}esh {G}enerator + and {D}elaunay {T}riangulator, in {\it Applied Computational Geometry: + Towards Geometric Engineering\/}, edited by M.~C. Lin and D.~Manocha, vol. + 1148 of {\it Lecture Notes in Computer Science\/}, pp. 203--222, + Springer-Verlag, 1996, from the First ACM Workshop on Applied Computational + Geometry. + +\bibitem[{{\it Shewchuk\/}(2001)}]{Shewchuk2002} +Shewchuk, J.~R., Delaunay refinement algorithms for triangular mesh generation, + {\it Computational Geometry: Theory and Applications\/}, {\it 22\/}, 1--3, + 2001. + +\bibitem[{{\it Stroustrup\/}(1997)}]{Stroustrup1997} +Stroustrup, B., {\it The {C}++ Programming Language\/}, 3rd ed., Addison-Wesley + Longman Publishing Co., Inc., Boston, MA, USA, 1997. + +\bibitem[{{\it van~den Broeke\/}(2008)}]{Broeke2008} +van~den Broeke, M., Depth and density of the {A}ntarctic firn layer, {\it + Arctic Ant. Alp. Res.\/}, {\it 40\/}, 432--438, 2008. + +\bibitem[{{\it van~der Veen and Whillans\/}(1989)}]{Veen1989} +van~der Veen, C.~J., and I.~M. Whillans, Force budget: I. theory and numerical + methods, {\it J. Glaciol.\/}, {\it 35\/}, 53--60, 1989. + +\bibitem[{{\it Vaughan et~al.\/}(2006)}]{Vaughan2006} +Vaughan, D., et~al., New boundary conditions for the {W}est {A}ntarctic {I}ce + {S}heet: Subglacial topography beneath {P}ine {I}sland {G}lacier, {\it + Geophys. Res. Lett.\/}, {\it 33\/}, 1--4, 2006. + +\bibitem[{{\it Vieli and Payne\/}(2003)}]{Vieli2003} +Vieli, A., and A.~Payne, Application of control methods for modelling the flow + of {P}ine {I}sland {G}lacier, {W}est {A}ntarctica, {\it Ann. Glaciol.\/}, + {\it 36\/}, 197--204, 2003. + +\bibitem[{{\it Vogel\/}(2002)}]{Vogel2002} +Vogel, C.~R., {\it Computational Methods for Inverse Problems\/}, Society for + Industrial and Applied Mathematics, Philadelphia, PA, USA, 2002. + +\end{thebibliography} +%}}} +\end{article} +%Figures{{{1 +\begin{figure} %PigBamg{{{2 + %\includegraphics[width=25pc]{figures/PigBamg.eps} + \caption{Anisotropic static adaptive mesh refinement, applied to the Pine Island Glacier, Antarctica. Frame a: InSAR surface velocity + from \cite{Rignot2008}, overlayed over adapted mesh (in white). Frame b: same InSAR surface velocity overlayed over regular mesh. + Both meshes comprise 2700 elements.} + \label{fig:pigbamg} +\end{figure} +%}}} +\end{document}