A new vertically integrated MOno-Layer Higher-Order (MOLHO) ice flow model

Numerical simulations of ice sheets rely on the momentum balance to determine how ice velocities change as the geometry of the system evolves. Ice is generally assumed to follow a Stokes flow with a nonlinear viscosity. Several approximations have been proposed in order to lower the computational cost of a full-Stokes stress balance. A popular option is the Blatter–Pattyn or higher-order model (HO), which consists of a three-dimensional set of equations that solves the horizontal velocities only. However, it still remains computationally expensive for long transient simulations. Here we present a depth-integrated formulation of the HO model, which can be solved on a two-dimensional mesh in the horizontal plane. We employ a specific polynomial function to describe the vertical variation in the velocity, which allows us to integrate the vertical dimension using a semi-analytic integration. We assess the performance of this MOno-Layer Higher-Order (MOLHO) model to compute ice velocities and simulate grounding line dynamics on standard benchmarks (ISMIP-HOM and MISMIP3D). We compare MOLHO results to the ones obtained with the original three-dimensional HO model. We also compare the time performance of both models in time-dependent runs. Our results show that the ice velocities and grounding line positions obtained with MOLHO are in very good agreement with the ones from HO. In terms of computing time, MOLHO requires less than 10 % of the computational time of a typical HO model, for the same simulations. These results suggest that the MOno-Layer Higher-Order formulation provides improved computational time performance and a comparable accuracy compared to the HO formulation, which opens the door to higher-order paleo simulations.

1 Introduction

Projecting the future evolution of the ice sheets of Greenland and Antarctica and their potential contribution to sea level rise often relies on computer simulations carried out by numerical ice sheet models (e.g., Aschwanden et al., 2019; Goelzer et al., 2020; Seroussi et al., 2020; Edwards et al., 2021). These ice sheet models solve a set of flow equations based on the conservation of momentum to obtain the ice velocity field over the entire ice sheet. There exist several ice flow models in the field of ice sheet/glacier modeling. One of these flow models is the full-Stokes equations (FS), a set of equations that solve the three-dimensional ice velocity as well as the ice pressure. The FS model is applicable to a wide range of systems: from small glaciers to continent-sized ice sheets. However, this model is computationally demanding, especially in time-dependent numerical simulations, which restricts its use to short-term projections or regional applications.

Over the last decades, simplified ice flow models were developed in order to decrease the computational cost of the momentum balance (e.g, Hutter, 1983; Morland, 1987; MacAyeal, 1989; Blatter, 1995; Pattyn, 2003; Goldberg, 2011). All of them are derived from the FS equations and rely on different approximations in the stress tensor or velocity field. A commonly applied approximation is the Blatter–Pattyn model (Blatter, 1995; Pattyn, 2003), also known as a higher-order model (HO). The HO model assumes a hydrostatic pressure and negligible contribution of horizontal gradient of vertical velocities to the velocity computation. The model is three-dimensional and solves the horizontal velocities over the entire ice sheet. However, for long transient runs and/or sensitivity analyses that require running the model for large ensembles, the HO model still demands relatively high computational resources. Recently, numerical schemes were applied to improve the computational efficiency of the HO model without compromising numerical accuracy (e.g., Langdon and Raymond, 1978; Bassis, 2010; Brinkerhoff and Johnson, 2015; Cuzzone et al., 2018; Shapero et al., 2021). These schemes rely on the finite element method (FEM), which allows the employment of polynomial functions of a degree equal to or greater than 2 to model the vertical variation in the horizontal velocities. These approaches decrease the number of vertical elements (layers) in the mesh, reducing the overall computational cost. A natural extension of such schemes is to employ a single layer of three-dimensional elements (e.g., triangular prisms) with a higher-order polynomial along the vertical axis (Brinkerhoff and Johnson, 2015; Shapero et al., 2021). This mono-layer, higher-order flow model presents a reasonable numerical accuracy in comparison to the original three-dimensional HO model, at a relatively lower computational cost (Brinkerhoff and Johnson, 2015). It has been used for simulations of glaciers and ice sheets (Brinkerhoff and Johnson, 2015; Brinkerhoff et al., 2017; Shapero et al., 2021). However, some technical aspects, such as the vertical integration of the viscosity, have not been fully resolved and the performance of this mono-layer higher-order model in marine ice sheet simulations has not been tested.

Here we present a finite element formulation of a new vertically integrated MOno-Layer Higher-Order (MOLHO) ice flow model inspired by the scheme employed by Brinkerhoff and Johnson (2015). Previous works have employed numerical integration (e.g., Gauss–Legendre quadrature) over the vertical dimension in triangles (Brinkerhoff and Johnson, 2015) or triangular prisms (Shapero et al., 2021) to account for the vertical shear and ice viscosity. Here, we pre-compute a depth-averaged ice viscosity, which allows us to perform an analytical vertical integration. Thus, the new formulation is written over the horizontal plane and relies on a two-dimensional mesh. The vertical variation in the horizontal velocities is described by a specific higher-order polynomial. This specific polynomial adds 2 degrees of freedom per node and allows us to integrate over the vertical axis. The vertical integration of the ice viscosity is evaluated by a semi-analytic scheme of low computational cost. The formulation is implemented in the Ice-sheet and Sea-level System Model (ISSM) v4.18 (Larour et al., 2012, 2020). We first run the ISMIP-HOM diagnostic experiments (Pattyn et al., 2008) to evaluate the velocity field for different wavelengths of bedrock elevation and compare the results to the solutions obtained using the three-dimensional HO model. We then run the MISMIP3D benchmark (Pattyn et al., 2013) to evaluate the behavior of grounding line dynamics in a marine ice sheet setup. With the MISMIP3D experiments we assess the grounding line position at steady state for different mesh resolutions as well as grounding line reversibility after imposed friction perturbations. We also run a time performance analysis to compare the simulation time spent by each flow model, MOLHO and HO, and conclude on the applicability of this model in the future.

The paper is organized as follows. We first present the description of the MOLHO formulation (Sect. 2), highlighting the main steps necessary to perform the vertical integration of the weak formulation (e.g, Sect. 2.3), especially the integration of the ice viscosity (Sect. 2.4). Then, we describe the numerical experiments performed to evaluate the new formulation (Sect. 3), followed by the results (Sect. 4). We finish the paper with a discussion (Sect. 5) and a summary of the conclusions of this work (Sect. 6). Additional information regarding the pre-computation of the ice viscosity as well as the resulting element stiffness matrix and load vector is provided in the Appendix sections.

2 Description of the MOno-Layer Higher-Order (MOLHO) model

2.1 Polynomial function for vertical shear

The specific polynomial employed to describe the vertical variation in the horizontal velocities is based on the shallow ice approximation (SIA, Hutter, 1983). In the SIA, the horizontal velocities ($\mathit{v}=[v_x,v_y{]}^{\mathit{⊺}}$) of an isothermal ice sheet following Glen's law are given by

$$\begin{array}{}\text{(1)}& \mathit{v}={\mathit{v}}^{\mathrm{b}}-\mathrm{2}(\mathit{\rho }g{)}^n\Vert \mathrm{\nabla }s{\Vert }^{n-\mathrm{1}}{\displaystyle \frac{A}{n+\mathrm{1}}}{H}^{n+\mathrm{1}}\left(\mathrm{1}-{\left({\displaystyle \frac{s-z}{H}}\right)}^{n+\mathrm{1}}\right)\mathrm{\nabla }s,\end{array}$$

where v^b ($=[v_x^{\mathrm{b}},v_y^{\mathrm{b}}{]}^{\mathit{⊺}}$) represent the basal horizontal velocities, H is the ice thickness, s is the ice surface elevation, A is the ice rate factor, n is the Glen law exponent, ρ is the ice density, and g is the gravitational acceleration. We write Eq. (1) as

$$\begin{array}{}\text{(2)}& \mathit{v}(x,y,z)={\mathit{v}}^{\mathrm{b}}(x,y)+{\mathit{v}}^{\mathrm{sh}}(x,y)\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right),\end{array}$$

where ${\mathit{v}}^{\mathrm{sh}}=[v_x^{\mathrm{sh}},v_y^{\mathrm{sh}}{]}^{\mathit{⊺}}$ is the internal deformation contribution to the surface velocities, and ζ is an auxiliary variable defined as

$$\begin{array}{}\text{(3)}& \mathit{\zeta }={\displaystyle \frac{s-z}{H}}=\mathrm{1}-{\displaystyle \frac{z-b}{H}},\end{array}$$

with b being the ice base elevation (note that $H=s-b$).

Equation (2) has been widely employed by the glaciological community in ice sheet models to set up boundary conditions (e.g., Raymond, 1983) and describe vertical deformation within the glacier body (e.g., Langdon and Raymond, 1978; Bueler and Brown, 2009; Bassis, 2010; Brinkerhoff and Johnson, 2015; Shapero et al., 2021).

In the SIA, v^sh values are defined using the local surface gradient only (Eq. 1). In the MOno-Layer Higher-Order (MOLHO) model, we employ Eq. (2) in a higher-order (HO) weak formulation such that v^sh values are evaluated only after solving a nonlinear system resulting from the HO stress-balance equations. Details of the computation of v^sh are shown in Sect. 2.3.

