Articles | Volume 16, issue 1
Research article
20 Jan 2022
Research article |  | 20 Jan 2022

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

Thiago Dias dos Santos, Mathieu Morlighem, and Douglas Brinkerhoff

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, Hutter1983; Morland1987; MacAyeal1989; Blatter1995; Pattyn2003; Goldberg2011). 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 (Blatter1995; Pattyn2003), 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 Raymond1978; Bassis2010; Brinkerhoff and Johnson2015; 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 Johnson2015; 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 Johnson2015). It has been used for simulations of glaciers and ice sheets (Brinkerhoff and Johnson2015; 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 Johnson2015) 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, Hutter1983). In the SIA, the horizontal velocities (v=[vx,vy]) of an isothermal ice sheet following Glen's law are given by

(1) v = v b - 2 ( ρ g ) n s n - 1 A n + 1 H n + 1 1 - s - z H n + 1 s ,

where vb (=[vxb,vyb]) 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

(2) v ( x , y , z ) = v b ( x , y ) + v sh ( x , y ) 1 - ζ n + 1 ,

where vsh=[vxsh,vysh] is the internal deformation contribution to the surface velocities, and ζ is an auxiliary variable defined as

(3) ζ = s - z H = 1 - z - b H ,

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., Raymond1983) and describe vertical deformation within the glacier body (e.g., Langdon and Raymond1978; Bueler and Brown2009; Bassis2010; Brinkerhoff and Johnson2015; Shapero et al.2021).

In the SIA, vsh 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 vsh values are evaluated only after solving a nonlinear system resulting from the HO stress-balance equations. Details of the computation of vsh 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 1-ζn+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 (Blatter1995; Pattyn2003)


(4) x μ v x y + v y x + y 4 μ v y y + 2 μ v x x + z μ v y z = ρ g s y .

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

(5) τ b = - α 2 v b ,

where vb=[vxb,vyb] is the horizontal velocity at the glacier base and α=α(x,y) is the friction coefficient that, in general, is a (nonlinear) function of the basal velocities vb. 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

(6) Ω 4 μ v x x + 2 μ v y y ϑ x x + μ v x y + v y x ϑ x y + μ v x z ϑ x z d Ω + Ω μ v x y + v y x ϑ y x + 4 μ v y y + 2 μ v x x ϑ y y + μ v y z ϑ y z d Ω + Γ b α 2 v x b ϑ x d Γ b + Γ b α 2 v y b ϑ y d Γ b = - Ω ρ g s x ϑ x d Ω - Ω ρ g s y ϑ y d Ω + Γ w 4 μ v x x + 2 μ v y y n x + μ v x y + v y x n y ϑ x d Γ w + Γ w 4 μ v y y + 2 μ v x x n y + μ v x y + v y x n x ϑ y d Γ w , v x = v x D , v y = v y D on Γ D , ϑ x , ϑ y V ,

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), nx and ny are the components of the unit vector pointing outward edge Γw, and ΓD is the glacier boundary where Dirichlet boundary conditions with values vxD,vyD are imposed.

2.3 Finite element discretization

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

(7) v x ( x , y , z ) = v x b ( x , y ) + v x sh ( x , y ) 1 - ζ n + 1 , v y ( x , y , z ) = v y b ( x , y ) + v y sh ( x , y ) 1 - ζ n + 1 .

We note that the decomposition employed here differs from the approach adopted in previous work (Brinkerhoff and Johnson2015). 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 vb and vsh are computed, the vertically averaged velocities, v, are obtained by

(8) v ( x , y ) = v b ( x , y ) + v sh ( x , y ) ( n + 1 ) ( n + 2 ) .

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

(9) v x ( x , y , z ) j = 1 n f v x , j b ϕ j ( x , y ) + j = 1 n f v x , j sh ϕ j ( x , y ) 1 - ζ n + 1 , v y ( x , y , z ) j = 1 n f v y , j b ϕ j ( x , y ) + j = 1 n f v y , j sh ϕ j ( x , y ) 1 - ζ n + 1 ,

with vx,jb,vy,jb and vx,jsh,vy,jsh being the values of the base velocities and surface shear velocities evaluated at nodal points j. In Eq. (9), nf 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., ϕ=ϕ(x,y).

