Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM)

[ 1 ] Ice flow models used to project the mass balance of ice sheets in Greenland and Antarctica usually rely on the Shallow Ice Approximation (SIA) and the Shallow-Shelf Approximation (SSA), sometimes combined into so-called “ hybrid ” models. Such models, while computationally efficient, are based on a simplified set of physical assumptions about the mechanical regime of the ice flow, which does not uniformly apply everywhere on the ice sheet/ice shelf system, especially near grounding lines, where rapid changes are taking place at present. Here, we present a new thermomechanical finite element model of ice flow named ISSM (Ice Sheet System Model) that includes higher-order stresses, high spatial resolution capability and data assimilation techniques to better capture ice dynamics and produce realistic simulations of ice sheet flow at the continental scale. ISSM includes several approximations of the momentum balance equations, ranging from the two-dimensional SSA to the three-dimensional full-Stokes formulation. It also relies on a massively parallelized architecture and state-of-the-art scalable tools. ISSM employs data assimilation techniques, at all levels of approximation of the momentum balance equations, to infer basal drag at the ice-bed interface from satellite radar interferometry-derived observations of ice motion. Following a validation of ISSM with standard benchmarks, we present a demonstration of its capability in the case of the Greenland Ice Sheet. We show ISSM is able to simulate the ice flow of an entire ice sheet realistically at a high spatial resolution, with higher-order physics, thereby providing a pathway for improving projections of ice sheet evolution in a warming climate.


Introduction
[2] Detailed and realistic modeling of the evolution of the Antarctic and Greenland Ice Sheets is needed to improve projections of sea level rise in a warming climate [Intergovernmental Panel on Climate Change (IPCC), 2007]. Such a modeling effort requires an accurate representation of the physics of ice motion, well-constrained boundary conditions at varying geometrical scales and computationally scalable software packages. In the past, large-scale ice sheet models were based on simplified models for ice flow, e.g., the Shallow Ice Approximation (SIA), which assumes that the basal shear stress of grounded ice is completely balanced by the gravitational driving stress [Hutter, 1982]. The SIA is computationally efficient and practical, especially for long-term reconstructions, but it is not relevant for short-term projections of the rapid evolution of ice streams, glaciers and ice shelves in response to climate change.
[3] On floating ice shelves, the SIA breaks down as vertical shear becomes negligible compared to lateral shear. The Shallow Shelf or Shelfy Stream Approximation (SSA) provides a viable, alternative solution that neglects vertical shear stresses and assumes that horizontal velocity is depth independent [MacAyeal, 1989;Morland and Zainuddin, 1987]. The SSA was initially designed for floating ice shelves and was extended to the case of fast-flowing ice streams, where most of the displacement is due to sliding. The SSA is a twodimensional model that solves for the depth-averaged velocity and therefore cannot be applied in areas where vertical variations in speed are not negligible like in the vicinity of grounding lines (where ice detaches from the bed and starts floating), ice stream margins, or the complex flow near an ice divide.
[4] Higher-order models, e.g., the three-dimensional (3D) Blatter-Pattyn approximation (BP) [Blatter, 1995;Pattyn, 2003], are more desirable in these situations. These models capture longitudinal stresses that are significant in fast-flowing ice streams, as well as components of the vertical shear stress that dominate in areas where flow is slow, and sliding at the base is small. These models are especially relevant for areas where both regimes of ice flow are of equal importance, as shown by Pattyn [1996] for the Shirase Glacier, East Antarctica. These models are usable over most of the ice sheet and ice shelf areas, but their assumptions are not strictly valid near grounding lines where a full-Stokes model (FS) is significantly more accurate [Lestringant, 1994;Nowicki and Wingham, 2008;Durand et al., 2009aDurand et al., , 2009bGagliardini et al., 2010;Morlighem et al., 2010]. FS is the most physically relevant model, but it is computationally expensive, therefore challenging to use at the continental scale at a high spatial resolution. Furthermore, FS may not be needed on large parts of the ice sheet where the resulting improvement in model accuracy will not be significant.
[5] Here we present a new finite element, thermomechanical numerical model of ice flow named ISSM (Ice Sheet System Model) that includes higher-order stress components, high spatial resolution capability and relies on a massively parallelized architecture. In the first part, we describe the equations adopted in ISSM to model ice flow, including its thermal regime, stress regime, boundary conditions and ice rheology. In the second part, we describe how we solve these equations using the Finite Element Method (FEM), static adaptive mesh refinement, scalable parallel iterative and direct solvers, and a specifically tailored software architecture. We then discuss our implementation of inverse methods, and what was done in ISSM to make it applicable at the continental scale. In the third part, we discuss the evaluation of ISSM against ice flow benchmarks. In the fourth part, we present an application of ISSM to the modeling of the entire Greenland Ice Sheet (GIS) to illustrate the large scale applicability of ISSM. In the fifth part, we discuss the results and various steps for improving the model. We conclude on the utilization of ISSM to improve projections of ice sheet contribution to sea level rise.

Model
[6] ISSM is a coupled thermomechanical ice flow model. It relies on the classical conservation laws for momentum balance, mass balance and energy balance, combined with constitutive material laws and boundary conditions. These laws have been described in numerous prior publications but are repeated here for the sake of completeness and to explain precisely how ISSM is implemented. In the following subsections, we describe the mechanical model, the thermal model, boundary conditions, calving front and grounding line dynamics. All notations and values for physical constants are presented in Table 1 for ease of use. Figure 1 can be used as reference for the geometrical information and notations used in the model.

