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 with 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 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.

1 Introduction

1.1 Ice sheet dynamics, sea-level rise, and grounding line migration

Numerical simulation of ice sheet flow is necessary to assess the future sea-level rise (SLR) due to melting of continental ice sheets and glaciers (Hanna et al., 2013) and to reconstruct the ice sheets of the past (Stokes et al., 2015; DeConto and Pollard, 2016) for comparison with measurements and validation of the models. Ice sheet model predictions are particularly sensitive to the numerical treatment of the grounding line (GL) (Durand and Pattyn, 2015; Konrad et al., 2018), the line where the ice sheet leaves the solid bedrock and becomes an ice shelf floating on water driven by buoyancy.

The distance that the GL moves may be long over paleo-timescales. In Kingslake et al. (2018) it is shown that the GL retreated several hundred kilometers in West Antarctica during the last 11 500 years and then advanced again after the isostatic rebound of the bed. The sensitivity, long time intervals, and long distances of the GL migration require a careful treatment of the GL and its neighborhood in the numerical method used to discretize the equations modeling the ice sheet dynamics. In this paper, we develop an accurate and efficient method for such problems.

1.2 Model equations

When the ice rests on the ground and is affected by large frictional forces on the bed, the ice flow is dominated by vertical shear stresses. On the other hand, when the ice is floating on water, it is the longitudinal stress gradient that controls the flow of the ice. The GL is in the transition zone between these two types of flow with a gradual change of the stress field (Schoof, 2011).

The most accurate ice model in theory is based on the full Stokes (FS) equations. A simplification of the FS equations by integrating the depth of the ice is the shallow shelf (or shelfy stream) approximation (SSA) (MacAyeal, 1989), which is often used for simulation of the coupling between a grounded ice sheet and a marine ice shelf. In the zone between the grounded ice and the floating ice, it is necessary to use the FS equations (Wilchinsky and Chugunov, 2000; Schoof and Hindmarsh, 2010; Docquier et al., 2011; Schoof, 2011) unless the ice is moving rapidly on the ground with low basal friction, when the SSA equations are accurate both upstream and downstream of the GL.

The evolution of the GL in simulations is sensitive to the model equations and the basal friction law. In the Marine Ice Sheet Model Intercomparison Project (MISMIP) (Pattyn et al., 2012, 2013), different ice models and implementations solve the same ice flow problems, and the predicted GL steady state and transient GL motion are compared. The results show that the position of the GL depends on the model equations (Pattyn et al., 2013). Predictions of the GL position and SLR are different for different ice models such as FS and SSA (Pattyn and Durand, 2013). Including equations with vertical shear stress at the GL such as the FS equations is crucial to accurately resolve GL dynamics in a wide range of circumstances.

The friction laws at the ice base depend on the effective pressure, the basal velocity, and the distance to the GL in different combinations in Leguy et al. (2014), Gagliardini et al. (2015), Brondex et al. (2017), and Gladstone et al. (2017). The GL position and SLR vary considerably depending on the choice of friction law. Given the friction law, the results are sensitive to its model parameters too (Gong et al., 2017).

1.3 Numerical methods

Parameters in the numerical methods used to simulate ice sheet flow influence the GL migration. Durand et al. (2009b) find that the mesh resolution along the ice bed has to be fine to obtain reliable solutions with FS in GL simulations. The GL is then located in a node of the fixed or static mesh. A mesh size below 1 km is necessary in Larour et al. (2019) to resolve the features at the GL. Adaptive meshes for a finite volume discretization of an approximation of the FS equations are employed in Cornford et al. (2013) to study the GL retreat and loss of ice in West Antarctica. The FS solutions of benchmark problems in Pattyn et al. (2013) computed by an implementation of the finite element method (FEM) in Elmer/ICE (Gagliardini et al., 2013) and FELIX-S (Leng et al., 2012) are compared in Zhang et al. (2017). The differences between these implementations are attributed to different treatment of a friction parameter at the GL and different assignment of grounded and floating nodes and element faces.

A subgrid scheme introduces an inner structure in the discretization element or mesh volume where the GL is located. Such schemes have been developed for simplifications of the FS equations. A subgrid model for the GL is tested in Gladstone et al. (2010b) for the one-dimensional (1D) SSA equation where the flotation condition for the ice defines the position of the GL. The GL migration is determined by the two-dimensional (2D) SSA equations discretized by the FEM in Seroussi et al. (2014). Subgrid models at the GL are compared to a model without an internal structure in the element. The conclusion is that sub-element parameterization is necessary to obtain accurate results at reasonable computational expense. 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 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 paleo-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.

1.4 Proposed method and outline of the paper

From the above we conclude that it seems crucial that the ice model include equations with vertical shear stress in the neighborhood 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 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 (Gagliardini et al., 2013). A subgrid discretization is proposed and tested for the element where the GL is located. The boundary conditions are 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 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.

2 Ice model

2.1 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 $\mathit{x}=(x,z{)}^T$ (Hutter, 1983). The nonlinear partial differential equations (PDEs) in the interior of the ice domain Ω are given by

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0},\\ -\mathrm{\nabla }\cdot \mathit{\sigma }=\mathit{\rho }\mathit{g},\end{array}\right.\end{array}$$

where the stress tensor is $\mathit{\sigma }=\mathit{\tau }\left(\mathit{u}\right)-p\mathbf{I}$ and the deviatoric stress tensor is $\mathit{\tau }\left(\mathit{u}\right)=\mathrm{2}\mathit{\eta }\left(\mathit{u}\right)\dot{\mathit{ϵ}}\left(\mathit{u}\right)$. The strain rate tensor is defined by

$$\begin{array}{}\text{(2)}& \dot{\mathit{ϵ}}\left(\mathit{u}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^T)=\left(\begin{array}{cc}{\dot{\mathit{ϵ}}}_{\mathrm{11}}& {\dot{\mathit{ϵ}}}_{\mathrm{12}}\\ {\dot{\mathit{ϵ}}}_{\mathrm{12}}& {\dot{\mathit{ϵ}}}_{\mathrm{22}}\end{array}\right).\end{array}$$

