A full Stokes subgrid scheme in two dimensions for simulation of grounding line migration in ice sheets using Elmer/ICE (v8.3)

The flow of large ice sheets and glaciers can be simulated by solving the full Stokes equations using a finite element method. The simulation is particularly sensitive to the discretization of the grounding line, which separates ice resting on bedrock and ice floating on water, and is moving in time. The boundary conditions at the ice base are enforced by Nitsche’s method and a subgrid treatment of the grounding line element. Simulations with the method in two dimensions for an advancing and a retreating grounding line illustrate the performance of the method. The computed grounding line position is compared 5 to previously published data with a fine mesh, showing that similar accuracy is obtained using subgrid modeling with more than 20 times coarser meshes. This subgrid scheme is implemented in the two dimensional version of the open source code Elmer/ICE.

A shallow approximation to FS with a subgrid scheme on coarse meshes is compared to FS in Feldmann et al. (2014) with similar results for the GL migration. Subgrid modeling and adaptivity are compared in Cornford et al. (2016) for a vertically integrated model. The thickness of the ice above flotation determines if the ice is grounded or floating. A fine mesh resolution 60 is necessary for converged GL positions with FS in Durand et al. (2009a, b). A dynamic mesh refinement and coarsening of the mesh following the GL would solve the problem in palaeo simulations when the GL moves long distances. An alternative is to introduce a subgrid scheme in the mesh elements where the GL is located in a static mesh and keep the mesh size coarser everywhere else in the ice sheet.

Proposed method and outline of the paper
From the above we conclude that, it seems crucial that the ice model includes equations with vertical shear stress in the neighbourhood of the GL, and one way to avoid the fine meshes that are otherwise needed close to the GL is to introduce a subgrid scheme in the discretization element where the GL is located. In this study, we develop such a numerical method for the FS equations in two dimensions introducing a subgrid scheme in the mesh element where the GL is located. Since the subgrid scheme is restricted to one element in a 2D vertical ice this is computationally inexpensive. In an extension to 3D, the 70 subgrid scheme would be applied along a line of elements in 3D. The results with numerical modeling will always depend on the mesh resolution but can be more or less sensitive to the mesh spacing and time steps. It depends on the equation, the mesh size, the mesh quality, and the finite element spaces in the approximation.
We solve the FS equations in a 2D vertical ice with the Galerkin method implemented in Elmer/ICE . A subgrid discretization is proposed and tested for the element where the GL is located. The boundary conditions are 75 imposed by Nitsche's method at the ice base in the weak formulation of the equations (Nitsche, 1971;Urquiza et al., 2014;Reusken et al., 2017). The linear Stokes equations are solved in Chouly et al. (2017a) with Nitsche's treatment of the boundary conditions. They solve the equations for the displacement but here we solve for the velocity using similar numerical techniques to weakly impose the Dirichlet boundary conditions on the normal velocity at the base. The frictional force in the tangential direction is applied on part of the element with the GL. The position of the GL within the element is determined in agreement 80 with theory developed for the linearized FS in Schoof (2011).
The paper is organized as follows. Section 2 is devoted to the presentation of the mathematical model of the ice sheet dynamics. In Sect. 3, the numerical discretization with FEM is given while the subgrid scheme around the GL is found in Sect.
4. The numerical results for a MISMIP problem are presented in Sect. 5. The extension to three dimensions (3D) is discussed in Sect. 6 and finally some conclusions are drawn in Sect. 7.

The full Stokes (FS) equations
To simulate flow in a 2D vertical cross-section of an ice sheet, we use the FS equations with coordinates x = (x, z) T (Hutter, 1983). The nonlinear partial differential equations (PDEs) in the interior of the ice domain Ω are given by 90 where the stress tensor is σ = τ (u) − pI and the deviatoric stress tensor is τ (u) = 2η(u)˙ (u). The strain rate tensor is defined bẏ I is the identity matrix, and the viscosity is defined by Glen's flow law η(u) = 1 2 (A(T )) − 1 n˙ 1−n n e ,˙ e = 1 2 tr(˙ (u)˙ (u)). (3)