To derive the MOLHO formulation, we use some identities related to vertical derivatives and vertical integration of the function $\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}$. Appendix B presents some useful identities employed in this work. Note that this polynomial function is derived for an isothermal ice sheet, which could be a source of differences in real applications when compared to a three-dimensional HO model. Nevertheless, one can modify the function to account for the vertical variation in ice temperature.

2.2 Weak formulation

The three-dimensional Blatter–Pattyn or higher-order (HO) ice flow model is defined as (Blatter, 1995; Pattyn, 2003)

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial x}}& {\displaystyle }\left(\mathrm{4}\mathit{\mu }{\displaystyle \frac{\partial v_x}{\partial x}}+\mathrm{2}\mathit{\mu }{\displaystyle \frac{\partial v_y}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(\mathit{\mu }\left({\displaystyle \frac{\partial v_x}{\partial y}}+{\displaystyle \frac{\partial v_y}{\partial x}}\right)\right)\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\mu }{\displaystyle \frac{\partial v_x}{\partial z}}\right)=\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial x}},\end{array}$$

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial x}}& \left(\mathit{\mu }\left({\displaystyle \frac{\partial v_x}{\partial y}}+{\displaystyle \frac{\partial v_y}{\partial x}}\right)\right)+{\displaystyle \frac{\partial }{\partial y}}\left(\mathrm{4}\mathit{\mu }{\displaystyle \frac{\partial v_y}{\partial y}}+\mathrm{2}\mathit{\mu }{\displaystyle \frac{\partial v_x}{\partial x}}\right)\\ & +{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\mu }{\displaystyle \frac{\partial v_y}{\partial z}}\right)=\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial y}}.\end{array}\end{array}$$

For the boundary conditions, we assume a negligible stress at the ice surface (i.e., τ^s=0) and a viscous friction law at the ice base defined as

$$\begin{array}{}\text{(5)}& {\mathit{\tau }}^{\mathrm{b}}=-{\mathit{\alpha }}^{\mathrm{2}}{\mathit{v}}^{\mathrm{b}},\end{array}$$

where ${\mathit{v}}^{\mathrm{b}}=[v_x^{\mathrm{b}},v_y^{\mathrm{b}}{]}^{\mathit{⊺}}$ is the horizontal velocity at the glacier base and $\mathit{\alpha }=\mathit{\alpha }(x,y)$ is the friction coefficient that, in general, is a (nonlinear) function of the basal velocities v^b. Note that other friction laws could be employed in this formulation. At the ice–ocean interface (i.e., calving front), a Neumann boundary condition based on ocean water pressure is applied.

Let 𝒱 be the space of kinematically admissible fields that satisfy the Dirichlet boundary conditions and whose first derivatives are square integrable on the glacier domain. The weak formulation (assuming non-homogeneous Dirichlet conditions on the model boundary and viscous sliding at the base) of HO is

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\underset{\mathrm{\Omega }}{\int }& \left(\mathrm{4}\mathit{\mu }{\displaystyle \frac{\partial v_x}{\partial x}}+\mathrm{2}\mathit{\mu }{\displaystyle \frac{\partial v_y}{\partial y}}\right){\displaystyle \frac{\partial {\mathit{\vartheta }}_x}{\partial x}}+\mathit{\mu }\left({\displaystyle \frac{\partial v_x}{\partial y}}+{\displaystyle \frac{\partial v_y}{\partial x}}\right){\displaystyle \frac{\partial {\mathit{\vartheta }}_x}{\partial y}}\\ & +\mathit{\mu }{\displaystyle \frac{\partial v_x}{\partial z}}{\displaystyle \frac{\partial {\mathit{\vartheta }}_x}{\partial z}}\mathrm{d}\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }\mathit{\mu }\left({\displaystyle \frac{\partial v_x}{\partial y}}+{\displaystyle \frac{\partial v_y}{\partial x}}\right){\displaystyle \frac{\partial {\mathit{\vartheta }}_y}{\partial x}}\\ & +\left(\mathrm{4}\mathit{\mu }{\displaystyle \frac{\partial v_y}{\partial y}}+\mathrm{2}\mathit{\mu }{\displaystyle \frac{\partial v_x}{\partial x}}\right){\displaystyle \frac{\partial {\mathit{\vartheta }}_y}{\partial y}}+\mathit{\mu }{\displaystyle \frac{\partial v_y}{\partial z}}{\displaystyle \frac{\partial {\mathit{\vartheta }}_y}{\partial z}}\mathrm{d}\mathrm{\Omega }\\ & +\underset{{\mathrm{\Gamma }}^{\mathrm{b}}}{\int }{\mathit{\alpha }}^{\mathrm{2}}{v}_x^{\mathrm{b}}{\mathit{\vartheta }}_x\mathrm{d}{\mathrm{\Gamma }}^{\mathrm{b}}+\underset{{\mathrm{\Gamma }}^{\mathrm{b}}}{\int }{\mathit{\alpha }}^{\mathrm{2}}{v}_y^{\mathrm{b}}{\mathit{\vartheta }}_y\mathrm{d}{\mathrm{\Gamma }}^{\mathrm{b}}\\ & =-\underset{\mathrm{\Omega }}{\int }\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial x}}{\mathit{\vartheta }}_x\mathrm{d}\mathrm{\Omega }-\underset{\mathrm{\Omega }}{\int }\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial y}}{\mathit{\vartheta }}_y\mathrm{d}\mathrm{\Omega }\\ & +\underset{{\mathrm{\Gamma }}^{\mathrm{w}}}{\int }\left[\left(\mathrm{4}\mathit{\mu }{\displaystyle \frac{\partial v_x}{\partial x}}+\mathrm{2}\mathit{\mu }{\displaystyle \frac{\partial v_y}{\partial y}}\right)n_x\right.\\ & \left.+\mathit{\mu }\left({\displaystyle \frac{\partial v_x}{\partial y}}+{\displaystyle \frac{\partial v_y}{\partial x}}\right)n_y\right]{\mathit{\vartheta }}_x\mathrm{d}{\mathrm{\Gamma }}^{\mathrm{w}}\\ & +\underset{{\mathrm{\Gamma }}^{\mathrm{w}}}{\int }\left[\left(\mathrm{4}\mathit{\mu }{\displaystyle \frac{\partial v_y}{\partial y}}+\mathrm{2}\mathit{\mu }{\displaystyle \frac{\partial v_x}{\partial x}}\right)n_y\right.\\ & \left.+\mathit{\mu }\left({\displaystyle \frac{\partial v_x}{\partial y}}+{\displaystyle \frac{\partial v_y}{\partial x}}\right)n_x\right]{\mathit{\vartheta }}_y\mathrm{d}{\mathrm{\Gamma }}^{\mathrm{w}},\\ & v_x=v_x^{\mathrm{D}},v_y=v_y^{\mathrm{D}}\mathrm{on}{\mathrm{\Gamma }}^{\mathrm{D}},\\ & \forall {\left[{\mathit{\vartheta }}_x,{\mathit{\vartheta }}_y\right]}^{\mathit{⊺}}\in \mathcal{V},\end{array}\end{array}$$

where Ω is the three-dimensional domain of the glacier, Γ^b is the ice base, Γ^w is the vertical face of the glacier boundary (e.g., calving front), n_x and n_y are the components of the unit vector pointing outward edge Γ^w, and Γ^D is the glacier boundary where Dirichlet boundary conditions with values $v_x^{\mathrm{D}},v_y^{\mathrm{D}}$ are imposed.

2.3 Finite element discretization

Based on the SIA, we decompose the horizontal velocities in MOLHO following Eq. (2):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{v}_x(x,y,z)& =v_x^{\mathrm{b}}(x,y)+v_x^{\mathrm{sh}}(x,y)\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right),\\ v_y(x,y,z)& =v_y^{\mathrm{b}}(x,y)+v_y^{\mathrm{sh}}(x,y)\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right).\end{array}\end{array}$$

We note that the decomposition employed here differs from the approach adopted in previous work (Brinkerhoff and Johnson, 2015). In the latter, the horizontal velocities were parameterized as the sum of depth-averaged velocities and the deviation from it (rather than basal and shear velocities) such that the first component of the stress balance solution (i.e., depth-averaged velocities) could be used directly in the continuity equation that controls the advection of the ice thickness. Here, once v^b and v^sh are computed, the vertically averaged velocities, $\stackrel{\mathrm{‾}}{\mathit{v}}$, are obtained by

$$\begin{array}{}\text{(8)}& \stackrel{\mathrm{‾}}{\mathit{v}}(x,y)={\mathit{v}}^{\mathrm{b}}(x,y)+{\mathit{v}}^{\mathrm{sh}}(x,y){\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}.\end{array}$$

To apply the finite element discretization to Eq. (6), we approximate the horizontal velocities by employing two-dimensional basis functions ϕ_j as follows:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{v}_x(x,y,z)& \approx \sum _{j=\mathrm{1}}^{n_{\mathrm{f}}}{v}_{x,j}^{\mathrm{b}}{\mathit{\varphi }}_j(x,y)+\sum _{j=\mathrm{1}}^{n_{\mathrm{f}}}{v}_{x,j}^{\mathrm{sh}}{\mathit{\varphi }}_j(x,y)\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right),\\ v_y(x,y,z)& \approx \sum _{j=\mathrm{1}}^{n_{\mathrm{f}}}{v}_{y,j}^{\mathrm{b}}{\mathit{\varphi }}_j(x,y)+\sum _{j=\mathrm{1}}^{n_{\mathrm{f}}}{v}_{y,j}^{\mathrm{sh}}{\mathit{\varphi }}_j(x,y)\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right),\end{array}\end{array}$$