I is the identity matrix, and the viscosity is defined by Glen's flow law:

$$\begin{array}{}\text{(3)}& \mathit{\eta }\left(\mathit{u}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathcal{A}\left(T^{\prime }\right)\right)}^{-\frac{\mathrm{1}}{n}}{\dot{\mathit{ϵ}}}_e^{\frac{\mathrm{1}-n}{n}},{\dot{\mathit{ϵ}}}_e=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\text{tr}\left(\dot{\mathit{ϵ}}\right(\mathit{u}\left)\dot{\mathit{ϵ}}\right(\mathit{u}\left)\right)}.\end{array}$$

Here $\mathit{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 𝒜(T^′), where T^′ is the ice temperature. For isothermal flow assumed here, the rate factor 𝒜 is constant. Finally, n is usually taken to be 3.

Figure 1 A two-dimensional schematic view of a marine ice sheet.

2.2 Boundary conditions

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 ${\mathrm{\Gamma }}_{\mathrm{b}}={\mathrm{\Gamma }}_{\mathrm{bg}}\cup {\mathrm{\Gamma }}_{\mathrm{bf}}$. At Γ_s, and Γ_bf, the floating part of Γ_b, we have

$$\begin{array}{}\text{(4)}& \begin{array}{rc}\mathit{\sigma }\mathit{n}={\mathit{f}}_{\mathrm{s}},& \mathit{\sigma }\mathit{n}={\mathit{f}}_{\mathrm{bf}},\end{array}\end{array}$$

respectively. The ice is stress free at Γ_s, f_s=0, and ${\mathit{f}}_{\mathrm{bf}}=-p_{\mathrm{w}}\mathit{n}$ at the ice–ocean interface Γ_bf, where p_w is the water pressure. Let

$${\mathit{\sigma }}_{\mathrm{nt}}=\mathit{t}\cdot \mathit{\sigma }\mathit{n},{\mathit{\sigma }}_{\mathrm{nn}}=\mathit{n}\cdot \mathit{\sigma }\mathit{n},u_{\mathrm{t}}=\mathit{t}\cdot \mathit{u},$$

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

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}_{\mathrm{nt}}+\mathit{\beta }(\mathit{u},\mathit{x})u_{\mathrm{t}}=\mathrm{0},u_{\mathrm{n}}=\mathit{n}\cdot \mathit{u}=\mathrm{0},-{\mathit{\sigma }}_{\mathrm{nn}}\ge p_{\mathrm{w}},\end{array}$$

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

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{\mathrm{nt}}=\mathrm{0},-{\mathit{\sigma }}_{\mathrm{nn}}=p_{\mathrm{w}}.\end{array}$$

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 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_{\mathrm{w}}=-{\mathit{\rho }}_{\mathrm{w}}g{z}_{\mathrm{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.

2.3 The free-surface equations

The boundaries Γ_s and Γ_b are time dependent and move according to two free-surface equations. The boundary Γ_bg follows 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

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial z_{\mathrm{s}}}{\partial t}}+u_{\mathrm{s}}{\displaystyle \frac{\partial z_{\mathrm{s}}}{\partial x}}-w_{\mathrm{s}}=a_{\mathrm{s}},\end{array}$$

where a_s denotes the surface mass balance and ${\mathit{u}}_{\mathrm{s}}=(u_{\mathrm{s}},w_{\mathrm{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

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial t}}+u_{\mathrm{b}}{\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial x}}-w_{\mathrm{b}}=a_{\mathrm{b}},\end{array}$$

where a_b is the basal mass balance and ${\mathit{u}}_{\mathrm{b}}=(u_{\mathrm{b}},w_{\mathrm{b}}{)}^T$ the velocity of the ice at Γ_bf. On Γ_bg, z_b=b(x); on Γ_bf, z_b>b(x).

The thickness of the ice is denoted by $H=z_{\mathrm{s}}-z_{\mathrm{b}}$ and depends on x and t.

2.4 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 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 τ₁₁ 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

