A comparison of two Stokes ice sheet models applied to the Marine Ice Sheet Model Intercomparison Project for plan view models (MISMIP3d)

We present a comparison of the numerics and simulation results for two “full" Stokes ice sheet models, FELIX-S (Leng et al., 2012) and Elmer/Ice (Gagliardini et al., 2013). The models are applied to the Marine Ice Sheet Model Intercomparison Project for planview models (MISMIP3D). For the diagnostic experiment (P75D) the two models give very similar results (<2% difference with respect to along-flow velocities) when using identical geometries and computational meshes, which we 5 interpret as an indication of inherent consistencies and similarities between the two models. For the Stnd, P75S, and P75R prognostic experiments, we find that FELIX-S (Elmer/Ice) grounding lines are relatively more retreated (advanced), results that are consistent with minor differences observed in the diagnostic experiment results and largely due to slightly different choices in the implementation of basal boundary conditions used by the two models. Based on current understanding, neither set of implementations can be argued to be more or less favorable. In this case, we propose that, at a particular resolution, 10 the span of grounding line positions from these two models provides one estimate for model error when treating the results from full Stokes models as a metric for accuracy in model intercomparison experiments. More importantly, we show that as grid resolution increases the grounding line positions for FELIX-S and Elmer/Ice appear to converge to within the estimated truncation error for Elmer/Ice.

