Impact of mesh resolution for MISMIP and MISMIP3d experiments using Elmer/Ice

The dynamical contribution of marine ice sheets to sea level rise is largely controlled by grounding line (GL) dynamics. Two marine ice sheet model intercomparison exercises , namely MISMIP and MISMIP3d, have been proposed to the community to test and compare the ability of models to capture the GL dynamics. Both exercises are known to present a discontinuity of the friction at the GL, which is believed to increase the model sensitivity to mesh resolution. Here, using Elmer/Ice, the only Stokes model which completed both intercomparisons, the sensitivity to the mesh resolution is studied from an extended MISMIP experiment in which the friction continuously decreases over a transition distance and equals zero at the GL. Using this MISMIP-like setup, it is shown that the sensitivity to the mesh resolution is not improved for a vanishing friction at the GL. For the original MISMIP experiment, i.e. for a discontinuous friction at the GL, we further show that the results are moreover very sensitive to the way the friction is interpolated in the close vicinity of the GL. In the light of these new insights, and thanks to increased computing resources, new results for the MISMIP3d experiments obtained for higher resolutions than previously published are made available for future comparisons as the Supplement.


Introduction
Marine terminating glaciers in Antarctica and Greenland control the dynamical contribution of these ice sheets to sea level rise.Among the processes at play, the retreat of the grounding line (GL) has a major impact on this dynamical contribution.Accurate modelling of GL dynamics is therefore a precondition for prognostic simulations of the future of ice sheets in a warming climate (Durand and Pattyn, 2015).Previous works have emphasised the importance of the mesh resolution around the GL (Vieli and Payne, 2005;Durand et al., 2009a, b;Pattyn et al., 2012;Durand and Pattyn, 2015) and how the friction is interpolated in the vicinity of the GL (Gladstone et al., 2012;Seroussi et al., 2014;Leguy et al., 2014).Two recent intercomparison exercises were designed to compare and test the ability of ice-sheet models to resolve the advance and retreat of the GL based on different perturbations.MISMIP was dedicated to two-dimensional flow-line geometry (Pattyn et al., 2012) and used an analytical solution (Schoof, 2007), whereas MISMIP3d was a fully threedimensional setup (Pattyn et al., 2013).
Elmer/Ice was the only Stokes model to complete the MISMIP experiment 3a (Pattyn et al., 2012), and it was one of only two Stokes models to perform the whole MIS-MIP3d experiments (Pattyn et al., 2013).Moreover, in the latter intercomparison exercise, the diagnostic experiment P75D was directly built from the geometry obtained with Elmer/Ice after the 100-year perturbation experiment.As the only Stokes model to perform the two intercomparison exercises, Elmer/Ice results are currently used as references for comparison with other models based on lower-order Stokes equations (e.g.Feldmann et al., 2014).The results of the MISMIP and MISMIP3d intercomparisons obtained with Elmer/Ice are also used as benchmarks to test Stokes models during their development.
Both MISMIP and MISMIP3d intercomparisons have confirmed that, except the heuristic approach prescribing the O. Gagliardini et al.: Friction at GL in Elmer/Ice boundary layer flux at the grounding line (Schoof, 2007), all other approaches require a fine resolution close to the grounding line to accurately describe its dynamics.One common feature of both MISMIP and MISMIP3d is the use of a constant sliding parameter over all the grounded part.Doing so, the friction at the GL presents a discontinuity, which is believed to increase the model sensitivity to the mesh size at the GL.This raises the following questions.Is the sensitivity of models to mesh resolution specific to the discontinuous friction imposed in both MISMIP and MISMIP3d?Are there alternative numerical methods that would decrease the sensitivity to the mesh resolution for a given setup?
Two recent contributions started answering these questions: the first by adopting a smoothed friction upstream the GL (Leguy et al., 2014) and the second by introducing a subgrid evaluation of the GL position (Seroussi et al., 2014).From a modified MISMIP setup and using the shallow shelf approximation (SSA) implemented on a fixed grid, Leguy et al. (2014) have shown that introducing a smooth transition between finite basal friction in the ice sheet and zero basal friction in the ice shelf significantly improves the numerical accuracy of the model.In other words, the sensitivity of the GL dynamics to the grid size is shown to be significantly reduced when the friction continuously decreases to zero upstream the GL.Importantly, by smoothing the friction, the physical problem is modified and will result in a more retreated steady-state GL position than the original MISMIP one.However, a smooth friction vanishing at the GL is certainly more realistic than a discontinuous one since one expects that the effective pressure is null at the GL.Using the MISMIP3d experiments, Seroussi et al. (2014) compared various parameterisations of the GL position for a finite element (FE) SSA model.Using the SSA, the GL position is directly evaluated from the floatation criterion and can therefore be located at any point of the domain and not only at the element nodes.In this way, the basal friction can be evaluated with a subgrid resolution.Their results, for a discontinuous friction at the GL (MISMIP3d), showed that subelement parameterisation of the GL significantly reduces the sensitivity of the results to the mesh size at the GL.The proposed methods, by estimating the GL position at a subgrid scale, acts similarly to an increased mesh resolution around the GL, but without the numerical cost associated with remeshing when the GL is moving.
For a Stokes model, the solution proposed by Leguy et al. (2014) might be an alternative as, unfortunately, the subelement parameterisation implemented by Seroussi et al. (2014) in their SSA model cannot be applied to solve the contact problem between the ice and its bed.Indeed, the contact condition can only be evaluated at the element nodes.In other words, for a Stokes model, the two alternatives are either to solve a modified problem which would be less sensitive to mesh resolution or to improve the accuracy of the model by increasing the mesh resolution.Obviously, the former solu-tion cannot be applied if one wants to solve the original MIS-MIP and MISMIP3d experiments.
The aim of this brief communication is to study, for the Elmer/Ice Stokes model, the impacts on the accuracy of a smooth transition of the friction at the GL and of the way the friction is implemented at the GL.It is first shown that for the Stokes solution, contrary to what is found by Leguy et al. (2014) for SSA, introducing a smooth transition of the friction at the GL has no significant effect on the sensitivity of the model to the grid size.In the case of a discontinuous friction at the GL, we then present three possible FE implementations of the friction at the GL and show that these different implementations result in significant differences in terms of GL dynamics for the well-defined MISMIP and MISMIP3d experiments.All the newly obtained MISMIP3d results are made available in the Supplement for future model comparisons.