We choose a similar decomposition for the test functions ϑx and ϑy:

(10) ϑ x = ϑ x b + ϑ x sh 1 - ζ n + 1 , ϑ y = ϑ y b + ϑ y sh 1 - ζ n + 1 ,

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

(11) ϑ x ( x , y , z ) = i = 1 n f ϑ x , i b ϕ i ( x , y ) + i = 1 n f ϑ x , i sh ϕ i ( x , y ) 1 - ζ n + 1 , ϑ y ( x , y , z ) = i = 1 n f ϑ y , i b ϕ i ( x , y ) + i = 1 n f ϑ y , i sh ϕ i ( x , y ) 1 - ζ n + 1 ,

where ϑx,ib,ϑy,ib and ϑx,ish,ϑy,ish 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 μ (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

(12) K 11 K 12 K 13 K 14 K 21 K 22 K 23 K 24 K 31 K 32 K 33 K 34 K 41 K 42 K 43 K 44 v x b v x sh v y b v y sh = F 1 F 2 F 3 F 4 ,

where Kij are matrices of size nf,e×nf,e, with nf,e being the total number of basis functions defined on each element. Here, we employ triangular elements and P1 Lagrangian basis functions and, therefore, nf,e=3. In the element matrix system (Eq. 12), vxb, vxsh, vyb, and vysh are vectors evaluated on each element's node (e.g., vxb=[vx,1b,vx,2b,,vx,nf,eb]). The loading vectors Fi are also evaluated on all nodes of each element. Appendices DG present the detailed technical steps of the vertical integration and the resulting expressions of each term (Kij and Fi) 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 μ=μ(x,y,z) is assumed to follow Glen's flow law (Glen1955):

(13) μ = 1 2 B ε ˙ e ( n - 1 ) / n ,

where B=A-1/n is the associated rate factor, and ε˙e is the effective strain rate tensor, which for HO is defined as

(14) ε ˙ e = ε ˙ x x 2 + ε ˙ y y 2 + ε ˙ x y 2 + ε ˙ x z 2 + ε ˙ y z 2 + ε ˙ x x ε ˙ y y ,

with each component defined as

(15) ε ˙ x x = v x x , ε ˙ y y = v y y , ε ˙ x y = 1 2 v x y + v y x , ε ˙ x z = 1 2 v x z , ε ˙ y z = 1 2 v y z .

Note that in Eq. (15), as assumed in the HO model, vz/x and vz/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 vx and vy (Eq. 9) are depth-independent. The expression of f(z) is a function of the polynomial function 1-ζn+1 and varies according to the definitions of matrices Kij 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:

(16) μ ( x , y ) = b s μ ( x , y , z ) f ( z ) d z ,

where μ(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 μ=μ(x,y) as follows:

(17) μ = b s μ ( z ) f ( z ) d z = μ b s f ( z ) d z = μ F ,

where F is the primitive function of f(z) evaluated between b and s (i.e., F=[F(z)]bs).

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

(18) μ = b s μ ( z ) f ( z ) d z H 2 ω i μ z i f z i ,

with ωi being the appropriate weights and zi the vertical coordinate evaluated at

(19) z i = H 2 ξ i + s + b 2 , ξ i [ - 1 , 1 ] .

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

(20) μ H 2 F ω i μ z i f z i .

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 μ(x,y) varies within the triangular elements. In ISSM, we compute the value of μ(x,y) through Eq. (20) at each quadrature point used during the two-dimensional numerical integration of the triangular-element matrices Kij of Eq. (12). We also note that all the components of the effective strain rate tensor ε˙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 1Constants used in the ISMIP-HOM experiment.

Download Print Version | Download XLSX

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:

(21) τ b = - C v b m - 1 v b ,

where τb is the basal stress, C is the friction coefficient, vb 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 (Schoof2007a, 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

(22) r ( x , y ) = 100 - x ,

with x in kilometers. The ice sheet covers an area of 800×50 km2. A no-flow condition (vx=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 (vy=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

(23) C * = C 1 - 0.75 × exp - x - x g 2 2 x c 2 - y - y g 2 2 y c 2 ,

where xg is the steady-state grounding line position at y=0, and yg=0. The spatial extent of the friction change along the x and y directions is given by xc=150 and yc=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 (Hecht2006), 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 2Mesh resolutions and the associated number of elements used in the MISMIP3D experiment.

* Used in the phase Stnd only.

Download Print Version | Download XLSX

Table 3Constants used in the MISMIP3D experiment.

Download Print Version | Download XLSX

Figure 1ISMIP-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., vx2+vy2). We employ 20 vertical layers for the HO model.


Figure 2ISMIP-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., vx2+vy2). We employ 20 vertical layers for the HO model.


4 Results


Figure 1 shows the Euclidean norm of the horizontal velocities (vx,vy) 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 Hindmarsh2010; Perego et al.2012) and the shelfy stream approximation (SSA; Morland1987; MacAyeal1989). 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 %1 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−1. 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 %2. 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−1.

Figure 3ISMIP-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.



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 xgh obtained with a mesh resolution h, defined as follows:

(24) ε x g h = x g 2 h - x g h x g 2 h ,

where xg2h 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 4MISMIP3D 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 (Schoof2007b; 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 5Grounding 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 6Time-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 4Displacements of the grounding line (ΔGL) at the end of the perturbation phase (P75S) obtained with MOLHO and HO models.

Download Print Version | Download XLSX

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−1, which is less than 1.5 %.

Figure 7Comparison 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  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 5Computational 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.

Download Print Version | Download XLSX

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  more3 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 μ 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 Johnson2015). 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  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 A1ISMIP-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:

Appendix C: Vertically integrated ice viscosity

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

(C1) f ( z ) = f 1 ( z ) = 1 , f 2 ( z ) = 1 - ζ n + 1 , f 3 ( z ) = 1 - ζ n + 1 2 , f 4 ( z ) = z 1 - ζ n + 1 2 ,

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

(C2) F = F 1 = H , F 2 = H ( n + 1 ) ( n + 2 ) , F 3 = 2 H ( n + 1 ) 2 ( 2 n + 3 ) ( n + 2 ) , F 4 = ( n + 1 ) 2 H ( 2 n + 1 ) .

Therefore, we have four expressions for the vertically integrated ad hoc ice viscosity μ that are used in different parts of the element stiffness matrix (see Appendix D):

(C3) μ = μ 1 = H 2 F 1 ω i μ z i f 1 z i , μ 2 = H 2 F 2 ω i μ z i f 2 z i , μ 3 = H 2 F 3 ω i μ z i f 3 z i , μ 4 = H 2 F 4 ω i μ z i f 4 z i .
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).



Appendix E: Stiffness matrix – basal friction

In ISSM, basal friction is written as a function of the basal velocities (vxb,vyb). Therefore, the matrices Kiib(i=1,3) below should be added to the respective stiffness matrices. The parameter α2=α2(x,y,vxb,vyb) depends on the friction law employed.

Appendix F: Loading vector – driving stress

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


Appendix G: Loading vector – calving front

If an element edge Γew is at the calving front, the vectors Fiw(i=1,,4) below should be added to the element loading vector. In the following integrals, nx and ny are the components of the unit vector pointing outward edge Γew.

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, ε˙(v), as


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


where ψ1 and ψ2 are the vertical basis functions defined for the MOLHO model as


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


where v1,i=vib and v2,i=vish are the basal and shear velocities computed over the horizontal mesh nodes i, respectively. Likewise, the test function ϑh∈𝒱h is defined as


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 vh∈𝒱h such that

Code availability

The MOno-Layer Higher-Order (MOLHO) model evaluated here is currently implemented in ISSM. The code can be downloaded, compiled, and executed following the instructions available on the ISSM website: (NASA2021). The public SVN repository for the ISSM code can also be found directly at (Larour et al.2020). The version of the code for this study, corresponding to ISSM release 4.18, is SVN version tag number 26413. The documentation of the code version used here is available at (Morlighem et al.2020).

Author contributions

DB proposed the idea of a single layer of a higher-order ice flow model. TDS and MM worked on the finite element formulation and implemented the MOno-Layer Higher-Order (MOLHO) model in ISSM. TDS designed the experimental setup and performed the simulations. TDS and MM led the analysis of the results. TDS led the initial writing of the paper. All authors contributed to writing the final version of the paper.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This work is from the PROPHET project, a component of the International Thwaites Glacier Collaboration (ITGC). This is ITGC contribution no. ITGC-054. We thank the editor, Benjamin Smith, and the two reviewers, Daniel Shapero and the one anonymous referee, for their positive and constructive comments, which improved the clarity of the paper.

Financial support

This research has been supported by the National Science Foundation (grant no. 1739031) and the Heising–Simons Foundation (grant no. 2021-3059).

Review statement

This paper was edited by Benjamin Smith and reviewed by Daniel Shapero and one anonymous referee.


Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M. E., and Wells, G. N.: The FEniCS Project Version 1.5, Arch. Numer. Soft., 3, 9–23,, 2015. a

Aschwanden, A., Fahnestock, M. A., Truffer, M., Brinkerhoff, D. J., Hock, R., Khroulev, C., Mottram, R., and Khan, S. A.: Contribution of the Greenland Ice Sheet to sea level over the next millennium, Sci. Adv., 5, eaav9396,, 2019. a

Bassis, J.: Hamilton-type principles applied to ice-sheet dynamics: new approximations for large-scale ice-sheet flow, J. Glaciol., 56, 497–513,, 2010. a, b

Blatter, H.: Velocity and stress-fields in grounded glaciers: A simple algorithm for including deviatoric stress gradients, J. Glaciol., 41, 333–344, 1995. a, b, c

Brinkerhoff, D., Truffer, M., and Aschwanden, A.: Sediment transport drives tidewater glacier periodicity, Nat. Commun., 8, 90,, 2017. a

Brinkerhoff, D. J. and Johnson, J. V.: Dynamics of thermally induced ice streams simulated with a higher-order flow model, J. Geophys. Res.-Earth, 120, 1743–1770,, 2015. a, b, c, d, e, f, g, h, i, j

Bueler, E. and Brown, J.: Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model, J. Geophys. Res.-Earth, 114, F03008,, 2009. a

Cornford, S. L., Martin, D. F., Graves, D. T., Ranken, D. F., Brocq, A. M. L., Gladstone, R. M., Payne, A. J., Ng, E. G., and Lipscomb, W. H.: Adaptive mesh, finite volume modeling of marine ice sheets, J. Comput. Phys., 232, 529–549,, 2013. a, b, c

Cuzzone, J. K., Morlighem, M., Larour, E., Schlegel, N., and Seroussi, H.: Implementation of higher-order vertical finite elements in ISSM v4.13 for improved ice sheet flow modeling over paleoclimate timescales, Geosci. Model Dev., 11, 1683–1694,, 2018. a

Edwards, T. L., Nowicki, S., Marzeion, B., Hock, R., Goelzer, H., Seroussi, H., Jourdain, N. C., Slater, D. A., Turner, F. E., Smith, C. J., McKenna, C. M., Simon, E., Abe-Ouchi, A., Gregory, J. M., Larour, E., Lipscomb, W. H., Payne, A. J., Shepherd, A., Agosta, C., Alexander, P., Albrecht, T., Anderson, B., Asay-Davis, X., Aschwanden, A., Barthel, A., Bliss, A., Calov, R., Chambers, C., Champollion, N., Choi, Y., Cullather, R., Cuzzone, J., Dumas, C., Felikson, D., Fettweis, X., Fujita, K., Galton-Fenzi, B. K., Gladstone, R., Golledge, N. R., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huss, M., Huybrechts, P., Immerzeel, W., Kleiner, T., Kraaijenbrink, P., Le clec'h, S., Lee, V., Leguy, G. R., Little, C. M., Lowry, D. P., Malles, J.-H., Martin, D. F., Maussion, F., Morlighem, M., O'Neill, J. F., Nias, I., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Radić, V., Reese, R., Rounce, D. R., Rückamp, M., Sakai, A., Shafer, C., Schlegel, N.-J., Shannon, S., Smith, R. S., Straneo, F., Sun, S., Tarasov, L., Trusel, L. D., Van Breedam, J., van de Wal, R., van den Broeke, M., Winkelmann, R., Zekollari, H., Zhao, C., Zhang, T., and Zwinger, T.: Projected land ice contributions to twenty-first-century sea level rise, Nature, 593, 74–82,, 2021. a

Feldmann, J., Albrecht, T., Khroulev, C., Pattyn, F., and Levermann, A.: Resolution-dependent performance of grounding line motion in a shallow model compared with a full-Stokes model according to the MISMIP3d intercomparison, J. Glaciol., 60, 353–360,, 2014. a, b, c

Gagliardini, O., Brondex, J., Gillet-Chaulet, F., Tavard, L., Peyaud, V., and Durand, G.: Brief communication: Impact of mesh resolution for MISMIP and MISMIP3d experiments using Elmer/Ice, The Cryosphere, 10, 307–312,, 2016. a

Glen, J. W.: The creep of polycrystalline ice, P. Roy. Soc. Lond. A, 228, 519–538,, 1955. a

Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., Abe-Ouchi, A., Shepherd, A., Simon, E., Agosta, C., Alexander, P., Aschwanden, A., Barthel, A., Calov, R., Chambers, C., Choi, Y., Cuzzone, J., Dumas, C., Edwards, T., Felikson, D., Fettweis, X., Golledge, N. ., Greve, R., Humbert, A., Huybrechts, P., Le clec'h, S., Lee, V., Leguy, G., Little, C., Lowry, D. P., Morlighem, M., Nias, I., Quiquet, A., Rückamp, M., Schlegel, N.-J., Slater, D. A., Smith, R. S., Straneo, F., Tarasov, L., van de Wal, R., and van den Broeke, M.: The future sea-level contribution of the Greenland ice sheet: a multi-model ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096,, 2020. a

Goldberg, D. N.: A variationally derived, depth-integrated approximation to a higher-order glaciological flow model, J. Glaciol., 57, 157–170,, 2011. a

Hecht, F.: BAMG: Bidimensional Anisotropic Mesh Generator, Tech. rep., FreeFem++, available at: (last access: 11 January 2022), 2006. a

Hutter, K.: Theoretical Glaciology, Mathematical Approaches to Geophysics, D. Reidel Publishing Company, Dordrecht, the Netherlands, ISBN 978-94-015-1167-4,, 1983. a, b

Langdon, J. and Raymond, C. F.: Numerical calculation of adjustment of a glacier surface to perturbations of ice thickness, Mater. Glyatsiol. Issled. Khron. Obsuzhdeniya, 32, 233–239, 1978. a, b

Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res.-Earth, 117, F01022,, 2012. a

Larour, E., Morlighem, M., and Seroussi, H.: Ice-Sheet and Sea-Level System Model, svn repository [code],, last access: 20 November 2020. a, b

MacAyeal, D.: Large-scale ice flow over a viscous basal sediment: Theory and application to ice stream B, Antarctica, J. Geophys. Res.-Solid, 94, 4071–4087,, 1989. a, b

Morland, L. W.: Unconfined ice shelf flow, in: Dynamics of the West Antarctic Ice Sheet, vol. 4 of Glaciology and Quaternary Geology, edited by: van der Veen, C. and Oerlemans, J., Springer, Dordrecht, the Netherlands, 99–116, ISBN 978-94-009-3745-1,, 1987. a, b

Morlighem, M., Seroussi, H., Larour, E., Schlegel, N., Borstad, C., de Fleurian, B., Adhikari, S., Bondzio, J., Sommers, A., McCormack, F., and dos Santos, T. D.: Ice Sheet System Model 2020 (4.18), User Guide, available at: (last access: 26 August 2021), 2020. a

NASA: Ice-sheet and Sea-level System Model, available at:, last access: 26 August 2021. a

Pattyn, F.: A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes, J. Geophys. Res.-Solid, 108, 2382,, 2003. a, b, c

Pattyn, F., Perichon, L., Aschwanden, A., Breuer, B., de Smedt, B., Gagliardini, O., Gudmundsson, G. H., Hindmarsh, R. C. A., Hubbard, A., Johnson, J. V., Kleiner, T., Konovalov, Y., Martin, C., Payne, A. J., Pollard, D., Price, S., Rückamp, M., Saito, F., Souček, O., Sugiyama, S., and Zwinger, T.: Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP–HOM), The Cryosphere, 2, 95–108,, 2008. a, b, c, d, e, f, g