The paper proceeds as follows. First, we give a brief overview of the governing Stokes-flow equations for ice flow, which are discretized and solved by Elmer/Ice and FELIX-S. We then discuss in some detail the implementation of boundary conditions -some specific to the problem of simulating marine ice sheets -and how they are implemented in the two models. A brief introduction to the MISMIP3d model setup is then given, followed by a presentation of experimental results for the two models.
We then give an in-depth discussion of the similarities and differences between results from the two models, our interpretation 5 of where these differences come from, and an assessment of their significance. We close with summary and concluding remarks. in expressed by the balance between the stress-tensor divergence and the gravitational body force, with σ representing the full stress tensor, ρ i the density of ice, and g the acceleration due to gravity.
The incompressibility of glacier ice is expressed as where u = (u, v, w) denotes the ice velocity vector. For glacier ice, the constitutive relation can be expressed as for a Newtonian 15 fluid, where τ is the deviatoric stress tensor, p is the isotropic ice pressure, I is the identity tensor,˙ is the strain rate tensor, and η is the "effective" viscosity, defined by Nye's generalization of Glen's flow law (Cuffey and Paterson, 2010), where n = (x,ŷ,ẑ) is the surface normal vector in a Cartesian reference frame. The lower-ice surface consists of two different boundary conditions (Durand et al., 2009a). For the "grounded" part of the flow where ice is in contact with bedrock (i.e., the ice-bedrock interface), the normal stress exerted by the ice body is larger than ocean water pressure. Here we apply the non-linear friction sliding law prescribed for the MISMIP3d experiments , 5 where C is a friction coefficient that is non-zero for grounded ice only, t is the bedrock tangent vector, and m is a friction-law exponent. For the floating part of the flow, where ice is detached from the bedrock (i.e., the ice-ocean interface, or the icebedrock interface at minimal floatation), the normal stress exerted by the ice body is smaller than or equivalent to the ocean water pressure, we apply a stress balance condition; normal stress, σ nn , is balanced by the pressure due to buoyancy, P w , 10 where z w is the sea level. For the case of z(x, y) = b(x, y), we also need to consider a "contact problem" (Durand et al., 2009a) to decide the actual location of the GL. We discuss the contact problem and its implementation in Elmer/Ice and FELIX in more detail below.
The evolution of both the upper and lower free surfaces are determined by a kinematic boundary condition, Both FELIX-S and Elmer/Ice discretize the Stokes-flow momentum balance equations using the Finite Element Method (FEM).
Both models have undergone extensive formal verification (see Gagliardini et al., 2013;Leng et al., 2013), have been subject to formal convergence studies (see Gagliardini et al., 2013;Leng et al., 2012Leng et al., , 2013, and have been shown to compare very favorably to one another when applied to the Ice Sheet Model Intercomparison for Higher-Order Models (ISMIP-HOM) (Pattyn 5 et al., 2008) experiments (see, Figures 6, 7, 10, 11, 13, and14 in Leng et al. (2012)). Additional details (and references) for Elmer/Ice are given in Gagliardini and Zwinger (2008) and Gagliardini et al. (2013), and for FELIX-S in Leng et al. (2012Leng et al. ( , 2013Leng et al. ( , 2014. Here, we provide a summary of several important similarities and differences between the numerical implementations used by Elmer/Ice and FELIX-S, noting that we view the differences as arbitrary. That is, there are not clear arguments for why one choice is superior to another and, in that sense, we view both methods as equally valid.

10
The first significant difference between Elmer/Ice and FELIX-S is in the choice of finite elements; Elmer/Ice uses hexahedral elements with P1-P1 basis functions (linear in velocity and pressure) and "bubble" function stabilization, whereas FELIX-S uses tetrahedral, Taylor-Hood elements with P2-P1 basis functions (quadratic in velocity, linear in pressure). The second important difference is that Elmer/Ice and FELIX-S use slightly different "masking" schemes for identifying grounded versus floating regions of the lower surface; Elmer/Ice marks the nodes bounding each element whereas FELIX-S marks the element 15 faces. The third important difference, which is a generic FEM implementation issue and not specific to the Stokes-flow problem, is in how the value of the basal friction coefficient, C, is applied at the Gaussian quadrature points. FELIX-S calculates the values of C at quadrature points directly (with the accuracy of integration increasing with the number of integration points), whereas Elmer/Ice interpolates the values of C at Gaussian quadrature points from nodal values (Gagliardini et al., 2016) (with the number of integration points needed for a given degree of accuracy determined by the order of the basis function).

20
In terms of similarities, Elmer/Ice and FELIX-S use the same scheme for evolving the free surfaces, based on an FEM discretization of the kinematic boundary equation (8) . The two models also use nearly identical implementations of the contact problem. For FELIX-S, the ocean water buoyancy pressure is compared to the normal stress of the ice on the bed and for Elmer/Ice, the ocean water buoyancy pressure is first integrated and then compared to the normal force of the ice on the bed (Durand et al., 2009a). While both models solve the contact problem at nodes, the information is used 25 slightly differently; Elmer/Ice uses it to decide if nodes in contact with the bed are floating or grounded whereas FELIX-S uses it to decide if nodes in contact with the bed constitute an element face that is floating or grounded. We return to the discussion of these different schemes and their impact on model output in Section 6.
Lastly, of the three potentially different ways for defining how the basal friction coefficient C varies over the area of a grounded-to-floating element -"Last Grounded" (LG), "Discontinuous" (DI), and "First Floating" (FF) (discussed in more simulation. We also return to this discussion in more detail in Section 6.

Experimental Setup
We provide a brief review of the MISMIP3d experimental setup, referring the reader to  for additional details. Three experiments are conducted and reported on; the "standard" prognostic experiment (Stnd), the prognostic, basal 5 sliding perturbation experiments (P75S and P75R), and the diagnostic experiment (P75D). The Stnd experiment is similar to that conducted in the original, two-dimensional MISMIP experiment for flowline models (Pattyn et al., 2012), where steadystate ice sheet GL positions are examined for a uniform, downward sloping (non-retrograde) bed in the along-flow (x) direction, with a uniform basal friction coefficient and uniform bed properties in the across-flow (y) direction. The goal is to compare three-dimensional model results to those from the two-dimensional test case, for which analytic solutions are available (Schoof,10 2007a). The prognostic P75S experiment starts from the steady-state geometry of the Stnd experiment and introduces a twodimensional, Gaussian perturbation (a slippery patch) to the basal sliding coefficient field, C(x,y), which introduces changes to the model state (velocity and geometry fields). The ice sheet geometry and GL are then allowed to advance for 100 years. The P75R experiment, which starts from the final state of the ice sheet at the end of the P75S experiment, returns the C(x,y) field to its original, uniform distribution, inducing GL retreat. The model is then integrated forward in time for another 100 years.

15
The P75D experiment compares the diagnostic model state when using the P75S geometry calculated by the Elmer/Ice model.
Below, we first report on the comparison between Elmer/Ice and FELIX-S for the P75D experiment. We then follow with a comparison for the Stnd, P75S and P75R experiments. is 50 m and across-flow resolution is varied from 2500 to 625 m. In this case the Elmer/Ice and FELIX-S nodal coordinates are very similar but not identical, as discussed further below in Section 5.3 (Note that we distinguish identical nodal coordinates as distinct from identical meshes, because the mesh can also be considered a function of element type, which are different for Elmer/Ice and FELIX-S).

The diagnostic experiment, P75D
We first compare the two models for the diagnostic experiment, P75D ( Figure 2). Both models use the same parameters (e.g., A, C, and m; see also Table 1) and, despite the different element types discussed above, have identical nodal coordinates over the entire model domain. From Figure 2, it is clear that the three velocity components (u, v and w) for Elmer/Ice and FELIX-

30
S are in very good agreement for both the upper and lower surfaces, an indication of inherent consistencies between the two models. For this experiment, the most direct comparison between Elmer/Ice and FELIX-S is afforded by the DI results as, prior to determining C, we directly interpolate the nodal basal boundary condition mask from the Elmer/Ice diagnostic solution onto the element-face mask used by FELIX-S. In general, for the horizontal velocity, u, the differences are relatively small over the entire model domain (<2%), relatively less near the ice divide and increase continuously from the GL to the ice shelf portion of the domain (Figure 3). For the v and w velocity components, we observe relatively larger discrepancies in the region of the 5 GL (around km 535 -555), but still very small differences (<5%) over the majority of the domain (Figure 3).
Despite efforts to make mesh, initial and boundary conditions, and parameter settings identical between the two models, several non-negligible differences discussed above are likely responsible for the small differences in velocities shown in Figure   3. The first likely cause for the small differences is the different boundary masking schemes; as noted above, FELIX-S marks the basal boundary faces in an element-wise manner versus the node-wise manner used by Elmer/Ice. To apply as similar as ice, as discussed further below). This may lead to small differences when assembling the element stiffness matrices and the 15 right hand side vectors (for the Dirichlet boundary conditions) as part of the FEM discretization of the Stokes system. Another likely cause for the minor velocity differences in the P75D experiment is the specification of the sliding coefficient C at Gauss quadrature points, as discussed above in Section 3. Finally, despite identical mesh coordinates, Elmer/Ice and FELIX-S use different element types, basis functions, and interpolation schemes as discussed above.
Overall, for the P75D experiment FELIX-S results in slightly larger horizontal velocities (u) at the GL than does Elmer/Ice. 20 As a result, FELIX-S exhibits a slightly (1%) larger ice flux through the GL than does Elmer/Ice. This systematic difference between the two models is likely a combination of the slightly different numerical choices discussed above. Again, as these choices appear arbitrary with respect to our current level of understanding, it is not clear that the implementation and results from one model can be distinguished as being superior to the other.

Stnd prognostic experiment 25
The comparison of Elmer/Ice and FELIX-S diagnostic experiment results demonstrate that model solutions are very close to one another when using identical nodal-mesh coordinates, but that slightly different numerics and/or implementations of boundary conditions result in non-zero differences in the model solutions. In turn, the prognostic experiments demonstrate how those biases accumulate and affect the time-integrated model solutions.
For the Stnd prognostic experiment, FELIX-S uses the same initial ice sheet geometry (based on the boundary-layer theory 30 solution of Schoof (2007b)), the same along-and across-flow resolution (50 m in the vicinity of the GL and 2500 m, respectively), and the same nodal mesh coordinates as Elmer/Ice. From this initial condition, the forward model is integrated for ∼1300 years, by which time the GL position is in equilibrium (based on the criteria that the relative rate of volume change is <10 −5 , the same criteria used by Elmer/Ice ). Both models demonstrate a continuous advance of the GL, with FELIX-S reaching a steady state GL position (x g ) of 519.85 km (Table 2) and Elmer/Ice reaching steady state positions of x g = 529.55, 526.80 and 522.35 km, for LG, DI and FF, respectively (Gagliardini et al., 2016). Apparently, FELIX-S produces a slightly smaller equlibrium-sized ice sheet with a GL position that is slightly more retreated than that of Elmer/Ice.
We attribute the different equilibrium GL locations to differences in the numerical schemes already discussed above. While the overall retreated grounding line of FELIX-S relative to Elmer/Ice is consistent with the minor velocity differences observed 5 -FELIX-S produces slightly higher along-flow velocities (and hence flux) upstream from, at, and downstream from the GL, with the time-integrated result of slightly thinner ice (and hence floatation) occurring slightly farther inland relative to Elmer/Ice -we note that Elmer/Ice velocities when using the FF scheme are significantly faster than for FELIX-S (Figures 1 and 2), and yet the Elmer/Ice GL when using the FF scheme is still advanced relative to that of FELIX-S (Table 2). Hence, other differences in the two numerical schemes must be more important in contributing to the observed steady-state GL location differences (we 10 return to the discussion of these differences in greater detail in Section 6). Regardless of the reasons, we note that the differences between the equilibrium positions for the FELIX-S and Elmer/Ice GL locations, for both DI and FF, are very close to or within the range of the estimated truncation error for Elmer/Ice at an along-flow resolution of 50 m (see Durand et al. (2009a), Figure   6, and Gagliardini et al. (2016), Figure 1c, and related discussions therein).
We repeat the Stnd prognostic experiment with FELIX-S but starting from an initially over-sized configuration, allowing the 15 ice sheet to shrink over time and the GL to retreat to its equilibrium position (as opposed to starting from an under-sized initial configuration with an advancing GL). In this case, an equilibrium GL position is reached after a forward model integration time of ∼1500 years and we find x g = 524.50 km, approximately a 5 km difference relative to the case with an advancing GL.
This difference in equilibrium GL positions under advanced versus retreated initial configurations is consistent with that found by Durand et al. (2009a) and Gagliardini et al. (2016) and is consistent with a model truncation error of ∼5 km at an along- Lastly, we conduct a convergence study for the Stnd experiment by comparing solution error against mesh resolution. In 25 order to control computational costs, the mesh is modified slightly relative to that discussed above. The number of across flow elements is unchanged but the number of vertical layers is reduced from 10 to 5. Also, we do not double the along-flow mesh resolution everywhere in the domain at each step in the convergence study but rather over a particular region within the vicinity of the grounding line (a similar approach has been taken in previous work, e.g. Durand et al. (2009a)). Simulations are conducted with along-flow resolution of 1600, 800, 400, 200, 100, and 50 m in this refined region. For the highest along-flow 30 resolution, which coincides with that of the Stnd experiment discussed above, the equilibrium GL position is 519.55 km (a difference of 0.30 km relative to when using 10 vertical layers). Figure 4 shows the Richardson estimate for the solution error versus the along-flow mesh resolution. Slight irregularities in the GL position as a function of increasing resolution result from doubling the mesh resolution in the along-flow direction and in the region of the GL, rather than over the entire mesh.
Regardless of these minor irregularities, the GL position is clearly convergent as a function of resolution with a convergence rate between linear and quadratic.

P75S and P75R prognostic experiments
In the P75S and P75R prognostic experiments, we investigate advance and retreat of the GL following a step-change perturbation in the basal friction distribution, for 100 years, and a return to the initial basal friction distribution, for a further 100 years 5 (the P75S and P75R experiments, respectively), as discussed above in Section 4. The initial condition for the P75S experiment is the steady-state GL position of the Stnd prognostic experiment discussed above. To manage computational costs, especially in experiments where sensitivity to mesh resolution is explored, both models employ regional refinement near the GL. Initial mesh resolution in this region is 50 m along flow and 2500 m across flow for both models but, because of the different equilibrium GL positions for the Stnd experiment, the area of refined mesh in FELIX-S and Elmer/Ice is located in slightly different decreases from 2500 m to 625 m), the "reversibility" -i.e. the return to the initial position -of the GL improves (note that we expect >>100 yrs to demonstrate full reversibility (Gagliardini et al., 2016)). More importantly, we also find that as the number of elements in the y direction increases from 20 to 80, the agreement between FELIX-S and Elmer/Ice increases for all of Elmer/Ice GL implementations (i.e., LG, DI, and FF; Figures 5-7). For the highest across-flow grid resolution, differences in the FELIX-S and Elmer/Ice DI and FF grounding line position changes are close to or below the published truncation error 20 for Elmer/Ice, and differences relative to Elmer/Ice LG are converging to that same value (Figure 8).

Discussion
As noted above, some fraction of the differences in the prognostic model simulation results can likely be attributed to the small differences in the model velocity fields, as seen in the P75D experiment. In turn, these differences are likely related to the different type of finite elements and basis functions used by the two models. However, we attribute the bulk of the prognostic 25 model simulation differences to differences in the treatment of the contact problem, and more importantly, to the different masking schemes used for the basal boundary conditions.
There are small differences in the way the contact problem is implemented in FELIX-S versus Elmer/Ice; while following the same physical basis for the contact problem, FELIX-S compares the normal stress and the sea water pressure acting at nodes, whereas Elmer/Ice compares the normal and sea water force acting at nodes (Durand et al., 2009a). The result may 30 be that, effectively, Elmer/Ice and FELIX-S "feel" slightly different normal forces (or pressures) at basal nodes of the ice-bed interface, resulting in slight differences when assessing whether a node (Elmer/Ice) or element (FELIX-S) is grounded or not.
Unfortunately, the different element types used by FELIX-S and Elmer/Ice do not allow for a definitive confirmation of this hypothesis.
Of greater importance, however, are the different treatments of the basal boundary condition masking schemes discussed in Section 3. Figure 1 provides a schematic summary of the differences in the Elmer/Ice and FELIX-S basal boundary masking schemes and demonstrates how those differences would impact the GL location in the two models for a particular "edge" case.

5
In the upper part of Figure 1, the nodes marked A and C are unambiguously floating (i.e., z(x, y, t) > b(x, y), so that no contact problem needs to be considered for those nodes). Because FELIX-S considers any element with one or more floating nodes to be floating, elements 3, 4, and 8 are all marked as floating, with the resulting FELIX-S GL position shown by the blue line in Figure 1. For the same geometric configuration, the node-based scheme used by Elmer/Ice defines a slightly different position for the GL, shown by the red line in Figure 1.

10
In addition to the slightly different grounding line locations, the different basal boundary masking schemes will lead to different profiles for C, as shown schematically in the lower part of Figure 1 where we plot approximate nodal C profiles for the two models. These differences come about because, for FELIX-S, the nodal matrix coefficients contain the contributions of C (and other variables) from the surrounding elements. As an example, consider profiles 1 and 3 in Figure 1. The C = 0 contributions from elements 3 and 8, and assuming additional floating elements to the north of element 3 and to the south 15 of element 8, reduce the matrix coefficients associated with the node along the centers of profiles 1 and 3 by a factor of approximately 2/6=1/3 (for these nodes, 2 of the 6 surrounding elements are floating), relative to Elmer/Ice. This estimate is only approximate because, in reality, the nodal coefficients contain additional terms related to the ice velocity and the basis functions, which are not uniform for all elements surrounding a node. Similarly, assuming that additional elements downstream of the FELIX-S GL are also floating, the coefficient at node B along profile 2 will be reduced by ∼5/6 relative to the Elmer/Ice 20 value, since 5 of the 6 surrounding elements are floating.
We attribute the majority of the differences observed in prognostic model simulations to these slight differences in GL position, and more importantly to these slight differences in the value of C. If we again consider profile 2 in Figure 1 (and to a lesser extent profiles 1 and 3), the relatively larger value of C for Elmer/Ice will lead to relatively less basal sliding there and, eventually, relatively thicker ice. This in turn will make it more likely that neighboring nodes may also eventually ground.

25
The overall, time-integrated result will be that, all other things being equal, the Elmer/Ice masking scheme will slightly favor grounding and/or grounding line advance relative to the FELIX-S scheme. This proposed difference in model behavior is consistent with the differences observed when the two models are applied to the prognostic experiments.
We further note that the differences between the nodal C profiles for Elmer/Ice and FELIX-S shown in Figure 1 are broadly similar to the differences between the DI and FF implementations in Elmer/Ice; despite the DI-like implementation of C in 30 FELIX-S, the different masking scheme results in C values at nodes that effectively "look" more similar to the FF implementation of Elmer/Ice (dashed C profile lines in Figure 1). Simulations using Elmer/Ice with these two different implementation demonstrate differences that are broadly similar to the FELIX-S and Elmer/Ice differences observed here for prognostic simulations; the FELIX-S equilibrium grounding line for the Stnd experiment is closest to that for Elmer/Ice when using the FF implementation (Table 2), the change in FELIX-S GL in the P75S experiment is closest to that observed for Elmer/Ice when 35 using the FF implementation (Table 2), and, at all across-flow resolutions, the advance and retreat curves for FELIX-S in the P75S and P75R experiments most closely resemble those for Elmer/Ice when using FF (Figures 5-7).
Based on our understanding of these model-to-model differences and their hypothesized impact on model simulations, we have a strong expectation that the differences in model outputs will decrease as model resolution increases. As the element size decreases, the gradients in ice sheet geometry that are the primary cause for differences in the nodal-versus element-5 based masking schemes will also decrease, and the two sets of model results should converge. Indeed this is exactly what we see for the P75S and P75R experiments. Similar to the observation of Gagliardini et al. (2016) that the LG, DI, and FF implementations in Elmer/Ice all converge to a similar solution with increasing resolution, we demonstrate here that the FELIX-S results also appear to converge to that same solution with increasing grid resolution (Figures 5-7). When considering the two most comparable implementations of the basal boundary condition masking schemes (FELIX-S and Elmer/Ice FF), the 10 two models agree for the P75S and P75R experiments to within the estimated truncation error for Elmer/Ice at all across-flow grid resolutions explored here (Figure 8). For the Elmer/Ice DI and LG implementations, the differences with FELIX-S as a function of grid resolution are also clearly converging ( Figure 8).

Conclusions
We have conducted a first, detailed comparison of two full Stokes ice sheet models, FELIX-S and Elmer/Ice, applied to the 15 MISMIP3d benchmark experiments. While previous informal comparisons have suggested very close agreement between the two models (Leng et al., 2012), here we explore the model similarities and differences much more carefully, focussing on how differences in model numerics lead to differences in model outputs when using identical meshes and forcing, and in particular on differences important for the simulation of marine ice sheet dynamics.
Overall, we find very close agreement between the two model outputs for cases where the impact of rather arbitrary choices 20 in the implementation of basal boundary conditions can be minimized; for the P75D experiment, diagnostic solutions (e.g., velocity fields) agree to within ∼2-5%. While it is difficult to attribute those small differences to particular numerical choices made by the two models, it is likely that different element types and basis functions and slightly different implementations of the contact problem for floating ice play a role. More significant differences between the two sets of model results are found for prognostic problems. Overall, we find that equilibrium grounding lines for FELIX-S are relatively more retreated than those 25 for Elmer/Ice (as demonstrated by the Stnd experiment) and that FELIX-S is slightly less inclined to ground, and hence less inclined to show grounding line advance then Elmer/Ice (as demonstrated by the Stnd and P75S and P75R experiments).
A detailed look at the two models strongly argues that differences in the basal boundary masking schemes and in the implementation of the basal friction coefficient are the source of these differences. As we are currently unable to judge whether or not one scheme is superior to the other, our results urge caution when interpreting the results from full Stokes models as does not prove) that, in the limit of high grid resolution, multiple full Stokes models can be shown to agree on a particular test case solution, despite small differences in their numerics.
Future efforts could improve on the work presented here by confirming the truncation error for the FELIX-S model, in order to understand if different numerics might be a means for further reducing model truncation error. Also, by running simulations at even finer grid resolutions, future efforts could definitively confirm that the results from multiple Stokes models converge in 5 the limit of very fine grid resolution.   the steady state GL position in the Stnd experiment. The rows for ∆xGLc and ∆xGLm denote the differences between xG 0 and the GL position at year 100 in the P75S experiment, at the centerline and margin, respectively. As it is invariant in the across-flow direction, we do not explore the sensitivity of the Stnd experiment to across-flow resolution. All GL positions and differences are given in km.