Sensitivity to mesh resolution and friction implementation
This section presents results on the sensitivity to the mesh resolution using a flow-line configuration.For that purpose, the GL dynamics is studied using a setup adapted from experiment 3a of the MISMIP intercomparison exercise (Pattyn et al., 2012).Experiment 3a assumes an overdeepened bedrock, a non-linear Weertman friction law and that the GL is evolved by step changes of the ice fluidity parameter.
Previous works have shown that steady-state position of GL could differ slightly depending on whether it is obtained from advancing or retreating GL but that this difference decreased with an increase in mesh resolution (Durand et al., 2009a).For a given mesh discretization, the accuracy of the model is therefore assessed as the difference between the retreat and advance steady positions.
Basal friction in the experiment 3a of MISMIP is imposed on the form of a non-linear Weertman sliding law, linking the basal shear stress and the sliding velocity: The original MISMIP 3a setup assumes a constant friction parameter C where the ice is grounded, i.e. for x ≤ x G , and perfect sliding at the interface between the ice and the ocean, i.e. for x > x G , x G being the GL position and assuming the horizontal velocity to be positive.
In order to smooth the friction upstream the GL, Leguy et al. (2014) have proposed a simple parameterisation of the effective pressure, the overburden pressure minus the water pressure, coupled with a Coulomb-type friction law.Here, following their idea, but assuming a simpler formulation, the friction parameter C of the original MISMIP experiment is modified as follows: Doing so, the friction linearly decreases over a distance L from C to 0 at the GL.Note that the physical problem is then modified and the steady solution for a given L > 0, as well as the transient phases, is expected to be different than those of the original MISMIP.When L = 0, the problem is equivalent to the original MISMIP and the friction presents a discontinuity at the GL.Because C * is estimated at the mesh nodes, and then interpolated on the element using the FE basis function, the same solution is expected for any L lower than or equal to the grid size.
The same type of mesh as the one used for producing the Elmer/Ice MISMIP results is used, with an evolving resolution along the flow direction (see Durand et al., 2009a, for more details).The discretization therefore refers to the minimum horizontal mesh size in the close vicinity of the GL.The model accuracy is studied for four mesh sizes, from 200 to 25 m, and L = 0, L = 60 and L = 500 m.Starting from the ice-sheet geometry for step 1 and step 5 of experiment 3a (see Pattyn et al., 2012, for more details), the ice fluidity for step 4 is then applied and the geometry is evolved until a steady state is obtained: one in advance (from step 1 to step 4) and one in retreat (from step 5 to step 4).
From Figs. 1 and 2a, one can clearly see that, for L > 0 (red and black curves), the problem is modified and so are the GL steady positions.The longer the length of the decreased friction (i.e. the larger L is), the less advanced the GL steady position.Simulations for L = 1000 m were even found to have their steady positions upstream the initial step 1 position and cannot be used therefore to test the model accuracy as both steady solutions are obtained in retreat mode.As shown in Fig. 2b, and contrary to what was found by Leguy et al. (2014), no improvement of the model accuracy is found when L is increased.For these simulations, the largest errors are even found for L = 500 m but with no significant differences from the other simulations.The reasons that might explain this different behaviour are multiple but most probably result from the two different flow approximations (SSA versus Stokes) and/or the adopted formulation to smooth the friction upstream the GL (form of the smoothing function and/or the typical length for the friction decay).However, these results seem to be in line with the ones obtained by Cornford with BISICLES (Cornford, personal communication; see the review material of this paper).
Moreover, in the case of a discontinuous friction at the GL (L = 0), three different numerical implementations of the friction in the close vicinity of the GL have been tested.The three implementations are presented in detail in the Supplement.The first is assuming that the GL defines the last grounded (LG) nodes and that friction is applied up to the nodes belonging to the GL.In the second, the nodes belonging to the GL are assumed to be the first floating (FF) nodes and are already freely slipping.The third one reproduces exactly the discontinuity (DI) of the friction at the nodes belonging to the GL.For the DI implementation, the friction at these nodes is only applied if integrating over an element where all other nodes are also in contact with the bedrock, but a free slip condition is applied if the node belongs to an element where at least one node is in contact with the ocean.The three implementations are illustrated in a twodimensional flow-line configuration in Fig. S1 of the Supplement.Note that as long as L > 0, all three implementations are equivalent and give the same results.Despite the DI implementation being certainly the most physical, up to now, all the published Elmer/Ice results were obtained using the LG method (Durand et al., 2009a(Durand et al., , b, 2011;;Gagliardini et al., 2010Gagliardini et al., , 2013;;Favier et al., 2012Favier et al., , 2014;;Drouet et al., 2013;Gudmundsson et al., 2012;Pattyn et al., 2012Pattyn et al., , 2013;;Krug et al., 2014).Note that other possible implementations, such as a constant friction value per element, would certainly yield other results.For L = 0, the three friction implementations (LG, DI and FF) converge to the same, most advanced, steady-state position when the mesh size is decreased.Nevertheless, as shown in Figs. 1 and 2a, for a given mesh size, differences on the steady GL positions from the three methods are of the same order as differences from advance to retreat for a given method.The LG method leads to the most advanced GL, the FF method to the least advanced GL and the DI method to an intermediate GL position.For a 200 m discretization, the difference between the LG and FF methods is larger than 15 km in both advance and retreat.The DI position is almost exactly half way between the LG and FF positions.With a 25 m resolution at the GL, these differences are reduced to less than 2 km in both advance and retreat.For the purpose of comparison, with a given method, the difference between advance and retreat is around ≈ 25 km at the resolution of 200 m and is decreased to less than 3 km at a resolution of 25 m.
Finally, Fig. 2a also shows the published Elmer/Ice GL position obtained in advance from step 3 to step 4 in Pattyn et al. (2012).This solution was produced using the same discretisation of 200 m at the GL, but not exactly the same mesh.Despite the same discretisation at the GL, there is a 3 km difference with the new LG solution for L = 0.In line with Durand et al. (2009b), these differences illustrate the sensitivity of the GL position not only to the mesh resolution at the GL, but also to the other mesh characteristics, and more specifically how strongly the mesh resolution is reduced downstream and upstream the GL.
As expected theoretically, the MISMIP flow-line study confirms that, despite a high jump in friction at the GL, all three implementations of the friction converge to an identical solution as the mesh resolution is improved, but they can lead to significantly different solutions for a too coarse mesh.
In the light of these significant differences between the three friction implementations for the MISMIP 3a experiment, the following section aims to quantify these differences for the MISMIP3d experiments.

Sensitivity to the lateral discretisation of MISMIP3d experiments
In this section, the three numerical implementations of the friction are compared using the prognostic experiments of MISMIP3d.New results for the diagnostic experiment P75D of MISMIP3d are also presented in Sect.S2 of the Supplement.The prognostic experiment in MISMIP3d is decomposed in three steps.First, assuming no lateral variation in y, a steady-state geometry is obtained for each model.In the second step, P75S, a Gaussian sliding perturbation is introduced precisely at the grounding line and centred on the axis of symmetry at y = 0 km.This constant perturbation is applied for the next 100 years.Finally, during the last step, P75R, the perturbation is removed and the GL moves back to its initial steady position.Only the first 100 years of the removal are studied.Note that for the grounding line to get back to its initial steady-state position might take much longer than 100 years as the behaviour in advance and retreat is not symmetrical.First, the steady GL positions for the three friction implementations are compared using meshes with the same resolution at the GL as the one used to obtain the MISMIP3d Elmer/Ice results (labeled LFA in (Pattyn et al., 2013)).
As expected from the previous section, the three methods result in three different GL positions x G 0 , the LG solution being more advanced by ≈ 7 km in comparison to the FF one (see Table S1 in the Supplement).It should be noticed that this distance is similar to the one obtained between the LG solution and the LFA solution published in Pattyn et al. (2013), using the same discretisation at the GL but not exactly the same mesh.This gives again an indication of how the results are sensitive to the mesh, and not only in the vicinity of the GL.It should also be noted that these differences stay much smaller than the differences obtained between the Stokes and SSA solutions (x G 0 ≈ 525 km for the Stokes against x G 0 ≈ 605 km for the SSA Pattyn et al., 2013;Seroussi et al., 2014;Feldmann et al., 2014).In what follows, the transient response is discussed relative to the steady GL position x G 0 obtained for each friction implementation.
The displacement of the GL relative to its initial steady position is found to be substantially different for the three friction implementations, for both the perturbation experiment P75S and the reversal of the perturbation experiment P75R (see Fig. S4).Such large differences for the transient response of the three methods can only be explained by a too coarse mesh.The steady solution being reasonably close, and independent of the lateral discretisation of the mesh (no transverse variation of any field so that the steady GL is a straight line perpendicular to the x direction), the source of discrepancy for the transient response certainly arises from the lateral discretisation.The number of lateral elements N y is only 20 for the previous simulations.The sensitivity of the transient response to the lateral discretisation is investigated by running the same experiment with two finer lateral mesh resolutions, everything else being the same.The results for lateral resolutions with N y = 20, N y = 40 and N y = 80 elements in the lateral direction are presented in the Figs.S4, S5 and S6, respectively.Figure 3 shows the differences from such lateral resolution visualised relative to the highest resolution N y = 80.As can be seen from Fig. 3, differences in the transient response of the three methods are significantly decreased when the lateral mesh refinement is increased.Nevertheless, even with the finest mesh (N y = 80), the difference between the methods stays relatively important (≈ 5 km between LG and FF at the end of the perturbation experiment, but to be compared to 17 km for N y = 20).Figure 3 indicates that the difference for the three methods between the higher resolution (N y = 80) and the two other mesh refinements (N y = 40 and N y = 20) is smaller for the DI method than the two others.In other words, the DI method seems to be less sensitive to the mesh refinement than the two other methods, certainly because it gives an intermediate solution whatever the mesh resolution.This is one more reason that justifies that the DI method should be preferentially adopted for future works.Note however that the decrease in mesh sensitivity is not as high as for the subgrid methods proposed for the SSA (Seroussi et al., 2014).
Higher lateral discretisation was not further explored for computing resource reasons, but this study clearly indicates that, as expected theoretically and shown in the previous section using the flow-line setup MISMIP, the difference between the three implementations is decreased as the mesh resolution is increased.Published LFA results (Pattyn et al., 2013) were obtained with a lateral discretisation of N y = 20 elements, which was certainly insufficient as shown by these new results using 40 and 80 lateral elements.For further comparisons, we recommend to use the more accurate results presented in Fig. S6 and provided as the Supplement.

Conclusions
In this paper, the sensitivity to the mesh resolution of the dynamical response of the GL is studied for different friction transition schemes upstream the GL.Contrary to Leguy et al. (2014), a smoother friction vanishing at the GL is not found to improve model sensitivity to mesh resolution.Explaining the reasons for such different behaviour is beyond the scope of this paper, but we encourage further works in that direction with various models and various smoothing functions for the friction upstream the GL.Having the friction smoothly decreasing to zero at the GL is certainly more realistic, as one expects the effective pressure to vanish at the GL.Therefore, even if it might present no advantage in terms of mesh sensitivity, such more realistic friction distribution should be preferred for future model intercomparisons.
In the case of a discontinuous friction, as in the MISMIP and MISMIP3d experiments, we have presented three possible implementations of the friction at the GL for a finite element formulation of the Stokes equations.So far, in all the applications using Elmer/Ice, it was assumed that the friction is applied up to the GL using the LG method.In doing so, the first elements immediately downstream from the GL undergo a little friction even if in contact with the ocean.We have shown that the treatment of the friction at the GL has a strong influence on both the velocity field and on the resulting GL dynamics for the mesh resolutions that were used to produce the MISMIP and MISMIP3d results.As expected theoretically, differences between the three implementations are shown to decrease as the mesh resolution is increased, but these differences remain substantial when using mesh resolutions that are numerically affordable for usual 3-D applications.Even for the smallest refinements accessed for the three-dimensional test case, differences are still observed.However, these differences are much smaller than those between Stokes and lower-order models.This gives an indication of the model error to be expected when performing GL dynamics simulations with a Stokes model.Moreover, using the MISMIP3d experiment, the lateral refinement is shown to have also a significant influence on the transient behaviour.
In the case of a discontinuous friction at the GL, we finally recommend to use the discontinuous DI implementation, which is certainly the most realistic and the less sensitive to the mesh refinement of the three.We also recommend to use these newly published results with finer mesh resolutions for future model comparison.
The Supplement related to this article is available online at doi:10.5194/tc-10-307-2016-supplement.

Figure 2 .
Figure 2. Experiment MISMIP 3a step 3: (a) grounding line positions as a function of resolution in advance (stars) and retreat (dots) for L = 0 m and the three GL implementations LG (brown), DI (purple) and FF (blue), L = 60 m (red) and L = 500 m (black), (b) model accuracy estimated from the difference between the retreat and advance GL steady positions (same colour legend).In (a), the large white star corresponds to the published GL position for step 4 of experience 3a in Pattyn et al. (2012) and the dot-dashed line is the Schoof (2007) solution.

Figure 3 .
Figure 3. Influence of lateral resolution from experiment MIS-MIP3d P75S and P75R: evolution of the absolute differences in kilometres between the highest resolution (N y = 80) and the two others (N y = 40 continuous line and N y = 20 dashed line) for the three different methods: LG (brown), DI (purple) and FF (blue), on the symmetry axis (y = 0; thick curves) and on the free-slip boundary (y = 50 km; thin curves).The initial results used to plot this figure are presented in Figs.S4, S5 and S6 for lateral resolutions N y = 20, N y = 40 and N y = 80, respectively.