with $v_{x,j}^{\mathrm{b}},v_{y,j}^{\mathrm{b}}$ and $v_{x,j}^{\mathrm{sh}},v_{y,j}^{\mathrm{sh}}$ being the values of the base velocities and surface shear velocities evaluated at nodal points j. In Eq. (9), n_f is the total number of basis functions in the entire mesh. By decomposing the velocities into two independent components, we can define the basis functions on a horizontal two-dimensional mesh only, i.e., $\mathit{\varphi }=\mathit{\varphi }(x,y)$.

We choose a similar decomposition for the test functions ϑ_x and ϑ_y:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\mathit{\vartheta }}_x& ={\mathit{\vartheta }}_x^{\mathrm{b}}+{\mathit{\vartheta }}_x^{\mathrm{sh}}\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right),\\ {\mathit{\vartheta }}_y& ={\mathit{\vartheta }}_y^{\mathrm{b}}+{\mathit{\vartheta }}_y^{\mathrm{sh}}\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right),\end{array}\end{array}$$

and they are defined using the same basis functions ϕ_i, as follows:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\mathit{\vartheta }}_x(x,y,z)& =\sum _{i=\mathrm{1}}^{n_{\mathrm{f}}}{\mathit{\vartheta }}_{x,i}^{\mathrm{b}}{\mathit{\varphi }}_i(x,y)+\sum _{i=\mathrm{1}}^{n_{\mathrm{f}}}{\mathit{\vartheta }}_{x,i}^{\mathrm{sh}}{\mathit{\varphi }}_i(x,y)\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right),\\ {\mathit{\vartheta }}_y(x,y,z)& =\sum _{i=\mathrm{1}}^{n_{\mathrm{f}}}{\mathit{\vartheta }}_{y,i}^{\mathrm{b}}{\mathit{\varphi }}_i(x,y)+\sum _{i=\mathrm{1}}^{n_{\mathrm{f}}}{\mathit{\vartheta }}_{y,i}^{\mathrm{sh}}{\mathit{\varphi }}_i(x,y)\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right),\end{array}\end{array}$$

where ${\mathit{\vartheta }}_{x,i}^{\mathrm{b}},{\mathit{\vartheta }}_{y,i}^{\mathrm{b}}$ and ${\mathit{\vartheta }}_{x,i}^{\mathrm{sh}},{\mathit{\vartheta }}_{y,i}^{\mathrm{sh}}$ are arbitrary coefficients.

We insert Eqs. (9) and (11) into the weak formulation (Eq. 6), and then we replace the ice viscosity μ by an ad hoc depth-averaged viscosity $\stackrel{\mathrm{‾}}{\mathit{\mu }}$ (discussed in Sect. 2.4). Thus, we analytically integrate the resulting formulation along the vertical axis. This generates a set of equations defined over the horizontal xy plane (two-dimensional mesh). The resulting nonlinear element matrix system is

$$\begin{array}{}\text{(12)}& \left[\begin{array}{cccc}{\mathbf{K}}_{\mathrm{11}}& {\mathbf{K}}_{\mathrm{12}}& {\mathbf{K}}_{\mathrm{13}}& {\mathbf{K}}_{\mathrm{14}}\\ {\mathbf{K}}_{\mathrm{21}}& {\mathbf{K}}_{\mathrm{22}}& {\mathbf{K}}_{\mathrm{23}}& {\mathbf{K}}_{\mathrm{24}}\\ {\mathbf{K}}_{\mathrm{31}}& {\mathbf{K}}_{\mathrm{32}}& {\mathbf{K}}_{\mathrm{33}}& {\mathbf{K}}_{\mathrm{34}}\\ {\mathbf{K}}_{\mathrm{41}}& {\mathbf{K}}_{\mathrm{42}}& {\mathbf{K}}_{\mathrm{43}}& {\mathbf{K}}_{\mathrm{44}}\end{array}\right]\left[\begin{array}{l}{\mathit{v}}_x^{\mathrm{b}}\\ {\mathit{v}}_x^{\mathrm{sh}}\\ {\mathit{v}}_y^{\mathrm{b}}\\ {\mathit{v}}_y^{\mathrm{sh}}\end{array}\right]=\left[\begin{array}{c}{\mathit{F}}_{\mathrm{1}}\\ {\mathit{F}}_{\mathrm{2}}\\ {\mathit{F}}_{\mathrm{3}}\\ {\mathit{F}}_{\mathrm{4}}\end{array}\right],\end{array}$$

where K_ij are matrices of size $n_{\mathrm{f},\mathrm{e}}\times n_{\mathrm{f},\mathrm{e}}$, with n_f,e being the total number of basis functions defined on each element. Here, we employ triangular elements and P1 Lagrangian basis functions and, therefore, $n_{\mathrm{f},\mathrm{e}}=\mathrm{3}$. In the element matrix system (Eq. 12), ${\mathit{v}}_x^{\mathrm{b}}$, ${\mathit{v}}_x^{\mathrm{sh}}$, ${\mathit{v}}_y^{\mathrm{b}}$, and ${\mathit{v}}_y^{\mathrm{sh}}$ are vectors evaluated on each element's node (e.g., ${\mathit{v}}_x^{\mathrm{b}}=[v_{x,\mathrm{1}}^{\mathrm{b}},v_{x,\mathrm{2}}^{\mathrm{b}},\mathrm{\dots },v_{x,n_{\mathrm{f},\mathrm{e}}}^{\mathrm{b}}{]}^{\mathit{⊺}}$). The loading vectors F_i are also evaluated on all nodes of each element. Appendices D–G present the detailed technical steps of the vertical integration and the resulting expressions of each term (K_ij and F_i) of Eq. (12). We note that we still rely on a numerical integration over the horizontal plane to evaluate the element matrices and loading vectors of Eq. (12).

2.4 Vertically integrated ice viscosity

A delicate aspect of the analytical depth integration of the HO weak formulation (Eq. 6) is how to handle the ice viscosity that is a function of ice velocity. The ice viscosity $\mathit{\mu }=\mathit{\mu }(x,y,z)$ is assumed to follow Glen's flow law (Glen, 1955):

$$\begin{array}{}\text{(13)}& \mathit{\mu }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{B}{{\dot{\mathit{\epsilon }}}_{\mathrm{e}}^{(n-\mathrm{1})/n}}},\end{array}$$

where $B=A^{-\mathrm{1}/n}$ is the associated rate factor, and ${\dot{\mathit{\epsilon }}}_{\mathrm{e}}$ is the effective strain rate tensor, which for HO is defined as

$$\begin{array}{}\text{(14)}& {\dot{\mathit{\epsilon }}}_{\mathrm{e}}=\sqrt{{\dot{\mathit{\epsilon }}}_{xx}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{yy}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{xy}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{xz}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{yz}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{xx}{\dot{\mathit{\epsilon }}}_{yy}},\end{array}$$

with each component defined as

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& {\dot{\mathit{\epsilon }}}_{xx}={\displaystyle \frac{\partial v_x}{\partial x}},\\ & {\dot{\mathit{\epsilon }}}_{yy}={\displaystyle \frac{\partial v_y}{\partial y}},\\ & {\dot{\mathit{\epsilon }}}_{xy}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial v_x}{\partial y}}+{\displaystyle \frac{\partial v_y}{\partial x}}\right),\\ & {\dot{\mathit{\epsilon }}}_{xz}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial v_x}{\partial z}},\\ & {\dot{\mathit{\epsilon }}}_{yz}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial v_y}{\partial z}}.\end{array}\end{array}$$

Note that in Eq. (15), as assumed in the HO model, $\partial v_z/\partial x$ and $\partial v_z/\partial y$ are neglected.

The vertical integration of the weak formulation, Eq. (6), requires the integration of the viscosity multiplied by a function of z only, f(z), since the basis functions ϕ(x,y) employed in the approximation of the ice velocities v_x and v_y (Eq. 9) are depth-independent. The expression of f(z) is a function of the polynomial function $\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}$ and varies according to the definitions of matrices K_ij in the element matrix system (Eq. 12). The four expressions of f(z) that appear in the element stiffness matrix are shown in Appendix C. Thus, we can write a generic expression for the vertical integration as follows:

$$\begin{array}{}\text{(16)}& {\mathit{\mu }}^{\prime }(x,y)=\underset{b}{\overset{s}{\int }}\mathit{\mu }(x,y,z)f\left(z\right)\mathrm{d}z,\end{array}$$

where ${\mathit{\mu }}^{\prime }(x,y)$ is the vertical integration of the ice viscosity multiplied by the function f(z).

In the following steps we drop the x and y dependency in the viscosity notation for simplicity. To achieve a vertically integrated formulation, let us assume that Eq. (16) can be rewritten by defining an ad hoc “vertical averaged” viscosity $\stackrel{\mathrm{‾}}{\mathit{\mu }}=\stackrel{\mathrm{‾}}{\mathit{\mu }}(x,y)$ as follows:

