www.the-cryosphere.net/2/67/2008/ © Author(s) 2008. This work is distributed under the Creative Commons Attribution 3.0 License. The Cryosphere

Abstract. The aim of this paper is to describe in detail how the benchmark tests ISMIP-HOM (Ice Sheet Model Intercomparison Project-Higher-Order ice-sheet Model) have been performed using the open source finite element (FE) code Elmer ( http://www.csc.fi/elmer ). The ISMIP-HOM setup consists of five diagnostic and one prognostic experiments, for both 2-D and 3-D geometries. For all the tests, the full-Stokes equations are solved. Some technical points concerning FE, such as mesh characteristics, stabilisation methods, numerical methods used to solve the linear system and parallel performance are discussed. For all these setups, the CPU time consumption in relation to the accuracy of the solution is analysed. Based on these findings, some general rules on optimising the computing time versus the accuracy of the results are deduced.


Introduction
Following the EISMINT (European Ice Sheet Modelling INiTiative, Payne et al., 2000) benchmark experiments, the Ice Sheet Model Intercomparison Project (ISMIP) aims to provide a comparison of the new generation of ice-sheet flow models.It is composed of three different tests.The comparison of the ice-dynamic response of Antarctic and Greenland ice sheet models applied to climatic warming (ISMIP-POLICE), coordinated by Philippe Huybrechts, is dedicated to global ice-sheet models.The simulation setup for Heinrich-type ice-sheet instabilities (ISMIP-HEINO), coordinated by Reinhard Calov and Ralf Greve, aims to compare the ability of 3-D ice-sheet flow models to simulate large scale surges induced by an activation and deactivation wave of the temperate basal area.Since the model domain represents a small synthetic circular ice-sheet, this Correspondence to: O. Gagliardini (gagliardini@lgge.obs.ujf-grenoble.fr)experiment is not limited to global ice-sheet flow models and can be performed by higher order flow models.The last experiment that compares Higher-Order Models (ISMIP-HOM, http://homepages.ulb.ac.be/ ∼ fpattyn/ismip/) is coordinated by Frank Pattyn and is addressed to all 2-D and 3-D ice flow models.
In the present paper we discuss in detail the results obtained by running the ISMIP-HOM experiments with the FE open source code Elmer.The complexity of the equations to be solved in the proposed tests is such that most of the standard FE software for fluid dynamics should be able to perform the experiments (at least the five diagnostic).Therefore, the main focus of such experiments lays more on the comparison of the results obtained for a various range of all the parameters that control the simulation, rather than in producing the results themselves.In this paper, we present a detailed study on the influence of the mesh parameters (e.g., number of nodes, type of elements, ratio of horizontal to vertical element densities), numerical methods used to solve the linear system (direct or iterative methods), and the influence of the tolerance criteria used to control the convergence of the non-linear iteration induced by the non-linear flow law of ice.In order to compare these results, the quality of the solution for given settings, relative to a reference solution, as well as the CPU time consumption have to be taken into account.Especially the latter, even if being machine dependent, should have been part of the model output for the ISMIP-HOM experiments, since the applicability of a model to a certain problem size is a fundamental information.All results presented in this paper have been obtained on a HP ProLiant Cluster of the Center for Scientific Computing (CSC, Finland).

Equations and numerical methods
We briefly summarise the equations to be solved and the numerical methods used to solve them by the finite element (FE) method.

Governing equations
For all these experiments, we are interested in solving the gravity-driven flow of ice over a rigid bedrock for various boundary conditions at the ice-bedrock boundary.The constituting law for the ice behaviour is given by a Norton-Hoff type law (Glen's law in Glaciology): where τ and D are the deviatoric stress and the strain-rate tensor, respectively.Both, linear and non-linear ice rheology are envisaged.For the latter, the effective viscosity η is strain-rate dependent and denotes as: where the second invariant of the strain-rate ˙ e is expressed by As the HOM experiments are defined to be isothermal, the fluidity parameter A, in Eq. ( 2), is a constant.With ice being a highly viscous and incompressible material, the mass conservation and the conservation of linear momentum reduce to the well-known set of Stokes equations: and where u is the velocity vector, ρ the ice density and g the gravity vector.The Cauchy stress tensor σ =τ + pI has been decomposed into its deviatoric part -replaced by the strainrate using Eq. ( 1) -and the isotropic pressure p.For both 2-D and 3-D experiments, the z-axis is aligned with the negative direction of the gravity vector and the x-axis points in the principal direction of the flow.
In the case of prognostic experiments, the surface elevation z=z s (x, y, t) is part of the solution.For this free surface, the following equation serves as the kinematic boundary condition: where u ⊥ =(u x , u y ), the operator ∇ ⊥ stands for the gradient evaluated in horizontal directions, i.e., ∇ ⊥ =(∂./∂x, ∂./∂y), and a is the accumulation-ablation function, considered as a vertical flux.For experiment F, vanishing accumulation is prescribed.

Boundary conditions
For each of these conservation equations, a various number of boundary conditions have to be applied, depending on the test case.Three different kinds of boundaries can be distinguished for all these tests.
2.2.1 Ice-bedrock interface z=z b (x, y, t) For the velocity field, either no sliding (u(x, y, z b )=0) -as for experiments A, B and E -or a friction law linking the sliding velocity to the basal shear stress -as for tests C and D -applies.Such friction law is of the form where n b is the unit normal vector pointing into the bedrock and t i the unit tangent vectors (i=1 in 2-D and i=1, 2 in 3-D).The sliding parameter A s is only a function of space (A s =1/β 2 in C, D and F).In contrary to the shallow ice approximation, the basal shear stress τ b is not equal to the basal driving stress ρg(z s −z b )∇ ⊥ z s , but rather part of the solution.

Upper surface z=z s (x, y, t)
For all these tests, the upper surface is a stress free surface, which implies that n s •(σ •n s )=p atm ≈ 0, where n s is the unit normal vector of the surface.Note that this condition does not imply that the isotropic pressure p vanishes at the surface.As a consequence of the variational formulation used in Elmer, this stress free condition is inherent (natural boundary condition), such that no explicit condition has to be applied for the Stokes problem.For the kinematic boundary condition at the free surface Eq. ( 6), the rate of accumulation, a, has to be prescribed.Within the FE formulation of Eq. ( 6), it acts as a source term.

Lateral boundaries of the model
For experiments A, B, C, D and F, periodicity of the vertically aligned in-and outflow boundaries is applied.In principle Elmer does support non-conformal periodical boundary mapping, nevertheless, in order to avoid inaccuracies introduced by interpolation, it is favourable to ensure conformal mapping, like we did in our applications.

Finite element formulation
We hereafter only discuss some technical aspects of interest in connection with the ISMIP-HOM experiments.For a detailed presentation of the FE implementation in Elmer, the reader can consult the Elmer documentation, which is to be found under http://www.csc.fi/elmer.

Stokes equation
Ice flow is governed by the Stokes problem for an incompressible fluid, corresponding to the solution of Eqs. ( 4) and (5).A numerical solution of this set requires stabilisation by either using the classical Taylor-Hood elements (quadratic interpolation of the velocity and linear interpolation of the pressure), the stabilised method (Franca and Frey, 1992) or the residual free bubbles method (Baiocchi et al., 1993).
The way how stress is computed from the velocity and how the isotropic pressure fields is obtained is of particular interest, as different methods can lead to sensitively different solutions.In this study, the deviatoric stress field is obtained by solving Eq. (1) using the following variational form with the scalar test functions : where D ij and η are computed from the nodal velocities using the derivative of the basis functions.

Free surface equation
The free surface elevation z s is discretized as where z i s is the (nodal) elevation at the i-th node inside the local element of the meshed ice sheet surface, and i stands for the nodal interpolation functions.
The discrete variational form of Eq. ( 6) is obtained by spatial integration over the ice sheet surface, using the test function .The free surface evolution is therefore solved on a two-dimensional space for a three-dimensional flow problem or on a one-dimensional space for a two-dimensional flow problem.At time step t + dt the following set of equations is then solved where z s (t) is the ice sheet surface at time step t, and u ⊥ =(u x , u y ) and u z are taken from the solution of the Stokes problem at time t.
Due to the hyperbolic nature of Eq. ( 10) the standard Galerkin method does not apply.Stabilisation is obtained by applying the stabilised method as presented by Donea and Huerta (2003) page 172, adding the element-wise defined term The stabilisation parameter includes the representative element diameter d, and the norm of the velocity ||u||.
At each time step, to avoid a distortion of the domain mesh caused by the moving free surface, the nodes of the domain mesh are re-distributed by solving a fictive elasticity problem introducing a mesh displacement vector d, where I the unity tensor.The fictive elastic properties of the mesh are given by the arbitrary Lame parameters µ and λ.The vertical component of the displacement at the free surface then is coupled to the free surface variable by the Dirichlet condition d • e z =z s (t).

Linear system
The discretisation and linearisation process of partial differential equations by the FE method leads to a linear system to be solved.There are two categories of methods to solve this linear system: direct and iterative methods.The first leads to an exact solution up to machine precision, but is not feasible on very large problems.The latter utilises a converging sequence of approximate solutions.
In our particular application two-dimensional problems were solved by the direct Unsymmetric MultiFrontal method (UMFPACK, Davis, 2004), whereas a BiConjugate Gradient Stabilised method (BiCGStab, Kelley, 1995) was applied to three-dimensional computations.After a renumbering in order to optimise the matrix bandwidth and a scaling of the resulting matrix, an incomplete LU (ILU) factorisation with the lowest fill-in order (ILU0) was used as pre-conditioner for the system.The convergence of the iterative method is obtained when the relative change of the residuals is lower than the criterion L .

Non-linear system
The non-Newtonian stress-strain relation introduces nonlinearities into the system.Linearisation is obtained by expressing the strain-rate in the effective viscosity (2) in terms of the velocity field taken from the previous iteration step, thus applying a fixed point iteration scheme.Convergence is checked upon the norm of the relative change of the field variables where |U n | stands for norm of the solution vector at the n-th non-linear iteration step.

Coupled solvers and time stepping scheme
Applying constant boundary conditions, the prognostic tests are integrated over time until a steady state is reached.For time-integration, an implicit scheme is applied.On a single time-level, the Stokes and free surface solvers are run iteratively until the convergence criterion C is reached.The same convergence estimate Eq. ( 14) as for the non-linear iteration is used but the reference value for the current solution U n stands now for the last converged value of the non-linear system before starting a new coupled iteration cycle.Using this implicit scheme, the time step can be increased with no decrease of the accuracy of the solution.If the time step is small enough, an explicit scheme is equivalent.
In order to guarantee good convergence, one should set

Mesh generation
Implied by the domain geometry, experiments A, B, C, D and F were performed using layered meshes consisting of N x horizontal and N z vertical layers in 2-D and N x and N y horizontal and N z vertical layers in 3-D.The mesh then consists of quadrilateral elements with 4 nodes in 2-D and hexahedral elements with 8 nodes in 3-D.
The geometry for experiment E was defined by a longitudinal profile of 5 km in length with a grid spacing of 100 m, based on the topography of the Haut Glacier d'Arolla from the year 1930.From the two sets of point coordinates describing the free surface as well as the bedrock, 50 sections (48 of quadrilateral and 2 of triangular topology) with linear connections between the supporting points at the surfaces were created (see Fig. 1).This brings along the advantage that variables for output required on these points do not have to be interpolated.
In order to compare the influence of longitudinal as well as vertical element density, four different computational meshes using layered structure and a fifth using unstructured meshing technique have been created using the pre-processor Gambit (see Table 1 and Fig. 2). 3 General settings using FEM 3.1 Mesh influence

Mesh density
Mesh density is one of the most important parameters that controls the quality of the solution.Higher resolution naturally results in higher quality, but the drawback is that CPU time consumption increases with mesh resolution.Due to the numerical methods used, the CPU time is found to be proportional to the number of degrees of freedom to the power 1.27 and 1.11 in 2-D and 3-D applications, respectively (see Fig. 3 and its caption).For the Stokes problem, the degrees of freedom are 3N in 2-D and 4N in 3-D, where N is the total number of nodes.In order to run one of the diagnostic experiments, the CPU cost is approximately 0.0024N 1.27 seconds for 2-D problems, whereas it is 0.061N 1.11 seconds in 3-D.This difference of regression between 2-D and 3-D cases certainly results from the different methods used to solve the linear system (direct in 2-D and iterative in 3-D).Note the discrepancy of the CPU time consumption for the different tests, indicating that the CPU time consumption depends on the geometry, i.e. domain length L, of the problem to be solved (and not only the plain amount of nodes).

Ratio of horizontal to vertical layers
In total six variations of experiments B and D were performed using a regular mesh consisting of N x horizontal and N z vertical layers.The ratio N x /N z was varied from 1 up to 10, while keeping an almost constant value for the overall amount of nodes, N .Since we constrained N x to be a multiple of 4 in order to have one node exactly at the minimum and maximum value of the sinusoidal functions (bedrock geometry for test B and sliding function for test D), N naturally is not exactly the same for all the meshes.The range of nodes was 10 278 ± 267.The comparison between the tests is done using the minimum and maximum values of the field variables computed at the surface and the bottom.In the whole range of the tested ratios N x /N z , the maximum relative difference for the velocity, stress and pressure fields is less than 1% for test B. For test D, this relative difference even is below 1% for all the field variables, except for the minimum value of shear stress at the bottom, which should vanish in the limit of A s →∞.Relative to the maximum value of the shear stress at the bottom (≈100 kPa for the 5 tests), this error is lower than 5% for all the meshes, but increases with a growing ratio N x /N z .
The aspect ratio H /L of the experiments, where H and L are respectively the height and the length of the domain, varies from 0.006 up to 0.20.As a common rule applied to the FE method, elements should not be too much elongated.For the six experiments B and D, one cannot con- clude from these results which ratio N x /N z between 1 to 10 is more suitable for these tests.Nevertheless, our results show that for these particular setups with large difference between horizontal and vertical gradients, elongated quadrilateral elements can be used.For L=160 km, this leads to an element aspect ration of 0.0125.Therefore, the same ratio N x /N z =240/120=2 was used to produce all the submitted results.Likewise, the 3-D meshes were also built with N x /N z =N y /N z =2 for experiments A and C.

Convergence criteria
The convergence of the solution is controlled by three criteria: -L for the convergence of the linear solution if an iterative method is used to solve the linear system (3-D experiments), -N L for the convergence of the non-linear iterations due to the non linear ice-rheology, -C for the convergence of the coupled problem for the prognostic experiments.
The two criteria L and N L can be increased in order to decrease the number of iterations and decrease the CPU time consumption.To quantify the influence of the convergence criteria, experiment A with domain length L=40 km has been performed for L =10 −7 up to L =10 −2 and O. Gagliardini and T. Zwinger: Elmer and the ISMIP-HOM experiments  is directly proportional to the two convergence criteria.Increasing the convergence criteria by one order of magnitude will decrease the accuracy of the solution by one order of magnitude.Most of the win in CPU is obtained by reduction of the needed number of non-linear iterations, as indicated in Fig. 4a.
In all experiments submitted L =10 −6 was set for the 3-D experiments (a direct method was used for the 2-D experiments) and NL =10 −5 was taken to be the non-linear convergence criterion for 2-D as well as 3-D non-linear problems.

Numerical stabilisation of the Stokes equation was achieved by different methods:
-Taylor-Hood elements using quadratic velocity basis function and linear pressure basis function, -linear elements stabilised using the residual free bubbles method (Baiocchi et al., 1993), -linear elements using the stabilisation method by Franca and Frey (1992).
For a given number of elements, the first two methods result in the same size of the linear system, whereas for the third, just half the size of the linear system is needed.Nevertheless, we established that in order to maintain the quality of the solution, the number of elements should be increased for the last method.The first two methods lead to the same results in the case of linear (or plane) boundaries (as for experiments A, B, C and D), but the number of mesh nodes for the bubbles is half compared to the Taylor-Hood method, such that tests A, B, C, D and F were performed using the bubble stabilisation method.Since the normal defined on the domain boundaries in experiment E shows discontinuities, linear elements stabilised with the bubbles method occurred to be less adapted to this problem than non-linear elements (see Sect. 4.3).

Parallel runs
Parallel runs were performed for experiment A with L=5 km.Elmer utilises domain decomposition to perform parallel runs.The mesh was partitioned such that the associated nodes of the periodic boundary conditions fall into the same partition.A good partition table should minimise the number of shared inter-partition nodes in order to decrease the information exchange between the processors and assign an equal amount of nodes to the partitions, such that load balance is guaranteed.It is found that the efficiency of the parallelisation is inverse proportional to the number of shared nodes, relative to the total number of nodes.A vertical partitioning was found not to converge, certainly because in that case the number of shared nodes was too large.As shown in Fig. 3, if using horizontal partitioning, the parallel efficiency is close to 1, i.e. the decrease of the CPU time consumption is approximately proportional to the number of partitions.

Specific characteristics of experiments
In the following, we want to highlight specific results obtained for each of the experiments.A complete description of all the ISMIP-HOM experiments can be found in Pattyn and Payne (2006) (http://homepages.ulb.ac.be/ ∼ fpattyn/ ismip/) and in the complementary material of Pattyn et al. (2008).3.
As shown in Fig. 2 of the complementary material, the shape of the horizontal surface velocity for L=5 km significantly differs from that for the five other cases.The surface velocity is larger over the bump, anti-correlated with the ice The Cryosphere, 2, [67][68][69][70][71][72][73][74][75][76]2008 www.the-cryosphere.net/2/67/2008/3.As for experiment B, the shape of the horizontal surface velocity for L=5 km differs from the five others.The same explanation as in experiment B still applies to this test.
The discontinuity at x=3L/4 observed on most of the p(z b ) curves is difficult to interpret (see Fig. 27 on the complementary material).It was found whatever the mesh discretisation was.A mesh with no node exactly at the position x=3L/4 revealed the same discontinuity.This discontinuity also appears for the 3-D experiments C (see Fig. 27 in the complementary material).More detailed study of the deviatoric stress fields over the whole domain indicates that A good method to verify the model results is to compare the component of the load caused by the gravity force parallel to the bed ρg1000L sin(0.1 • ) to the integral of the basal drag over the whole bedrock L 0 σ nt dx.Obviously, these two quantities should cancel each other.For all the six experiments, this balance was found to be fulfilled with a relative error lower than 0.6%.

Experiment E
Experiment E000 (Arolla flow line without slip section) was run on all five grids presented in Sect.2.4.The submitted experiments have been run using the settings summarised in Table 4. Results of experiments E are shown in Figs.28 and 29 of the complementary material.Comparison of the CPU time spent for read-in of data and solution of the flow problem are depicted in Fig. 6.
The corresponding results for the Cartesian components of the surface velocities as well as the Cartesian component τ xy (z b ) of the deviatoric stress tensor and the pressure difference p(z b ) at the bedrock are shown in Fig. 7.
The results show a quantitative equal behaviour of the surface velocity field between the five different meshes.Nevertheless, a higher frequency as well as amplitude of the oscillation in the stress and pressure curves obtained at the bedrock can be seen on the grids using more than just one mesh interval between the supporting points of the geometry.Comparison to a run with second order elements on the mesh E 20 20 2nd revealed that this is a numerical artifact caused by the use of first order elements.Figure 8 clearly shows, that the curve for p(z b ) for the second order elements is missing these oscillations.
The linear interpolation between the data points used to construct the mesh is certainly responsible for these oscillations.These findings lead to the conclusion that the solution is not only influenced by the numerical method used to solve the equations, but also the method used to generate a computational grid from a digital elevation model (DEM).In our case we stick to the instructions of the test coordinator to use linear connections along bedrock and free surface between the given DEM points.We agree that this linear interpolation was necessary for an objective comparison of all the models with different node numbers on the boundaries.For other applications, we would recommend to perform a non-linear interpolation between the given DEM points.In Zwinger et al. (2007), we chose Nonuniform Rational B-Splines (Les and Wayne, 1997) interpolation in order to get a smooth bedrock geometry, which enabled us to do the runs using linear elements.
As one can deduce from the horizontal position of the inserted markers from the result obtained with the mesh E 1 20 (where the basal mesh points coincide with the supporting points of the geometry), the peaks in the solution obtained with the linear element mesh E 20 20 occur exactly at these supporting points of the geometry.The most evident explanation hence is that the kinks in the bedrock topography occurring at mesh points (due to the interconnection using straight lines) are responsible for these oscillations.This phenomenon is omitted using the second order element functions, since the first derivative of the test functions does not show unsteadiness and the kinks are being smoothed out.
Experiment E001 (Arolla flow line with slip section) was run on three of the grids.The results are displayed in Fig. 9.Here the situation concerning the oscillations of the stress and pressure curves is even more pronounced, as singular points, where a sudden change from no-slip to slip takes place, are imposed on the bedrock boundary.Again, an additional run on the mesh E 20 20 2nd revealed that -apart from the peaks at the singular points -oscillations can be avoided by applying second order elements.The comparison between different meshes and element types is depicted in    5. Since the central area over the bump is the zone of interest, the horizontal element size is increased from the centre to the lateral boundaries with a factor 2. To ensure the mass conservation, periodic boundary conditions are applied be- tween the upstream and downstream boundaries.For the two lateral boundaries, a zero flux and the periodicity for the velocity within the plane of the boundary are prescribed.
The steady state solution is reached by running a time dependent simulation with constant boundary conditions using a time step size equivalent to one year.As expected, the surface elevation velocity is a decreasing function of time.As shown in Fig. 11b (thick lines), the CPU consumption needed for one time step is a decreasing function of time, indicating that the number of coupled iterations needed for a converged solution during one time step is also decreasing with time (from 6 to 1).The two tests were stopped when the norm of the relative difference between two time steps of all the variables (velocity and surface elevation) was lower than 10 −4 .As shown in Fig. 11b (thin lines), this corresponds to a maximal vertical velocity of the surface of few centimetres per year.It took a time-integration equivalent to 197 (exp.F without sliding) and 228 (exp.F with sliding) years to comply with these convergence criteria.
Since there is no accumulation on the surface, the initial volume V 0 =9.9685×10 12 m 3 should be conserved.As shown in Fig. 11a, for both experiments, the relative change of volume is of the order of the precision (<10 −6 ).

Conclusions
All ISMIP-HOM benchmark tests have been performed using the finite element code Elmer.Some technical points linked to FE, such as element types and densities, numerical solution methods and parallel performance have been discussed.reported together with the parameters of the model setup.
From our sensitivity experiments, we have established that horizontally elongated elements did not cause any drop in solution quality for experiments A to D. For experiment E, the oscillations arising at the supporting points of the geometry were eliminated using second-order element functions.The few tests performed using parallel computing, with a parallel efficiency close to one, are very promising for future 3-D applications.
As shown in Pattyn et al. (2008), results obtained with Elmer are very close to the ones obtained with the other full-Stokes model.Consequently, a relatively high level of confidence can be attributed to these results, even if no analytical solution is available to validate the accuracy of all these results.In Pattyn et al. (2008), it is also shown that the observed deviations between the models are the consequence of the model hypothesis (ranging from SIA to full-Stokes), and not caused by the choice of the numerical methods and their implementation.In other words, models built on the same hypothesis give similar results.For the full-Stokes models, this clearly indicates that the results are very close to the exact solution of the problem.
Output figures for all the tests, as well as the original output files in the format defined by the ISMIP-HOM conveners, are attached as complementary material to this paper.One purpose of this article is to provide fully documented results of all ISMIP-HOM runs obtained with Elmer.Such material can be used in the future in helping the developments of new higher-order models, by allowing an easy comparison with some existing well documented benchmark tests.

Fig. 1 .
Fig. 1.Three sections (horizontal length 100 m) of the lower end of Arolla glacier geometry flow-line (experiment E) model with mesh E 20 20 applied.

Fig. 3 .
Fig. 3. CPU time consumption (s) as a function of the total degrees of freedom, for 2-D and 3-D experiments A, B, C and D. Power regressions y=0.0006x 1.27 and y=0.013x 1.11 for 2-D and 3-D experiments are plotted using straight and dotted lines, respectively.Dashed and dash-dot lines are parallel to the 3-D regression line with a power law factor divided by 4 and 16, respectively.The stars associated with the number indicate results for parallel runs for experiment A with L=5 km and the corresponding number of processors.

Fig. 4 .
Fig. 4. (a) CPU time consumption (circle) and number of non-linear iterations N NL (cross) relative to the reference run ( L =10 −7 ) and (b) maximum of the relative difference with the reference run ( L =10 −7 ) on the velocities as a function of L = N L /10, for experiment A with length L=40 km.

4. 1
Experiments A and B Results of experiments A are shown in Figs. 1 to 7 of the complementary material and the settings for the submitted experiments are given in Table 2. Results of experiments B are shown in Figs. 8 to 12 of the complementary material, http://www.the-cryosphere.net/2/67/2008/tc-2-67-2008-supplement.zip, as well as the results for Experiment A in the plane y=L/4.Settings for the submitted experiments B are given in Table can be explained by mass conservation: variation of horizontal flux cannot be balanced anymore by vertical flux at the free surface since the vertical velocity would be too large for the given flow depth.Therefore, the horizontal flux is more or less constant inducing larger velocity for smaller depth, and vice-versa.This feature is not observed in 3-D certainly because parts of the ice can flow around the obstacle by its lateral sides.The maximum of the absolute value of the difference between the isotropic and hydrostatic pressure at the bed, p(z b )=p(z b )−ρ g H , decreases with increasing L (see Fig.1of the complementary material).For L=160 km, the isotropic pressure at the base is very close to the hydrostatic pressure (| p(z b )|<4 kPa, to be compared to the mean basal isotropic pressure p≈9000 kPa).For the 3-D tests A, the output variables plotted in the plane y=L/4 are very close to that of the 2-D tests B (see Figs. 8 to 12 of the complementary material).This indicates that the third direction y, perpendicular to the mean flow direction, does not play a significant role for the flow.Even if the bedrock gradients in y and x directions are of the same order of magnitude, the 3-D flow conditions are not far from the plane-strain assumption made for the 2-D flow.The results get closer with increasing length of the domain.4.2 Experiments C and D Results of experiments C are shown in Figs. 13 to 21 of the complementary material file and the settings for the submitted experiments are given in Table 2. Results of experiments D and experiments C in the plane y=L/4 are shown in Figs.22 to 27 of the complementary material.Settings for the submitted experiments D are given in Table

Fig. 8 .
Fig. 8. Experiment E000: comparison of p(z b ) along the bedrock between linear and second-order element runs.

Fig. 10 .
Fig. 10.Experiment E001: comparison of p(z b ) along the bedrock between linear and second-order element runs.

Fig. 11 .
Fig. 11.Change of total ice volume relative to the initial volume (left) and maximum value of the surface change dz s / dt (thin lines, right) and CPU time consumption per one time step (thick lines, right) as a function of the time t for experiments F000 (solid lines) and F001 (dashed lines).

Table 1 .
Mesh parameters for experiment E.

Table 2 .
Settings for submitted experiments A and C.

Table 3 .
Settings for submitted experiments B and D.

Table 4 .
Settings for submitted experiment E.

Table 5 .
Settings for submitted experiment F.