Capabilities and performance of Elmer / Ice , a new-generation ice sheet model

The Fourth IPCC Assessment Report concluded that ice sheet flow models, in their current state, were unable to provide accurate forecast for the increase of polar ice sheet discharge and the associated contribution to sea level rise. Since then, the glaciological community has undertaken a huge effort to develop and improve a new generation of ice flow models, and as a result a significant number of new ice sheet models have emerged. Among them is the parallel finite-element model Elmer/Ice, based on the opensource multi-physics code Elmer. It was one of the first fullStokes models used to make projections for the evolution of the whole Greenland ice sheet for the coming two centuries. Originally developed to solve local ice flow problems of high mechanical and physical complexity, Elmer/Ice has today reached the maturity to solve larger-scale problems, earning the status of an ice sheet model. Here, we summarise almost 10 yr of development performed by different groups. Elmer/Ice solves the full-Stokes equations, for isotropic but also anisotropic ice rheology, resolves the grounding line dynamics as a contact problem, and contains various basal friction laws. Derived fields, like the age of the ice, the strain rate or stress, can also be computed. Elmer/Ice includes two recently proposed inverse methods to infer badly known parameters. Elmer is a highly parallelised code thanks to recent developments and the implementation of a block preconditioned solver for the Stokes system. In this paper, all these components are presented in detail, as well as the numerical performance of the Stokes solver and developments planned for the future.


Introduction
Since the 2007 IPCC report (Solomon et al., 2007), theoretical glaciology has taken a big leap towards improved ice sheet flow models, in order to provide reliable future estimates of the dynamical contribution of ice sheets to sea level rise.These models were originally designed to reconstruct the evolution of ice sheets over past glaciological cycles, neglecting short-term responses and local features.The new challenge of running ice sheet models to provide estimates of future sea level rise has created the need for a new generation of ice sheet models (Vaughan and Arthern, 2007;Gillet-Chaulet and Durand, 2010;Blatter et al., 2011;Kirchner et al., 2011;Alley and Joughin, 2012).This new generation of ice sheet models includes a set of requisites that are essential to provide a sufficiently accurate description of the ice flow dynamics.