Mechanical Model
[7] ISSM comprises four different ice flow models, based on a gradual simplification of the full-Stokes (FS) momentum balance equations. The first model is FS: Tr where r Á s is the divergence vector of the stress tensor, s, r the ice density, g the acceleration due to gravity, _ ɛ the strain rate tensor and Tr the trace operator.
[8] Equation (1) expresses the balance of stresses, and equation (2) expresses the incompressibility of flow. In equation (1), acceleration is neglected, following a scale analysis by Reist [2005] that shows ice flow acceleration is negligible even in the most extreme surges or streamflow that may occur in a glacier.
[9] The material constitutive law describes the deformation of ice under stress. For incompressible viscous fluids, the constitutive law has the form: where s′ = s + pI is the deviatoric stress tensor, I the identity matrix, p the ice pressure and m the ice effective viscosity. The ice effective viscosity is assumed to be nonlinear and follow a Norton-Hoff law [Glen, 1955]: where B is the ice hardness, n Glen's law exponent and _ ɛ e the effective strain rate (defined as the second invariant of the strain rate tensor). B is temperature dependent and follows an Arrhenius law: where A 0 is the flow factor, Q the activation energy for ice creep, R the universal gas constant and T* = T À bp the absolute temperature corrected for the dependence of melting point on pressure, b being the rate of change of melting point with pressure. The values of these constants follow Payne et al. [2000] (see Table 1). Alternatively, ice hardness, B, could be chosen to follow the temperature-dependent relationship discussed by Cuffey and Paterson [2010].
[10] Using equation (3) and Glen's flow law (equation (4)), Equation (1) can be rewritten in terms of velocity components and pressure: where (u, v, w) are the x, y and z components of the velocity vector v, in the (x, y, z) Cartesian coordinate system, with z in the vertical direction. Equations (6), (7), (8), and (9) are at the core of FS. This model includes four unknowns (u, v, w, p) and it is therefore computationally demanding to solve for.
[11] The second model is a simplified three-dimensional (3D) model from Blatter [1995] and Pattyn [2003] (BP) derived from FS by making two assumptions: (1) the horizontal gradients of vertical velocities are negligible compared to the vertical gradients of horizontal velocities: and (2) the bridging effects [van der Veen and Whillans, 1989] are negligible, which reduces the third equation of the momentum balance (equation (8)) to: For a more exhaustive presentation of this derivation, we refer the reader to Schoof and Hindmarsh [2010]. The first two FS equations (equations (6) and (7)) are then reduced to: [12] These two equations are at the core of the BP model, which has only two unknowns (u and v) instead of four. Vertical velocity, w, is recovered through incompressibility, and is hence decoupled from the initial system of equations.
[13] The third model, the Shallow-Shelf or Shelfy-Stream Approximation (SSA) [MacAyeal, 1989] results from further assuming that vertical shear is negligible: This assumption makes equations (12) and (13) independent of z. Depth averaging yields: where m is the depth-averaged ice viscosity, H the local ice thickness and s the upper surface elevation. Vertical velocity, w, is again recovered through incompressibility (equation (2)). SSA is a two-dimensional (2D) model with two unknowns (u and v). It is therefore considerably faster to solve than BP because it is a 2D model whereas BP is a 3D model.
[14] The last model implemented in ISSM is the Shallow Ice Approximation (SIA) [Hutter , 1983]. In this model, only the deviatoric stress components s′ xz and s′ yz are included. The horizontal gradients of vertical velocity are also neglected compared to the vertical gradients of horizontal velocities, so the FS equations are reduced to: [15] This simple model is the basis of many ice sheet models used in long term reconstruction of ice sheet dynamics [Payne and Baldwin, 2000;Ritz et al., 1997;Huybrechts et al., 2003]. It is computationally efficient, as each equation (for components u and v) can be solved independently, involving only one degree of freedom.

Mechanical Boundary Conditions
[16] For all four models (FS, BP, SSA and SIA), boundary conditions are required to ensure a solution can be computed. At the surface of the ice sheet, we assume a stress-free surface: where n is the unit outward-pointing normal vector. Friction is applied at the ice-bedrock interface. The basal drag follows an empirical relationship presented by Paterson [1994] and first calibrated by Budd et al. [1979], written in a viscous-type law of friction: where v b is the basal velocity vector tangential to the glacier base plane, t b the tangential component of the external force s Á n, and a 2 a positive constant (i.e., stress opposes ice motion). Paterson [1994], usually includes the effective pressure in the basal friction law, to account for the presence of water lubricating the ice-bedrock interface. However, ISSM is not currently capable of modeling the distribution of water under an ice sheet.
[17] At the ice-seawater interface, water pressure is imposed according to: where r w is the water density and z the vertical coordinate equal to zero at sea level (positive upwards). For the BP model, the ice pressure is taken equal to the lithostatic pressure so that this boundary condition may be expressed in terms of the deviatoric stress: [18] We assume this approximation to be valid as long as the term ∂ ∂z (2m ∂w ∂z ) can be neglected in equation (11). A more complete implementation of this boundary condition assuming non lithostatic pressure is currently being implemented for such cases.
[19] For SSA, this boundary condition is depth-integrated: where b is the elevation of the ice lower surface.
[20] On all other boundaries where stresses are not specified, velocity constraints are applied. This is necessary in order for the formulations to be complete and for non-singular solutions to exist. Usually, observed surface velocities are used to constrain such boundaries.

Thermal Model
[21] The ice hardness, B, in equation (4) is temperature dependent. The thermal equation is derived from the energy balance equation and includes conduction-advection in three directions. We do not use s-coordinates [Pattyn, 2003] but rely on the Arbitrary Lagrangian-Eulerian (ALE) method [Hughes et al., 1981;Donea et al., 2004]. The mesh velocity, w, is therefore included in the equation to account for mesh displacements: where T is the ice temperature, k th the ice thermal conductivity, c the ice heat capacity, assumed to be both constant, D the Laplace operator, and F the heat production term (deformational heating): where _ ɛ e is the effective strain rate.
[22] Ice temperature needs to be kept below the pressure melting point. Our strategy for implementing this will be presented in the Numerical Implementation Section.

Thermal Boundary Conditions
[23] The thermal model presented above is constrained using the following boundary conditions. At the surface, air temperature is prescribed as a Dirichlet boundary condition equal to the mean annual air temperature, T s : On grounded ice, we impose the following relationship between the geothermal heat flux, G, and the frictional heat flux: Under floating ice, we use the simple parameterization from Holland and Jenkins [1999]: with c pM the mixed layer specific heat capacity, g the thermal exchange velocity, and T f = 273.15 À bp the temperature of freezing of seawater. This parameterization does not explicitly account for melting at the base, which has been shown to be a very important process. Further work is indeed required to couple ISSM and ocean circulation models in order to accurately capture such melting rates. For all remaining boundaries, a zero flux boundary condition is applied.