95
Here u = (u, w) T is the vector of velocities, ρ is the density of the ice, p denotes the pressure, and the gravitational vector is denoted by g. The viscosity η is a function of the rate factor A(T ) where T is the ice temperature. For isothermal flow assumed here, the rate factor A is constant. Finally, n is usually taken to be 3.

100
At the boundary Γ of the ice domain Ω we define the normal outgoing vector n and tangential vector t (see Fig. 1). In the 2D vertical case considered here, the ice sheet geometry is constant in y. The ice surface is denoted by Γ s and the ice base is At Γ s and Γ bf , the floating part of Γ b , we have that respectively. The ice is stress-free at Γ s , f s = 0, and f bf = −p w n at the ice/ocean interface Γ bf where p w is the water pressure.

105
Let where σ nn and σ nt are the normal and tangential components of the stress and u t is the tangential component of the ice velocity at the ice base. Then for the slip boundary Γ bg , the grounded part of Γ b where the ice rests on the bedrock, we have a friction law for the sliding ice n t a s * * * * where u n is the normal component of the ice velocity. The type of friction law is determined by the friction coefficient β (≥ 0).
At Γ bf , there is a balance between σ nn and p w and the contact is friction-free, β = 0. Then At the GL, the boundary condition switches from β > 0 and u n = 0 on Γ bg to β = 0 and a free u n on Γ bf . In a 2D vertical 115 cross-section of ice, the GL is the point (x GL , z GL ) shared between Γ bg and Γ bf .
The ocean surface is at z = 0, and p w = −ρ w gz b . The density of sea water is denoted byρ w , z b is the z-coordinate of Γ b , and g is the vertical component of the gravitational force.

The free surface equations
The boundaries Γ s and Γ b are time-dependent and move according to two free surface equations. The boundary Γ bg follows 120 the fixed bedrock with coordinates (x, b(x)).
The z-coordinate of the ice surface position z s (x, t) at Γ s (see Fig. 1) is the solution of an advection equation where a s denotes the surface mass balance and u s = (u s , w s ) T the velocity at the ice surface in contact with the atmosphere.
Similarly, the z-coordinate for the ice base z b of the floating ice at Γ bf satisfies where a b is the basal mass balance and u b = (u b , w b ) T the velocity of the ice at Γ bf . On Γ bg , The thickness of the ice is denoted by H = z s − z b and depends on x and t.

A first order solution close to the grounding line
The 2D vertical solution of the FS equations in Eq. (1) with a constant viscosity, n = 1 in Eq. (3), is expanded in small 130 parameters in Schoof (2011). The solutions in different regions around the GL are connected by matched asymptotics. Upstream of the GL at the grounded part, x < x GL , the leading terms in the expansion satisfy a simple relation in scaled variables close to the GL. Across the GL, the ice velocity u, the flux of ice uH, and the depth-integrated normal or longitudinal stress τ 11 in Eq. (2) are continuous. By including higher order terms in the expansion in small parameters, it is shown in Schoof (2011, Sect. 4.7) that the ice surface slope is continuous and Archimedes' flotation condition is not satisfied immediately downstream of the GL. A rapid variation in the vertical velocity w in a short distance interval at the GL causes oscillations in the ice surface in the analysis as also observed in FS simulations in Durand et al. (2009a). The flotation condition in (9) determines where the GL is in SSA in Docquier et al. (2011); Drouet et al. (2013).
In Schoof (2011, Sect. 4.3), the solution to the FS in a 2D vertical cross-section of ice is expanded in two parameters, ν and .