$$\begin{array}{}\text{(17)}& {\mathit{\mu }}^{\prime }=\underset{b}{\overset{s}{\int }}\mathit{\mu }\left(z\right)f\left(z\right)\mathrm{d}z=\stackrel{\mathrm{‾}}{\mathit{\mu }}\underset{b}{\overset{s}{\int }}f\left(z\right)\mathrm{d}z=\stackrel{\mathrm{‾}}{\mathit{\mu }}F,\end{array}$$

where F is the primitive function of f(z) evaluated between b and s (i.e., $F=\left[F\right(z){]}_b^s$).

If Eq. (16) was evaluated by employing a Gaussian quadrature (between −1 and 1), the integral would be approximated by

$$\begin{array}{}\text{(18)}& {\mathit{\mu }}^{\prime }=\underset{b}{\overset{s}{\int }}\mathit{\mu }\left(z\right)f\left(z\right)\mathrm{d}z\approx {\displaystyle \frac{H}{\mathrm{2}}}\sum {\mathit{\omega }}_i\mathit{\mu }\left(z_i\right)f\left(z_i\right),\end{array}$$

with ω_i being the appropriate weights and z_i the vertical coordinate evaluated at

$$\begin{array}{}\text{(19)}& z_i={\displaystyle \frac{H}{\mathrm{2}}}{\mathit{\xi }}_i+{\displaystyle \frac{s+b}{\mathrm{2}}},{\mathit{\xi }}_i\in [-\mathrm{1},\mathrm{1}].\end{array}$$

By combining Eqs. (17) and (18), we identify the ad hoc averaged viscosity as follows:

$$\begin{array}{}\text{(20)}& \stackrel{\mathrm{‾}}{\mathit{\mu }}\approx {\displaystyle \frac{H}{\mathrm{2}F}}\sum {\mathit{\omega }}_i\mathit{\mu }\left(z_i\right)f\left(z_i\right).\end{array}$$

The accuracy of the assumption employed in Eq. (17) depends on the order of the numerical integration of Eq. (20). We test a different number of integration points in Appendix A. Based on these numerical tests, we employ an integration order equal to 5 to perform the ISMIP-HOM and MISMIP3D experiments. We note that the ad hoc viscosity $\stackrel{\mathrm{‾}}{\mathit{\mu }}(x,y)$ varies within the triangular elements. In ISSM, we compute the value of $\stackrel{\mathrm{‾}}{\mathit{\mu }}(x,y)$ through Eq. (20) at each quadrature point used during the two-dimensional numerical integration of the triangular-element matrices K_ij of Eq. (12). We also note that all the components of the effective strain rate tensor ${\dot{\mathit{\epsilon }}}_{\mathrm{e}}$ (Eq. 15) are evaluated using the horizontal velocities defined by Eq. (7).

3 Numerical experiments

In all numerical experiments performed here, a Picard iteration scheme is used to solve the nonlinear stress balance equations. We employ the bi-conjugate gradient and block–Jacobi precondition to solve the linear system at each Picard iteration.

3.1 ISMIP-HOM setup

To assess the performance of MOLHO, we first compare the velocities from HO and MOLHO based on the ISMIP-HOM benchmark (Pattyn et al., 2008). This benchmark aims to test the response of ice flow models in diagnostic runs using different scales of glacier geometries and flow regimes. Here, we perform experiments A, C, and E, for which only the ice velocities are computed (i.e., transient simulations are not considered). Table 1 shows the parameters employed in these experiments. We describe here the main characteristics of the experiments, and the details are found in Pattyn et al. (2008).

Table 1 Constants used in the ISMIP-HOM experiment.

Experiment A consists of a flow over a sinusoidal bumpy bed with different wavelengths. The wavelengths (domain size) vary from 5 to 160 km. The amplitude of the bumps is 500 m. The bed is inclined in the x direction. Basal velocities are set to zero, and, therefore, this experiment is suitable to assess the performance of models in simulating internal deformation only (vertical shear). Periodic boundary conditions are imposed over the lateral boundaries of the model domain.

Experiment C is similar to experiment A. The bedrock is parallel to ice surface and has no bump. The ice base is allowed to slide over the bed following a Weertman-type sliding law:

$$\begin{array}{}\text{(21)}& {\mathit{\tau }}^{\mathrm{b}}=-C\parallel {\mathit{v}}^{\mathrm{b}}{\parallel }^{m-\mathrm{1}}{\mathit{v}}^{\mathrm{b}},\end{array}$$

where τ^b is the basal stress, C is the friction coefficient, v^b is the ice base velocity, and m is the sliding law exponent. In experiment C, m=1 and the friction coefficient C is defined by a sinusoidal function with different wavelengths that also vary from 5 to 160 km. This experiment assesses the performance of the model in stream-type flows.