Pattyn, F., Perichon, L., Durand, G., Favier, L., Gagliardini, O., Hindmarsh, R. C. A., Zwinger, T., Albrecht, T., Cornford, S., Docquier, D., Fürst, J. J., Goldberg, D., Gudmundsson, G. H., Humbert, A., Hütten, M., Huybrechts, P., Jouvet, G., Kleiner, T., Larour, E., Martin, D., Morlighem, M., Payne, A. J., Pollard, D., Rückamp, M., Rybak, O., Seroussi, H., Thoma, M., and Wilkens, N.: Grounding-line migration in plan-view marine ice-sheet models: results of the ice2sea MISMIP3d intercomparison, J. Glaciol., 59, 410–422,, 2013. a, b, c, d, e

Perego, M., Gunzburger, M., and Burkardt, J.: Parallel finite-element implementation for higher-order ice-sheet models, J. Glaciol., 58, 76–88,, 2012. a

Rathgeber, F., Ham, D. A., Mitchell, L., Lange, M., Luporini, F., Mcrae, A. T. T., Bercea, G.-T., Markall, G. R., and Kelly, P. H. J.: Firedrake: Automating the Finite Element Method by Composing Abstractions, ACM Trans. Math. Softw., 43, 24,, 2016. a

Raymond, C. F.: Deformation in the Vicinity of Ice Divides, J. Glaciol., 29, 357–373,, 1983. a