140
The aspect ratio of the ice ν is the quotient between a typical scale of the thickness of the ice H and a horizontal length scale L, ν = H/L, and is ν times the quotient between the longitudinal and the shear stresses τ 11 and τ 12 in Eq.
(2). If ν 5/2 1 then in a boundary layer close to the GL and x < x GL it follows from the equations that the leading terms in the solution in scaled variables satisfy

145
On floating ice τ 22 − p + p w = 0 and the hydrostatic flotation criterion Eq. (9) is fulfilled. This is a first order approximation of the second relation in Eq. (6). On the grounded ice domain, we have τ 22 − p + p w < 0.
Introducing the notation and letting H bw = −z b be the thickness of the ice below the sea level yields If x < x GL then χ a < 0 in the neighborhood of x GL on Γ bg and if x > x GL then χ a = 0 and Eq. (9) holds true on Γ bf . Suppose that z s and z b are linear in x. Then χ a is also linear in x. In numerical experiments with the linear FS (n = 1) in Nowicki and In Sect. 4, we take this same approach but use an indicator χ(x) orχ(x) derived from the solutions of the nonlinear FS 155 equations to estimate the GL position. These indicators are approximated by χ a (x, z b ).
In this section we state the weak form of Eq. (1), introduce the spatial FEM discretization used for Eq. (1), and give the time-discretization of Eq. (7) and (8). 160 We start by defining the mixed weak form of the FS equations. Introduce k = 1 + 1/n, k * = 1 + n with n from Glen's flow law and the spaces

The weak form of the FS equations
see, e.g. Jouvet and Rappaz (2011, Eq. (3.7)), Chen et al. (2013, Sect. 3.1), Martin and Monnier (2014, Eq. (21)). The weak is satisfied, where The last term in B N is added in the weak form in Nitsche's method (Nitsche, 1971) to impose the Dirichlet condition u n = 0 weakly on Γ bg . It can be considered as a penalty term. Since u = u n n + u t t, the contribution of the tangential force can also 170 be written βu · v when u n = 0. The value of the positive parameter γ 0 depends on the physical problem and h is a measure of the mesh size on Γ b . The sensitivity of the GL positions for different values of γ 0 is shown in Sect. 5. The first term in B N symmetrizes the boundary term B Γ + B N on Γ bg and vanishes when u n = 0. The boundary term F Γ (v) is from the buoyancy force at the ice/ocean interface in (6) where p w depends on z b on Γ bf .

The discretized FS equations
We employ linear Lagrange elements with Galerkin Least Square (GLS) stabilization (Franca and Frey, 1992;Helanow and Ahlkrona, 2018) to avoid spurious oscillations in the pressure using the standard MISMIP setting in Elmer/ICE (Durand et al., 2009a;Gagliardini et al., 2013) approximating solutions in the spaces V k and Q k * in Eq. (13).
The mesh is constructed from a footprint mesh on the ice base and then extruded with the same number of layers equidistantly in the vertical direction according to the thickness of the ice sheet. To simplify the implementation in 2D, the footprint mesh 180 on the ice base consists of N + 1 nodes at In general, the GL is somewhere in the interior of an interval [x i−1 , x i ] and it crosses the interval boundaries as it moves forward in the advance phase and backward in the retreat phase of the ice. The advantage with Nitsche's method of formulating the boundary conditions is that if x GL ∈ [x i−1 , x i ] then the boundary integral over the interval can be split into two parts in Eq.
In the GL element, we have  (14). A large γ 0 will force u n to be small. The penalty term may change inside an element as in (15) where it is = 0 only in the grounded part.
The resulting system of nonlinear equations form a nonlinear complementarity problem (Christensen et al., 1998). The 195 distance d between the base of the ice and the bedrock at time t and at x is If d > 0 on Γ bf then the ice is not in contact with the bedrock and σ nn + p w = 0 and if σ nn + p w < 0 on Γ bg then the ice and the bedrock are in contact and d = 0. Hence, the complementarity relation in the vertical direction is