The domain of experiment E is based on the 5 km center line of a real alpine glacier (Haut Glacier d'Arolla, Switzerland). The ice flows only in the x direction, and the width of the glacier (y direction) is set to be constant and equal to 1 m over the entire glacier extent. The ice flow rate is uniform over the entire domain (see Table 1). Two basal boundary conditions are employed in this experiment: frozen bed and slip bed. In the former, basal velocities are set to zero, and in the latter, a small region of perfect slip (i.e., C=0) is imposed between x=2.2 and x=2.5 km.

We employ the same horizontal mesh resolution in both HO and MOLHO models in all the experiments. To generate the three-dimensional mesh for HO, we extrude the horizontal mesh using 20 vertical layers. Due to convergence issues in the nonlinear solver, we use only eight vertical layers to perform experiment E with HO.

3.2 MISMIP3D setup

The MISMIP3D benchmark (Pattyn et al., 2013) allows us to assess the convergence of the grounding line position at steady state and the reversibility of the grounding line in the MOLHO flow model. In this experiment, an analytical solution for the steady-state grounding line only exists for the shelfy stream approximation that does not account for vertical shear (Schoof, 2007a, b). Thus, for comparison purposes, we also run this benchmark with the three-dimensional HO model.

The setup consists of a marine ice sheet flowing along the x axis over a uniformly sloping bedrock whose elevation, r(x,y) (in meters, negative if below sea level), is described by

$$\begin{array}{}\text{(22)}& r(x,y)=\mathrm{100}-x,\end{array}$$

with x in kilometers. The ice sheet covers an area of 800×50 km². A no-flow condition (v_x=0) is imposed at the ice divide (x=0). A Neumann boundary condition based on ocean water pressure is applied at the calving front (x=800 km). We impose a free-slip condition (v_y=0) at the lateral boundaries of the domain (y=0 and y=50 km).

The basal sliding is described by a Weertman friction law as given by Eq. (21). The values of C and m as well as other parameters used in the experiments are shown in Table 3.

The MISMIP3D experiment is divided into three phases: (i) steady state (Stnd), (ii) basal friction perturbation (P75S), and (iii) basal friction restoration (P75R).

In phase Stnd, the model starts from a 100 m thick ice shelf and runs forward in time until a steady state is reached. Here, the steady state is reached after 30 000 years of simulation.

The phase P75S starts from the end of phase Stnd. The friction coefficient C is replaced by C^* at the beginning of phase P75S. The new coefficient C^* is defined as

$$\begin{array}{}\text{(23)}& C^{*}=C\left[\mathrm{1}-\mathrm{0.75}\times \exp \left(-{\displaystyle \frac{{\left(x-x_{\mathrm{g}}\right)}^{\mathrm{2}}}{\mathrm{2}{x}_{\mathrm{c}}^{\mathrm{2}}}}-{\displaystyle \frac{{\left(y-y_{\mathrm{g}}\right)}^{\mathrm{2}}}{\mathrm{2}{y}_{\mathrm{c}}^{\mathrm{2}}}}\right)\right],\end{array}$$

where x_g is the steady-state grounding line position at y=0, and y_g=0. The spatial extent of the friction change along the x and y directions is given by x_c=150 and y_c=10 km, respectively. The model runs for 100 years with C^*. This friction perturbation forces the grounding line to migrate asymmetrically: it advances along y=0 and retreats along y=50 km.

At the end of phase P75S, the friction coefficient is restored to its initial value, C. Then, the model runs for 30 000 years to reach steady state again. This is the phase P75R. The grounding line reversibility is achieved if the positions of the grounding line at the end of phases Stnd and P75R are sufficiently close to each other.

Previous intercomparison model results showed a strong mesh-resolution dependency to achieve grounding line reversibility (Pattyn et al., 2013; Feldmann et al., 2014). This dependency decreases if a sub-element friction parameterization is employed (Cornford et al., 2013; Feldmann et al., 2014; Seroussi et al., 2014). Here, we use sub-element parameterization type I (SEP1) as presented in Seroussi et al. (2014). To assess the robustness of the implementation of the MOLHO flow model, we run the three phases with different mesh resolutions (Table 2). The meshes are generated using the Delaunay triangulation (Hecht, 2006), and a spatially uniform resolution is applied to most of them. To save computational resources, we generate the 250 m resolution mesh by refining the 500 m resolution mesh between x=400 and x=650 km. The HO model mesh is generated by extruding the two-dimensional mesh employed in the MOLHO model, except for the two finest resolutions (500 and 250 m) for which we coarsen all triangular elements located at x<500 km and at x>620 km before the mesh extrusion. We employ 10 equally spaced vertical layers in all HO model meshes. The number of elements for each model and mesh resolutions is shown in Table 2. We note that the finest mesh (250 m) is used to analyze the convergence of the steady-state grounding line (phase Stnd) only, since grounding line reversibility is already achieved with the coarser meshes (see Sect. 4.2).

Table 2 Mesh resolutions and the associated number of elements used in the MISMIP3D experiment.

^* Used in the phase Stnd only.

Table 3 Constants used in the MISMIP3D experiment.

Figure 1 ISMIP-HOM experiment A: surface velocity at y=0.25L for different wavelengths, L, obtained with the three-dimensional higher-order model (HO) and with the two-dimensional MOno-Layer Higher-Order (MOLHO) model. The surface velocity is computed by taking the Euclidean norm of the horizontal velocities at the ice surface (i.e., $\sqrt{v_x^{\mathrm{2}}+v_y^{\mathrm{2}}}$). We employ 20 vertical layers for the HO model.

Figure 2 ISMIP-HOM experiment C: surface velocity at y=0.25L for different wavelengths, L, obtained with the three-dimensional higher-order model (HO) and with the two-dimensional MOno-Layer Higher-Order (MOLHO) model. Results obtained with L1L2 and shelfy stream approximation (SSA) models are also shown. The surface velocity is computed by taking the Euclidean norm of the horizontal velocities at the ice surface (i.e., $\sqrt{v_x^{\mathrm{2}}+v_y^{\mathrm{2}}}$). We employ 20 vertical layers for the HO model.

4 Results

4.1 ISMIP-HOM

Figure 1 shows the Euclidean norm of the horizontal velocities (v_x,v_y) at the ice surface and at y=0.25L for several wavelengths (L) of bedrock bumps for experiment A. MOLHO and HO produce virtually the same surface speeds for wavelengths equal to or higher than 40 km. The difference between the speeds is about 4 % for L=40 km and about 2 % for L=160 km. For shorter bump lengths (≤20 km), MOLHO overestimates the ice speeds in comparison to HO. The differences are up to 60 % for a wavelength equal to 5 km. The differences are ≤11 % for L=20 km.

In the case where the ice base slides over the bed (experiment C), MOLHO and HO yield similar surface velocities for all wavelengths (Fig. 2). The differences between both models vary from 0.05 % for L=5 km to 1.2 % for a L=160 km. For comparison purposes, we also run experiment C with two other vertically integrated models that are also commonly employed for ice streams: a version of the L1L2 model (Schoof and Hindmarsh, 2010; Perego et al., 2012) and the shelfy stream approximation (SSA; Morland, 1987; MacAyeal, 1989). The L1L2 considers vertical shear in the ice viscosity computation, while the SSA model does not account for vertical shear stresses. For larger wavelengths of friction coefficient (≥40 km), all four models yield similar velocities, as expected (Pattyn et al., 2008). The differences in speed between L1L2 and SSA in comparison to HO are up to 6 % for L=40 km and <5 % for L=160 km. For shorter wavelengths (<40 km), vertical shear starts playing a role on the ice flow, which may compromise the results obtained with L1L2 and SSA models in comparison to HO. The maximum speed differences vary from 3 % to 5 % for L=5 and L=20 km, respectively, for both L1L2 and SSA models. On the other hand, MOLHO “approximates” the HO results even for the shortest L.

In the real glacier setup (experiment E) and considering a frozen bed, MOLHO overestimates the surface velocities up to 40 %¹ over about the first 60 % of the glacier flow line (Fig. 3). This region presents steeper bedrock than the region close to the glacier snout (see the glacier profile in Pattyn et al., 2008). Considering the whole domain, the RMSE (root-mean-square error) is 5.5 m yr⁻¹. Imposing a small region of perfect sliding in the model domain causes the glacier to speed up. The difference in speeds achieved with MOLHO and HO is <40 %². The differences are higher over the first 3 km of the glacier, as already noted in the frozen bed setup. The RMSE is 7.6 m yr⁻¹.

Figure 3 ISMIP-HOM experiment E: surface velocity in the ice flow direction obtained with the three-dimensional higher-order model (HO) and with the two-dimensional Mono-Layer Higher-Order (MOLHO) model. Two basal conditions are presented: frozen bed (a) and slip bed (b). We employ eight vertical layers for the HO model.

4.2 MISMIP3D

The steady-state positions of the grounding line (GL) at y=0 for different mesh resolutions are shown in Fig. 4. We also present the grounding line positions obtained with the HO model. To estimate the convergence error in GL positions as we increase the mesh resolution, we compute the relative error of the steady-state GL position $x_{\mathrm{g}}^h$ obtained with a mesh resolution h, defined as follows:

$$\begin{array}{}\text{(24)}& {\mathit{\epsilon }}_{x_{\mathrm{g}}}^h=\left|{\displaystyle \frac{x_{\mathrm{g}}^{\mathrm{2}h}-x_{\mathrm{g}}^h}{x_{\mathrm{g}}^{\mathrm{2}h}}}\right|,\end{array}$$

where $x_{\mathrm{g}}^{\mathrm{2}h}$ is the position of the grounding line obtained with the closest coarser mesh resolution. The estimated convergence errors for both MOLHO and HO models are also shown in Fig. 4.

Figure 4 MISMIP3D experiment: steady-state positions (a) and estimated error convergence (b) of the grounding line (GL) at y=0 obtained with different mesh resolutions using HO and MOLHO. The error convergence is estimated using Eq. (24). For comparison, the steady-state grounding line position predicted by the boundary layer (BL) theory that is based on the SSA model is also shown.

Using a coarser mesh resolution (5 km), the grounding lines obtained with both HO and MOLHO extend a few kilometers beyond 600 km, reaching the steady-state GL location predicted by the boundary layer analysis and achieved by numerical simulations using SSA flow models (Schoof, 2007b; Cornford et al., 2013; Feldmann et al., 2014; Seroussi et al., 2014). However, as obtained by other ice sheet models employing full-Stokes, hybrid, or higher-order flow equations, the grounding line at steady state tends to be located between 520 and 580 km (e.g., Pattyn et al., 2013; Cornford et al., 2013; Gagliardini et al., 2016). As seen in Fig. 4, the grounding lines obtained with MOLHO and HO approach 560 km as the mesh gets finer. With the finest mesh used here (250 m), the steady-state GL positions computed with MOLHO and HO models are 562.100 and 559.916 km, respectively, a difference of less than 0.5 % between them.

As also shown in Fig. 4, both HO and MOLHO flow models show similar convergence of the relative error in GL position. The relative errors decrease with mesh resolution, reaching about 1 % with a mesh resolution of 250 m.

At the end of the restoration phase (P75R), the GL position is similar to the position obtained in the steady-state phase (Stnd, Fig. 5), which indicates that grounding line reversibility is fully achieved with MOLHO. Overall, the simulated displacement of the grounding line in the perturbation phase (P75S) using MOLHO is comparable to the ones simulated with the HO model, as seen in Fig. 5. At y=0, the GL advances about 11–13 km, while at y=50 km the GL retreats 4–6 km (Table 4). As shown in Table 4, the displacements of the grounding line in the P75S experiment using MOLHO are in good agreement with the ones obtained with HO.

Figure 5 Grounding line positions of the MISMIP3D experiments obtained with different mesh resolutions. Black lines are the steady-state grounding line positions (Stnd). Red lines are the positions obtained after 100 years of basal friction perturbation (P75S). Dashed blue lines are the new steady-state grounding line positions obtained after resetting the friction coefficient to its original value (P75R).

Figure 6 Time-dependent grounding line positions of the MISMIP3D experiments obtained with different mesh resolutions. Solid lines are the grounding line positions at y=0. Dashed lines are the grounding line positions at y=50 km. Red lines are the positions obtained during the basal friction perturbation experiment (P75S). Blue lines are the positions obtained during the first 100 years after restoring the friction coefficient to its original value (P75R). Note that the x axes are reversed for the P75R experiment results.

Table 4 Displacements of the grounding line (ΔGL) at the end of the perturbation phase (P75S) obtained with MOLHO and HO models.

The transient response of the grounding line at y=0 and y=50 km during the P75S experiment and the first 100 years of the restoration phase (P75R) is shown in Fig. 6. As shown in the figure, the perturbation and restoration responses of the grounding obtained with MOLHO and HO are similar.

To analyze the differences in ice velocity in MOLHO and HO, we perform diagnostic runs employing the same MISMIP3D-type ice sheet profile. For this run, we use the ice geometry (surface, base, thickness, and grounding line) obtained with HO in the phase Stnd and with a mesh resolution of 1 km. In Fig. 7 we show the surface and basal velocities around the grounding line, where there are more pronounced differences between the model results. We also show the depth-averaged velocity (Eq. 8 for MOLHO model), which controls the advection of ice downstream. Interestingly, for this ice sheet profile, we find that MOLHO produces slightly lower velocities compared to HO (Fig. 7), which may explain why the steady-state grounding lines are slightly advanced in comparison to the HO model simulations (Fig. 5). The differences over the region shown in Fig. 7 are up to 4.5 m yr⁻¹, which is less than 1.5 %.

Figure 7 Comparison of ice velocities in the x direction obtained with MOLHO and HO using the same glacier geometry (thickness and grounding line). The figure shows the velocities at the surface and base of the glacier as well as the depth-averaged velocities. The velocities are computed using the steady-state geometry obtained with the phase Stnd with the HO model and horizontal mesh resolution of 1 km. Dashed lines are the grounding line positions.

4.3 Time performance

To assess and illustrate the computational cost of MOLHO in time-dependent ice sheet simulations, we run the first 500 years of the spin-up phase of the MISMIP3D setup (Stnd). This phase is long enough to simulate a large grounding line migration (advance, in this case), therefore requiring several iterations of the stress balance solver.

We use four different mesh resolutions (5, 2, 1, and 0.5 km) and three different ice flow approximations: MOLHO, HO, and SSA. For this time performance comparison, we use the same MOLHO mesh for both HO and SSA. For HO, we extrude the horizontal mesh using 10 vertical layers. The number of elements of each MOLHO mesh is shown in Table 2, while Table 5 shows the total number of degrees of freedom for each mesh and flow model. The experiments are performed in parallel using 28 computational cores (CPUs) in a 2× Intel Xeon E5-2680v4 2.8 GHz with 128 GB of memory. The time step is 0.125 year, which fulfills the Courant–Friedrichs–Lewy (CFL) condition for all meshes.

Table 5 Computational time to run the first 500 years of the MISMIP3D spin-up phase (Stnd). For comparison purposes, we show the simulation times obtained by three different ice flow models: MOLHO, HO, and SSA. SSA and HO employ the same mesh as MOLHO. For HO, we use 10 vertical layers. The number of elements for each mesh resolution employed in MOLHO and SSA is shown in Table 2.

^* Estimated.

As shown in Table 5, the computation time using MOLHO is, overall, 1 order of magnitude lower in comparison to the time spent in HO. The ratio between the simulation times obtained with MOLHO and HO varies with mesh resolution. For coarser meshes, MOLHO is 10× faster than HO. For the finer mesh resolutions, MOLHO is almost 20× faster than HO. The number of degrees of freedom for the HO model is 5.5× higher than MOLHO, for all mesh resolutions employed here.

In comparison to SSA, MOLHO is between 1.5 and 3.5× slower, depending on the mesh resolution. This is probably due to different rates of convergence of the (iterative) linear solver (during the Picard iterations) that vary with the number of degrees of freedom of each mesh employed here (Table 5). MOLHO has, by construction, 2× more degrees of freedom than SSA for any given mesh, which demands theoretically 2× more³ computational time.

5 Discussion

Higher-order models were developed to approximate the Stokes solution at a lower computational cost. The ISMIP-HOM initiative (Pattyn et al., 2008) was the first intercomparison effort to compare higher-order and full-Stokes models in experiments where the SIA model is not valid. In that intercomparison exercise, results from models based on the Blatter–Pattyn approximation (the “LMLa models” in the ISMIP-HOM nomenclature; HO in this paper) differed from Stokes results in cases where all stress components play an equally important role in ice dynamics (i.e., experiment A for ≤20 km). In these cases, the higher-order assumptions may not be appropriate. The results from LMLa and full-Stokes models showed a better agreement for lower aspect ratios in which longitudinal stresses become dominant (experiment A for ≥40 km and experiment C). Also, the LMLa and full-Stokes results showed relatively good agreement for the Haut Glacier d’Arolla geometry (experiment E) with the no basal slip condition; abrupt changes in basal condition (slip–no-slip jumps) led to a poor agreement between all models, however.

Overall, MOLHO's results (ice velocities and grounding line positions) are in close agreement with the results obtained with a full three-dimensional HO model. This is especially true for the marine ice sheet with fast sliding (e.g., MISMIP3D). Similar results are also observed in the ISMIP-HOM (Pattyn et al., 2008) experiment C, where basal sliding is present (see Fig. 2). In experiment A, MOLHO overestimates ice speeds compared to HO for deformation-dominated cases (shorter wavelengths of bedrock elevation). Our computation of the vertically averaged ice viscosity $\stackrel{\mathrm{‾}}{\mathit{\mu }}$ might be the reason behind some of the velocity overestimation, since the longitudinal stresses present in the upper layers of the ice are being “watered down” by the averaging with the softer ice at the bed. This suggests that the choice of the vertical quadrature scheme (e.g., order and/or vertical distribution of Gaussian points) plays an important role in the results of the MOLHO formulation (Brinkerhoff and Johnson, 2015). Also, the determinant of the Jacobian in the HO model is not constant within the prismatic elements (due to the mapping between the “real” and reference elements), which might be another source to explain these differences.

By running the three phases of the MISMIP3D setup, we obtain grounding line reversibility for all mesh resolutions employed here. The evolution of the grounding line during the friction perturbation and restoration phases (P75S and P75R, respectively) is consistent with previous intercomparison results (Pattyn et al., 2013). These results confirm the robustness of our implementation and the ability of MOLHO to simulate grounding line dynamics in marine ice sheets at a comparable accuracy as HO. Also, using the MOLHO flow model, the convergence of the steady-state grounding line positions with increasing mesh resolutions follows the ones achieved with the three-dimensional higher-order model. This suggests that the ice velocity computation and the friction parameterizations employed in MOLHO are consistent with the HO model implemented in ISSM.

In terms of time performance, the computational cost of MOLHO is relatively higher than the SSA model. This is expected since MOLHO has 2× more degrees of freedom than SSA and also includes a vertical integration of the ice viscosity at every Gauss point when we compute the element stiffness matrix. In comparison to HO with 10 vertical layers employed in the time performance analysis, the computational cost is at least 10× lower. While it is possible to decrease the simulation time of HO by decreasing the number of vertical layers (at the risk of reducing numerical accuracy) and by improving the scalability (i.e., increasing the number of CPUs and/or computational nodes), there will always be bottlenecks related to three-dimensional meshes that will increase computational cost in comparison to two-dimensional meshes.

In MOLHO, the vertical variation in the horizontal velocities is prescribed by a SIA-derived polynomial function of the order equal to n+1. In the experiments performed here, we assumed isothermal ice flow. However, the vertical integration of the ice viscosity can accommodate vertical variation in ice temperature by changing the polynomial function accordingly, e.g., solving for the analytical SIA solution and using a scaled version of that as the vertical basis function.

6 Conclusions

In this work, we present the formulation of a vertically integrated MOno-Layer Higher-Order (MOLHO) model and compare its performance with the three-dimensional Blatter–Pattyn model (HO) using two benchmarks: ISMIP-HOM and MISMIP3D. In the experiments with no basal sliding of the ISMIP-HOM setup (experiment A), MOLHO produces surface velocities close to the ones obtained by the HO model for wavelengths of bedrock elevation equal to or higher than 40 km. For shorter wavelengths (<40 km) where the approximation in vertical shear may break down, MOLHO tends to overestimate surface velocities. When the ice base is allowed to slide (experiment C), MOLHO and HO produce similar surface velocities for all wavelengths of friction coefficient. These results suggest that MOLHO is more suitable to simulate larger ice masses such as marine ice sheets than small ice caps. Yet, for the Haut Glacier d'Arolla experiment (experiment E), MOLHO overestimates upstream velocities up to 40 % but produces similar speed over the downstream part of the glacier profile. In the MISMIP3D experiments, the steady-state grounding line positions obtained with MOLHO are very close to the positions obtained with the HO model. The grounding line migrations have a similar temporal behavior. This suggests that the grounding line dynamics achieved with MOLHO are consistent with the original HO model. By carrying out a time performance analysis in time-dependent runs, we found that the computation cost of MOLHO is at least 10× lower in comparison to the HO model. MOLHO is therefore an excellent alternative to HO for long simulations, such as paleo ice sheet modeling.

Appendix A: Integration order of the ice viscosity

The accuracy of the semi-analytic integration of the vertical axis in MOLHO formulation depends on the integration order employed in the ice viscosity computation (Eq. 20). Figure A1 shows the surface speeds for the ISMIP-HOM experiments obtained with four orders tested here: 3, 5, 10, and 15. We only present the cases for the shortest wavelength (5 km) for experiments A and C and the frozen bed case for experiment E.

For experiments A and C, the differences between the velocities obtained with an integration order equal to 3 and equal to 15 are smaller than 1 %. The difference reduces to 0.1 % when we compare orders 5 and 15. Integration orders 10 and 15 yield similar surface velocities: the differences are smaller than 0.0005 %.

For experiment E, the speed differences are up to 3 % comparing integration orders 3 and 15; between 5 and 15, the differences reduce to 0.5 %. Comparing orders 10 and 15, the speed differences are smaller than 0.005 %.

Figure A1 ISMIP-HOM experiments: surface velocities for experiments A, C, and E obtained with the MOLHO model using different integration orders for the ice viscosity computation. Only the short-wavelength (5 km) case is shown for experiments A and C. For experiment E, only the frozen bed case is shown.

Appendix B: Useful identities

Using the definition of ζ given by Eq. (3), we use the following identities:

$$\begin{array}{}\text{(B1)}& {\displaystyle }{\displaystyle \frac{\partial }{\partial z}}\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right)={\displaystyle \frac{n+\mathrm{1}}{H}}{\mathit{\zeta }}^n,\text{(B2)}& {\displaystyle }\underset{b}{\overset{s}{\int }}{\mathit{\zeta }}^{n+\mathrm{1}}\mathrm{d}z={\displaystyle \frac{H}{(n+\mathrm{2})}},\text{(B3)}& {\displaystyle }\underset{b}{\overset{s}{\int }}\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right)\mathrm{d}z=H{\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}},\text{(B4)}& {\displaystyle }\underset{b}{\overset{s}{\int }}{\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right)}^{\mathrm{2}}\mathrm{d}z=\mathrm{2}H{\displaystyle \frac{(n+\mathrm{1}{)}^{\mathrm{2}}}{(\mathrm{2}n+\mathrm{3})(n+\mathrm{2})}},\text{(B5)}& {\displaystyle }\underset{b}{\overset{s}{\int }}{\left({\displaystyle \frac{\partial }{\partial z}}\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right)\right)}^{\mathrm{2}}\mathrm{d}z={\displaystyle \frac{(n+\mathrm{1}{)}^{\mathrm{2}}}{H(\mathrm{2}n+\mathrm{1})}},\text{(B6)}& {\displaystyle }\underset{b}{\overset{s}{\int }}\mathit{\zeta }\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right)\mathrm{d}z={\displaystyle \frac{H}{\mathrm{2}}}{\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{3})}},\text{(B7)}& {\displaystyle }\underset{b}{\overset{\mathrm{0}}{\int }}{\mathit{\zeta }}^{n+\mathrm{1}}\mathrm{d}z={\displaystyle \frac{H}{n+\mathrm{2}}}\left[\mathrm{1}-{\left({\displaystyle \frac{s}{H}}\right)}^{n+\mathrm{2}}\right],\text{(B8)}& {\displaystyle }\underset{b}{\overset{\mathrm{0}}{\int }}{\mathit{\zeta }}^{n+\mathrm{2}}\mathrm{d}z={\displaystyle \frac{H}{(n+\mathrm{3})}}\left[\mathrm{1}-{\left({\displaystyle \frac{s}{H}}\right)}^{n+\mathrm{3}}\right],\text{(B9)}& {\displaystyle }\begin{array}{rl}\underset{b}{\overset{\mathrm{0}}{\int }}z\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right)\mathrm{d}z& =-{\displaystyle \frac{b^{\mathrm{2}}}{\mathrm{2}}}-{\displaystyle \frac{sH}{n+\mathrm{2}}}\left[\mathrm{1}-{\left({\displaystyle \frac{s}{H}}\right)}^{n+\mathrm{2}}\right]\\ & +{\displaystyle \frac{H^{\mathrm{2}}}{(n+\mathrm{3})}}\left[\mathrm{1}-{\left({\displaystyle \frac{s}{H}}\right)}^{n+\mathrm{3}}\right].\end{array}\end{array}$$