Mass Transport Model
[24] The mass conservation equation relates ice thickness to ice flux divergence, surface mass balance and basal mass balance. This equation is used to calculate the rate of ice thickness change: where H is ice thickness, v = ( u , v ) the depth-averaged velocity, _ M s the surface mass balance (m a À1 in ice equivalent, positive for accumulation, negative for ablation) and _ M b the basal melting rate (m a À1 in ice equivalent, positive when melting, negative when freezing). This equation is used to change the glacier thickness at each time step. In order to update surface and bed position using the new thickness, we implement the following strategy: for grounded ice, the bedrock is assumed to stay constant in time, and the entire thickness is added to the bedrock position to infer the new surface; for floating ice, the thickness increment is distributed between the bed and surface using a hydrostatic assumption. This strategy is not ideal for floating ice in the immediate vicinity of the grounding line, where it has been shown that hydrostatic equilibrium is not always fulfilled [Lestringant, 1994]. This is specifically an issue for the FS model which does not assume hydrostatic equilibrium, and which can be very sensitive to variations in the ice geometry near the grounding line. Further work is therefore required to implement a formulation of the transport equations using the local surface velocity to transport the upper free surface, similar to what has been implemented by Nowicki and Wingham [2008] and Durand et al. [2009aDurand et al. [ , 2009b.

Mass Transport Boundary Conditions
[25] The following boundary conditions are used for the mass transport model in case of regional scale models, where boundaries define a closed domain. Ice thickness is then constrained on the inflow boundary, or at the ice divide (where mass advection is negligible): where H obs is the measured ice thickness. At outflow boundaries, we apply a free-flux boundary condition so that the mass flux is left unconstrained. For large scale ice sheet models, where the only boundary is the ice front, no constraint is needed at the ice divide, only free-flux boundary conditions for the entire ice front.

Ice Front Dynamics
[26] The present-day ice front position is imposed as a mesh boundary, which remains fixed in time. We apply a free-flux boundary condition at the ice front where the melt rate and calving rate are assumed to exactly compensate for the ice speed at the front, i.e., we assume a stable, fixed ice front position. In future developments, we will allow this boundary to migrate with time, with a migration controlled by a calving law. At this time, ISSM includes the simplest form of calving front dynamics.
[27] Indeed, this criterion is also used for the entire ice margin of the ice sheet. In case the margin retreats, ISSM assumes a minimum ice thickness of 1 m, which allows for the retreating ice to become dynamically decoupled from the rest of the ice sheet, without introducing instabilities in the ice flow dynamics, and without the need for actively imposing velocity constraints when the ice is fully retreated. This scenario will also be covered more efficiently when dynamic boundary migration is implemented.

Grounding Line Dynamics
[28] Grounding line migration is a key control on ice flow dynamics [Schoof, 2007a[Schoof, , 2007bNowicki, 2007;Durand et al., 2009bDurand et al., , 2009a, which is important to capture in order for transient ice flow models to be realistic. In ISSM, we use a simple 3D grounding line migration criterion based on the hydrostatic equilibrium (see Figure 1). At each time step of the transient ice flow solution, we check the following for every vertex of the mesh: 1. The criteria b ≤ b a where b a is the depth of the glacier bed or seafloor. For most ice sheet/ice shelf configurations, b is negative. If this condition is verified for a floating vertex (i.e., on an ice shelf), we ground the vertex and force b = b a .
2. The criteria b > b HE where b HE is the depth of the bottom of the ice in hydrostatic equilibrium: b HE = ÀHr/r w . If this condition is verified for a grounded vertex (i.e., on the ice sheet), we unground the vertex and force b = b HE .
[29] These two criteria 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. This means that we do not allow vertices that are not connected to the grounding line to unground. Of course, the presence of hydrological networks connected to the open ocean near the grounding line is contrary to this assumption. However ISSM is not currently capable of modeling such water distribution under an ice sheet, especially when interactions with the open ocean are considered. Further work is therefore required to include better hydrological/ocean circulation interactions into our migration criteria.
[30] Furthermore, our criteria for grounding line migration need to be further refined, using a full-Stokes computation of the stress-balance of the ice sheet/ice shelf system, and migration based on contact and pressure conditions at the grounding line [Nowicki and Wingham, 2008;Durand et al., 2009aDurand et al., , 2009b. Work is currently in progress toward implementing such criteria.

Numerical Implementation
[31] The previous section dealt with the classic equations for modeling ice flow. Solving these equations is a challenge when the goal is to capture ice flow at the continental scale, e.g., in Greenland or Antarctica. In this section, we describe how we approach this challenge in ISSM using the FEM, static adaptive mesh refinement, parallel iterative solvers, and inverse methods for the different approximations of ISSM.

FEM Discretization and Treatment of Nonlinearities
[32] ISSM employs the Continuous Galerkin Finite Element Method for discretization of the ice flow equations. This method allows for the use of unstructured meshes, which are most efficient to capture varying spatial resolutions across an ice sheet. For the set of diagnostic equations, triangular Lagrange P1 elements (piecewise linear) are used in 2D and prismatic P1 elements in 3D, except for FS. These elements are simple to implement and provide stable discretizations. For FS, they become unstable, and condensed MINI elements [Gresho and Sani , 2000b] have to be used to fulfill the compatibility Ladyzhenskaya-Babuška-Brezzi (LBB) condition. In order to treat the material nonlinearity introduced by the viscosity dependence on the effective strain rate (equation (4)), we rely on a Picard iteration scheme [Hindmarsh and Payne, 1996;Reist, 2005;De Smedt et al., 2010]. Several stopping criteria are implemented: 1. The criteria, res < 1 , where res is defined as kKU old À Fk/ kFk where K is the current iteration stiffness matrix, F the current load vector, U old the solution vector at the previous iteration, kÁk the Euclidean norm, and 1 a threshold (usually taken a 10 À4 , see Table 1). This criterion ensures the convergence of the mechanical stress equilibrium.
2. The criteria, rel < 2 , where rel is defined as k(U old À U)/ U old k ∞ , where U is the solution vector of the current iteration, kÁk ∞ the infinity norm, and 2 a threshold (usually taken as 10 À2 , see Table 1). This criterion ensures the convergence of the solution vector in relative terms.
3. The criteria, abs < 3 , where abs is defined as kU old À Uk ∞ and 3 a threshold (usually taken as 10 m/a, see Table 1). This criterion ensures convergence of the solution vector in absolute terms.
[33] All three criteria can be combined, although ISSM ensures that at least the first criterion on the convergence of the mechanical stress equilibrium is always met.
[34] For the thermal model, P1 elements are used, combined with the Streamline Upwind Petrov-Galerkin (SUPG) [Gresho and Sani, 2000a] formulation to prevent potential numerical oscillations due to dominant advection terms. Temperature is kept below the pressure melting point using an iterative penalty-based scheme, similar to those used for contact problems [Courant, 1943]. This method is similar in nature to Zwinger et al. [2007], with the difference that penalties are used instead of Dirichlet boundary conditions. Once the number of applied penalties does not vary anymore for additional iterations, we consider that convergence has been reached.
[35] For the mass transport model, we use an implicit finite difference approximation in time, which is more stable than explicit schemes. For this type of hyperbolic problems, instabilities still develop when relying on the standard Galerkin FEM, and for this case, we use either streamline upwinding or discontinuous Galerkin finite elements [Johnson et al., 1984;Brezzi et al., 2004].
[36] Using the FEM for the thermal, mechanical and mass transport equations allows us to significantly improve numerical performance, when combined with a judicious selection of the underlying mesh. Static adaptive mesh refinement is the method we have implemented to optimize the number of degrees of freedom solved in each set of equations.

Mesh Refinement
[37] At the continental scale, solving for ice flow velocity involves a large number of degrees of freedom (hereafter referred to as dof), which can prove computationally challenging, even using parallel technologies and computer clusters. In order to minimize the number of dofs, ISSM relies on static anisotropic adaptive mesh refinement. This technique distorts the mesh and reduces discretization errors, while minimizing the number of dofs required. The mesh refinement is done once during model setup (static), as opposed to during the solution run (dynamic). This approach is physically realistic provided ice flow does not evolve drastically, which restricts this approach to short term transient runs. For longer transients, further work is required to ensure dynamic adaptation of the mesh to capture the evolving ice flow.
[38] 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 [Habashi et al., 2000] provided that the solution is regular enough. For each element of a mesh, E, the difference between a P1 interpolated field, u h , and the exact field, u, is bounded as follows: where C is a constant that depends only on the space dimension, h is the characteristic length of the element, H u (x) is the Hessian matrix of u (x), and |H u (x)| its spectral norm.
[39] According to this equation, if we compute the Hessian of the observed surface velocities, we can optimize the mesh distribution accordingly. Here, we implement an edge-based anisotropic mesh optimization methodology inspired by the YAMS and BAMG libraries [Frey, 2001;Hecht, 2006] to equi-distribute the interpolation error in each direction of each element to control the global interpolation error. An example of static adaptive anisotropic mesh refinement is shown in Figure 2 for Jakobshavn Isbrae, West Greenland.
[40] Figure 2a shows a uniform 1,500 element mesh comprised of triangles of approximately equal area. In comparison, Figure 2b shows a 1,500 element mesh adapted using InSAR surface velocities. Shear margins are well captured by the algorithm, while the interior of the ice sheet, where ice flow deformation is weak, is captured by a coarser mesh.
Comparison between the two meshes shows that for an equal number of elements, anisotropic mesh refinement captures ice flow far more efficiently, resulting in tremendous computational gains.

Solvers and Performance
[41] ISSM relies on the suite of solvers provided by the Portable, Extensible Toolkit for Scientific Computation package (PETSc) [Balay et al., 1997[Balay et al., , 2008; see also Portable, extensible toolkit for scientific computation, 2011, http://www.mcs.anl.gov/petsc] to compute solution vectors for the diagnostic, mass transport and thermal models. These solvers are used on a per-case basis, according to how well conditioned the equations are.
[42] For the SSA, we rely on the iterative GMRES solver [Saad and Schultz, 1986] combined with a Jacobi preconditioner. The SSA is sufficiently well conditioned for this solver to converge consistently. For the BP, the conditioning is not as good and we have to switch to a preconditioning based on the Multigrid Algebraic method [Briggs et al., 2000]. This solver converges consistently for this type of model. For FS, the system of equations is a saddle point problem, for which special preconditioning is needed, based on block factorization of the stiffness matrix, and preconditioning of the Schur complement [Benzi et al., 2005]. Work is currently in progress to implement this solver. As an alternative, ISSM uses the Multifrontal Massively Parallel Sparse direct solver (MUMPS) [Amestoy et al., 2001[Amestoy et al., , 2006. Although poorly scalable, this solver does not suffer from convergence issues, as it relies on a direct solving method.
[43] The performance of ISSM is entirely dependent on the solver phase, where 95% of the computational time is spent, irrespective of the ice flow approximation used. Scalability is improved for the SSA and BP, for which iterative solvers are used. Further work is required to implement a scalable FS solver, similar to what was implemented by Leng et al. [2010], but the performance of the MUMPS solver still ensures that reasonable computational times are reached, as shown later on.

Inverse Methods
[44] Ice flow models are difficult to constrain when parameters such as ice hardness or basal drag coefficient are unknown. For example, the basal drag coefficient, a in equation (19), cannot be measured directly, and it is one of the most critical parameters controlling ice flow. This is why we need inverse methods.
[46] We use a partial differential equation constrained optimization algorithm similar to MacAyeal [1992], 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 calculates the gradient of the cost function with respect to the unknown parameter. This cost function is generally taken as: where G s is defined as the ice upper surface domain.
[47] This cost function works better in areas of highvelocity than in slow moving regions because the adjoint state (Lagrange multipliers vector) is larger where the velocity misfit is high, which occurs in regions of high speed. To minimize this effect, we use a cost function that measures the logarithm of the misfit instead [Morlighem et al., 2010]: where ɛ is a minimum velocity used to avoid zero velocity, and log is the natural logarithm. This cost function enables a robust estimation of the basal drag coefficient or depthaveraged ice hardness after only a few iterations over the entire model domain.
[48] To derive the adjoint, ice viscosity is assumed to be independent of the velocity. This is convenient because for a linear material law, the system of ice flow equations becomes self-adjoint. This assumption is not rigorous because viscosity depends on the strain rate, but this simplification allows an easy calculation of the adjoint state for all three ice flow models (SSA, BP and FS), is widely used in the literature [MacAyeal, 1992[MacAyeal, , 1993 and has a limited impact on the control method convergence [Goldberg and Sergienko, 2011].
[49] This data assimilation technique has been successfully extended to BP and FS [Morlighem et al., 2010]. The major difference with SSA is that only the surface horizontal velocities are taken into account in the cost function evaluation, while the gradient with respect to the basal drag coefficient, a, or the depth-averaged ice hardness, B, are computed at the base only. In 3D models, we have yet to implement ice hardness inversions across the entire thickness of the ice sheet. We are so far limited to depth-averaged hardness.
[50] A Tikhonov regularization term which penalizes the oscillations of the unknown parameter, resulting from the noise in the observations, can be added to the misfit to stabilize the inversion [Vogel, 2002]. This term is defined as: where G b is the bed domain, a is the inverted parameter, and l the Tikhonov parameter. In order to choose l, an L-curve analysis could be employed, as done by Jay-Allemand et al.
[2011], however, for the GIS, the cost of such an analysis would be prohibitive. The choice of l also depends heavily on the type of basin, cost-function and model. For our application to the GIS, our strategy was to increase l to improve the parameter regularity, while retaining significant similarity with the l = 0 misfit. This strategy led to a value of l for the GIS inversion of the basal drag coefficient of 8.1 Â 10 À15 for both logarithmic (m 3 /s/Pa) and absolute (m 5 /s 3 /Pa) cost functions. The approach is qualitatively similar to the L-curve analysis, while avoiding the need for repeated inversions of the entire ice sheet, which is prohibitive. In terms of convergence, given the fact that we rely on a succession of logarithmic and absolute cost functions in each iteration step of the inversion, an objective convergence criterion is difficult to assess. Our strategy to decide on the convergence of the inversion is essentially qualitative, and is based on an estimate of how flat the misfit becomes after a certain amount of iterations, for both logarithmic and absolute cases.
[51] For the remainder of this study, we limit ourselves to the case of a basal drag coefficient inversion, where ice hardness is computed using the thermal model and equation (5). For each step of the optimization, one theoretically needs to re-calculate a thermomechanical equilibrium solution and update ice hardness, B, accordingly on grounded ice, to ensure consistency between ice flow velocity and ice viscosity. This can be done at the regional scale [Morlighem et al., 2010], on models for which computation time is small enough to allow for simultaneous inversion of the basal drag coefficient and update of the thermal regime. Practically, for large scale models such as Antarctica or Greenland, computation times are such that the inversion cannot include updates to the thermal regime. This implies that the temperatures for such large models is computed initially, before the inversion is started, and then kept constant throughout the inversion. This allows computation times for each inversion step to remain manageable.

Software Architecture and Management
[52] ISSM relies on a suite of languages and libraries to achieve high performance computing and scalability. The primary software language is C/C++ [Stroustrup, 1997], which allows for flexible and object oriented development. This kernel is then integrated within the MATLAB environment (MathWorks, software, 2008) to facilitate ease of use, and linked to the Message Passing Interface MPI  and the PETSc libraries to achieve parallelization. Tight integration of these capabilities ensures high performance computing and a significant amount of scalability. For more details on this software integration, and additional packages relied upon to ensure among others accurate verification and validation as well as easy installation, we refer the reader to Appendix A.

Verification of ISSM With Benchmarks
[53] Here we describe the evaluation of ISSM using the ISMIP-HOM (Ice Sheet Model Intercomparison Project for Higher-Order Models) [Pattyn et al., 2008] benchmark. This benchmark targets validation for three dimensional ice flow models. Our results show an excellent agreement with the participating models, for all A to F tests. Here, we present some of these results, specifically for tests A, C and F. We do not show results for tests B, D and E as these tests are specific to flowband models. In experiments A, C and F, we rely on triangular elements and a regular mesh with 50 intervals in both horizontal directions and 20 extruded layers vertically.
[54] Figure 3 shows test A results of the ISMIP-HOM benchmark using the BP and FS. This test is a diagnostic experiment that involves an ice slab flowing over a sinusoidal bumpy-bed with zero velocity at the base. We use the penalty method to impose periodic boundary conditions on the sides. The test is performed for six domain lengths, L, ranging from 5 km to 160 km. Our modeled surface velocity agrees well with the benchmark velocity computed by other software in the benchmark (Figure 3). Results from both higher-order and full-Stokes models are within the interval described by the benchmark models [Pattyn et al., 2008].
[55] Experiment C of the ISMIP-HOM benchmark is similar to A, except for a flat bedrock and a bed that is not frozen. The basal friction coefficient is prescribed by a sinusoidal periodic law, so periodic boundary conditions are also applied here. Our results using both FS and BP agree quite well with the benchmark results (Figure 4), and are also within the intervals described by other models. We note that the results for both BP and FS are very similar at all domain lengths.
[56] The last ISMIP-HOM benchmark presented here is experiment F, which is a transient experiment for a slab of ice that flows over a sloping bed. The initial bedrock and surface are parallel and ice thickness is prescribed at 1,000 m. A Gaussian bump is introduced at the center of the basal topography. As in the previous two experiments, penalties are applied to simulate the periodic boundary conditions. They are used for both velocity and thickness. The free surface and velocity evolve until a steady state solution is reached. This experiment is run with basal sliding (slip bed) and without basal sliding (no-slip bed). Results are presented in Figure 5 for BP and FS. The steady state surface found by both models and for both basal conditions agrees very well with the results of the ISMIP-HOM benchmarks. The steady state surface velocity found for the non-slip bed is slightly higher than the results of other models. Steady state surface velocity solved by ISSM is higher for the BP solution and lower for the FS solution by about 1 m a À1 compared to the case of the slip bed experiment. Only a handful of ice sheet models have successfully performed these tests. In addition, this test exhibits challenges in terms of mass conservation, especially for finite element based treatments of mass transport [Gagliardini and Zwinger, 2008]. In order to check whether our implementation conserves mass throughout the transient run, we computed the total mass of the slab of ice at the beginning of the run and at the end. Both masses were identical within 8 digits. This demonstrates that our numerical implementation efficiently conserves mass within single float precision.

Example Modeling of the Greenland Ice Sheet
[57] Running a transient ice flow model at the continental scale and at high spatial resolution is a significant challenge given the sheer size of the discretized system of equations involved. Here, we discuss the application of ISSM to infer the basal drag coefficient of the Greenland Ice Sheet using three different formulations: (1) SSA, (2) BP, and (3) FS. We also discuss the application of ISSM to model transient ice flow using the higher-order BP, over long time periods (500 years) using inversion methods to initialize unknown boundary conditions. A detailed scientific interpretation of the results is out of the scope of this paper, and will be addressed in a companion paper. Our intent here is to demonstrate that ISSM can be used in diagnostic and transient mode, and efficiently rely on data assimilation methods to replicate the ice flow velocity of an entire ice sheet with high accuracy.

Input Data and Model Setup 5.1.1. Input Data
[58] The following data sets are used to constrain our inversion. InSAR surface velocities are from [Rignot et al.,  [Allen, 2009]. Surface elevation is deduced from bed elevation and ice thickness data from CReSIS. The spatial extent of the ice sheet is from the fuzzy ground cover mask of Box and Decker [2010]. We treat the entire ice sheet as grounded, hence floating ice shelves in the north are ignored to simplify the problem. Boundary conditions for the thermal model are surface temperature from Ettema et al. [2009] and geothermal heat flux from Shapiro and Ritzwoller [2004].
[59] For the transient ice flow simulations, we rely on the model setup presented on the SeaRISE website (http://websrv. cs.umt.edu/isis/index.php/SeaRISE_Assessment), which is fairly similar to our inversion setup. The datasets include mean annual surface temperature from Fausto et al. [2009], precipitation rates from Burgess et al. [2010], basal heat flux from Shapiro and Ritzwoller [2004], and bedrock topography, ice thickness and surface elevation from Bamber et al. [2001]. Observed surface velocities are also provided by Joughin et al. [2010] and filled with balance velocities [see, e.g., Bamber et al., 2000] in the areas where observations are not available. Discontinuities introduced by this patching remain small enough at the level of mesh resolution implemented for the inversions not to focus on the velocity jumps at the boundaries. At high resolution though, this could be an issue, which should be investigated further.

Model Initialization
[60] As explained earlier, running large-scale, fullycoupled, thermal-mechanical inversions is prohibitive given our current computational capabilities. Here, we infer ice viscosity at the beginning of the inversion and keep it constant throughout the initialization. To calculate ice viscosity, we need to calculate the strain heating and basal shear heating components from the thermal model based on a flow simulation that best matches the observed surface velocities. To solve this recursive problem, we model an initial flow using an ad-hoc basal drag coefficient (a in equation (19)) of 1.5 Â 10 5 (Pa s/m) 1/2 for surface velocity greater than 60 m a À1 and 9 Â 10 5 (Pa s/m) 1/2 for surface velocity less than 60 m a À1 . In addition, we use a Dirichlet boundary condition at the surface to impose observed surface velocities. The resulting flow field is used to calculate a thermal regime and ice viscosity that best match the observations. By keeping ice viscosity constant throughout the computation, we improve computational efficiency to the point of making the inversion possible at the continental scale. However, because ice viscosity is evaluated using an ad-hoc basal drag law, the contribution from basal friction induced heating is probably not accurately constrained, especially in fast-moving ice streams where the ice-bed interface is not frozen. This should impact the basal drag coefficient inversion most in fast-moving areas. However, for these areas, the basal drag coefficient will be small overall, which should mitigate the issue.
[61] Once ice viscosity is computed, we then relax the surface Dirichlet boundary conditions, so that we can carry out a basal drag inversion with a fully unconstrained forward model. However, we still need to impose ice velocity at one location to avoid the situation of a singular forward model. Here, we impose a zero ice velocity on the Geikie Plateau, East Greenland, because this is an area where few quality thickness data are available, and where therefore it is almost impossible to model ice flow realistically.
[62] For the transient simulation, an identical process is carried out to initialize the inversion of the basal friction using a BP model.

Mesh
[63] For the inversion, the Greenland Ice Sheet is meshed using static adaptation based on the InSAR surface velocities. We use 115,000 elements in the horizontal plane (62,000 vertices). Resolution is 2 km at the coast, where most of the deformation occurs, and 30 km in the interior, where ice velocity is less than a few m a À1 . For the 2D SSA, this mesh results in a system of 124,000 dofs, which is solvable on any small-scale cluster. For BP and FS, the 2D mesh is extruded vertically with 20 layers. This results in a system of 2,500,000 dofs for the BP model and 5,000,000 dofs for the 3D FS. For such large systems, we need large clusters, i.e. with at least a few hundred CPUs. For the transient run, the SeaRISE experiment relies on a different mesh, with 10 vertical layers, and 446,000 elements overall (for 270,000 vertices), resulting in a system of 540,000 dofs.
[64] Because we keep the number of layers constant over the entire GIS, large variations in size from one element to another will occur, especially vertically. Indeed, ice thickness can vary from meters to thousands of meters across the entire GIS, which represents a variation of 3 orders of magnitude. To assess the impact of bad aspect-ratios in our elements, we carried out a comprehensive study, presented in Table 2, where we varied the number of vertical layers and the vertical aspect ratio (between the smallest and longest element) for the GIS. The results show that between 2 and 20 vertical layers, the aspect ratio drops by 1 order of magnitude. The condition number for the discretized system of equations (also defined at the radius of the eigenspectrum of the discretized system of equations) correspondingly increases by almost 2 orders of magnitude for the BP model and FS models and almost 3 orders of magnitude for the 3D thermal model. While this tends to show a larger impact on the thermal model, the magnitude of the condition number  [Pattyn et al., 2008]. Linearized results from Gudmundsson [2003] are also shown in black for the no-slip bed case. for the thermal model is still 8 orders of magnitude smaller than for BP and FS models. Therefore, the impact of aspectratios on the 3D thermal model remains small. For the BP and FS models, aspect-ratios impact FS models more, which explains why iterative solvers fare better for BP models than for FS models. These results also show how increasingly difficult solving for BP and FS is with increasing vertical resolution, as the conditioning of the system matrices degrades rapidly. A compromise must therefore be found between the need for increased vertical resolution and the difficulty of solving such vertically refined models.

Performance
[65] To solve the large system of equations for the Greenland model, we used the NAS's Pleiades cluster with 30 Intel Westmere nodes, 6 CPUs per node. Each Westmere node carries 24 Gb of memory, which is necessary for the MUMPS solver to operate efficiently. Each CPU accesses 4 Gb of RAM, which is necessary for the MUMPS Symbolic LU factorization during the direct solve phase. Not all three models converged equally fast, but we stop the inversion in each case after 50 gradient descents. The run time statistics are listed in Table 3. The results demonstrate that ISSM is able to perform higher-order data assimilation at large scale under 2 hours CPU and FS diagnostic modeling at the continental scale, with data assimilation on an entire ice sheet completed in less than 18 hours CPU.
[66] In terms of scalability, Figure 6 and Table 4 present results on the parallel efficiency and computation times for a series of models of increasing size (from 100 dofs to 1,000,000 dofs) and complexity (2D SSA, 3D Thermal, 3D BP and 3D FS), using an increasing number of CPUs (1 to 180, identical to the CPUs used in the inversion of basal friction on the GIS). These results demonstrate that at low number of CPUs (<10), parallel efficiency runs at 50% or more (except for small models with less than 1,000 dofs). Here, parallel efficiency is defined as the ratio t 1 /(N t N ), where N is the number of CPUs for the run, t N the amount of computational time required for the run and t 1 the amount of computational time required for the same run using only 1 CPU. At higher number of CPUs, efficiency degrades quickly, but still remains above 10% for most models. Such evolution is to be expected for parallel direct solvers such as MUMPS. A comparative study based on a parallel iterative solver for the 2D SSA model (in red on Figure 6) shows much improved parallel efficiency, albeit with a similar decrease at higher number of CPUs. Table 4 also shows that for large sized models (100,000 dofs), run times are on the order of minutes, which demonstrates ISSM's efficiency in handling large and complex models.
[67] For our transient run, scalability studies were not carried out, given the prohibitive amount of time and resources this would require, given that we are running using a higher-order BP model. However, we can rely on the SeaRISE transient run to provide some statistics that demonstrate the applicability of ISSM to continental scale transient ice flow simulations. The model comprises 540,000 dofs which are run for 500 years, for a total of 3,000 time steps. The runs were computed on the NAS's Pleiades cluster in 40 hours, on 40 Harpertown nodes (6 CPUs per node, for a total of 240 CPUs). The statistics for the run are also presented in Table 3. This experiment proves that (1) ISSM is capable of handling higher-order transient ice flow simulations at the continental scale, using high-resolution anisotropic meshes; and (2) ISSM can evolve the geometry of an ice sheet in time. Of course, many issues remain regarding this simulation, especially in the initial setup, but such analysis of the SeaRISE experiment are outside the scope of the present study.

Results
[68] The basal drag coefficient for all three models is shown in Figure 7, with corresponding modeled and observed surface velocities in Figure 8. Differences between all three modeled velocities and observed surface velocity are shown in Figures 9a-9c. Differences between all three modeled velocities are shown in Figures 9e-9g.
[69] All three inversions show similar patterns of basal drag coefficient, however with different amplitudes. For the 2D SSA model (Figure 7a), the basal drag coefficient is the lowest, because this model does not include internal vertical deformation and thereby results in stiffer ice. Explained otherwise, 2D SSA only treats depth-averaged ice velocity that are best-fitted to surface velocities. This implies that ice velocity near the bedrock is faster, hence resulting in lower friction. In the BP model (Figure 7b), the basal drag coefficient is 50% higher in the interior as a result of increased vertical internal deformation. In 3D FS, the basal drag coefficient is another 30% higher than for the BP because the model includes all modes of vertical shearing. This is probably due to the fact that FS captures the entire stress tensor, which is equivalent to making ice softer on average for the same flow parameters. However, other reasons might be at play, such as the influence of bridging effects [Morlighem et al., 2010], or the fact that convergence might not have been fully reached for the whole ice sheet, or in certain specific areas.
[70] All three velocity calculations are in excellent agreement with observations (Figures 9a-9c). Most of the discrepancies between models and observations are near the coastline, where surface velocities are highest, reaching up to 11 km/a on the Jakobshavn Glacier. Intrinsic differences between models (Figures 9d, 9e, and 9f) are on average less than 50 m/a, demonstrating the ability of even the simplest Aspect Ratio, A r (characteristic length ratios from smallest to longest element) for vertical extrusion layers 2, 3, 4, 5, 6, 10, 15, and 20. Three models are explored: the 3D BP, 3D FS, and 3D Thermal models. The last row of the table presents the ratio of 20 vertical layers to 2 vertical layers. models to match velocities on the GIS. Largest differences are found near the coast, which is expected, because this is where we attribute most variations between lower-order and higher-order models. It is however difficult to explain variations between all three velocities to either differences between the models or variable rates of convergence of the solutions. Both factors are indeed intrinsically combined and difficult to separate. Given enough CPU time, all solutions might eventually converge to similar velocity fields, but this has not been verified given the computational challenges it would entail. This is however a conclusion for the diagnostic  Figure 6. Scalability of ISSM diagnostic (2D SSA, 3D BP and 3D FS models) and 3D thermal solutions for the Greenland Ice Sheet, expressed in terms of parallel efficiency (in % of linear scaling) vs number of processing elements (1 to 180 Intel Westmere cores on the NAS Pleiades cluster). For each type of solution, several sized models were implemented, ranging from 100 to 1 million dofs (see corresponding legend in the BP Diagnostic frame). All runs in blue rely on the MUMPS direct solver. All runs in red rely on a parallel iterative GMRES solver (preconditioned using the Additive Schwartz Method and an overlap of 2). model only. The conclusions are probably very different for transient models.
[71] The presence of fast flow features along the coast and their spatial extent inland is well captured in all three models, as well as the presence of outlet glaciers within terminal valleys only a few km wide. Of course, because our inversion is essentially an optimization to best-fit observed surface velocities, such fit is expected from our method. However, for a model such as 2D SSA, because all components of the stress tensor are not accounted for, we would expect the model to be incapable of best-fitting the observations, even when forced through an optimization algorithm. The interesting feature here is that this does not appear to be the case, and 2D SSA appears to be very efficient at capturing surface ice flow even for ice flow regimes where it is at the limit of its applicability.
[72] Figure 10 shows results for the transient simulation over 500 years based on the SeaRISE control experiment, for a higher-order BP model. Results include ice thickness, surface velocity and depth-averaged temperatures at times t = 0, t = 100 yr and t = 500 yr. They show a stable evolution of the ice sheet, with a slight increase in the ice volume of approximately 1% over 500 years. The ice sheet thins slightly in the interior, and thickens near the coastline, especially near outlet glaciers. Temperature variations can be observed on 79 North Glacier (slight decrease) and along the southwest coast, where temperature increase by about 5°C. Such variations are probably due to the initial setup of the basal friction, which takes into account a steady state instead of a transient temperature field. The SeaRISE control run relies on a constant precipitation rate, we therefore expect our transient run to be very stable. This does appear to be approximately the case for the thickness and the velocity field, but temperature variations indeed show the difficulty of correctly spinning up such transients, and the extensive work that remains in addressing such difficulties.

Discussion
[73] In the example case of the Greenland Ice Sheet, we run the flow models at an intermediate level of spatial resolution, i.e., 30 km in the interior down to 2 km at the coast. Ultimately, the models should be run at a resolution closer to one ice thickness at the coast ($600 m) and in the interior (3 km), which is also closer to the requirements in spatial resolution for grounding line dynamics [Nowicki and Wingham, 2008;Morlighem et al., 2010;Durand et al., 2009b]. The main problem for running such simulations is both the computing power needed and the lack of detailed ice thickness data at that spatial scale Morlighem et al., 2011]. In addition, if spatial resolution is decreased below one ice thickness, inverse methods may not improve model spin-up Table 4. Computation Times (in s) for a 100,000 Dofs Large Model (Based on the Greenland Ice Sheet), Using 1,2,5,10,20,50,100,150,and 180 CPUs and Running the 2D SSA, 3D BP, 3D FS, and 3D Thermal Models  2D SSA  67  39  23  17  13  12  12  12  12  3D BP  540  290  150  90  67  45  43  40  41  3D FS  1200  640  300  170  100  60  47  45  45  3D Thermal  360  200  100  71  50  39  38  any further, because at that resolution, relevant information is smeared by the ice flow model itself [Gudmundsson, 2003;Hindmarsh, 2004;Truffer, 2004].
[74] The example case of Greenland demonstrates that inverse methods are readily applicable at the continental scale with ISSM, even in the computationally demanding case of the FS. Inverse methods make ice flow models more realistic by initializing unknown boundary conditions at the ice-bed interface, at present time, without any spin up. The state of mechanical equilibrium of the ice sheet is therefore better constrained. In terms of thermal regime, our method does not ensure strict consistency between the mechanical regime (which forces the frictional heating term in the thermal model) and thermal regime. To consider the impact of this discrepancy, we evaluate the difference between the temperature field after inversion of the basal friction and the initial model setup (relying on an ad-hoc basal drag coefficient distribution). The temperature field after inversion is on average 2°C warmer over the entire GIS, with peaks at 5°C in the North-Western area (Petermann Glacier and Humboldt Glacier in particular). This amounts to a difference of about 1% over the entire GIS, which is attributable to the increased frictional heating captured by the inversion.
[75] Nevertheless, the computational cost of inversion is much lower than long (10,000 yr) transient simulations required for model spin up. This means that ISSM can be used to help and constrain long transient runs [Huybrechts et al., 1996;Ritz et al., 1997;Pollard and DeConto, 2009;Greve, 2005] so as to converge to accurate present time configurations. It can also be used to run short term (500 yr) transient flow simulations of the entire ice sheet, assuming that the pattern of basal friction shown in Figure 7 remains constant with time. As shown in our results for the SeaRISE control run, many issues remain regarding the interaction between inversion methods and spin-up of transient models. One of the most important ones is that the inversion assumes that the ice sheet is in thermal steady state, which is clearly not the case. Differences between the real temperature field and our model occur at the beginning of the run, which makes it difficult to assess how much of the subsequent evolution of the temperature field is real or model induced. Another issue is that our transient model assumes that basal friction remains constant through time. However, the hydrological network under the ice sheet is certain to evolve through time, which will impact basal friction significantly. The absence of a hydrological model in ISSM currently precludes us from capturing such changes. Finally, issues regarding the availability of datasets at the continental scale and at high resolution remain. Such datasets are currently being improved upon, notably by Operation IceBridge, but entire regions remain that cannot be constrained reliably, among them southwest Greenland and the Geikie Plateau, where reliable thickness datasets are not yet available.
[76] Data assimilation techniques are however likely to become an integral and growing part of ice sheet models, especially as more observations are gathered about ice sheets from a variety of airborne and spaceborne instruments. A similar approach has been developed for data assimilation in ocean models [Marshall et al., 1997]. In this study, data assimilation only employed ice surface velocity. In the future, similar capabilities will be developed to incorporate other variables, e.g., changes in surface elevation measured from altimeters or changes in ice mass measured from time-variable gravity.
[77] The ability of ISSM to model an entire ice sheet on a reasonable size computer cluster is not the result of a single numerical breakthrough but the result of combining multiple, complex numerical capabilities. These unique capabilities are as follow. First, ISSM uses the FEM, which makes it possible to rely on anisotropic meshes and thereby limit the number of dofs. A regular grid of Greenland at 1 km requires 7 million dofs. The same level of precision can be obtained with an unstructured mesh comprising one tenth of the elements. The gain is at least a factor 10 in solution time. Second, the architecture of ISSM is optimized for large scale data assimilation. ISSM is coded in C/C++, fully object oriented, which enables harmonious model development. ISSM relies on MPICH as a parallel communication layer, which enables it to run on large clusters, e.g., the NAS's Pleiades Cluster (more than 100,000 Intel cores). To reach full scale capability, partitioning of the elements at the earliest stages of every ice flow model is key and is achieved successfully in ISSM using METIS. Similarly, ISSM uses iterative solvers from PETSc to achieve scalability. At present and pending further development, ISSM is already capable of tackling systems of several Figure 9. (a-c) Differences between modeled and observed surface velocities (in m a À1 ) for the 2D Shelfy-Stream, 3D Blatter/Pattyn, and 3D full-Stokes models, respectively. (d-f) Differences between modeled surface velocities for the 3D Blatter/Pattyn vs 2D Shelfy-Stream, 3D full-Stokes vs 3D Blatter/ Pattyn, and 3D full-Stokes vs 2D Shelfy-Stream, respectively. millions dofs, with computation times on the order of hours for continental scale simulations, including for the FS model.
[78] Several modules are currently being developed to improve ISSM. One module will include a calving law to constrain ice front positions and their time evolution, following the work of Alley et al. [2008], Walter et al. [2010], and Bassis [2011]. Second, we will include a more sophisticated, higher-order scheme of grounding line dynamics as by Nowicki [2007] to capture ice flow instabilities, nonlinearities and hysteretic effects, which may not be well captured with our current hydrostatic criterion. Third, we will include a model of subglacial hydrology that will constrain the representation of basal friction. This would enable more realistic descriptions of the time evolution of basal drag as by Johnson [2002]. By computing the water distribution evolution in time, it would be possible to assess the value of the effective water pressure, which is an integral part of any basal friction law [Johnson, 2002;Le Brocq et al., 2009]. Indeed, water at the ice-bed interface increases lubrication of bedrock, which in turn increases basal sliding [Weertman, 1957]. Finally, ISSM will include a capability to automatically compute the adjoint state of a model using automatic differentiation compilers, as pioneered by ocean models [Marshall et al., 1997] and compilers such as OpenAD [Utke et al., 2008]. These compilers calculate the adjoint for any model for which a simple analytical solution might not exist. This capability would enable ISSM to assimilate data sets from a variety of sensors and platforms to reduce model uncertainties. We will rely on ADIC2 [Narayanan et al., 2010], to try and address issues that compilers such as OpenAD currently encounter with C++ object orientation and integration of complex solver libraries such as PETSc.

Conclusions
[79] The capability to run large scale, high spatial resolution, higher-order ice flow models, constrained using inverse methods, is critical to improving reconstructions or Figure 10. Transient ice flow simulation of the Greenland Ice Sheet, using the higher-order BP model, from present-time to 500 years into the future. (a-c) Thickness, H; velocity, V; and depth-averaged temperature, T avg at time t = 0, as initialized by the basal friction inversion. (d-f) Thickness, H; velocity, V; and depth-averaged temperature, T avg at time t = 100 years. (g-i) Thickness, H; velocity, V; and depthaveraged temperature, T avg at the end of the transient run, t = 500 years.
projections of the mass balance of large ice sheets. The Ice Sheet System Model (ISSM) offers a solution to this problem by using finite element methods, anisotropic static mesh adaptation, scalable solvers and parallel technologies. Unknown model parameters such as basal drag are constrained using inverse methods applied to satellite-derived observations of surface motion. ISSM offers varying orders of complexity ranging from simple 2D models to the full-Stokes 3D solution. ISSM is user friendly, highly modular, and may be run on massive computer clusters or a single CPU within MATLAB. The implementation of ISSM is validated with standard benchmarks and its applicability at the large scale is demonstrated in the case of Greenland using higher-order and full-Stokes models. ISSM is now public domain, released under a three clause BSD License, available at http://issm.jpl.nasa. gov. This will provide a pathway for the development of more advanced, more realistic numerical models of ice sheet flow that will improve our capability to project ice sheet evolution.

Appendix A: Software Architecture and Management
[80] ISSM relies on the C language [Kernighan, 1988] for the numerical implementation of finite elements. For management of all objects, the C++ language [Stroustrup, 1997] is used, which strongly relies on polymorphic capabilities. This allows ISSM to be both flexible (new finite elements are 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).
[81] To improve ease of use, ISSM is hosted in MATLAB (MathWorks, software, 2008), a common scientific platform. The C/C++ core is interfaced to the MATLAB environment using the External API. This ensures that pre and post-processing can be carried out using MATLAB modules (also called MEX functions). These modules behave as standard routines, but they encapsulate intrinsic ISSM capabilities. The result is a seamless integration of ISSM's ice flow model within MATLAB.
[82] ISSM is first and foremost a parallel architecture. While it can be used in serial mode within MATLAB, it is designed to run massively parallel computations on large clusters, e.g., the NASA Advanced Supercomputing (NAS) Pleiades cluster. When running on these clusters, ISSM relies on its C/C++ core compiled as a stand-alone executable. Parallelism is achieved through the Message Passing Interface . This library is flexible enough to allow runs on distributed as well as shared memory clusters.
[83] The numerics are implemented with PETSc. 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 is used in serial or parallel mode the same way. PETSc also provides access to a wide array of direct and iterative solvers, proprietary as well as external, along with corresponding preconditioners.
[84] In addition to MATLAB and PETSc, ISSM relies on a wide array of external packages, of which we only list the most important ones.
2. VALGRIND: instrumentation framework for building dynamic analysis tools [Nethercote and Seward, 2007; see also Valgrind, 2011, http://valgrind.org]. This package is used to find memory leaks in the code, as well as for optimizing memory management.
3. METIS: Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices [Karypis and Kumar, 1998]. This package is used to partition objects such as elements and vertices across a cluster. This partitioning scheme results in partitions that have equal numbers of elements on each cluster node. This also ensures a good computational load balance, because resulting stiffness matrices are wellpartitioned, which results in tremendous gains during the solve phase.
[85] ISSM also provides a suite of tools to automatically carry out nightly runs. These nightly runs consist of approximately one hundred test cases based on a variety of configurations that range from synthetic square ice shelves and ice sheets to simplified versions of Pine Island and Nioghalvfjerdsfjorden glaciers. More than 2000 output fields are tested against archived files and an HTML report is automatically generated and sent by email every night to all developers, for verification and validation. This ensures that no error is introduced in the software during capability development.
[86] To ensure easiness of installation of our suite of tools, we rely on Autotools [Vaughan et al., 2000], a software package for automatic installation of makefiles. ISSM also uses a versioning system for hosting the code [Pilato et al., 2008], which allows multiple developers to synchronize development of the code and regular maintenance. Finally, ISSM comes with a full-fledged documentation, which is compiled at each new installation and is synchronized on the ISSM website (http://issm.jpl.nasa.gov).