200
The contact friction law is such that β > 0 when x < x GL and β = 0 when x > x GL . The complementarity relation along the ice base at x is then the non-negativity of d and In particular, these relations are valid at the nodes x = x j , j = 0, 1, . . . , N .
The complementarity condition also holds for u n and σ nn such that without any sign constraint on u n except for the retreat phase when the ice leaves the ground and u n < 0.
Similar implementations for contact problems using Nitsche's method are found in Chouly et al. (2017a, b), where the unknowns in the PDEs are the displacement fields instead of the velocity in Eq. (1). Analysis in Chouly et al. (2017a) suggests that Nitsche's method for the contact problem can provide a stable numerical solution with an optimal convergence rate.

210
The nonlinear equations, Eq. (14), for the nodal values of u and p are solved by Picard iterations. The system of linear equations in every Picard iteration is solved directly by using the MUMPS linear solver in Elmer/ICE. The condition on is used to decide if the node x j is geometrically grounded or floating. It is computed at each time step and is not changed during the nonlinear iterations (Picard). The procedure for solution of the nonlinear FS equations is outlined in Algorithm 1. In two dimensions, the GL will be located in one element.

215
Algorithm 1 Solve the FS equations For a given mesh, compute dj, j = 0, 1, ..., N, for all the nodes xj at the ice base.
Mark node j as geometrically grounded if dj < 10 −3 , otherwise floating.
Find the element which contains both geometrically grounded and floating nodes, and mark the grounded node in this element as 'GL node'.
Compute the residual of the FS equations with the initial guess of the solution.
while the residual is larger than the tolerance do Assemble the FEM matrix for the interior of the domain Ω.
for the boundary elements on Γ b do if has 'GL node' then Mark the current element as a 'potential GL element'.
Use the subgrid scheme in Algorithm 3 of Sect. 4 for the assembly. Compute the solution and the residual. The advection equations for the moving ice boundary in Eq. (7) and (8) are discretized in time by a finite difference method and in space by FEM with linear Lagrange elements for z s and z b . An artificial diffusion stabilization term is added, making the spatial discretization behave like an upwind scheme in the direction of the velocity as implemented in Elmer/ICE.
The advection equations Eq. (7) and Eq. (8) are integrated in time by a semi-implicit method of first order accuracy. Let 220 c = s or b. Then the solution is advanced from time t to t +1 = t + ∆t with the time step ∆t by The spatial derivative of z c is approximated by FEM as described above. A system of linear equations is solved at t +1 for z +1 c . This time discretization and its properties are discussed in Cheng et al. (2017) and summarized in Algorithm 2.
Algorithm 2 Time scheme of the GL migration problem Start from an initial geometry Ω 0 defined by z 0 b , z 0 s .
Solve the FS equations on Ω with Algorithm 1, to get the solution u .
Solve for z +1 b and z +1 s with u by the semi-implicit Euler method.

