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

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 5 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 of 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 model (MOLHO) 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 10 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 MOnoLayer 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. 15


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 20 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

Polynomial function for vertical shear
The specific polynomial employed to describe the vertical variation of the horizontal velocities is based on the Shallow Ice Approximation (SIA, Hutter, 1983). In SIA, the horizontal velocities (v = [v x , v y ] ) of an isothermal ice sheet following Glen's law are given by: 70 where v b (= v b x , v b y ) are 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's law exponent, ρ is the ice density, and g is the gravitational acceleration. We write Equation 1 as: where v sh = v sh x , v sh y is the internal deformation contribution to the surface velocities, and ζ is an auxiliary variable defined as: Equation 2 has been widely employed by the glaciological community in ice sheet models to setup boundary conditions (e.g., Raymond, 1983) and describe vertical deformation within the glacier body (e.g., Bueler and Brown, 2009;Bassis, 2010;Brinkerhoff and Johnson, 2015;Shapero et al., 2021).
In SIA, v sh are defined using the local surface gradient only (Eq. 1). In the MOno Layer Higher-Order Model (MOLHO), 80 we employ Eq. 2 into a Higher-Order (HO) weak formulation such that v sh 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 1 − ζ n+1 . Appendix B presents some useful identities employed in this work.

Weak formulation 85
The three-dimensional Blatter-Pattyn or Higher-Order (HO) ice flow model is defined as (Blatter, 1995;Pattyn, 2003): 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: 90 where 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 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 V 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 95 model boundary and viscous sliding at the base) of HO is: 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), and Γ D is the glacier boundary where Dirichlet boundary conditions with values v D x , v D y are imposed.

100
Based on SIA, we decompose the horizontal velocities in MOLHO following Eq. 2: 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 105 velocities) could be used directly into the continuity equation that controls the advection of the ice thickness. Here, once v b and v sh are computed, the vertically-averaged velocities,v, are obtained by: To apply the finite element discretization to Eq. 6, we approximate the horizontal velocities by employing two-dimensional basis functions φ j as follows: with v b x,j , v b y,j and v sh x,j , v sh y,j 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., φ = φ (x, y).
We choose a similar decomposition for the test functions ϑ x and ϑ y : and they are defined using the same basis functions φ i , as follows: where ϑ b x,i , ϑ b y,i and ϑ sh x,i , ϑ sh y,i are arbitrary coefficients. We insert Eq. 9 and Eq. 11 into the weak formulation (Eq. 6), and then we replace the ice viscosity µ by an ad hoc depth-120 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: where K ij are matrices of size n f,e ×n f,e , with n f,e being the total number of basis functions defined on each element. Here, we 125 employ triangular elements and P1-Lagrangian basis functions and, therefore, n f,e = 3. In the element matrix system (Eq. 12), and v sh y are vectors evaluated on each element's node (e.g. The loading vectors F i are also evaluated on all nodes of each element. Appendix D presents 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.

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 (Glen, 1955): where B = A −1/n is the associated rate factor, andε e is the effective strain rate tensor, which for HO is defined as: 135ε e = ε 2 xx +ε 2 yy +ε 2 xy +ε 2 xz +ε 2 yz +ε xxεyy , with each component defined as: 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-140 independent. The expression of f (z) is a function of the polynomial function 1 − ζ n+1 and varies according to the definitions of matrices K ij into 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: In the following steps we drop the x and y dependency in the viscosity notation for simplicity. To achieve a vertically 145 integrated formulation, let's assume that Eq. 16 can be rewritten by defining an ad hoc 'vertical averaged' viscosityμ =μ (x, y) as follows: where F is the primitive function of f (z) evaluated between b and s (i.e., F = [F (z)] s b ). If Eq. 16 was evaluated by employing a Gaussian quadrature (between -1 and 1), the integral would be approximated by: with ω i being the appropriate weights and z i the vertical coordinate evaluated at: By combining Eq. 17 and Eq. 18, we identify the ad hoc averaged viscosity as follows: 155 The accuracy of the assumption employed in Eq. 17 depends on the order of the numerical integration of Eq. 20. We test different number of integration points in Sect. A. Based on these numerical tests, we employ 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 K ij of Eq. 12. We also note that all the components of the effective 160 strain rate tensorε e (Eq. 15) are evaluated using the horizontal velocities defined by Eq. 7.

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