Appendix C: Vertically integrated ice viscosity

The function f(z) can have the following expressions:

$$\begin{array}{}\text{(C1)}& f\left(z\right)=\left\{\begin{array}{cl}{f}_{\mathrm{1}}\left(z\right)& =\mathrm{1},\\ f_{\mathrm{2}}\left(z\right)& =\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right),\\ f_{\mathrm{3}}\left(z\right)& ={\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right)}^{\mathrm{2}},\\ f_{\mathrm{4}}\left(z\right)& ={\left(\frac{\partial }{\partial z}\left(\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}\right)\right)}^{\mathrm{2}},\end{array}\right.\end{array}$$

where we keep ζ (Eq. 3) to simplify the notation. The vertical integration of the primitive of f, F, reads

$$\begin{array}{}\text{(C2)}& F=\left\{\begin{array}{cl}{F}_{\mathrm{1}}& =H,\\ F_{\mathrm{2}}& =H\frac{(n+\mathrm{1})}{(n+\mathrm{2})},\\ F_{\mathrm{3}}& =\mathrm{2}H\frac{(n+\mathrm{1}{)}^{\mathrm{2}}}{(\mathrm{2}n+\mathrm{3})(n+\mathrm{2})},\\ F_{\mathrm{4}}& =\frac{(n+\mathrm{1}{)}^{\mathrm{2}}}{H(\mathrm{2}n+\mathrm{1})}\end{array}\right..\end{array}$$

Therefore, we have four expressions for the vertically integrated ad hoc ice viscosity $\stackrel{\mathrm{‾}}{\mathit{\mu }}$ that are used in different parts of the element stiffness matrix (see Appendix D):

$$\begin{array}{}\text{(C3)}& \stackrel{\mathrm{‾}}{\mathit{\mu }}=\left\{\begin{array}{cl}{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}& =\frac{H}{\mathrm{2}{F}_{\mathrm{1}}}\sum {\mathit{\omega }}_i\mathit{\mu }\left(z_i\right)f_{\mathrm{1}}\left(z_i\right),\\ {\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}& =\frac{H}{\mathrm{2}{F}_{\mathrm{2}}}\sum {\mathit{\omega }}_i\mathit{\mu }\left(z_i\right)f_{\mathrm{2}}\left(z_i\right),\\ {\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}& =\frac{H}{\mathrm{2}{F}_{\mathrm{3}}}\sum {\mathit{\omega }}_i\mathit{\mu }\left(z_i\right)f_{\mathrm{3}}\left(z_i\right),\\ {\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{4}}& =\frac{H}{\mathrm{2}{F}_{\mathrm{4}}}\sum {\mathit{\omega }}_i\mathit{\mu }\left(z_i\right)f_{\mathrm{4}}\left(z_i\right)\end{array}\right..\end{array}$$

Appendix D: Stiffness matrix

By inserting Eqs. (9) and (11) into the weak formulation Eq. (6), and by employing the useful identities shown in Appendix B and the four expressions of the ad hoc ice viscosity defined in Appendix C, we integrate the terms of the element stiffness matrix along the vertical axis. In the following integrals, Ω_e is the two-dimensional element domain defined over the horizontal xy plane (e.g., triangle, for Delaunay triangulation-based meshes).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{11}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\mathrm{4}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{12}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{4}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}\right){\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{13}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\mathrm{2}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{14}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{2}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}\right){\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{21}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{4}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}\right){\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{22}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{4}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.{\displaystyle \frac{\mathrm{2}(n+\mathrm{1}{)}^{\mathrm{2}}}{(\mathrm{2}n+\mathrm{3})(n+\mathrm{2})}}+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{4}}{\mathit{\varphi }}_j{\mathit{\varphi }}_i{\displaystyle \frac{(n+\mathrm{1}{)}^{\mathrm{2}}}{H(\mathrm{2}n+\mathrm{1})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{23}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{2}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}\right){\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{24}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{2}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.{\displaystyle \frac{\mathrm{2}(n+\mathrm{1}{)}^{\mathrm{2}}}{(\mathrm{2}n+\mathrm{3})(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{31}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\mathrm{2}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{32}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{2}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right){\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{33}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\mathrm{4}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{34}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{4}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right){\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{41}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{2}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right){\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{42}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{2}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.{\displaystyle \frac{\mathrm{2}(n+\mathrm{1}{)}^{\mathrm{2}}}{(\mathrm{2}n+\mathrm{3})(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{43}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{4}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right){\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{44}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\mathrm{4}H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial y}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial y}}+H{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{3}}{\displaystyle \frac{\partial {\mathit{\varphi }}_j}{\partial x}}{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.{\displaystyle \frac{\mathrm{2}(n+\mathrm{1}{)}^{\mathrm{2}}}{(\mathrm{2}n+\mathrm{3})(n+\mathrm{2})}}+{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{4}}{\mathit{\varphi }}_j{\mathit{\varphi }}_i{\displaystyle \frac{(n+\mathrm{1}{)}^{\mathrm{2}}}{H(\mathrm{2}n+\mathrm{1})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\end{array}$$

Appendix E: Stiffness matrix – basal friction

In ISSM, basal friction is written as a function of the basal velocities ($v_x^{\mathrm{b}},v_y^{\mathrm{b}}$). Therefore, the matrices ${\mathbf{K}}_{ii}^{\mathrm{b}}(i=\mathrm{1},\mathrm{3})$ below should be added to the respective stiffness matrices. The parameter ${\mathit{\alpha }}^{\mathrm{2}}={\mathit{\alpha }}^{\mathrm{2}}(x,y,v_x^{\mathrm{b}},v_y^{\mathrm{b}})$ depends on the friction law employed.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{11}}^{\mathrm{b}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\mathit{\alpha }}^{\mathrm{2}}{\mathit{\varphi }}_j{\mathit{\varphi }}_i\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_{\mathrm{33}}^{\mathrm{b}}=\left[\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\mathit{\alpha }}^{\mathrm{2}}{\mathit{\varphi }}_j{\mathit{\varphi }}_i\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right]\end{array}$$

Appendix F: Loading vector – driving stress

The driving stress is integrated over the two-dimensional element domain, Ω_e. The resulting loading vectors are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{F}}_{\mathrm{1}}=\left[-\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\mathit{\rho }gH{\displaystyle \frac{\partial s}{\partial x}}{\mathit{\varphi }}_i\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right],\\ {\displaystyle }& {\displaystyle }{\mathit{F}}_{\mathrm{2}}=\left[-\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\mathit{\rho }gH{\displaystyle \frac{\partial s}{\partial x}}{\mathit{\varphi }}_i{\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right],\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{F}}_{\mathrm{3}}=\left[-\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\mathit{\rho }gH{\displaystyle \frac{\partial s}{\partial y}}{\mathit{\varphi }}_i\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right],\\ {\displaystyle }& {\displaystyle }{\mathit{F}}_{\mathrm{4}}=\left[-\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\mathit{\rho }gH{\displaystyle \frac{\partial s}{\partial y}}{\mathit{\varphi }}_i{\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{2})}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right].\end{array}$$

Appendix G: Loading vector – calving front

If an element edge ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{w}}$ is at the calving front, the vectors ${\mathit{F}}_i^{\mathrm{w}}(i=\mathrm{1},\mathrm{\dots },\mathrm{4})$ below should be added to the element loading vector. In the following integrals, n_x and n_y are the components of the unit vector pointing outward edge ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{w}}$.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{F}}_{\mathrm{1}}^{\mathrm{w}}=\left[\underset{{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{w}}}{\int }\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }g{H}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{w}}g{b}^{\mathrm{2}}\right)n_x{\mathit{\varphi }}_i\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{w}}\right]\\ {\displaystyle }& {\displaystyle }{\mathit{F}}_{\mathrm{2}}^{\mathrm{w}}=\left[\underset{{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{w}}}{\int }{\mathit{\rho }}_{\mathrm{w}}g\left\{-{\displaystyle \frac{b^{\mathrm{2}}}{\mathrm{2}}}-{\displaystyle \frac{sH}{n+\mathrm{2}}}\left[\mathrm{1}-{\left({\displaystyle \frac{s}{H}}\right)}^{n+\mathrm{2}}\right]\right.\right.\\ {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{H^{\mathrm{2}}}{(n+\mathrm{3})}}\left[\mathrm{1}-{\left({\displaystyle \frac{s}{H}}\right)}^{n+\mathrm{3}}\right]\right\}{n}_x{\mathit{\varphi }}_i\\ {\displaystyle }& {\displaystyle }\left.+\mathit{\rho }gH\left[{\displaystyle \frac{H}{\mathrm{2}}}{\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{3})}}\right]n_x{\mathit{\varphi }}_i\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{w}}\right]\\ {\displaystyle }& {\displaystyle }{\mathit{F}}_{\mathrm{3}}^{\mathrm{w}}=\left[\underset{{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{w}}}{\int }\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }g{H}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{w}}g{b}^{\mathrm{2}}\right)n_y{\mathit{\varphi }}_i\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{w}}\right]\\ {\displaystyle }& {\displaystyle }{\mathit{F}}_{\mathrm{4}}^{\mathrm{w}}=\left[\underset{{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{w}}}{\int }{\mathit{\rho }}_{\mathrm{w}}g\left\{-{\displaystyle \frac{b^{\mathrm{2}}}{\mathrm{2}}}-{\displaystyle \frac{sH}{n+\mathrm{2}}}\left[\mathrm{1}-{\left({\displaystyle \frac{s}{H}}\right)}^{n+\mathrm{2}}\right]\right.\right.\\ {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{H^{\mathrm{2}}}{(n+\mathrm{3})}}\left[\mathrm{1}-{\left({\displaystyle \frac{s}{H}}\right)}^{n+\mathrm{3}}\right]\right\}{n}_y{\mathit{\varphi }}_i\\ {\displaystyle }& {\displaystyle }\left.+\mathit{\rho }gH\left[{\displaystyle \frac{H}{\mathrm{2}}}{\displaystyle \frac{(n+\mathrm{1})}{(n+\mathrm{3})}}\right]n_y{\mathit{\varphi }}_i\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{w}}\right]\end{array}$$