end for
A numerical stability problem in z b is encountered in the boundary condition at Γ bf when the FS equations are solved in 225 Durand et al. (2009a). It is resolved by expressing z b in p w at Γ bf with a damping term. An alternative interpretation of the idea in Durand et al. (2009a) and an explanation follow below.
The relation between u n and u t at Γ bf and u where z bx denotes ∂z b /∂x. Inserting u b and w b from Eq. (21) into Eq. (8) yields Instead of discretizing Eq. (22) explicitly at t +1 with u n to determine p +1 w , the base coordinate is updated implicitly Assuming that z bx is small, the time step restriction in Eq. (23) is estimated by considering a 2D slab of the floating ice of 235 width ∆x and thickness H. Newton's law of motion yields where M = ∆x(z s − z b )ρ is the mass of the slab. Dividing by M , integrating in time for u n (t m ), letting m = + 1 or , and approximating the integral by the trapezoidal rule for the quadrature yields with the parameters Then insert u m n into Eq. (23). All terms in u m n from time steps i < m are collected in the sum ∆tF m−1 . Then Eq. (23) can be written 245 For small changes in z b in Eq. (24), the explicit method with m = is stable when ∆t is so small that When H = 100 m on the ice shelf, ∆t < 6.1 s which is far smaller than the stable steps for Eq. (20). Choosing the implicit scheme with m = + 1, the bound on ∆t is x j is denoted by E j . The sequence of E j , j = 1, . . . , N, approximates Γ b . The grounding line element containing the GL is E i . Depending on how the mesh is created from the initial geometry and updated during the simulation, the first floating node at x i , as well as the GL element, can be either on the bedrock (as in Fig. 2) or at the ice base above the bedrock (as in Fig. 3), even though the corresponding analytical solutions are identical. Denote the situation in Fig. 2 as case i, and the one in Fig.   275 3 as case ii. The physical boundary conditions of the two cases are different only at the GL element. More precisely, in case i, the net force in the vertical direction on the node x i is pointing inward, namely χ(x i ) = σ nn (x i ) + p w (x i ) > 0, whereas in case ii, the floating condition σ nn (x i ) + p w (x i ) = 0 is satisfied in the node x i . The directions of the vertical net force at the nodes x i−1 and x i are shown by the arrows in the upper panels of Fig. 2 and 3. Consequently, the external forces and boundary conditions imposed on the GL element are different in the two cases. For instance, in case i, the GL element is considered These two cases are similar to the LG and FF cases in Gagliardini et al. (2016) implying that the numerical solutions in the two cases are different, especially on a coarse mesh (mesh size at about 100 m or larger). Thus, we propose a subgrid scheme to reduce these differences in the spatial discretization and to capture the GL migration without using a fine mesh resolution 285 (< 100 m). The schematic drawing of the subgrid scheme for the two cases is shown in the middle panels of Fig. 2 and 3. The GL element is divided into the grounded (red) and floating (blue) parts by the estimated GL positionx GL on E i , which is the numerical approximation of the analytical GL position x GL .
The GL moves toward the ocean in the advance phase and away from the ocean in the retreat phase. First, we consider case i in the advance phase and define the indicator by 290 which vanishes on the floating ice and is negative and approximately equal to χ a = τ 22 − p + p w in Eq. (11) on the ground since the slope of the bedrock is small and n ≈ (0, −1) T . Because of the poor spatial resolution of the coarse mesh, χ(x i ) is positive.
The water pressure p w (x) is a linear function of x on the GL element and the numerical solution of σ nn (x) is also piecewise linear on every element with the standard Lagrange elements in Elmer/ICE . Hence, it makes sense to approximate the analytical GL position x GL byx GL by linear interpolation in the current framework. This approach fits well 300 with case i since the indicator χ(x) has opposite signs at x i−1 and x i , see the middle panel of Fig. 2 wherex GL is marked by a red ' * '. It guarantees the existence and uniqueness ofx GL on the GL element.
Another situation in the advance phase is case ii shown in Fig. 3. As the elements on both sides of the node x i are geometrically floating, the boundary condition imposed on x i becomes χ(x i ) = σ nn (x i )+p w (x i ) = 0. However, the implicit treatment of the ice base moves the z-coordinate of the node x i towards the bedrock with u n > 0 in Eq. (23) as discussed in Sect. 3.3.