170
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 175 slide over the bed following a Weertman-type sliding law: 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.

180
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 km and x=2.5 km. 185 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 8 vertical layers to perform experiment E with HO.

MISMIP3D setup
The MISMIP3D benchmark (Pattyn et al., 2013) allows us to assess the convergence of the grounding line position at steady 190 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 m, negative if below sea level), is described by: with x in km. The ice sheet covers an area of 800×50 km 2 . 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 200 parameters used in the experiments are shown in Table 3.
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.

205
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: where x g is the 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 km and y c = 10 km, respectively. The model runs for 100 years with C * . This friction 210 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 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 reversibil-215 ity (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 Delaunay triangulation (Hecht, 2006), and a spatially uniform resolution is applied to most of them. To save computational resources, we generate the 250-m   which we coarsen all triangular elements located at x < 500 km and at x > 620 km before the mesh extrusion. We employ ten 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 analyse the convergence of the steady state grounding line (phase for 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 235 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 240 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 km 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 245 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 induces 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 .

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 h g obtained with a mesh resolution h, defined as follows: where x 2h g 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.
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 260 simulations using SSA flow models (Schoof, 2007b;Cornford et al., 2013;Feldmann et al., 2014;Seroussi et al., 2014).    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 to 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 265 km 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 270 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 To analyze the differences in ice velocity in MOLHO and HO, we perform diagnostic runs employing the same MISMIP3Dtype 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 Figure 7 we show the surface and basal velocities around the grounding 280 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 km yr −1 , which is less than 1.5%.

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 290 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 ten 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 CFL condition for all meshes.

295
As shown in Table 5, the computation time using MOLHO is, overall, one 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. 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 employs the same mesh as MOLHO.
For HO, we use ten vertical layers. The number of elements for each mesh resolution employed in MOLHO/SSA is shown in Table 2 In comparison to SSA, MOLHO is between 1.5 to 3.5× slower, depending on the mesh resolution. This is probably due to 300 different rates of convergence of the linear solver during the Picard iterations. MOLHO has, by construction, 2× more degrees of freedom than SSA for any given mesh, which demands theoretically 4× more 3 computational time.

Discussion
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, respec-305 tively) are 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 similar accuracy as HO.
Using the MOLHO flow model, the convergence of the steady state grounding line positions with mesh resolutions follows the ones achieved with the three-dimensional Higher-Order model. This suggests that the ice velocity computation and the 310 friction parameterizations employed in MOLHO are consistent with the HO model implemented in ISSM. Yet, there is still room to improve the convergence rate of both MOLHO and HO models with mesh resolution.
Overall, MOLHO's results (ice velocities and grounding line positions) are in close agreement with the results obtained with a full three-dimensional 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 speed relative to the 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 Johnson, 320 2015). Also, the determinant of the Jacobian in 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.
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 ten vertical layers employed in the time 325 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 of the horizontal velocities is prescribed by a SIA-derived polynomial function of order 330 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.

Conclusions
In this work, we presented the formulation of a vertically integrated,

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 355 employed in the ice viscosity computation (Eq. 20). Figure A1 shows the surface speeds for the ISMIP-HOM experiments obtained with four orders tested: 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. Integrating orders 10 and 15 yield 360 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%.
The function f (z) can have the following expressions (see Appendix D): where we keep ζ (Eq. 3) to simplify the notation. The vertical integration of the primitive of f , F , reads: Therefore, we have four expressions for the vertically integrated ad hoc ice viscosityμ that are used in different parts of the 380 element stiffness matrix (see Appendix D): Appendix D: Stiffness matrix By inserting Eq. 9 and Eq. 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 385 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). In ISSM, basal friction is written as a function of the basal velocities (v b x , v b y ). Therefore, the matrices K b ii (i = 1, 3) below 405 should be added to the respective stiffness matrices. The parameter α 2 = α 2 x, y, v b x , v b y depends on the friction law employed. If an element edge Γ w e is at the calving front, the vectors F w i (i = 1, . . . , 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 Γ w e .