O. Gagliardini et al.: Elmer/Ice model
As a first requisite, these models must be able to describe the ice flow heterogeneity, and particularly the major contribution of individual ice streams to the total ice discharge.This requires the use of an unstructured mesh in the horizontal plane (e.g.Gillet-Chaulet et al., 2012;Larour et al., 2012;Seddik et al., 2012) or of adaptive multi-grid methods (Cornford et al., 2013b).These mesh techniques are essential to produce hundred-metre-scale grid sizes in areas of interest, especially near the coast, while for the interior regions where variations in velocity gradients are small, classic grid sizes can be kept to save computing resources.Grid refinement is even more essential when considering the dynamics of the grounding line, i.e. the boundary between the grounded ice sheet and the floating ice shelf, because a grid size that is too large gives inconsistent grounding line dynamics (Durand et al., 2009;Pattyn et al., 2013).
The second important requisite is to have an accurate description of the complex state of stress prevailing in ice streams to solve the full-Stokes system, or at least to adopt a higher-order asymptotic formulation.As shown by the ISMIP-HOM inter-comparison exercise (Pattyn et al., 2008), higher-order models are needed to describe the ice flow in areas where the basal topography and slipperiness vary greatly, which are generally the most dynamic regions within ice sheets.Higher-order models are also necessary to properly describe the dynamics of the grounding line.The MISMIP inter-comparison (Pattyn et al., 2012) indicated the need to solve the full-Stokes equations near the grounding line to obtain fully accurate results.
The consequences of these first two requisites, i.e. high numerical resolution at places of interest and higher-order formulations, are a high computing cost and the necessity to develop parallel codes, able to run on hundreds of CPUs.Recent studies (Larour et al., 2012;Gillet-Chaulet et al., 2012;Seddik et al., 2012;Cornford et al., 2013b) have fulfilled these requirements and have shown that by deploying highperformance computing (HPC) techniques this challenge can be successfully taken on.In this context, Elmer/Ice takes advantage of being backed by a large open-source community that also develops new numerical and HPC techniques for the code (e.g.Malinen, 2007).
The third requisite, and from the physical point of view the most challenging, is to implement physically founded boundary conditions.These improvements are far more complex and it will take more time to fully address them in the ice sheet flow models.The recently observed changes in coastal glacier dynamics (e.g.Moon et al., 2012) are certainly driven by changes in ice sheet and ice shelf boundary conditions, and consequently linked to changes in the ocean and atmosphere components of the climatic system.In the simplest cases, changes in the climatic components directly drive the changes at the boundaries of the ice mass.This is the case for surface air temperature or ocean temperature which directly drive the temperature boundary condition of the upper surface or the bottom ice/ocean interface, respectively.In other more complex cases, the link between changes in the ocean and/or atmosphere and changes in the ice flow is indirect.Intermediate processes (often not observable) are involved, as in the case for example of the link between surface runoff and basal sliding or ocean temperature and calving rate.Thus, a dedicated model is required to describe the processes responsible for the transfer of these changes to the ice mass.Driving this dedicated transfer model might require coupling the ice sheet model with an atmosphere or an ocean model.
The last important requisite for a forecast model is to be able to simulate present-day observations with as much fidelity as possible (Aschwanden et al., 2013).This point must be addressed clearly using data assimilation techniques and specific inverse methods to estimate the less well known parameters of the model (e.g.Heimbach and Bugnion, 2009;Arthern and Gudmundsson, 2010;Morlighem et al., 2010).
Recent ice sheet model developments have started to fulfil some of these priority requisites, leading the way toward the new generation of ice sheet models (Bueler and Brown, 2009;Pollard and DeConto, 2009;Rutt et al., 2009;Larour et al., 2012;Leng et al., 2012;Winkelmann et al., 2011;Favier et al., 2012;Gillet-Chaulet et al., 2012).Among them, the Elmer/Ice model already includes many of these requisites.Elmer/Ice is the glaciological extension of Elmer, the opensource finite element (FE) software developed by CSC in Finland (http://www.csc.fi/elmer/).Elmer is a multi-physics code base from which it was possible to develop new specialised modules for computational glaciology while maintaining the compatibility with the main Elmer distribution.Thus, Elmer/Ice still benefits from the developments of the standard Elmer distribution.In this paper, for simplicity we refer to Elmer/Ice even if some of the features described belong to the main Elmer distribution.Elmer/Ice was not originally designed as an ice sheet model since the first applications were restricted to local areas of interest or glaciers (Le Meur et al., 2004;Zwinger et al., 2007;Zwinger and Moore, 2009).Elmer/Ice was primarily developed to solve the flow of anisotropic polar ice and the evolution of its strain-induced fabric (Gillet-Chaulet et al., 2006;Durand et al., 2007;Seddik et al., 2008Seddik et al., , 2011;;Ma et al., 2010;Martín and Gudmundsson, 2012).It has since then been used to model the flow of a cold firn-covered glacier using a dedicated snow/firn rheological law (Zwinger et al., 2007).Elmer/Ice has been the only full-Stokes model to perform the whole set of the ISMIP-HOM experiments (Gagliardini and Zwinger, 2008;Pattyn et al., 2008) and is still the only full-Stokes model to participate in the grounding line experiments MISMIP (Pattyn et al., 2012).Elmer/Ice was further used as a reference for the later MISMIP3d experiments (Pattyn et al., 2013).Recently, data assimilation was implemented within Elmer/Ice (Jay-Allemand et al., 2011;Schäfer et al., 2012;Gillet-Chaulet et al., 2012) to infer poorly known parameters such as basal drag.Today, Elmer/Ice is the only three-dimensional full-Stokes model that solves the grounding line dynamics (Favier et al., 2012), and it will be the only full-Stokes model able to run forecast simulations for the whole Greenland ice sheet for the coming AR5 IPPC report, in the framework of both SeaRISE (Seddik et al., 2012) and ice2sea (Gillet-Chaulet et al., 2012;Shannon et al., 2013;Edwards et al., 2013) programmes.
In this paper, we summarise ten years of consistent developments and present the current state of the new-generation ice sheet model Elmer/Ice (Elmer library version 7.0 SVN revision 5955).We only focus on the past developments that are relevant for simulations of three-dimensional ice sheets.Specific developments regarding two-dimensional flow line or glacier applications are not presented here, but one can consult previous publications on these types of applications (the complete list of Elmer/Ice publications can be found on http://elmerice.elmerfem.org/).Section 2 presents the governing equations implemented in Elmer/Ice.The associated boundary conditions are discussed in Sect.3. Other useful equations, such as the equation to evaluate the age of the ice, are presented in Sect. 4. Section 5 is dedicated to the inverse methods implemented in Elmer/Ice.Some technical aspects related to the resolution of these equations in the framework of the FE method are discussed in Sect.6.The efficiency of Elmer/Ice was verified by standard convergence and scalability tests described in Sect.7. Finally, we provide some insights into the future planned developments in Sect.8.

Ice flow equations
Ice is a fluid with an extremely high viscosity that flows very slowly so that inertia and acceleration terms entering the momentum equation can be neglected.Therefore, the three-dimensional velocity field and the pressure field of an ice mass flowing under gravity are obtained by solving the Stokes equations over the ice volume .The Stokes equations express the conservation of linear momentum divσ and the mass conservation In these equations, ρ is the ice density, g = (0, 0, −g) the gravity vector, u = (u, v, w) the ice velocity vector, σ = τ − pI the Cauchy stress tensor with p = −trσ /3 the isotropic pressure, τ the deviatoric stress tensor and I the identity matrix.This system of equations of unknowns u and p is closed by adopting one of the rheological laws presented in the next section.The conditions that are applied on the boundary of the volume are discussed in Sect.3.

Rheological laws for polar ice
Even if most ice sheet models assume an isotropic rheological law for ice, it is well known that the viscous response of polar ice can be strongly anisotropic, and that this response depends on the crystal orientation distribution, i.e. the ice fabric (e.g.Gagliardini et al., 2009).Elmer/Ice includes the classic isotropic Glen's flow law as well as two anisotropic flow laws.As shown in various applications, the anisotropy of polar ice has a strong influence on the overall flow (Zwinger et al., 2013) and will in turn modify the age-depth relationship (Gillet-Chaulet et al., 2006;Seddik et al., 2011).In central parts of ice sheets, ice anisotropy and the development of fabric can explain the observed heterogeneity of the ice deformation along a drilling (Durand et al., 2007).On the coastal area, due to the large contrast of the stress regimes for the grounded part and for the ice shelf, the ice anisotropy induces an apparent hardening of the ice up to a factor 10 when ice moves from grounded to floating (Ma et al., 2010).
When ice is assumed to behave as an isotropic material, its rheology is given by a Norton-Hoff power law, known as Glen's law in glaciology, which links the deviatoric stress τ with the strain rate ε: where the effective viscosity η is defined as In Eq. ( 4), ε2 e = tr(ε 2 )/2 is the square of the second invariant of the strain rate and A = A(T ) is a rheological parameter which depends on T , the ice temperature relative to the pressure melting point, via an Arrhenius law.The enhancement factor E in Eq. ( 4) is often used to account for anisotropy effects, by prescribing an ad hoc value depending on the ice age and/or type of flow.Due to the state of stress, E is expected to be greater than 1 for grounded ice of polar ice sheets, whereas a value lower than 1 should be used for floating ice shelves (Ma et al., 2010).A compressible form of Glen's law (Gagliardini and Meyssonnier, 1997), well adapted to describe the flow of firn, is also implemented in Elmer/Ice (Zwinger et al., 2007).
Both implemented anisotropic flow laws depend on the ice polycrystalline fabric, which is described by its second-and fourth-order orientation tensors a (2) and a (4) , respectively, defined as where c is the crystal c axis unit vector and denotes the average over all the grains that compose the polycrystal.By definition odd-order orientations tensors are null, and the higher the order of the orientation tensor the better the description of the fabric.However, it can be shown that with a linear flow law, knowing the second-and fourth-order orientation tensors is sufficient to uniquely define the macroscopic flow law (Gillet-Chaulet et al., 2005 (2) 33 = 1/3, for a single maximum fabric with its maximum in the third direction, a (2) 22 < 1/3, and for a girdle-type fabric in the plane (x 1 , x 2 ), a (2) 33 < 1/3 and a (2) 22 > 1/3.In addition to three eigenvalues, three Euler angles are necessary to uniquely define a (2) with respect to a general reference frame.It can be shown analytically with a linear flow that if the second-and fourth-order orientation tensors have the same eigenframe, the polycrystal behaviour will exhibit orthotropic symmetries (Gillet-Chaulet et al., 2006).The equations for the fabric evolution are presented in Sect.2.5.
The first anisotropic flow law implemented in Elmer/Ice is the non-linear General Orthotropic Flow Law (GOLF, Gillet-Chaulet et al., 2005;Ma et al., 2010).The GOLF provides a non-collinear and non-linear relation between strain rate and stress, using the concept of structure tensors.In its initial form, the ice was assumed to behave as a linearly viscous orthotropic material.In more recent works (Martín et al., 2009;Ma et al., 2010), the GOLF has been extended to a nonlinear form by adding an invariant in the anisotropic linear law.The simplest choice is either to add the second invariant of the strain rate εe (Martín et al., 2009) or the second invariant of the deviatoric stress τ e (with τ 2 e = tr(τ 2 )/2, Pettit et al., 2007;Ma et al., 2010).No theoretical or experimental results are available today to discard one of these two solutions, and other solutions based on anisotropic invariants of the deviatoric stress and/or the strain rate are also possible.In Elmer/Ice, both solutions are implemented.Using the second invariant of the deviatoric stress, for a given fabric and a given state of stress, the corresponding strain rate relative to the isotropic response is the same for the linear and nonlinear cases.Using the strain-rate invariant in the same way as Martín et al. (2009) leads to an opposite definition of the anisotropy ratios: for a given strain rate, the corresponding stress relative to the isotropic response is the same for the linear and non-linear cases.When using the stress second invariant, the GOLF reads The six dimensionless anisotropy viscosities η r (a (2) ) and η r+3 (a (2) ) (r = 1, 2, 3) are functions of eigenvalues of the second-order orientation tensor a (2) , which represent a measure of the anisotropy strength.The three structure tensors M r are given by the dyadic products of the three eigenvectors of a (2) , which then represent the material symmetry axes.In the method proposed by Gillet-Chaulet et al. (2006), the six dimensionless viscosities η r (a (2) ) are tabulated as a function of the fabric strength (i.e. the a (2) i ) using a micro-macro model.Various micro-macro models, from the assumption of uniform stress within the ice polycrystal to the assumption of uniform strain rate, as well as different crystal anisotropy can be used to tabulate the six viscosities η r .The most realistic polycrystalline response is obtained using the viscoplastic self-consistent model (VPSC, Castelnau et al., 1996, 1998), with the two crystal anisotropy parameters chosen so that the experimentally observed polycrystal anisotropy is reproduced (Gillet-Chaulet et al., 2006;Ma et al., 2010).When the ice is isotropic, η r = 0 and η r+3 = 1 (r = 1, 2, 3), then the GOLF (6) reduces to Glen's isotropic flow law (3) with E = 1.
The second anisotropic flow law implemented in Elmer/Ice is the Continuum-mechanical Anisotropic Flow model based on an anisotropic Flow Enhancement factor (CAFFE, Seddik et al., 2008;Placidi et al., 2010).The CAFFE model assumes collinearity between the strain rate and deviatoric stress tensors, so that the general form of Glen's law (3) is not modified, but the enhancement factor E is a function of the polycrystalline deformability D such that with The polycrystalline deformability D is a function of strain rate and fabric.When D = 0, the minimal enhancement factor E min is reached, which corresponds to an uni-axial compression on a single maximum fabric.For an isotropic fabric, D = 1 and the response is identical whatever the strain rate, whereas the maximal enhancement E max is obtained for D = 5/2, which corresponds to a single maximum fabric undergoing simple shearing.The adopted form for the polycrystalline deformability, which verifies the above criteria, reads (9)

Evolution of the surface boundaries
For transient simulations, the upper and lower boundaries of the domain are allowed to evolve, following an advection equation.Evolution of the upper surface z = z s (x, y, t) is given by where (u s , v s , w s ) are the surface velocities obtained from the Stokes solution and a s = a s (x, y, t) is the accumulation/ablation prescribed as a vertical component only.
Elmer/Ice provides a surface melting parameterisation based on positive degree-day (PDD) method (Reeh, 1991), supplemented by the semi-analytical solution for the PDD integral by Calov and Greve (2005) (Seddik et al., 2012).The accumulation/ablation distribution can also be inferred from a regional climate model either directly as in Gillet-Chaulet et al. (2012) and Shannon et al. (2013) or using a surface elevation parameterisation as in Edwards et al. (2013).
The lower surface of an ice sheet is either in contact with the bedrock or the ocean.The evolution of the lower surface z = z b (x, y, t) is given as where (u b , v b , w b ) are the basal velocities and a b⊥ = a b⊥ (x, y, t) is the melting/accretion function, taken perpendicular to the surface.Assuming a rigid, impenetrable bedrock z = b(x, y), the following topological conditions must be fulfilled by z s and z b : The weak formulation of Eq. ( 10) or Eq. ( 11), in combination with the constraints (12) forms a variational inequality.Technically, it is solved using a method of imposed Dirichlet conditions that are released by a criterion based on the residual, as described in Sect.6.5.In Gagliardini et al. (2010), melting below the ice shelf was prescribed using a parameterised expression following Walker et al. (2008).As discussed in Sect.8, a proper description of the basal melting below ice shelves will certainly require the coupling of Elmer/Ice with an ocean model or at least the implementation of a plume-type model.
The margin boundary of an ice sheet is either land-or marine-terminated, depending on whether the bedrock elevation at the ice front is located above or below sea level, respectively.In both cases, the front position evolves with time and its evolution is governed by the imbalance between ice flux and ablation/basal melting/calving processes.Land-terminated fronts can be treated classically by adopting a minimal ice thickness h min , so that the exact condition ( 12) is replaced by the less strict one z s (x, y, t) ≥ b(x, y) + h min (and Where the ice sheet is marine-terminated, this type of treatment cannot be applied because the sea water pressure and lateral buttressing forces would not be correctly taken into account.The front boundary of a marine-terminated ice sheet must therefore be allowed to move over time, as a function of the calving rate and ice flux at the margin. Assuming that the calving front is a vertical surface, it can be described by the implicit function F c (x, y, t) = 0 (Greve and Blatter, 2009).Denoting by grad F c = (∂F c /∂x, ∂F c /∂y, 0) its gradient, N c = | grad F c | the norm and n c = grad F c /N c the unit normal vector (assumed to point out of the ice), the calving front evolves as follows where c ⊥ is the calving rate.The latter is defined as the ice volume flux across the calving front, c ⊥ = (u − w c ) • n c , where w c is the kinematic velocity of the calving front (Greve and Blatter, 2009).Implementation of calving laws to evaluate the calving rate c ⊥ is part of the developments currently ongoing in Elmer/Ice, as discussed more in details in Sect.8. Moving the mesh both vertically (upper and lower surface) and horizontally (calving front) induces additional terms in the convection part of equations and in turn technical issues that are discussed in Sect.6.1.

Heat equation
The temperature within the ice is obtained from the general balance equation of internal energy and reads where κ = κ(T ) and c v = c v (T ) are the heat conductivity and specific heat of ice, respectively.The last term in the heat equation represents the amount of energy produced by viscous deformation.The ice temperature T is bounded by the pressure melting point T m , so that T ≤ T m , or equivalently T ≤ 0, with T = T − T m being the homologous temperature entering the Arrhenius law to estimate Glen's parameter in Eqs. ( 4) and ( 6).This inequality, as well as temperaturedependent material properties, make the solution of the heat transfer equation a non-linear problem which is solved using an iterative method as presented in Sect.6.5.

Fabric description and evolution
Assuming that recrystallisation processes do not occur and that the ice fabric is induced solely by deformation, the evolution of the second-order orientation tensor a (2) defined by Eq. ( 5) can be written as where W is the spin tensor defined as the antisymmetric part of the velocity gradient.The tensor C is defined as The interaction parameter α controls the relative weighting of the strain rate ε and the deviatoric stress τ in the fabric evolution Eq. ( 15).When α = 0, the fabric evolution is solely controlled by the state of strain rate, whereas in the case where α = 1 the fabric evolves under the influence of the deviatoric stress solely.In between, as for the VPSC, both the strain rate and deviatoric stress contribute to the fabric evolution.Assuming ι = 1, an interaction parameter α = 0.06 is in accordance with the crystal anisotropy and the VPSC model used to derive the polycrystalline behaviour (Gillet-Chaulet  et al., 2006).Seddik et al. (2008Seddik et al. ( , 2011) ) adopted instead α = 0 and a value of ι lower than 1.In Eq. ( 15), the fourth-order orientation tensor is evaluated assuming a closure approximation giving a (4) as a tensorial function of a (2)  (Chung and Kwon, 2002;Gillet-Chaulet et al., 2006).Theoretically, recrystallisation processes, such as continuous and migration recrystallisation, can be included by adding terms in Eq. ( 15) to parameterise on the polycrystalline scale the phenomena occurring at the grain scale (Seddik et al., 2011).Because experimental data are currently missing, these parameterisations have not yet been validated and are not presented here.

Boundary conditions
For all the equations presented above, classic Dirichlet, Neumann, Robin, symmetric and periodic boundary conditions can be applied on the boundary of the domain.In this section, we present the conditions to be applied on the different boundaries of an ice sheet for the main equations presented above, and we focus more specifically on the treatment of the basal boundary.

Ice/atmosphere boundary
The upper free surface z = z s (x, y, t), also denoted s , is in contact with the atmosphere and is therefore a stress-free surface, so that where n s is the normal outward-pointing unit vector to the free surface.For the dating equation, fabric equations and all other transport equations, Dirichlet conditions are applied on the upper surface only where the ice velocity enters the domain (mainly in the accumulation area).Where z = z s and u • n s ≤ 0, the temperature is equal to the imposed surface temperature, T (x, y, z s , t) = T s (x, y, t), and the fabric is assumed to be isotropic, a (2) (x, y, z s , t) = I/3.For the heat equation, a heat flux can be imposed at the upper surface to account for melt-water refreezing.

Ice/bedrock boundary
The lower interface z = z b (x, y, t), also denoted b , may be in contact with either the sea or the bedrock, so two kinds of boundary conditions coexist on a single surface.The conditions to be applied where the ice is in contact with the sea are presented in the next section.Where the ice is in contact with the bedrock (i.e.z b = b), the following conditions apply: where ) are the basal shear stresses and basal velocities, respectively, defined in terms of tangent vectors t i and normal outward-pointing unit vector to the bedrock n b .Note that the boundary condition Eq. ( 18) for the Stokes problem is equivalent to the free-surface Eq. ( 11).The effective pressure N is defined as the difference between the ice normal stress and the water pressure, such as N = −σ nn − p w with σ nn = n b • σ n b .Equation ( 18) is the no-penetration condition accounting for basal melting (a b⊥ < 0) or basal accretion (a ⊥b > 0), whereas Eq. ( 19) stands for the general form of a friction law.When f f = 0, the ice slides perfectly over the bedrock, whereas when f f → +∞ basal sliding is null.The three different friction laws implemented in Elmer/Ice are presented below.
The first friction law linearly relates the basal shear stress to the basal velocity, such as where β ≥ 0 is the basal friction parameter.As shown later, this simple law is used for data assimilation and in this case β is a control parameter.
The second law implemented in Elmer/Ice is a Weertmantype sliding law: where u b is the norm of the sliding velocity β m is a sliding parameter and m an exponent.When m = 1, the Weertman-type friction law Eq.( 21) reduces to the linear law Eq.( 20).Theoretically, in the case of ice sliding without cavitation over an undulating bed, m is equal to 1/n (Lliboutry, 1968), where n is Glen's law exponent.
The third friction was proposed by Schoof (2005) from mathematical expansions and by Gagliardini et al. (2007) from FE simulations.This law describes the flow of clean ice over a rigid bedrock when cavitation is likely to occur: where χ = 1/(C n N n A s ), α q = (q −1) q−1 /q q , A s is the sliding parameter in the absence of cavitation and n Glen's law exponent, resulting in a non-linear relation between the basal drag σ nt i and the basal sliding velocity u t i .The maximal value of σ nt i is C and the exponent q ≥ 1 controls the post-peak decrease.When the post-peak exponent q is equal to 1, the basal drag tends asymptotically to its maximum value C (no post-peak decrease).Note that in the limit case where N 0, the sliding parameter A s and the friction parameter β m are inversely proportional.As shown by Schoof (2005), the coefficient C should be chosen smaller than the maximum local positive slope of the bedrock topography at a decimetre to metre scale, so that the ratio σ nt i /N ≤ C fulfills Iken's bound (Iken, 1981).The friction law Eq.( 22) is strongly related to the water pressure p w through the effective pressure N. The law Eq.( 22) can then be used to couple the hydrology and the ice dynamics.The hydrological model and its implementation in Elmer/Ice are presented in de Fleurian et al. (2013).
For the heat equation, the geothermal heat flux q geo is imposed where the basal temperature is lower than the pressure melting point (T < T m or T < 0), and the following Neumann-type boundary condition applies: where |σ nt i u t i | is the heat energy induced by basal friction.
Where the temperature melting point is reached (T = T m ), the amount of melted water is estimated from the imbalance of heat fluxes and surface production: where L is the latent heat of ice.

Ice/sea boundary
At the bottom surface z = z b (x, y, t) where the ice is in contact with the ocean (i.e.z b > b) and at the front of the ice sheet, the normal stress is equal to the sea pressure p w (z, t), which evolves vertically as follows: where ρ w is the sea water density and l w the sea level.The Neumann condition applied on these ice/ocean interfaces is thus

Grounding line dynamics
The position of the grounding line is part of the solution and can evolve with time.Its position at each time step is determined by solving a contact problem.The contact is tested by comparing at each node where z b = b the normal force R n exerted by the ice on the bedrock and the equivalent water force F w .R n is directly evaluated from the residual of the Stokes system, whereas F w is obtained by integrating the water pressure over the boundary elements using the boundary element shape functions.Then, if R n > F w and z b = b, the boundary conditions Eqs. ( 18) and ( 19) apply; whereas if R n = F w and z b = b, or z b > b, the boundary condition Eq. ( 26) applies instead.

Auxiliary equations
The goal of an ice sheet simulation, usually, is to obtain information on either the geometry, the age/depth relationship or simply the exerted stresses and forces on a particular surface in contact with the ice.This section introduces the methods needed to obtain such information.

Age equation
The age A of the ice at each point of the ice sheet domain is obtained by solving the following equation: where z = z s and u • n s ≤ 0, the age of the ice is zero, i.e.A(x, y, z s,t ) = 0 (Zwinger and Moore, 2009).By solving the age equation we can compute isochrones and determine dating as a function of depth at an ice core (drilled or planned) location.Input parameters entering other equations might also be age-dependent, such as the enhancement factor for example.

Depth and elevation
It is often very useful to know the depth below the upper surface or the height above the bedrock at each point of the ice sheet domain.For example, it can be used to prescribe parameterisation of the temperature or the ice fabric fields as a function of depth.With the FE method, using unstructured meshes, the depth d(x, y, z, t) = z s − z or the height h(x, y, z, t) = z − z b at any point M(x, y, z) cannot be estimated directly because nodes are not necessarily vertically aligned.Therefore, we compute the depth d (or equivalently height h) field by solving the following equations: with the boundary conditions d = 0 on z = z s or h = 0 on z = z b .Effectively, we solve, here for the height h, the following system: with the boundary condition h| b = 0 and the unity vector e z in the vertical direction.The variational form is obtained after integrating Eq. ( 29) by parts and accounting for the boundary condition Eq. ( 30), leading to a degenerated Laplace equation of the form

Stress and strain rate
Elmer/Ice includes solvers to compute the Cauchy stress, deviatoric stress or strain rate fields from the Stokes solution, and also includes eigenvalues of these tensor variables.
This results in solving six independent equations, one for each of the six independent components of the stress tensor.In a similar manner, components εij of the nodal strain rate tensor are obtained from the following variational form: (33)

Inverse methods within Elmer/Ice
The ice effective viscosity η(x, y, z) in Eq. ( 3) and the basal friction coefficient β(x, y) in Eq. ( 20) are two particularly important input fields when modelling the flow of real glaciological systems.However, these two parameters are used to represent complex processes, and their values in situ are poorly constrained and can vary by several orders of magnitude with time and space.On the other hand, our knowledge of some of the outputs of the model (surface velocity, surface elevations) has considerably increased recently with data acquired by remote spatial observation.Two variational inverse methods have been implemented within Elmer/Ice to constrain η(x, y, z) and β(x, y) in diagnostic simulations from topography and surface horizontal velocity data.Both methods are based on minimising a cost function that measures the mismatch between the model and the observations.The two methods are briefly described below and their implementation in Elmer/Ice is verified in Sect.7.

Robin inverse method
This method, initially proposed by Arthern and Gudmundsson (2010), consists in solving alternatively the natural Neumann-type problem, defined by Eqs.(1) and ( 2) and the surface boundary conditions (17), and the associated Dirichlet-type problem, defined by the same equations except that the Neumann upper-surface condition Eq. ( 17) is replaced by a Dirichlet condition where observed surface horizontal velocities are imposed, such that u = u obs and v = v obs for z = z s . (34) The cost function that expresses the mismatch between the solutions of the two models is given by where superscripts N and D refer to the Neumann and Dirichlet problem solutions, respectively.The Gâteaux derivatives of the cost function J o with respect to the parameters η and β for perturbations η and β , respectively, are given by where the symbol ε2 e denotes the square of the second invariant of the strain rate as defined for Eq. ( 4) and | • | defines the norm of the velocity vector.Note that this derivative is exact only for a linear rheology and thus is only an approximation of the true derivative of the cost function when using Glen's flow law Eq.( 3) with n > 1 in Eq. (4).

Control inverse method
For a linear isotropic rheology (a scalar viscosity η independent of the velocity, i.e. n = 1 in Eq. 4), the Stokes system of equations is self-adjoint.Denoting by λ and q the adjoint variables corresponding to u and p, respectively, they are solutions of the following equations: where ελ is the equivalent of the strain rate tensor constructed with λ.For a non-linear rheology, the operator used by the forward solver (Stokes operator) remains self-adjoint when equipped with the Newton linearisation (Petra et al., 2012).
The cost function is chosen to measure the mismatch between the modelled and observed surface velocities where j is the mismatch measure function and u obs are the observed surface velocities.The choice of j can be casedependent and will affect the boundary condition terms of the adjoint system.For example, as the surface velocity direction is mainly governed by topography, we can discard the error on the velocity direction and express j as where subscript H refers to the horizontal component of the velocity vectors (Gillet-Chaulet et al., 2012).The Gâteaux derivatives of Eq. ( 40) with respect to η and β are obtained as follows:

Regularisation
When working with non-perfect (noisy) data, it is necessary to add a regularisation term in the cost function to improve the conditioning of the inverse problem and ensure the existence of a unique minimum.The regularisation term is based on a priori information on the solution either from measurements, from analytical solutions (Raymond Pralong and Gudmundsson, 2011), or from assumptions on the spatial variations of the variable.In Elmer/Ice, a smoothness constraint on a variable α can be imposed in the form of a Tikhonov regularisation penalising the first spatial derivatives of α as in Morlighem et al. (2010), Jay-Allemand et al. ( 2011), and Gillet-Chaulet et al. ( 2012): The Gâteaux derivative of J reg with respect to α for a perturbation α is obtained by The total cost function to minimise then reads where λ is a positive ad hoc parameter.The cost function minimum is therefore no longer the best possible fit to observations, but a compromise (through the tuning of λ) between fitting with observations and smoothness in α.

Minimisation
The Gâteaux derivatives of J o are given by a continuous scalar product represented by the integral terms in Eqs. ( 37) and ( 42).When discretized on the FE mesh, these equations are transformed into a discrete Euclidean product as follows: where γ represents η or β, ∇ γ J o is the continuous Fréchet derivative of J o , the expression of which is given by comparison with Eqs. ( 37) and ( 42), ∇ i γ J o is its value at mesh node i = (1, . . ., N p ) and W i is the nodal weight associated with node i and computed following the standard integration scheme.The sum of all weights is the volume (or area) of the FE mesh.The discrete gradients of J o at each mesh node used for the minimisation are then given by W i ∇ i γ J o and account for the volume or area surrounding each node.
The minimisation of the cost function J o with respect to η i or β i is done using the limited memory quasi-Newton routine M1QN3 (Gilbert and Lemaréchal, 1989) implemented in Elmer/Ice in reverse communication mode.This method uses an approximation of the second derivatives of the cost function and is therefore more efficient than a fixed-step gradient descent.
How we define the inner product used to compute the Gâteaux derivatives affects the definition of the Fréchet derivatives, and could affect the convergence of the minimisation, but does not affect the minimum we are seeking to achieve.As for glaciological applications, velocities and strain rates can vary by several orders of magnitude inside the domain, and we have observed that including the nodal weights in the definition of the Fréchet derivatives leads to good convergence properties when using an unstructured mesh where large elements correspond to low-velocity areas, and vice versa.Possible alternatives are, for the control inverse method (Morlighem et al., 2010), to use a cost function that measures the logarithm of the misfit or, for the Robin inverse method (Arthern and Gudmundsson, 2010), to use a spatially varying step size rather than a fixed step in the gradient descent algorithm, as proposed in Schäfer et al. (2012).

Mesh and deforming geometry
Ice sheets and ice caps have a very small aspect ratio, horizontal dimensions being much larger than the vertical dimensions, and therefore meshing requires special care.The strategy commonly adopted in Elmer/Ice for meshing glaciers, ice sheets and ice caps is to mesh first the horizontal 2-D footprint and then extrude it vertically.These meshes are then vertically structured with the same number of layers over the whole domain, whereas the horizontal dimension can be meshed using an unstructured mesh.This is one of the main advantage of a FE ice flow model in comparison to the classically used finite difference or volume methods for which the grid has the same size over all the domain, unless a mesh adaptive method is implemented (Cornford et al., 2013a).
The unstructured mesh of the footprint can be created using triangle-shaped elements of various sizes to account for the spatial heterogeneity of the variables gradient.The horizontal size of the elements can be controlled using, for example, a metric constructed from the Hessian matrix of observed surface speed (Gillet-Chaulet et al., 2012).Technically, optimising the mesh sizes according to this metric is done using the freely available anisotropic mesh adaptation software YAMS (Frey and Alauzet, 2005).Because of the overall size of ice sheets, the mesh is then partitioned and all partitions are solved in parallel using the Message Passing Interface (MPI).In Elmer/Ice, the mesh can be generated either by extrusion as a preprocessing step or by a built-in mesh extrusion feature which operates on the parallel level.This internal procedure efficiently removes some of the possible bottlenecks in preprocessing as the maximum mesh size is no longer constrained by serial operations.Also, in the case of an extruded mesh, certain operations become trivial, as for example modifying the geometry or computing the depth or elevation, which efficiently becomes a one-dimensional problem.
For transient simulation, the geometry of the ice sheet is evolving with time and the mesh has to be deformed to follow these changes.The common approach to deform geometries in Elmer/Ice, if dealing with unstructured meshes, is to rearrange the nodes by solving a pseudo linear elasticity problem.Any mesh displacement, x, in Elmer/Ice is relative to the initial mesh position, x 0 , i.e. x(t) = x 0 + x(t).A deformation of the surface, for instance, can be induced by a changing free-surface elevation, h.Hence, the prescribed vertical deformation here is x(t) • e z = h(t) − h(t = 0).Inside the bulk mesh, the corresponding deformation is then obtained by solving where Y and κ are respectively a pseudo Young's modulus and Poisson ratio, describing the resistance against the deformation and its directional ratio.Further, describes the symmetric strain tensor In consequence, the induced mesh velocity from the recomputation of x by Eq. ( 48) from one discrete time-level t to t + t then is given by The continuum equations, as presented above, strictly, are valid only in a fixed reference frame.In ice sheets or glaciers, however, the geometry by nature is not fixed.In a fixed reference frame, if u is the fluid velocity, the total change of a scalar property, , consists of the local change and of a convective part u • grad .
For instance, if we solve Eq. ( 14), we should take any induced mesh velocity Eq. ( 50) into account.This is done by the arbitrary Lagrangian-Eulerian (ALE) formulation, which is based on the Reynolds transport theorem (e.g.Greve and Blatter, 2009).Conservation simply demands that in the general reference frame of a moving mesh, Eq. ( 51) changes into A slight deviation from this is for the kinematic boundary condition Eq. ( 10), as the convection term is only in the horizontal plane, i.e.
The same is applied to Eq. ( 11) for the evolution of z b .A special case of Eq. ( 52) is when the surface is considered to move horizontally (it does vertically by definition) at the speed of the fluid particles, i.e. u m = u s .This, for instance, is needed if dealing with advancing fronts in marine-terminated glaciers.In this case, the new position of the surface is determined by x(t + t) = x(t) + u(t) t.In terms of the absolute mesh update x, this means that x(t + t) = x(t)+ u(t) t, which simply reflects Eq. ( 50) under u m = u s and a s = 0.

Variational formulations
Elmer/Ice relies on the FE method, and all the equations presented above are solved using a discretised variational form, leading in turn to solving a linear system for which the unknowns are the nodal value of the variable.The aim of this section is to present some technical details for the variational forms of both the Stokes and transport equations.

Stokes equations
The discrete variational form of the Stokes system Eqs.( 1) and ( 2) is obtained by integration over the ice domain using the vector-valued weight function and the scalar weight function , In the relation given above, the left-hand-side term in the momentum equation given in Eq. ( 1) has been integrated by parts.One part is re-formulated by applying Green's theorem, transforming it from an integral over the domain into one over the closed boundary of the domain , for which Neumann or Newton boundary conditions can be set (e.g.vanishing deviatoric surface stress components).The numerical solution of Eq. ( 54) is obtained by either using the stabilised method from Franca and Frey (1992) or the residual free bubbles method in Baiocchi et al. (1993).
For non-linear rheology, e.g.n = 3 in Glen's law (Eq.4) and Eq. ( 54) is non-linear and needs to be solved iteratively.For the (k + 1)-st iterations of the non-linear loop, the effective viscosity η k+1 is estimated from Eq. ( 4) using the previously computed velocity field u k (fixed-point method or Picard method) or a Newton linearisation such as Convergence is obtained much faster with the Newton iteration than with the Picard method, but the former algorithm is known to diverge when starting too far from the solution.In practice, the solution is initially approached by performing some Picard iterations, and then activating the Newton linearisation for the last iterations.Because the convergence is problem-dependent and depends on the initial solution, there is no general rule stating when to start the Newton iteration scheme.The efficiency of the Newton method is illustrated on a test case in Sect.7.2.

Transport equations
13 ), equations for the age of ice (Eq.27) or the ice fabric evolution (Eq.15) can be expressed using the generic form where i = 1, K = 0 and F = 1 for the dating equation and where i = 1, . . ., 5, and K and F are vector functions of C, W, a (2) and a (4) for the fabric evolution equations.For incompressible fluids (ice), div(A i u) simplifies to u• grad (A i ).Equation ( 57) is a first-order hyperbolic equation, and is nonlinear when solving for the fabric evolution.
The variational formulation of the transport equations is obtained by multiplying Eq. ( 57) by the test function and integrating over the ice volume .Because in the case of a vector solution A the set of equations is solved iteratively for each component independently, the variational formulation is presented for a scalar A, and reads The second term is then integrated by parts, so that Dirichlet conditions have to be applied on all the boundaries of where the flow velocity is directed inside the ice domain.Because of the missing diffusion terms in Eq. ( 59), the classic Galerkin method is unstable.This transport equation is either solved using the discontinuous Galerkin method proposed by Brezzi et al. (2004) or a semi-Lagrangian method (Staniforth and Côté, 1991;Martín et al., 2009).

Preconditioned linear solvers
The discretisation and linearisation of the varying viscosity Stokes system lead to linear systems which cannot be solved efficiently by using standard linear solvers.A special preconditioned version of the generalised conjugate residual (GCR) method has therefore been implemented recently into Elmer/Ice to obtain effective parallel solutions of these systems.This new preconditioner utilises the natural block structure of the associated linear algebra problem and is derived from approximating the associated pressure Schur complement matrix S as S ≈ (1/η)M, where η denotes the viscosity corresponding to the current non-linear iterate and M is the mass matrix corresponding to the pressure approximation (for similar solvers for varying viscosity flows see Grinevich and Olshanskii, 2009;Burstedde et al., 2009;Geenen et al., 2009;ur Rehman et al., 2011).Results of scalability tests done with this block preconditioned solver are presented in Sect.7.3.Note that in conjunction with this Stokes solver version, we employ a bubble stabilisation strategy based on utilising bubble basis functions corresponding to the high-order version of the finite element method.

Normal consistency
All boundary conditions involving vector (velocities) or tensor values (stress) are in need of a consistent description of the surface normals.Nodal normals, by nature of the discretisation applied in the FE method, are not uniquely defined, especially with linear elements.Thus, the representation of surfaces has non-continuous derivatives at nodes.For ice sheet boundary conditions, this is a problem occurring typically at the bedrock interface in the presence of sliding where a Dirichlet-type non-penetration condition has to be applied, i.e. u•n = 0. Depending on the nodal definition of n, this condition can lead to artificial source and sinks for mass and momentum in very uneven parts of the bedrock.Gillet-Chaulet et al. (2012) have shown that using the average of the normal to the elements sharing the node to estimate the nodal normal can lead to an artificial mass loss of up to 10 % of total ice discharge at the margin of Greenland.Recently, the mass-conserving way of deducing the nodal surface (Walkley et al., 2004) has been implemented in Elmer.For a node x j , with an element-correlation number N j , the surface normal is derived by The relation above constructs the nodal surface normal, n j , as a sum of the normals evaluated at the adjacent elements, k , using the same weighting functions, ϕ k , as for the momentum equation.

Accounting for inequality
As presented in Sect.2, the ice temperature T in Eq. ( 14) is bounded by the pressure melting point T m , and the free surface z b and z s in Eqs. ( 10) and ( 11) must fulfil the inequalities z s (x, y, t) ≥ b(x, y)+h min and z b (x, y, t) ≥ b(x, y).The variational inequality is solved using a method of imposed Dirichlet condition that are released by a criterion based on the residual.Let be the matrix equation of the unconstrained system.In FE method terminology, A is the system matrix, h the solution vector and a the body force.We then have to solve Eq. ( 61) under a constraint vector h min .Here, we choose a minimum value, hence a lower constraint, but the same method works also with upper bounds, like it is applied in the case for constraining the temperature with the local pressure melting point.This lower constraint reads In order to enforce Eq. ( 62) we apply the same method as in Zwinger et al. (2007): -A node i violating h i > h min i is set as "active".
-For each "active" node a Dirichlet condition h i = h min i is introduced into Eq.( 61).This is achieved by setting the ith row of the system matrix to A ij = δ ij , where δ ij is the Kronecker symbol, and the ith entry of the force vector to a i = h min i .Doing so for all active nodes results in an altered, constrained system -Instead of Eq. ( 61) we are now solving Eq. ( 63), obtaining a solution vector h .
h in turn is inserted into the unconstrained system Eq.( 61), defining the residual -If an earlier "active" node is found to comply with Eq. ( 62), it is taken of the list if, and only if, R i < 0.
This algorithm is repeated as long as there is no change in the active node set and the convergence criteria imposed for the solver are met.
In a converged state, the physical meaning of the residual can be interpreted.For the heat equation Eq. ( 14), the residual represents the additional cooling needed to comply with the inequality T < T m .For the free-surface equations Eqs. ( 10) and ( 11), the residual can be interpreted as the pernode additionally needed accumulation/ablation to meet the constraint Eq. ( 12).

Convergence tests
Convergence of the Stokes solver is tested by running the same problem with an increased mesh resolution.The purpose of this exercise is to verify the model and compare the efficiency of the various elements and associated stabilisation methods available within Elmer.For three-dimensional geometries, the Stokes equation can be solved using 8-node (linear) or 20-node (quadratic) hexahedron elements, or 6node (linear) wedge elements.Stabilisation of the Stokes equations is done either using the stabilised method (Franca and Frey, 1992) or the residual free bubbles method (Baiocchi et al., 1993).
The Stokes solver is verified using the manufactured analytical solution first proposed in Sargent and Fastook (2010) and subsequently modified and corrected in Leng et al. (2013).Here, we use exactly the same geometry and set of parameters as in Leng et al. (2013).All meshes are structured and defined by the number of elements in x, y and z directions.The mesh discretisation is made to vary from a very coarse mesh (20 × 20 × 5, 10 584 degrees of freedom for the linear element) up to a fine mesh (160 × 160 × 40, 4 251 044 degrees of freedom).The modelled domain is the bumpy bed of ISMIP-HOM experiment A (Pattyn et al., 2008) with the 80 km length.For this geometry, the finer mesh does correspond to horizontal and mean vertical resolutions of 500 and 25 m, respectively.The convergence rate for the various elements and stabilisation methods is obtained as the slope of the L 2 relative error norm function of the grid size refinement.The L 2 relative error norm of two arbitrary vectors u and v is defined as For simplicity, it is plotted here as a function of the cubic root of the inverse of the degrees of freedom, which is proportional to the grid size refinement.These curves for the three components of the velocity and isotropic pressure are plotted in Fig. 1.They show that whatever the element type and stabilisation method, the rates of convergence (slopes of the curves) are similar and close to 3 for velocity and pressure.A rate of convergence of 3 is greater than the theoretical Results of the convergence tests: L 2 relative error norm between Elmer/Ice and analytical solutions for the 3 components of the velocity (u, v, w) and the pressure, p, as a function of the grid size (which is proportional to the inverse of the cubic root of the degrees of freedom) for Franca and Frey (1992) stabilisation with (black) 6-node wedge element, (red) 8-node hexahedron element and with (blue) 20-node hexahedron element; and for the residual free bubbles stabilisation (Baiocchi et al., 1993) with (green) 8-node hexahedron element and with (magenta) 20-node hexahedron element.The black dashed line indicates a rate of convergence of 3.
value expected (e.g.Ern and Guermond, 2004), especially for the linear element and pressure.Surprisingly, the same rate of convergence is obtained for linear and quadratic elements, and for a given discretisation the quality of the solution is even better using linear elements, so that the use of quadratic elements is not recommendable, at least in this particular example.For this application and the quadratic 20-node hexahedron element, the residual free bubbles method is found to be less accurate than the stabilisation method of Franca and Frey (1992).

Picard versus Newton linearisation
Picard and Newton schemes for the non-linear solution of the Stokes equations are compared by performing the ISMIP-HOM experiment A005 (Pattyn et al., 2008;Gagliardini and Zwinger, 2008) for two different initial conditions.The first one assumes null velocity and pressure, whereas the second initial condition is equal to the SIA solution for this problem.
The switch from the Picard to the Newton iterative scheme is controlled by a criterion on δ u p ,u p+1 , the L 2 relative error norm, Eq. ( 65), between the previous p and current p + 1 velocity fields of the non-linear iteration loop.The same diagnostic simulation is repeated for switch criteria of 10 −6 The colour of the dot-dashed lines indicates the value of the switch criterion of the corr 36 Fig. 2. Evolution of the L 2 relative error norm between two consecutive solutions of the Stokes system, δ u p ,u p+1 , as a function of the non-linear iteration, for a switch criterion from the Picard to the Newton scheme of 10 −6 (Picard only, black curve), red 10 −2 , green 10 −1 , blue 1 and magenta 2 (Newton only), for the ISMIP-HOM experiment A005 with initial conditions (circle) assuming zero velocity and pressure and (triangle) estimated from the SIA solution.The colour of the dot-dashed lines indicates the value of the switch criterion of the corresponding colour curve.
The evolution of δ u p ,u p+1 as a function of the non-linear iteration indices is presented in Fig 2 .As expected, the Newton scheme is quadratically convergent, while Picard converges only linearly (Paniconi and Putti, 1994).When the initial condition is null velocity and pressure, it takes 40 Picard iterations to converge, whereas with Newton's method alone, it requires only 10 iterations.Surprisingly, even if it takes less Picard iterations to converge for the SIA initial condition, the convergence of the Newton solver is only obtained if Picard iterations are performed from δ u p ,u p+1 < 10 −1 .This example shows that Newton's method can diverge if the initial condition is too far from the converged solution.A switch criterion δ u p ,u p+1 < 10 −2 is found to work in most cases and it reduces the non-linear iterations by a factor about 2. Because the CPU consumption is almost proportional to the number of non-linear iterations performed within one time step, switching from Picard to Newton iterative schemes can reduce CPU time by the same factor.

Elmer/Ice scalability
The scalability of the new block preconditioned solver is tested and compared with the parallel sparse direct solver MUMPS (Amestoy et al., 1998).For this purpose, the  diagnostic Stokes solution is computed using the present-day Greenland geometry.The Greenland footprint is first meshed using regular triangle elements and then vertically extruded.Different meshes are constructed by varying the horizontal element size from 5 km to 10 km, but all have 20 vertical layers.The size of the tested meshes varies from 708 000 up to 4 580 000 nodes.Temperature and basal drag are imposed using the same fields as in Gillet-Chaulet et al. (2012).Newton iterations are used after the convergence criterion reaches 5 × 10 −2 .
The results obtained for strong scalability, i.e. a constant problem size with different partitionings, and for weak scalability, i.e. a constant load per CPU using different mesh sizes, are presented in Figs. 3 and 4, respectively.If the elapsed time is t n for a number of partitions n, then, for a strong scaling test, the scalability of a solver can be characterised either by its acceleration t ref /t n or its efficiency (n/n ref )t ref /t n , where ref stands for the reference simulation (often the one with the smallest mesh size).For a weak scaling test, the acceleration depicted in Fig. 4 is defined as The weak scalability experiment uses a constant number of 4200 nodes per partition in combination with an increasing number of partitions from 168 up to 1092.Weak scalability is found to be greater than 60 % even for the largest test case.Efficiency greater than 100 % is obtained with the new block preconditioned method for the strong scalability, whereas for a number of partitions greater than 100, MUMPS was always found to scale badly due to an increase of the required memory.This new solution strategy using the block preconditioned solver clearly opens the door to applications one order of magnitude larger in mesh size than what we were able to achieve so far using a direct solver like MUMPS.

Inverse methods validation
We test the two inverse methods previously presented in a 2-D example resembling a calving glacier.As our objective is to validate the numerical implementation, we use a linear rheology for which the two inverse methods implemented are exact.
Our domain is 20 km long, the bed elevation is constant and equal to −900 m.The free-surface elevation decreases linearly from 200 m at x = 0 km to 100 m at x = 20 km.The free surface is stress-free, we prescribe an homogeneous Dirichlet condition of 50 m a −1 for the horizontal velocity at x = 0 km, we apply a Neumann condition (hydrostatic sea pressure) at x = 20 km, and we apply a linear sliding law and a non-penetration condition at the bedrock.
We generate a reference solution with The surface velocities computed from this reference solution are then used as perfect synthetic observations (twin experiments).
The first step in the validation process is to assess the ability of each solver independently to reconstruct the synthetic observations.For the sliding coefficient β, we start  39 Fig. 5. Evolution as a function of the number of iterations of (top) the cost function relative to the initial cost function and (bottom) the norm of the gradient vector relative to the initial gradient norm for Robin (black curves) and control (red curves) inverse methods for the inversion of β (solid curve), E η (x, y) (dashed curve) and E η (x, y, z) (dotted curve).from the initial guess β i = 10 −3 .The viscosity η is expressed as η = E η η 0 , and the optimisation is done on the viscosity enhancement factor E η with initial guesses E η = 1 and η 0 = 15.It is possible a priori that several distributions of E η , especially in the vertical direction, can lead to the same surface velocities; therefore, we made it possible in the model to invert E η only in the horizontal plane (x, y) when using vertically extruded meshes.The gradient of J with respect to E η at a given position (x, y) is then obtained as the vertical sum of the nodal gradients at position (x, y).
To ensure that β and E η remain positive during the optimisation, β is expressed as β = 10 α 1 and E η is expressed as The optimisation is then done with respect to the α i (i = 1, 2).The evolution of the cost function and norm of the gradient obtained for each test is given in Fig. 5.Both the cost function and the gradient norm decrease and tend toward zero with the number of iterations.
The second step in the validation process is to verify the following approximation: For a given perturbation α i , the left-hand-side term is evaluated by computing J (α i + hα i ) with the direct model for several values of h, and the right-hand-side term is computed directly from the nodal gradients.This test is done  39 Fig. 6.Ratio (h) obtained with the Robin (black curves) and control (red curves) inverse methods for perturbations of the enhancement factor to the viscosity E η (x, y).
for the initial conditions of the previous twin experiments.We also test the implementation of the Tikhonov regularisation, Eq. ( 44), by choosing the cost function such as J = b 0.5(∂α i /∂x) 2 d where α i = 1.0 −9 (1.0+sin(2.0π× 4x/L)).
For each solver, we use 10 random perturbation fields α i where each nodal value of α i is a random number between −50 % and 50 % of the mean value of α i .The gradient computed from the two inverse methods is verified using the ratio (h) defined as An example of the evolution of (h) as a function of h is shown in Fig. 6 for the perturbation of the enhancement factor to the viscosity E η (x, y).For both inverse methods and all experiments, i.e. perturbation of β (not shown), E η (x, y, z) (not shown), E η (x, y) (Fig. 6) and the Tikhonov regularisation (not shown), the ratio (h) is found to decrease as h decreases and it reaches a value typically lower than 10 %.Such values are already satisfactory; nonetheless we could obtain even more accurate gradients by automatically deriving the code itself.
Figure 7 illustrates the difference of results obtained when inverting for E η only in the horizontal plane or in the whole ice volume.As can be seen, the two inferred fields E η (x, y) and E η (x, y, z) are significantly different even if surface velocity and observation are in comparable agreement.This indicates that a non-unique solution can be obtained when the number of control parameters starts to be larger than the number of observation points.

Outlook
A number of the requisites for an ice sheet model as discussed in the Introduction have already been implemented in Elmer/Ice, and especially those necessary to accurately describe the flow of polar ice.Nevertheless, as for other ice sheet models, the physical processes at the boundaries and their coupling with the other components of the climate system can still be improved in the near future.This is the prerequisite for running any forecast simulations of ice sheets, and not only sensitivity experiments based on more-or-less crude parameterisations that link changes in the atmosphere and the ocean to changes at the boundaries of the ice sheet.
Our efforts in the near future will be dedicated to improving the physical description of the ice/atmosphere, ice/ocean and ice/bedrock boundaries, as well as the models describing how the pertinent variables at these interfaces are distributed.
For the basal boundary condition, numerical modelling (Schoof, 2010) or direct measurements (Sole et al., 2011) seem to indicate a very complex relation, most certainly nonlinear, between changes in surface runoff and modulation of basal sliding.Two ingredients would then be required to fully account for the complexity of basal processes in relation with changes in surface runoff: (i) a proper basal friction law depending on the effective basal pressure (i.e.Eq. 22), and (ii) an associated hydrological model to describe the basal water pressure distribution.This hydrological model is currently under development and will be presented in de Fleurian et al. (2013).
Changes in the front position of marine-terminated glaciers seem to have a great influence on the upstream ice flow by modulating the buttressing force (Vieli and Nick, 2011).Determining the rate at which icebergs are calved for many different configurations remains an open question in glaciology.Submarine melting acting at the calving front of glaciers certainly increases calving rate by undercutting the ice (Rignot et al., 2010;O'Leary and Christoffersen, 2013).A general calving law, especially for 3-D configurations, still needs to be formulated (Benn et al., 2007).Better knowledge of the stress distribution at the front of glaciers and of the submarine melting distribution, as well as a reliable ice damage model (e.g.Pralong, 2005;Jouvet et al., 2011), are the required ingredients to describe calving at the front of marineterminated glaciers.In Elmer/Ice, the already implemented ALE formulation for the free surface accounts for moving ice sheet margin boundaries.Because Elmer/Ice solves the full-Stokes system, all components of the stress field are known and can therefore be used to evaluate ice damage.Work is in progress to implement an ice damage rheological law following Pralong (2005) with the aim of using damage iso-surface to locate the calving surface and move accordingly the front surface.
Melting from beneath the ice shelves is certainly one the most important triggers of the observed recent ice stream accelerations (e.g.Payne et al., 2004;Dupont and Alley, 2005).Not only is the total amount of basal melting important, but also its spatial distribution (Gagliardini et al., 2010).For numerical and technical reasons, coupling an ocean model with an ice sheet model is still a challenging issue.An intermediate approach we would like to explore as a preliminary step towards a complete coupling with an ocean model is the implementation within Elmer/Ice of a plume-type model (Holland and Feltham, 2006).

Conclusions
We have presented in detail the Elmer/Ice ice sheet flow model, from the equations implemented to the way they are solved using the FE method.Elmer/Ice contains a high mechanical description of ice flow: it solves the complete Stokes equations without any approximation, includes two complex anisotropic flow laws, resolves the grounding line dynamics as a contact problem and incorporates a basal friction law accounting for cavitation.Temperature and fabric fields within the ice sheet domain can be determined in a coupled manner with the flow solution.Other equations allowing for the derivation of secondary variables from the Stokes solution, such as the age of the ice, the stress or strain rate fields, are also implemented.Two recent inverse methods have been implemented in Elmer/Ice that make it possible to infer the poorly known parameters to construct the initial state of the ice sheet.From a technical point of view, Elmer/Ice reaps the benefits of the FE method, and provides an easy mesh adaptation method to focus on areas of interest.Elmer/Ice is a highly parallelised code and, as a recent important improvement, the new block preconditioned solver will in the near future lead to increase significantly the size of the solved problems.As listed in the previous section, there is still a need for future improvements and new developments, particularly by linking more tightly the pertinent variables controlling the flow at the boundaries, like the basal water pressure, with the other components of the climatic system.This next step is the requisite for driving ice sheet forecast simulations and furnishing reliable estimates of ice-sheet-induced sea level rise for the coming centuries.

Fig. 2 .
Fig.2.Evolution of the L 2 relative error norm between two consecutive solutions of th a function of the non-linear iteration, for a switch criterion from the Picard to the New only, black curve), red 10 −2 , green 10 −1 , blue 1 and magenta 2 (Newton only), for the IS with initial conditions (circle) assuming zero velocity and pressure and (triangle) estim

Fig. 4 .Fig. 4 .
Fig.4.Efficiency for a weak scalability experiment using the block precondi nodes in all meshes of 4200 and meshes from 0.708 × 10 6 nodes up to 4.58 × 1

Fig. 5 .
Fig. 5. Evolution as a function of the number of iterations of (top) the cost function relative to the initial cost function and (bottom) the norm of the gradient vector relative to the initial gradient norm for the Robin (black curves) and Control (red curves) inverse methods for the inversion of β (solid curve), Eη(x,y) (dashed curve) and Eη(x,y,z) (dotted curve).

Fig. 6 .
Fig.6.Ratio ∆(h) obtained with the Robin (black curves) and Control (red curves) inverse methods for perturbations of the enhancement factor to the viscosity Eη(x,y).

Fig. 5 .
Fig. 5. Evolution as a function of the number of iterations of (top) the cost function rela and (bottom) the norm of the gradient vector relative to the initial gradient norm for Control (red curves) inverse methods for the inversion of β (solid curve), Eη(x,y) ( (dotted curve).

Fig. 6 .
Fig. 6.Ratio ∆(h) obtained with the Robin (black curves) and Control (red curves) inv of the enhancement factor to the viscosity Eη(x,y).

Fig. 7 .Fig. 7 .
Fig. 7. Comparison of the obtained inverted enhancement factor field for the Control method assuming only changes in a horizontal plane (Eη(x,y), left) or changes in all 3 directions (Eη(x,y,z), right).The inverted field Eη(x,y) in the left is virtually identical to the one prescribed to obtain the observed velocities (Eq.(66)).

2013 1306 O. Gagliardini et al.: Elmer/Ice model
In addition, calculating of the stress from the velocity and isotropic pressure fields is a matter of interest because different methods can lead to noticeably different solutions.In Elmer/Ice, the components σ ij of the nodal Cauchy stress field are obtained from an existing Stokes solution (u, p) by writing the variational version of the constitutive law in a componentwise manner as www.geosci-model-dev.net/6/1299/2013/Geosci.Model Dev., 6, 1299-1318,