305
The result is that p w defined by the implicit z b in (23) satisfies σ nn + p w > 0 in (27) and χ(x i ) > 0.
The implicit treatment of the ice base has the consequence that only case ii occurs in the retreat phase. When the FS equations are solved, the implicit update of the ice base with u n < 0 in Eq. (23) implies that the last grounded node in the previous time step is leaving the bedrock when the ice is retreating and the GL moves back to the adjacent element. Case i will not appear in In this circumstance, χ(x i ) = 0 in the floating node and a correction of χ(x) is introduced 310 into case ii byχ iñ Here p b (x) = −ρ w gb(x) is the water pressure on the bedrock corresponding to linear extrapolation of the pressure for x > x GL along the element on the bedrock. Furthermore,χ(x) ≥ χ(x). Notice that p b (x i ) = p w (x i ) > p w (x i ), wherex i is a point on the bedrock with the same x coordinate of x i , as illustrated in the middle panel of Fig. 3. Both χ(x) in (27) andχ(x) in (29) are 315 nonlinear in x but the numerical approximation of them will vary linearly in x. A solutionx GL is found by linear interpolation ofχ(x) between the nodes x i−1 and x i as in Eq. (28). It follows from Eq. (28) thatx GL is located on the element boundary, see Figs. 2 and 3. If we compare with case i, this correction can be considered as using σ nn (x GL ) to approximate σ nn (x GL ) on a virtual element between x i−1 andx i , see Fig. 3. The positionx GL is a numerical approximation of the analytical GL position, although it is not geometrically in contact with the bedrock.