Appendix H: Compact weak formulation of MOLHO

Here we present a compact form of the weak formulation of MOLHO, which could be useful for implementation in software frameworks like FEniCS (Alnæs et al., 2015) and Firedrake (Rathgeber et al., 2016). It is based on the notation presented in Brinkerhoff and Johnson (2015).

First, we define the higher-order strain rate tensor, $\dot{\mathit{\epsilon }}\left(\mathit{v}\right)$, as

$$\dot{\mathit{\epsilon }}\left(\mathit{v}\right)=\left[\begin{array}{llc}\mathrm{2}\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}& \frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial v_x}{\partial y}+\frac{\partial v_y}{\partial x}\right)& \frac{\mathrm{1}}{\mathrm{2}}\frac{\partial v_x}{\partial z}\\ \frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial v_x}{\partial y}+\frac{\partial v_y}{\partial x}\right)& \mathrm{2}\frac{\partial v_y}{\partial y}+\frac{\partial v_x}{\partial x}& \frac{\mathrm{1}}{\mathrm{2}}\frac{\partial v_y}{\partial z}\end{array}\right].$$

Second, we define the space of trial and test functions, 𝒱^h⊂𝒱, using a tensor product as follows:

$$\begin{array}{ll}{\displaystyle }{\mathcal{V}}^h& {\displaystyle }=\left(\mathrm{span}\left\{\sum _{i=\mathrm{1}}^{n_{\mathrm{f}}}{\mathit{\varphi }}_i(x,y)\left[{\mathit{\psi }}_{\mathrm{1}}+{\mathit{\psi }}_{\mathrm{2}}\left(z\right)\right]\right\},\right.\\ {\displaystyle }& {\displaystyle }\left.\mathrm{span}\left\{\sum _{i=\mathrm{1}}^{n_{\mathrm{f}}}{\mathit{\varphi }}_i(x,y)\left[{\mathit{\psi }}_{\mathrm{1}}+{\mathit{\psi }}_{\mathrm{2}}\left(z\right)\right]\right\}\right),\end{array}$$