Schoof, C.: Marine ice-sheet dynamics. Part 1. The case of rapid sliding, J. Fluid Mech., 573, 27–55,, 2007a. a

Schoof, C.: Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res.-Earth, 112, 1–19,, 2007b.  a, b

Schoof, C. and Hindmarsh, R. C. A.: Thin-Film Flows with Wall Slip: An Asymptotic Analysis of Higher Order Glacier Flow Models, Q. J. Mech. Appl. Math., 63, 73–114,, 2010. a

Seroussi, H., Morlighem, M., Larour, E., Rignot, E., and Khazendar, A.: Hydrostatic grounding line parameterization in ice sheet models, The Cryosphere, 8, 2075–2087,, 2014. a, b, c

Seroussi, H., Nowicki, S., Payne, A. J., Goelzer, H., Lipscomb, W. H., Abe-Ouchi, A., Agosta, C., Albrecht, T., Asay-Davis, X., Barthel, A., Calov, R., Cullather, R., Dumas, C., Galton-Fenzi, B. K., Gladstone, R., Golledge, N. R., Gregory, J. M., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huybrechts, P., Jourdain, N. C., Kleiner, T., Larour, E., Leguy, G. R., Lowry, D. P., Little, C. M., Morlighem, M., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Reese, R., Schlegel, N.-J., Shepherd, A., Simon, E., Smith, R. S., Straneo, F., Sun, S., Trusel, L. D., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., Zhao, C., Zhang, T., and Zwinger, T.: ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070,, 2020. a

Shapero, D. R., Badgeley, J. A., Hoffman, A. O., and Joughin, I. R.: icepack: a new glacier flow modeling package in Python, version 1.0, Geosci. Model Dev., 14, 4593–4616,, 2021. a, b, c, d, e


Excluding the first 15 % of the glacier domain since the velocities are relatively small there.


See footnote 1.


Considering an optimal iterative linear solver with time complexity of 𝒪(n), where n is the number of degrees of freedom.

Short summary
Projecting the future evolution of Greenland and Antarctica and their potential contribution to sea level rise often relies on computer simulations carried out by numerical ice sheet models. Here we present a new vertically integrated ice sheet model and assess its performance using different benchmarks. The new model shows results comparable to a three-dimensional model at relatively lower computational cost, suggesting that it is an excellent alternative for long-term simulations.