320
Since we have p b (x) = p w (x) and χ(x) =χ(x) at the GL element in case i, we can simply useχ(x) to findx GL for the two cases by replacing χ in (28) byχ.
The domains Γ bg and Γ bf are separated atx GL as in Eq. (15) and the integrals on the GL element are calculated with a high-order integration scheme as in Seroussi et al. (2014). We introduce two step functions H N (x) and H β (x) to include and exclude quadrature points in the integration of Nitsche's term and the slip boundary condition, respectively. They are defined The penalty term in Nitsche's method restricts the motion of the element in the normal direction. It is only imposed on an element which is fully geometrically on the ground in case i. On the contrary in case ii, the GL element E i is not in contact with the bedrock, see Fig. 3. The normal velocity on the element should not be forced to zero and only the floating 330 boundary condition is then used on the GL element. Nitsche's penalty term should be imposed on all the fully geometrically grounded elements and partially on the GL element in the advance phase as in case i. The step function H N (x) indicates how Nitsche's method is implemented on the basal elements, see the lower panels of Fig. 2 and 3 for the two cases. The penalty term contributes to the integration only when H N (x) = 1.
The slip coefficient β is treated similarly with the step function H β (x), where H β (x) = 1 is on the fully geometrically 335 grounded elements and H β (x) = 0 on the floating elements. To further smooth the transition of β at the GL, the step function is set to be 1/2 in parts of the GL element before integrating using the high order scheme. The convergence of the nonlinear iterations is improved in this way. In case i, full friction is applied at the grounded part between x i−1 andx GL of the GL element since this part is also geometrically grounded in the analytical solution of the FS as in Fig. 2. Then, the friction is lower in the remaining part of E i . For the floating part betweenx GL and x i in case ii, there is no friction and H β (x) = 0 and 340 we have reduced friction between x i−1 andx GL , see the lower panel of Fig. 3. The boundary integral Eq. (15) on E i is now rewritten with the two step functions as A summary of the numerical treatment of the GL is: -Advance phase ⇒ indicator χ in (27), case i or case ii 345 -Retreat phase ⇒ indicatorχ in (29), case ii The case is determined by the geometry of the GL element and the sign of the indicator χ.
The algorithm for the GL element is: Algorithm 3 Subgrid modeling for the GL element Take all the 'potential GL elements' and solve χ(x) = 0 (advance phase) orχ(x) = 0 (retreat phase) to findxGL and the GL element.
Determine which case this GL element belongs to by checking the geometrical conditions at xi.
Specify HN (x) and H β (x) based onxGL depending on the case and the advance or retreat phase.
Equations (1), (7), and (8) form a system of coupled nonlinear equations. They are solved by Elmer/ICE v.8.3 in the same manner as Durand et al. (2009b); Gagliardini et al. (2013Gagliardini et al. ( , 2016. The detailed procedure is explained in Algorithms 1, 2, and with χ orχ and the solution σ nn (x) from the previous nonlinear iteration, and the step functions H N and H β are adjusted accordingly. The water pressure p b is fixed since the ice geometry is not changed during the nonlinear iterations.

355
The numerical experiments follow the MISMIP benchmark (Pattyn et al., 2012) and a comparison is made with the results in Gagliardini et al. (2016). Using the experiment MISMIP 3a, the setups are exactly the same as in the advancing and retreating simulations in Gagliardini et al. (2016). The experiments are run with spatial resolutions of ∆x = 4 km, 2 km, 1 km and 0.5 km.
The mesh at the base is extruded vertically in 20 layers with equidistantly placed nodes in each vertical column. The time step is ∆t = 0.125 year for all four resolutions to eliminate time discretization errors when comparing different spatial resolutions.

360
The dependence on γ 0 in (30) for the retreating ice is shown in Fig. 4  The distance between the steady state GL positions of the retreat and the advance phases is shown in Fig. 6 (b). The maximal distance is about 6 km at ∆x = 1 km with the subgrid model, whereas in Gagliardini et al. (2016), the resolution has to be below 50 m to achieve a similar result. We observe oscillations at the ice surface near the GL in all the experiments as expected from Durand et al. (2009a); Schoof 375 (2011). A zoom-in plot of the surface elevation computed with four different mesh sizes ∆x = 0.5, 1, 2, 4 km after an advance simulation to t = 10000 years is found to the left in Fig. 7. The abscissa is the distance from the steady state GL position for each mesh size. In general, the estimated GL position does not coincide with any nodes even at the steady state but it may be close to a node.
The ratio between the thickness below sea level H bw and the ice thickness H is shown to the right in Fig. 7. As in the left 380 panel, the ratio is plotted versus the distance from the GL achieved with the particular mesh size. The horizontal, purple, dashdotted line represents the ratio of ρ/ρ w . The solutions vary smoothly over the mesh with ∆x = 0.5 km which appears to be a sufficient resolution in both panels of Fig. 7. Moreover, the solutions converge regularly toward the solution with ∆x = 0.5 km when the mesh size decreases.
The result for ∆x = 0.5 km in the right panel confirms that the hydrostatic assumption Hρ = H bw ρ w in Eq. (9)  at the GL are almost discontinuous as analyzed in Schoof (2011). With the subgrid model, the rapid variation is captured with ∆x = 0.5 km. The convergence for decreasing mesh size behaves smoothly. w is zoomed-in close to the GL with the distance to the mesh dependent GL on the x-axis. Seroussi et al. (2014) describe four different subgrid models (NSEP, SEP1, SEP2 and SEP3) for the friction in SSA and evaluate them in a FEM discretization on a triangulated, planar domain. The flotation criterion is applied at the nodes of the triangles. In the NSEP, an element is floating or not depending on how many of the nodes of that element are floating. In the other three methods, an inner structure in the triangular element is introduced. One part of a triangle is floating and one part is grounded. The amount of friction in a triangle with the GL is determined by the flotation criterion. Either the friction 400 coefficient is reduced, the integration in the element only includes the grounded part, or a higher order polynomial integration (SEP3) is applied. Faster convergence as the mesh is refined is observed for the latter methods compared to the first method.
The discretization of the friction in Sect. 4 is similar to the SEP3 method but the FS equations also require a subgrid treatment of the normal velocity condition. In the method for the FS equations in Gagliardini et al. (2016), the GL position is in a node and the friction coefficient is approximated in three different ways. The coefficient is discontinuous at the node in one case (DI 405 in Gagliardini et al. (2016)). Our coefficient is also discontinuous but at the estimated location of the GL between the nodes.
The convergence of the steady state GL position toward the reference solutions in Gagliardini et al. (2016) is observed in the simulations in Fig. 5 and 6. However, as the meshes we used are at least 20 times larger than the 25 m finest resolution in Gagliardini et al. (2016), it has probably not reached the convergence asymptote. At the current resolutions, the discretization introduces a strong mesh effect such as the two different geometrical interpretations in the two cases mentioned in Sect. 4.

410
The subgrid scheme is able to provide a more accurate representation of the GL position and the boundary conditions, but the numerical solution of the velocity field, pressure as well as the two free surfaces are still computed on the coarse mesh, which are the main sources of the numerical errors. Additional uncertainty at the GL is introduced by the approximation of the bedrock geometry, the friction at the GL, and the modeling of the ice/ocean interaction. It is shown in Cheng and Lötstedt (2020) that the solution at the GL is particularly sensitive to variation in the geometry and friction at the ice base.

415
Our method can be extended to a triangular mesh covering Γ b in the following way (considering linear Lagrange functions).
The condition on χ in Eq. (27) orχ in Eq. (29) is applied on the edges of each triangle T in the mesh. If χ < 0 in all three nodes then T is grounded. If χ ≥ 0 in all nodes then T is floating. The GL passes inside T if χ has a different sign in one of the nodes. Then the GL crosses the two edges where χ < 0 in one node and χ ≥ 0 in the other node. In this way, a continuous reconstruction of a piecewise linear GL is possible on Γ b . The same tests are applied toχ. The FEM approximation is modified 420 in the same manner as in Sect. 4 using step functions in Nitsche's method.
An alternative to a subgrid scheme is to introduce static or dynamic adaptation of the mesh on Γ b with a refinement at the GL as in e.g. Gladstone et al. (2010a); Cornford et al. (2013); Drouet et al. (2013). In general, a fine mesh is needed at the GL and in an area surrounding it. Since the GL moves long distances in simulations of palaeo-ice sheets, the adaptation should be dynamic, permit refinement and coarsening of the mesh varying in time, and be based on some estimate of the numerical 425 error of the method. In shorter time intervals, a static adaptation may be sufficient since the GL will move a shorter distance.
Furthermore, shorter time steps are necessary for numerical stability in static and dynamic mesh adaptation schemes. A static adaptation is determined once before the simulation starts. Introducing a time dependent, dynamic mesh with adaptivity into an existing code requires a substantial coding effort and will increase the computational work considerably compared to a static mesh. Subgrid modeling is easier to implement and the increase in computing time is small. A combination of dynamic mesh 430 adaptation and subgrid discretization may be the ultimate solution. Then the mesh at the GL would be adapted to resolve the variation in the interior of the ice at the GL while the subgrid modeling would handle the discontinuity at the basal boundary.

Conclusions
A subgrid scheme at the GL has been developed and tested in the SSA model for 2D vertical ice flow in Gladstone et al. (2010b) and in Seroussi et al. (2014), for the friction in the vertically integrated model BISICLES (Cornford et al., 2013) for 2D flow 435 in Cornford et al. (2016), and for the PISM model mixing SIA with SSA in 3D in Feldmann et al. (2014). Here we propose a subgrid scheme for the FS equations for a 2D vertical cross-section of ice, implemented in Elmer/ICE, that can be extended to 3D. The mesh is static and the moving GL position within one element is determined by linear interpolation with an auxiliary function χ(x) orχ(x). Only in that element, the FEM discretization is modified to accommodate the discontinuities in the boundary conditions.

440
The numerical scheme is applied to the simulation of a 2D vertical ice sheet with an advancing GL and one with a retreating GL. The model setups for the tests are the same as in one of the MISMIP examples (Pattyn et al., 2012) and in Gagliardini et al. (2016). The solution converges smoothly in the neighborhood of the GL when the mesh size is reduced. Comparable results to Gagliardini et al. (2016) are obtained using the subgrid scheme with more than 20 times larger mesh sizes. A larger mesh size also allows a longer time step for the time integration.