$$\begin{array}{}\text{(9)}& H\mathit{\rho }=-z_{\mathrm{b}}{\mathit{\rho }}_{\mathrm{w}}\end{array}$$

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 Eq. (9) determines where the GL is in SSA in Docquier et al. (2011) and 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 ϵ. The aspect ratio of the ice ν is the quotient between a typical scale of the thickness of the ice ℋ and a horizontal length scale ℒ, $\mathit{\nu }=\mathcal{H}/\mathcal{L}$, and ϵ is ν times the quotient between the longitudinal and the shear stresses τ₁₁ and τ₁₂ in Eq. (2). If ${\mathit{\nu }}^{\mathrm{5}/\mathrm{2}}\ll \mathit{ϵ}\ll \mathrm{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

$$\begin{array}{}\text{(10)}& {\mathit{\tau }}_{\mathrm{22}}-p={\mathit{\sigma }}_{\mathrm{22}}=\mathit{\rho }g(z-z_{\mathrm{s}}).\end{array}$$

On floating ice ${\mathit{\tau }}_{\mathrm{22}}-p+p_{\mathrm{w}}=\mathrm{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 ${\mathit{\tau }}_{\mathrm{22}}-p+p_{\mathrm{w}}<\mathrm{0}$.

Introducing the notation

$$\begin{array}{}\text{(11)}& {\mathit{\chi }}_{\mathrm{a}}(x,z)={\mathit{\tau }}_{\mathrm{22}}-p+p_{\mathrm{w}}=\mathit{\rho }g(z-z_{\mathrm{s}}(x\left)\right)-{\mathit{\rho }}_{\mathrm{w}}g{z}_{\mathrm{b}}\left(x\right)\end{array}$$

and letting $H_{\mathrm{bw}}=-z_{\mathrm{b}}$ be the thickness of the ice below the sea level yields

$$\begin{array}{}\text{(12)}& {\mathit{\chi }}_{\mathrm{a}}(x,z_{\mathrm{b}})=-g(\mathit{\rho }H-{\mathit{\rho }}_{\mathrm{w}}{H}_{\mathrm{bw}}).\end{array}$$

If x<x_GL, then χ_a<0 in the neighborhood of x_GL on Γ_bg; 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 Wingham (2008), χ_a(x,z_b) varies linearly in x for x<x_GL.

In Sect. 4, we take this same approach but use an indicator χ(x) or $\stackrel{\mathrm{̃}}{\mathit{\chi }}\left(x\right)$ derived from the solutions of the nonlinear FS equations to estimate the GL position. These indicators are approximated by χ_a(x,z_b).

3 Discretization by FEM

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 Eqs. (7) and (8).

3.1 The weak form of the FS equations

We start by defining the mixed weak form of the FS equations. Introduce $k=\mathrm{1}+\mathrm{1}/n$, $k^{*}=\mathrm{1}+n$ with n from Glen's flow law and the spaces

$$\begin{array}{}\text{(13)}& {\mathbf{V}}_k=\mathit{\{}\mathit{v}:\mathit{v}\in (W^{\mathrm{1},k}\left(\mathrm{\Omega }\right){)}^{\mathrm{2}}\mathit{\}},Q_{k^{*}}=\mathit{\{}q:q\in L^{k^{*}}\left(\mathrm{\Omega }\right)\mathit{\}};\end{array}$$

see, for example, Jouvet and Rappaz (2011, Eq. 3.7), Chen et al. (2013, Sect. 3.1), and Martin and Monnier (2014, Eq. 21). The weak solution (u,p) of Eq. (1) is obtained as follows. Find $(\mathit{u},p)\in {\mathbf{V}}_k\times Q_{k^{*}}$ such that for all $(\mathit{v},q)\in {\mathbf{V}}_k\times Q_{k^{*}}$ the equation

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& A\left(\right(\mathit{u},p),(\mathit{v},q\left)\right)+B_{\mathrm{\Gamma }}(\mathit{u},\mathit{v},p)+B_{\mathcal{N}}(\mathit{u},\mathit{v},q)\\ & =F\left(\mathit{v}\right)+F_{\mathrm{\Gamma }}\left(\mathit{v}\right)\end{array}\end{array}$$

is satisfied, where

$$\begin{array}{rl}& A\left(\right(\mathit{u},p),(\mathit{v},q\left)\right)=\underset{\mathrm{\Omega }}{\int }\mathrm{2}\mathit{\eta }\left(\mathit{u}\right)\dot{\mathit{ϵ}}\left(\mathit{u}\right):\dot{\mathit{ϵ}}\left(\mathit{v}\right)\mathrm{d}\mathit{x}\\ & -b(\mathit{u},q)-b(\mathit{v},p),\\ & b(\mathit{u},q)=\underset{\mathrm{\Omega }}{\int }q\mathrm{\nabla }\cdot \mathit{u}\mathrm{d}\mathit{x},\\ & B_{\mathrm{\Gamma }}(\mathit{u},\mathit{v},p)=-\underset{{\mathrm{\Gamma }}_{\mathrm{bg}}}{\int }\left({\mathit{\sigma }}_{\mathrm{nn}}(\mathit{u},p)\mathit{n}\cdot \mathit{v}+{\mathit{\sigma }}_{\mathrm{nt}}(\mathit{u},p)\mathit{t}\cdot \mathit{v}\right),\\ & \mathrm{d}s=\underset{{\mathrm{\Gamma }}_{\mathrm{bg}}}{\int }\left(-{\mathit{\sigma }}_{\mathrm{nn}}(\mathit{u},p)\mathit{n}\cdot \mathit{v}+\mathit{\beta }(\mathit{t}\cdot \mathit{u})(\mathit{t}\cdot \mathit{v})\right)\mathrm{d}s,\\ & B_{\mathcal{N}}(\mathit{u},\mathit{v},q)=-\underset{{\mathrm{\Gamma }}_{\mathrm{bg}}}{\int }{\mathit{\sigma }}_{\mathrm{nn}}(\mathit{v},q)\mathit{n}\cdot \mathit{u}\mathrm{d}s\\ & +{\mathit{\gamma }}_{\mathrm{0}}\underset{{\mathrm{\Gamma }}_{\mathrm{bg}}}{\int }{\displaystyle \frac{\mathrm{1}}{h}}(\mathit{n}\cdot \mathit{u})(\mathit{n}\cdot \mathit{v})\mathrm{d}s,\\ & F\left(\mathit{v}\right)=\underset{\mathrm{\Omega }}{\int }\mathit{\rho }\mathit{g}\cdot \mathit{v}\mathrm{d}\mathit{x},\\ & F_{\mathrm{\Gamma }}\left(\mathit{v}\right)=-\underset{{\mathrm{\Gamma }}_{\mathrm{bf}}}{\int }{p}_{\mathrm{w}}\mathit{n}\cdot \mathit{v}\mathrm{d}s.\end{array}$$

The last term in B_𝒩 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 $\mathit{u}=u_{\mathrm{n}}\mathit{n}+u_{\mathrm{t}}\mathit{t}$, the contribution of the tangential force can also be written βu⋅v when u_n=0. The value of the positive parameter γ₀ 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 γ₀ is shown in Sect. 5. The first term in B_𝒩 symmetrizes the boundary term B_Γ+B_𝒩 on Γ_bg and vanishes when u_n=0. The boundary term F_Γ(v) is from the buoyancy force at the ice–ocean interface in EQ. (6) where p_w depends on z_b on Γ_bf.

3.2 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 on the ice base consists of N+1 nodes at ${\mathit{x}}_{\mathrm{i}}=(x_{\mathrm{i}},z_{\mathrm{b}}(x_{\mathrm{i}}\left)\right),i=\mathrm{0},\mathrm{\dots },N,$ with x coordinates x_i and a constant mesh size $\mathrm{\Delta }x=x_{\mathrm{i}}-x_{\mathrm{i}-\mathrm{1}}$.

In general, the GL is somewhere in the interior of an interval $[x_{\mathrm{i}-\mathrm{1}},x_{\mathrm{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_{\mathrm{GL}}\in [x_{\mathrm{i}-\mathrm{1}},x_{\mathrm{i}}]$ then the boundary integral over the interval can be split into two parts in Eq. (14) such that $(x,z_{\mathrm{b}}(x\left)\right)\in {\mathrm{\Gamma }}_{\mathrm{bg}}$ when $x\in [x_{\mathrm{i}-\mathrm{1}},x_{\mathrm{GL}}]$, and if $x\in [x_{\mathrm{GL}},x_{\mathrm{i}}]$ then $(x,z_{\mathrm{b}}(x\left)\right)\in {\mathrm{\Gamma }}_{\mathrm{bf}}$. In the GL element, we have

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{B}_{\mathrm{\Gamma }}+B_{\mathcal{N}}=& \underset{[x_{\mathrm{i}-\mathrm{1}},x_{\mathrm{GL}}]}{\int }-\left({\mathit{\sigma }}_{\mathrm{nn}}\right(\mathit{u},p)\mathit{n}\cdot \mathit{v}+{\mathit{\sigma }}_{\mathrm{nn}}(\mathit{v},q)\mathit{n}\cdot \mathit{u})\\ & +\mathit{\beta }(\mathit{t}\cdot \mathit{u})(\mathit{t}\cdot \mathit{v})+{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{0}}}{h}}(\mathit{n}\cdot \mathit{u})(\mathit{n}\cdot \mathit{v})\mathrm{d}s,\\ F_{\mathrm{\Gamma }}=& -\underset{[x_{\mathrm{GL}},x_{\mathrm{i}}]}{\int }{p}_{\mathrm{w}}\mathit{n}\cdot \mathit{v}\mathrm{d}s,\end{array}\end{array}$$

with the integration element  ds following Γ_b. There is a change of the boundary condition in the middle of the FEM element where the GL is located. With a strong formulation of the boundary condition u_n=0, the basis functions in V_k share this property, and the condition changes from the grounded node x_i−1 where the basis function satisfies u_n=0 to the floating node at x_i with a free u_n without taking the position of the GL inside $[x_{\mathrm{i}-\mathrm{1}},x_{\mathrm{i}}]$ into account. With the weak formulation in Nitsche's method, the standard basis functions we use do not satisfy u_n=0 strictly. The boundary condition is imposed on the solution by the additional penalty term multiplied by γ₀ in B_𝒩 in Eq. (14). A large γ₀ will force u_n to be small. The penalty term may change inside an element as in Eq. (15), where it is not equal to zero only in the grounded part.

The resulting system of nonlinear equations forms a nonlinear complementarity problem (Christensen et al., 1998). The distance d between the base of the ice and the bedrock at time t and at x is

$$\begin{array}{}\text{(16)}& d(x,t)=z_{\mathrm{b}}(x,t)-b\left(x\right)\ge \mathrm{0}.\end{array}$$

If d>0 on Γ_bf, then the ice is not in contact with the bedrock and ${\mathit{\sigma }}_{\mathrm{nn}}+p_{\mathrm{w}}=\mathrm{0}$; if ${\mathit{\sigma }}_{\mathrm{nn}}+p_{\mathrm{w}}<\mathrm{0}$ on Γ_bg, then the ice and the bedrock are in contact and d=0. Hence, the complementarity relation in the vertical direction is

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& d(x,t)\ge \mathrm{0},{\mathit{\sigma }}_{\mathrm{nn}}+p_{\mathrm{w}}\le \mathrm{0},\\ & d(x,t)({\mathit{\sigma }}_{\mathrm{nn}}+p_{\mathrm{w}})=\mathrm{0}\text{on}{\mathrm{\Gamma }}_{\mathrm{b}}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(18)}& \mathit{\beta }\ge \mathrm{0},\mathit{\beta }(x,t)d(x,t)=\mathrm{0}\end{array}$$

on Γ_b. In particular, these relations are valid at the nodes x=x_j, $j=\mathrm{0},\mathrm{1},\mathrm{\dots },N$.

The complementarity condition also holds for u_n and σ_nn such that

$$\begin{array}{}\text{(19)}& \begin{array}{l}{\mathit{\sigma }}_{\mathrm{nn}}+p_{\mathrm{w}}\le \mathrm{0},u_{\mathrm{n}}({\mathit{\sigma }}_{\mathrm{nn}}+p_{\mathrm{w}})=\mathrm{0}\text{on}{\mathrm{\Gamma }}_{\mathrm{b}},\end{array}\end{array}$$

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.

The nonlinear equation 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 d_j=d(x_j) 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.

3.3 Discretization of the advection equations

The advection equations for the moving ice boundary in Eqs. (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 Eqs. (7) and (8) are integrated in time by a semi-implicit method of first-order accuracy. Let c=s or b. Then the solution is advanced from time t^ℓ to $t^{\mathrm{\ell }+\mathrm{1}}=t^{\mathrm{\ell }}+\mathrm{\Delta }t$ with the time step Δt by

$$\begin{array}{}\text{(20)}& z_c^{\mathrm{\ell }+\mathrm{1}}=z_c^{\mathrm{\ell }}+\mathrm{\Delta }t\left(a_c^{\mathrm{\ell }}-u_c^{\mathrm{\ell }}{\displaystyle \frac{\partial z_c^{\mathrm{\ell }+\mathrm{1}}}{\partial x}}+w_c^{\mathrm{\ell }}\right).\end{array}$$

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_c^{\mathrm{\ell }+\mathrm{1}}$. This time discretization and its properties are discussed in Cheng et al. (2017) and summarized in Algorithm 2.

A numerical stability problem in z_b is encountered in the boundary condition at Γ_bf when the FS equations are solved in 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 ${\mathit{u}}_{\mathrm{b}}=\mathit{u}(x,z_{\mathrm{b}}(x\left)\right)$ is

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& {\mathit{u}}_{\mathrm{b}}=\left(\begin{array}{c}{u}_{\mathrm{b}}\\ w_{\mathrm{b}}\end{array}\right)=\left(\begin{array}{c}{z}_{\mathrm{b}x}\\ -\mathrm{1}\end{array}\right){\displaystyle \frac{u_{\mathrm{n}}}{\sqrt{\mathrm{1}+z_{\mathrm{b}x}^{\mathrm{2}}}}}\\ & +\left(\begin{array}{c}\mathrm{1}\\ z_{\mathrm{b}x}\end{array}\right){\displaystyle \frac{u_{\mathrm{t}}}{\sqrt{\mathrm{1}+z_{\mathrm{b}x}^{\mathrm{2}}}}},\end{array}\end{array}$$

where z_bx denotes $\partial z_{\mathrm{b}}/\partial x$. Inserting u_b and w_b from Eq. (21) into Eq. (8) yields

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial t}}=a_{\mathrm{b}}-u_{\mathrm{n}}\sqrt{\mathrm{1}+z_{\mathrm{b}x}^{\mathrm{2}}}.\end{array}$$

Instead of discretizing Eq. (22) explicitly at t^ℓ+1 with $u_{\mathrm{n}}^{\mathrm{\ell }}$ to determine $p_{\mathrm{w}}^{\mathrm{\ell }+\mathrm{1}}$, the base coordinate is updated implicitly,

$$\begin{array}{}\text{(23)}& z_{\mathrm{b}}^{\mathrm{\ell }+\mathrm{1}}=z_{\mathrm{b}}^{\mathrm{\ell }}+\mathrm{\Delta }t\left(a_{\mathrm{b}}^{\mathrm{\ell }+\mathrm{1}}-u_{\mathrm{n}}^{\mathrm{\ell }+\mathrm{1}}\sqrt{\mathrm{1}+(z_{\mathrm{b}x}^{\mathrm{\ell }+\mathrm{1}}{)}^{\mathrm{2}}}\right),\end{array}$$

in the evaluation of p_w in F_Γ(v) in Eq. (14).

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 width Δx and thickness H. Newton's law of motion yields

$$M{\dot{u}}_{\mathrm{n}}=Mg-\mathrm{\Delta }x{p}_{\mathrm{w}},$$

where $M=\mathrm{\Delta }x(z_{\mathrm{s}}-z_{\mathrm{b}})\mathit{\rho }$ is the mass of the slab. Dividing by M, integrating in time for u_n(t^m), letting $m=\mathrm{\ell }+\mathrm{1}$ or ℓ, and approximating the integral by the trapezoidal rule for the quadrature yields

$$\begin{array}{rl}& u_{\mathrm{n}}\left(t^m\right)=\underset{\mathrm{0}}{\overset{t^m}{\int }}g+{\displaystyle \frac{g{\mathit{\rho }}_{\mathrm{w}}}{\mathit{\rho }}}{\displaystyle \frac{z_{\mathrm{b}}}{z_{\mathrm{s}}-z_{\mathrm{b}}}}\mathrm{d}s\approx g{t}^m\\ & +{\displaystyle \frac{g{\mathit{\rho }}_{\mathrm{w}}}{\mathit{\rho }}}\sum _{i=\mathrm{0}}^m{\mathit{\alpha }}_i{\displaystyle \frac{z_{\mathrm{b}}^i}{z_{\mathrm{s}}^i-z_{\mathrm{b}}^i}}\mathrm{\Delta }t=u_{\mathrm{n}}^m,\end{array}$$

with the parameters

$${\mathit{\alpha }}_i=\mathrm{0.5},i=\mathrm{0},m,{\mathit{\alpha }}_i=\mathrm{1},i=\mathrm{1},\mathrm{\dots },m-\mathrm{1}.$$

Then insert $u_{\mathrm{n}}^m$ into Eq. (23). All terms in $u_{\mathrm{n}}^m$ from time steps i<m are collected in the sum ΔtF^m−1. Then Eq. (23) can be written

$$\begin{array}{}\text{(24)}& \begin{array}{rl}& z_{\mathrm{b}}^{\mathrm{\ell }+\mathrm{1}}=z_{\mathrm{b}}^{\mathrm{\ell }}-\mathrm{\Delta }{t}^{\mathrm{2}}{\displaystyle \frac{g{\mathit{\rho }}_{\mathrm{w}}}{\mathrm{2}\mathit{\rho }}}{\displaystyle \frac{z_{\mathrm{b}}^m}{z_{\mathrm{s}}^m-z_{\mathrm{b}}^m}}\\ & +\mathrm{\Delta }t\left(a_{\mathrm{b}}^{\mathrm{\ell }}-g{t}^m-\mathrm{\Delta }t{F}^{m-\mathrm{1}}\right).\end{array}\end{array}$$

For small changes in z_b in Eq. (24), the explicit method with m=ℓ is stable when Δt is so small that

$$\begin{array}{}\text{(25)}& |\mathrm{1}-\mathrm{\Delta }{t}^{\mathrm{2}}{\displaystyle \frac{g{\mathit{\rho }}_{\mathrm{w}}}{\mathrm{2}H\mathit{\rho }}}|\le \mathrm{1}.\end{array}$$

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=\mathrm{\ell }+\mathrm{1}$, the bound on Δt is

$$\begin{array}{}\text{(26)}& \mathrm{1}/|\mathrm{1}+\mathrm{\Delta }{t}^{\mathrm{2}}{\displaystyle \frac{g{\mathit{\rho }}_{\mathrm{w}}}{\mathrm{2}H\mathit{\rho }}}|\le \mathrm{1};\end{array}$$

i.e. there is no bound on positive Δt for stability, but accuracy will restrict Δt.

Much longer stable time steps are possible at the surface and the base of the ice with a semi-implicit method (Eq. 20) and a fully implicit method (Eq. 23) compared to an explicit method. For example, the time step for the problem in Eq. (20) with 1 km mesh size can be up to a couple of months. Therefore, we use the scheme in Eq. (20) for Eqs. (7) and (8) and the scheme in Eq. (23) for Eq. (22) and p_w as in Durand et al. (2009a). The difference between the approximations of z_b in Eqs. (20) and (23) is 𝒪(Δt²).

4 Subgrid scheme around the grounding line

The basic idea of the subgrid scheme for the FS equations in this paper follows the GL parameterization (SEP3) for SSA in Seroussi et al. (2014) and the analysis of FS in Schoof (2011). The GL is located at the position where the ice is on the ground and the flotation criterion is perfectly satisfied such that ${\mathit{\sigma }}_{\mathrm{nn}}=-p_{\mathrm{w}}$. In the FS equations, the hydrostatic assumption Eq. (9) may not be valid close to the GL. Therefore, the GL position can not be determined by simply checking the total thickness of the ice H against the depth below sea level H_bw. Instead, the flotation criterion is computed by comparing the water pressure with the numerical normal stress component orthogonal to the boundary inspired by the first-order analysis in Sect. 2.4.

The numerical solutions (e.g., Gagliardini et al., 2016; Gladstone et al., 2017) converge to the analytical solution of the FS PDE as the mesh size decreases. The analytical solution satisfies $z_{\mathrm{b}}(x,t)>b\left(x\right)$ with the boundary conditions in Eq. (6) at the base of the floating ice; where the ice is in contact with the bedrock $z_{\mathrm{b}}(x,t)=b\left(x\right)$, the boundary conditions are given by Eq. (5). Examples of the analytical solution are demonstrated by the thin light blue lines in Figs. 2 and 3 with a black “*” at the analytical GL position x_GL. The two figures share the same analytical solution. However, as illustrated in Figs. 2 and 3, the basal boundary of the ice z_b(x,t) does not conform with the mesh from the spatial discretization. In particular, the GL position x_GL of the analytical solution does not coincide with any of the nodes, but it usually stays on the bedrock b(x) between the last grounded (x_i−1) and the first floating (x_i) nodes; see Figs. 2 and 3. The linear element boundary between any x_j−1 and x_j is denoted by ℰ_j. The sequence of ${\mathcal{E}}_j,j=\mathrm{1},\mathrm{\dots },N$, approximates Γ_b. The grounding line element containing the GL is ℰ_i.

Figure 2 Schematic figure of the GL in case i, with the arrows indicating the direction of the net forces in the vertical direction. The light blue line is the analytical solution of the ice sheet with the analytical GL position x_GL. The red line is the grounded boundary Γ_bg, the dark blue line is the floating boundary Γ_bf, and the brown line is the bedrock topography b(x). (a) The last grounded and first floating nodes as defined in Elmer/ICE. (b) Linear interpolation to approximate the numerical GL position ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$. (c) The step functions ℋ_𝒩(x) and ℋ_β(x) indicate the area for Nitsche's penalty and slip boundary conditions.

Figure 3 Schematic figure of the GL in case ii, with the arrows indicating the direction of the net forces in the vertical direction. The colors of the lines follow those in Fig. 2. (a) The last grounded and first floating nodes as defined in Elmer/ICE. The node x_i is fully geometrically floating, and the net force is 0. (b) Linear interpolation to approximate the numerical GL position ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$. The point ${\widehat{\mathit{x}}}_{\mathrm{i}}$ on the bedrock has the same x coordinate as x_i. (c) The step functions ℋ_𝒩(x) and ℋ_β(x) indicate the area for Nitsche's penalty and slip boundary conditions.

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. 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 $\mathit{\chi }\left({\mathit{x}}_{\mathrm{i}}\right)={\mathit{\sigma }}_{\mathrm{nn}}\left({\mathit{x}}_{\mathrm{i}}\right)+p_{\mathrm{w}}\left({\mathit{x}}_{\mathrm{i}}\right)>\mathrm{0}$, whereas in case ii the floating condition ${\mathit{\sigma }}_{\mathrm{nn}}\left({\mathit{x}}_{\mathrm{i}}\right)+p_{\mathrm{w}}\left({\mathit{x}}_{\mathrm{i}}\right)=\mathrm{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 Figs. 2a and 3a. 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 as geometrically grounded (defined as in Algorithm 1), shown with in red in Fig. 2a. In case ii, the GL element is treated as geometrically floating and colored in blue in Fig. 3a.

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 of 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 (<100 m). The schematic drawing of the subgrid scheme for the two cases is shown in the panels b of Figs. 2 and 3. The GL element is divided into the grounded (red) and floating (blue) parts by the estimated GL position ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$ on ℰ_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 as

$$\begin{array}{}\text{(27)}& \mathit{\chi }\left(\mathit{x}\right)={\mathit{\sigma }}_{\mathrm{nn}}+p_{\mathrm{w}},\end{array}$$

which vanishes on the floating ice and is negative and approximately equal to ${\mathit{\chi }}_{\mathrm{a}}={\mathit{\tau }}_{\mathrm{22}}-p+p_{\mathrm{w}}$ in Eq. (11) on the ground since the slope of the bedrock is small and $\mathit{n}\approx (\mathrm{0},-\mathrm{1}{)}^T$. Because of the poor spatial resolution of the coarse mesh, χ(x_i) is positive.

To determine the position ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$, we solve $\mathit{\chi }\left({\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}\right)={\mathit{\sigma }}_{\mathrm{nn}}\left({\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}\right)+p_{\mathrm{w}}\left({\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}\right)=\mathrm{0}$ by linear interpolation between χ(x_i−1) and χ(x_i) such that

$$\begin{array}{}\text{(28)}& {\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}={\mathit{x}}_{i-\mathrm{1}}-{\displaystyle \frac{\mathit{\chi }\left({\mathit{x}}_{i-\mathrm{1}}\right)}{\mathit{\chi }\left({\mathit{x}}_{i-\mathrm{1}}\right)-\mathit{\chi }\left({\mathit{x}}_i\right)}}({\mathit{x}}_{i-\mathrm{1}}-{\mathit{x}}_{\mathrm{i}}).\end{array}$$

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 (Gagliardini et al., 2013). Hence, it makes sense to approximate the analytical GL position x_GL by ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$ by linear interpolation in the current framework. This approach fits well with case i since the indicator χ(x) has opposite signs at x_i−1 and x_i; see Fig. 2b, where ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$ is marked with a red “*”. It guarantees the existence and uniqueness of ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{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 $\mathit{\chi }\left({\mathit{x}}_{\mathrm{i}}\right)={\mathit{\sigma }}_{\mathrm{nn}}\left({\mathit{x}}_i\right)+p_{\mathrm{w}}\left({\mathit{x}}_i\right)=\mathrm{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. The result is that p_w defined by the implicit z_b in Eq. (23) satisfies ${\mathit{\sigma }}_{\mathrm{nn}}+p_{\mathrm{w}}>\mathrm{0}$ in Eq. (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 that situation since z_b(x_i)>b(x_i). In this circumstance, χ(x_i)=0 in the floating node and a correction of χ(x) is introduced into case ii by $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ in

$$\begin{array}{}\text{(29)}& \stackrel{\mathrm{̃}}{\mathit{\chi }}\left(\mathit{x}\right)={\mathit{\sigma }}_{\mathrm{nn}}\left(\mathit{x}\right)+p_{\mathrm{b}}\left(\mathit{x}\right).\end{array}$$

Here $p_{\mathrm{b}}\left(\mathit{x}\right)=-{\mathit{\rho }}_{\mathrm{w}}gb\left(x\right)$ 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, $\stackrel{\mathrm{̃}}{\mathit{\chi }}\left(\mathit{x}\right)\ge \mathit{\chi }\left(\mathit{x}\right)$. Notice that $p_{\mathrm{b}}\left({\mathit{x}}_{\mathrm{i}}\right)=p_{\mathrm{w}}\left({\widehat{\mathit{x}}}_{\mathrm{i}}\right)>p_{\mathrm{w}}\left({\mathit{x}}_{\mathrm{i}}\right)$, where ${\widehat{\mathit{x}}}_{\mathrm{i}}$ is a point on the bedrock with the same x coordinate of x_i, as illustrated in Fig. 3b. Both χ(x) in Eq. (27) and $\stackrel{\mathrm{̃}}{\mathit{\chi }}\left(\mathit{x}\right)$ in Eq. (29) are nonlinear in x, but the numerical approximation of them will vary linearly in x. A solution ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$ is found by linear interpolation of $\stackrel{\mathrm{̃}}{\mathit{\chi }}\left(\mathit{x}\right)$ between the nodes x_i−1 and x_i as in Eq. (28). It follows from Eq. (28) that ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{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 ${\mathit{\sigma }}_{\mathrm{nn}}\left({\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}\right)$ to approximate σ_nn(x_GL) on a virtual element between x_i−1 and ${\widehat{\mathit{x}}}_{\mathrm{i}}$; see Fig. 3. The position ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$ is a numerical approximation of the analytical GL position, although it is not geometrically in contact with the bedrock.

Since we have p_b(x)=p_w(x) and $\mathit{\chi }\left(\mathit{x}\right)=\stackrel{\mathrm{̃}}{\mathit{\chi }}\left(\mathit{x}\right)$ at the GL element in case i, we can simply use $\stackrel{\mathrm{̃}}{\mathit{\chi }}\left(\mathit{x}\right)$ to find ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$ for the two cases by replacing χ in Eq. (28) by $\stackrel{\mathrm{̃}}{\mathit{\chi }}$.

The domains Γ_bg and Γ_bf are separated at ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{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, ℋ_𝒩(x) and ℋ_β(x), to include and exclude quadrature points in the integration of Nitsche's term and the slip boundary condition, respectively. They are defined for case i in Fig. 2 and for case ii in Fig. 3. To achieve a reasonable numerical accuracy within the GL element, as suggested in Seroussi et al. (2014), at least 10th-order Gaussian quadrature is used.

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. In contrast, in case ii the GL element ℰ_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 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 ℋ_𝒩(x) indicates how Nitsche's method is implemented on the basal elements; see Figs. 2c and 3c for the two cases. The penalty term contributes to the integration only when ℋ_𝒩(x)=1.

The slip coefficient β is treated similarly with the step function ℋ_β(x), where ℋ_β(x)=1 is on the fully geometrically grounded elements and ℋ_β(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 and ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{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 ℰ_i. For the floating part between ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$ and x_i in case ii, there is no friction, ℋ_β(x)=0, we have reduced friction between x_i−1 and ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$; see Fig. 3c. The boundary integral Eq. (15) on ℰ_i is now rewritten with the two step functions as

$$\begin{array}{}\text{(30)}& \begin{array}{rl}& B_{\mathrm{\Gamma }}+B_{\mathcal{N}}=\underset{{\mathcal{E}}_i}{\int }-{\mathcal{H}}_{\mathcal{N}}\left({\mathit{\sigma }}_{\mathrm{nn}}\right(\mathit{u},p)\mathit{n}\cdot \mathit{v}+{\mathit{\sigma }}_{\mathrm{nn}}(\mathit{v},q)\mathit{n}\cdot \mathit{u})\\ & +{\mathcal{H}}_{\mathit{\beta }}\mathit{\beta }(\mathit{t}\cdot \mathit{u})(\mathit{t}\cdot \mathit{v})+{\mathcal{H}}_{\mathcal{N}}{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{0}}}{h}}(\mathit{n}\cdot \mathit{u})(\mathit{n}\cdot \mathit{v})\mathrm{d}s,\\ & F_{\mathrm{\Gamma }}=\underset{{\mathcal{E}}_i}{\int }(\mathrm{1}-{\mathcal{H}}_{\mathcal{N}})p_{\mathrm{w}}\mathit{n}\cdot \mathit{v}\mathrm{d}s.\end{array}\end{array}$$

A summary of the numerical treatment of the GL is as follows:

• advance phase ⇒ indicator χ in Eq. (27), case i or case ii;

• retreat phase ⇒ indicator $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ in Eq. (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 as follows:

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) and Gagliardini et al. (2013, 2016). The detailed procedure is explained in Algorithms 1, 2, and 3. The solution to the nonlinear FS system is computed with Picard iterations to a 10⁻⁵ relative error with a limit of maximal 25 nonlinear iterations. The ${\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{GL}}$ position is determined dynamically during each fixed-point iteration by solving Eq. (28) with χ or $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ and the solution σ_nn(x) from the previous nonlinear iteration, and the step functions ℋ_𝒩 and ℋ_β are adjusted accordingly. The water pressure p_b is fixed since the ice geometry is not changed during the nonlinear iterations.

5 Results

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, 2, 1, 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.

The dependence on γ₀ in Eq. (30) for the retreating ice is shown in Fig. 4 with γ₀ between 10⁴ and 10⁹. The estimated GL positions do not vary with different choices of γ₀ from 10⁵ to 10⁸, which suggests a suitable range of γ₀. If γ₀ is too small (γ₀≪10⁴), oscillations appear in the estimated GL positions. If γ₀ is too large (γ₀≫10⁸), then more nonlinear iterations in Algorithm 1 are needed in each time step. The same dependency of γ₀ is observed for the advancing experiments and for different mesh resolutions as well. The results are not very sensitive to γ₀ and for the remaining experiments we choose γ₀=10⁶.

Figure 4 The MISMIP 3a retreat experiment with Δx=1 km for different choices of γ₀ in the time interval [0,10 000] years.

The GL position during the transient simulations in the advance and retreat phases is displayed in Fig. 5, and the steady-state results (at t=10 000) are shown in Fig. 6 for different mesh resolutions. The range of the steady-state solutions from Gagliardini et al. (2016) with mesh resolution from 25 to 200 m are shown as background shaded regions in red. We achieve similar GL migration results for both the advance and retreat experiments with at least 20-times-larger mesh resolutions. The GL position is insensitive to the variation in mesh size between 0.5 and 4 km.

The distance between the steady-state GL positions of the retreat and the advance phases is shown in Fig. 6b. 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.

Figure 5 The MISMIP 3a experiments for the GL position when $t\in [\mathrm{0},\mathrm{10}\mathrm{000}]$ with Δx=4, 2, 1, and 0.5 km for the advance (solid) and retreat (dashed) phases.

Figure 6 The MISMIP 3a experiments at the final time t=10 000 with the resolutions at Δx=4, 2, 1, and 0.5 km. (a) The GL positions in the advance (★) and retreat (•) phases. (b) The distance between the retreat and the advance x_GL at the steady states. The shaded regions indicate the range of the results in Gagliardini et al. (2016) with 20-times-smaller mesh resolutions from 25 to 200 m, with the axis scale shown in red at the top of the plot.

We observe oscillations at the ice surface near the GL in all the experiments as expected from Durand et al. (2009a) and Schoof (2011). A zoomed-in plot of the surface elevation computed with four different mesh sizes (Δx=0.5, 1, 2, and 4 km) after an advance simulation to t=10 000 years is found in Fig. 7a. 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.

Figure 7 Details of the solutions for the advance experiment with $\mathrm{\Delta }x=\mathrm{0.5},\mathrm{1},\mathrm{2},\mathrm{4}$ km after 10 000 years. The solid dots represent the nodes of the elements, and the vertical, dashed red lines indicate the GL position. (a) The oscillations in z_s at the ice surface near GL. (b) The flotation criterion is evaluated by H_bw∕H. The ratio ρ∕ρ_w is drawn in a horizontal, dash-dotted purple line.

The ratio between the thickness below sea level H_bw and the ice thickness H is shown in Fig. 7b. As in Fig. 7a, the ratio is plotted versus the distance from the GL achieved with the particular mesh size. The horizontal, dash-dotted purple 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 panel b confirms that the hydrostatic assumption Hρ=H_bwρ_w in Eq. (9) is not valid in the FS equations for x>x_GL close to the GL and at the GL position (cf. Durand et al., 2009a; Schoof, 2011). For x<x_GL we have $H_{\mathrm{bw}}/H<\mathit{\rho }/{\mathit{\rho }}_{\mathrm{w}}$ since H_bw decreases and H increases. The conclusion from numerical experiments in van Dongen et al. (2018) is that the hydrostatic assumption and the SSA equations approximate the FS equations well for the floating ice beginning at a short distance away from the GL.

The surface and the base velocity solutions from the advance experiment are displayed in Fig. 8 with Δx=0.5, 1, 2, and 4 km after 10 000 years. The horizontal velocities on the two surfaces are on top of each other for all Δx with negligibly small differences on the floating ice as expected. The vertical velocities w on the surface (dotted curves) and the base (solid curves) 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.

Figure 8 The velocities u (a) and w (b) of the ice in the advance experiment with Δx=0.5, 1, 2, and 4 km after 10 000 years. The solutions at the upper surface are the solid curves, and the solutions at the lower surface are the dotted curves. The vertical velocity w is zoomed in close to the GL, with the distance to the mesh-dependent GL on the x axis.

6 Discussion

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. The friction 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 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 Figs. 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. 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, and 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.

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 $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ in Eq. (29) is applied on the edges of each triangle 𝒯 in the mesh. If χ<0 in all three nodes, then 𝒯 is grounded. If χ≥0 in all nodes, then 𝒯 is floating. The GL passes inside 𝒯 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 $\stackrel{\mathrm{̃}}{\mathit{\chi }}$. The FEM approximation is modified 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, for example, Gladstone et al. (2010a), Cornford et al. (2013), and 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 paleo-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 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 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.

7 Conclusions

A subgrid scheme at the GL was 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 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 $\stackrel{\mathrm{̃}}{\mathit{\chi }}\left(\mathit{x}\right)$. Only in that element is the FEM discretization modified to accommodate the discontinuities in the boundary conditions. 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.