where ψ₁ and ψ₂ are the vertical basis functions defined for the MOLHO model as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{1}}=\mathrm{1},\\ {\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{2}}=\mathrm{1}-{\mathit{\zeta }}^{n+\mathrm{1}}.\end{array}$$

Note that ψ₁ and ψ₂ are related to the basal velocities and vertical shear deformation, respectively. Thus, the horizontal ice velocities, v, are approximated by v^h∈𝒱^h:

$${\mathit{v}}^h(x,y,z)=\sum _{i=\mathrm{1}}^{n_{\mathrm{f}}}{\mathit{\varphi }}_i(x,y)\left[{\mathit{v}}_{\mathrm{1},i}{\mathit{\psi }}_{\mathrm{1}}+{\mathit{v}}_{\mathrm{2},i}{\mathit{\psi }}_{\mathrm{2}}\left(z\right)\right],$$

where ${\mathit{v}}_{\mathrm{1},i}={\mathit{v}}_i^{\mathrm{b}}$ and ${\mathit{v}}_{\mathrm{2},i}={\mathit{v}}_i^{\mathrm{sh}}$ are the basal and shear velocities computed over the horizontal mesh nodes i, respectively. Likewise, the test function ϑ^h∈𝒱^h is defined as

$${\mathit{\vartheta }}^h(x,y,z)=\sum _{i=\mathrm{1}}^{n_{\mathrm{f}}}{\mathit{\varphi }}_i(x,y)\left[{\mathit{\vartheta }}_{\mathrm{1},i}{\mathit{\psi }}_{\mathrm{1}}+{\mathit{\vartheta }}_{\mathrm{2},i}{\mathit{\psi }}_{\mathrm{2}}\left(z\right)\right],$$

with ϑ_1,i and ϑ_2,i being arbitrary.

In the following, Ω is the horizontal xy plane of the model domain, n is the unit vector pointing outward edge Γ^w, and the gradient operator (∇) applies to the entire three-dimensional domain (i.e., the x, y, and z directions). Therefore, the weak formulation is written as follows.

Find v^h∈𝒱^h such that

$$\begin{array}{ll}{\displaystyle }\underset{\mathrm{\Omega }}{\int }\underset{b}{\overset{s}{\int }}\mathrm{\nabla }{\mathit{\vartheta }}^h& {\displaystyle }\cdot \left(\mathrm{2}\mathit{\mu }\dot{\mathit{\epsilon }}\left({\mathit{v}}^h\right)\right)\mathrm{d}z\mathrm{d}\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}^{\mathrm{b}}}{\int }{\mathit{\vartheta }}^h\cdot \left({\mathit{\alpha }}^{\mathrm{2}}{\mathit{v}}^h\right)\mathrm{d}{\mathrm{\Gamma }}^{\mathrm{b}}\\ {\displaystyle }& {\displaystyle }=-\underset{\mathrm{\Omega }}{\int }\underset{b}{\overset{s}{\int }}{\mathit{\vartheta }}^h\cdot \left(\mathit{\rho }g\mathrm{\nabla }s\right)\mathrm{d}z\mathrm{d}\mathrm{\Omega }\\ {\displaystyle }& {\displaystyle }+\underset{{\mathrm{\Gamma }}^{\mathrm{w}}}{\int }{\mathit{\vartheta }}^h\cdot \left(\mathrm{2}\mathit{\mu }\dot{\mathit{\epsilon }}\left({\mathit{v}}^h\right)\mathit{n}\right)\mathrm{d}{\mathrm{\Gamma }}^{\mathrm{w}},\\ {\displaystyle }& {\displaystyle }{\mathit{v}}^h={\mathit{v}}^D\mathrm{on}{\mathrm{\Gamma }}^D,\\ {\displaystyle }& {\displaystyle }\forall {\mathit{\vartheta }}^h\in {\mathcal{V}}^h.\end{array}$$