Comparison of ice dynamics using full-Stokes and Blatter–Pattyn approximation: application to the Northeast Greenland Ice Stream

Full-Stokes (FS) ice sheet models provide the most sophisticated formulation of ice sheet flow. However, their applicability is often limited due to the high computational demand and numerical challenges. To balance computational demand and accuracy, the so-called Blatter–Pattyn (BP) stress regime is frequently used. Here, we explore the dynamic consequences of using simplified approaches by solving FS and the BP stress regime applied to the Northeast Greenland Ice Stream. To ensure a consistent comparison, we use one single ice sheet model to run the simulations under identical numerical conditions. A sensitivity study to the horizontal grid resolution (from 12.8 to a resolution of 0.1 km) reveals that velocity differences between the FS and BP solution emerge below ∼ 1 km horizontal resolution and continuously increase with resolution. Over the majority of the modelling domain both models reveal similar surface velocity patterns. At the grounding line of the 79^∘ North Glacier the simulations show considerable differences whereby the BP model overestimates ice discharge of up to 50 % compared to FS. A sensitivity study to the friction type reveals that differences are stronger for a power-law friction than a linear friction law. Model differences are attributed to topographic variability and the basal drag, in which neglected stress terms in BP become important.

1 Introduction

The most comprehensive description of ice sheet flow is given by the full-Stokes (FS) equations. Such a model is considered the most accurate available, capable of describing highly dynamic ice sheets, including ice streams, ice shelves and grounding line migration. But this formulation is also the most computationally expensive, both from the numerical perspective and because FS effects only occur at higher resolutions.

Due to increasing computer power, the emergence of new generation ice sheet models (ISMs) in the last decades allowed for ice flow models to directly rely on FS equations or approaches with only a few simplifications. Although ice dynamical processes require an accurate representation of the stresses, e.g. in regions of high variability in basal topography and/or basal slipperiness (Gudmundsson, 2003; Hindmarsh, 2004), simplifications to FS are still necessary for long time integrations or for large ensemble modelling to keep the computing time affordable. For this reason, ice flow modelling has frequently relied on simplified mathematical models to balance computational costs and accuracy, such as the three-dimensional and computationally efficient Blatter–Pattyn (BP) formulation (Blatter, 1995; Pattyn, 2003). Although BP neglects several components of the full system (e.g. bridging stresses), it is considered valid in the majority of an ice sheet.

Current ice sheet projections within the ISMIP6 effort (Ice Sheet Model Intercomparison Project for CMIP6; Nowicki et al., 2016) cover a large diversity in model approaches and approximations to the Stokes flow (Goelzer et al., 2018; Seroussi et al., 2019; Goelzer et al., 2020; Seroussi et al., 2020), but only one contribution solved the FS problem (Seroussi et al., 2020, model UTAS_ElmerIce with a variable resolution between 4 and 40 km). Up to now, there is no clear understanding of whether FS simulations have the potential to narrow uncertainties in current sea-level predictions from the Greenland and Antarctic ice sheets. However, the ability to accurately simulate the behaviour of outlet glaciers has been cited in the reports of the Intergovernmental Panel on Climate Change as a prerequisite for reliable ice sheet projections (IPCC, 2013; Meredith et al., 2019; Oppenheimer et al., 2019).

In earlier ISMIPs (Ice Sheet Model Intercomparison Projects) the performance in simulating idealized ice stream and marine ice sheet dynamics by ISMs was explored: e.g. in ISMIP-HOM (Ice Sheet Model Intercomparison Project for Higher-Order Models; Pattyn et al., 2008), MISMIP2d (Marine Ice Sheet Model Intercomparison Project;  Pattyn et al., 2012), MISMIP3d (Pattyn et al., 2013), and MISMIP+ (Cornford et al., 2020). The ISMIP-HOM exercise focused on so-called higher-order (HO) and FS models and in particular reveals that the FS models show a smaller spread and hence are in better agreement with one another and with available analytical solutions. In the earlier MISMIP experiments, any systematic study of differences in solutions from Stokes models was not possible due to the limited number of FS models¹. The recent MISMIP+ effort reflects that the choice of the friction law plays a larger role compared to, for example, ice flow approximations (two FS models contribute to this exercise). From these synthetic scenarios, it yet needs to be assessed if discrepancies between FS and simpler models are minor compared to other uncertainties in ice sheet modelling given the huge additional computational cost when running FS models (about 10 times slower than BP; Larour et al., 2012, Table 3 therein).

In contrast to the idealized ISMIP experiments, the influence of FS on ice dynamical behaviour compared to simpler models was tested on realistic problems (e.g. Leysinger Vieli and Gudmundsson, 2004; Le Meur et al., 2004; Morlighem et al., 2010; Seddik et al., 2012, 2017; Favier et al., 2014). These exercises provide important insights into the transient response and indicate that FS produces different results compared to simpler models. For example, in most of the projection experiments conducted, FS tends to contribute less to sea-level rise than simpler models under identical projection scenarios (Seddik et al., 2012, 2017; Favier et al., 2014). Morlighem et al. (2010) concluded that treating ice flow of the fast-flowing Pine Island Glacier (PIG, Antarctica) with a steeply rising bed near the grounding line with FS is essential. However, assessing whether FS is urgently needed compared to simpler models is complicated because of many interacting processes. (i) Comparability is often limited by the use of different ISMs as it is not entirely clear how much of the differences is due to numerical treatments as, for example, different horizontal grids resolve the bed differently or to discretization schemes. (ii) Differences are also subject to, for example, assumptions on basal flow conditions or the initialization approach. (iii) Moreover, when transient simulations are forced with climate scenarios, the response involves numerous processes and interactions that make it difficult to separate causes.

Here, we intend to overcome model intercomparison shortcomings by performing consistent numerical experiments by solving the FS system or using the BP approximation. Our comparison focuses on FS and BP as the BP is used as a compromise for FS to balance computational costs and accuracy. The employed ISM includes the option to run the FS model or to activate the BP mode, allowing us to compare the results consistently as they are computed under the same numerical conditions. Here, we strive to fill the gap of FS and BP intercomparisons, stepping away from synthetic scenarios and performing a high-resolution real system application. We select two regions of the Northeast Greenland Ice Stream (NEGIS) as an investigation area. Both regions cover slow- and fast-flowing ice, as well as smooth and highly variable bed topography. The first domain focuses on the ice stream, while the second domain covers the NEGIS outlets at the 79^∘ North Glacier (79NG) and Zacharias Isbræ (ZI). Our study does not treat the transient behaviour; we rather aim to address ice dynamical differences between the stress regimes. As differences between FS and BP are expected to emerge in regions with steep bed gradients or high aspect ratios (height$/$length ratio ϵ), we step-wise increase spatial resolution from 12.8 km to a resolution of 0.1 km.

2 Theoretical background

2.1 The full-Stokes equations

In the case of ice sheets, the Reynolds number is very low (e.g. Fowler and Larson, 1978), and therefore the inertial term in the Navier–Stokes equations can be neglected. This approximation is often referred to as Stokes flow. The core of the FS equation system builds the mass balance and momentum equation for incompressible ice and is written as

$$\begin{array}{}\text{(1)}& {\displaystyle }\text{div}\mathit{v}=\mathrm{0},\text{(2)}& {\displaystyle }\text{div}\mathbf{t}={\mathit{\varrho }}_{\mathrm{i}}\mathit{g},\end{array}$$

with the density of ice ${\mathit{\varrho }}_{\mathrm{i}}=\mathrm{917}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, the three-dimensional velocity field $\mathit{v}=(v_x,v_y,v_z)$ in Cartesian coordinates, the gravitational acceleration vector pointing downward $\mathit{g}=-\mathit{g}{\mathit{e}}_z$ $(\mathit{g}=\mathrm{9.81}\mathrm{m}{\mathrm{s}}^{-\mathrm{2}})$ and the Cauchy stress tensor t. We split the Cauchy stress into a deviatoric part t^D and an isotropic pressure p:

$$\begin{array}{}\text{(3)}& \mathbf{t}={\mathbf{t}}^{\mathrm{D}}+p\mathbf{I},\end{array}$$

with $p=-\frac{\mathrm{1}}{\mathrm{3}}\mathrm{tr}\left(\mathbf{t}\right)$ and I being the identity tensor. The constitutive equation for the non-Newtonian fluid links strain rates to velocity gradients:

$$\begin{array}{}\text{(4)}& {\mathbf{t}}^{\mathrm{D}}=\mathrm{2}\mathit{\eta }\mathbf{D},\end{array}$$

where D is the strain rate tensor ($\mathbf{D}=\frac{\mathrm{1}}{\mathrm{2}}(\text{grad}\mathit{v}+(\text{grad}\mathit{v}{)}^T$)). The viscosity is given by the Glen–Steinemann flow law (Glen, 1955; Steinemann, 1954):

$$\begin{array}{}\text{(5)}& \mathit{\eta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}EA(T,W{)}^{-\mathrm{1}/n}{\dot{\mathit{\epsilon }}}_e^{(\mathrm{1}-n)/n},\end{array}$$

with the flow law exponent n=3, the enhancement factor E, the rate factor A depending on temperature T and microscopic water content W, and the effective strain rate ${\dot{\mathit{\epsilon }}}_e$ being the second invariant of the strain-rate tensor. To the latter, a small value of ${\dot{\mathit{\epsilon }}}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{30}}$ s⁻¹ is added to keep the term non-zero.

2.2 Boundary conditions

The upper surface is assumed to be traction free. The ice sheet base is subject to basal sliding according to a friction law for the basal shear stress ${\mathit{t}}_{\mathrm{b}}^{\mathrm{D}}$ in the tangential plane:

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}{\mathit{t}}_{\mathrm{b}}^{\mathrm{D}}& ={\mathit{\tau }}_{\mathrm{b}}=\left[\right(\mathbf{I}-\mathit{n}{\mathit{n}}^T\left){\mathbf{t}}^{\mathrm{D}}\right]\cdot \mathit{n}\\ & =-{\mathit{\beta }}^{\mathrm{2}}(\mathbf{I}-\mathit{n}{\mathit{n}}^T)\cdot \mathit{v}=-{\mathit{\beta }}^{\mathrm{2}}{\mathit{v}}_{\mathrm{b}},\end{array}\text{(7)}& {\displaystyle }\mathit{v}\cdot \mathit{n}=\mathrm{0},\end{array}$$

with the unit normal vector n pointing out of the ice and v_b being the velocity in the tangential plane at the base. Basal refreezing or melting is neglected. The basal drag parameter β² often includes the effective pressure to account for the presence of water lubricating the ice–bed interface, as well as bed properties such as roughness. For floating ice we use the same boundary condition but prescribe free slip (β²=0). The parameterization for β² and the remaining lateral boundary condition are given below for each experiment set-up.

2.3 Model classification

The mass and momentum balance (Eqs. 1 and 2) can be written in terms of velocity components and pressure p:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\displaystyle \frac{\partial v_x}{\partial x}}+{\displaystyle \frac{\partial v_y}{\partial y}}+{\displaystyle \frac{\partial v_z}{\partial z}}=\mathrm{0},\\ & {\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{2}\mathit{\eta }{\displaystyle \frac{\partial v_x}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(\mathit{\eta }{\displaystyle \frac{\partial v_x}{\partial y}}+\mathit{\eta }{\displaystyle \frac{\partial v_y}{\partial x}}\right)\\ & +{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial v_x}{\partial z}}+\mathit{\eta }\underset{\mathrm{¯}}{\underset{\mathrm{¯}}{{\displaystyle \frac{\partial v_z}{\partial x}}}}\right)-{\displaystyle \frac{\partial p}{\partial x}}=\mathrm{0},\\ & {\displaystyle \frac{\partial }{\partial x}}\left(\mathit{\eta }{\displaystyle \frac{\partial v_x}{\partial y}}+\mathit{\eta }{\displaystyle \frac{\partial v_y}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(\mathrm{2}\mathit{\eta }{\displaystyle \frac{\partial v_y}{\partial y}}\right)\\ & +{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial v_y}{\partial z}}+\mathit{\eta }\underset{\mathrm{¯}}{\underset{\mathrm{¯}}{{\displaystyle \frac{\partial v_z}{\partial y}}}}\right)-{\displaystyle \frac{\partial p}{\partial y}}=\mathrm{0},\\ & \underset{\mathrm{¯}}{\underset{\mathrm{¯}}{{\displaystyle \frac{\partial }{\partial x}}}}\left(\mathit{\eta }\underset{\mathrm{¯}}{\underset{\mathrm{¯}}{{\displaystyle \frac{\partial v_x}{\partial z}}}}+\mathit{\eta }\underset{\mathrm{¯}}{\underset{\mathrm{¯}}{{\displaystyle \frac{\partial v_z}{\partial x}}}}\right)+\underset{\mathrm{¯}}{\underset{\mathrm{¯}}{{\displaystyle \frac{\partial }{\partial y}}}}\left(\mathit{\eta }\underset{\mathrm{¯}}{\underset{\mathrm{¯}}{{\displaystyle \frac{\partial v_y}{\partial z}}}}+\mathit{\eta }\underset{\mathrm{¯}}{\underset{\mathrm{¯}}{{\displaystyle \frac{\partial v_z}{\partial y}}}}\right)\\ & +{\displaystyle \frac{\partial }{\partial z}}\left(\mathrm{2}\mathit{\eta }{\displaystyle \frac{\partial v_z}{\partial z}}\right)-{\displaystyle \frac{\partial p}{\partial z}}=-{\mathit{\varrho }}_{\mathrm{i}}g,\end{array}\end{array}$$

as well as the effective strain rate occurring in Eq. (5):

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\dot{\mathit{\epsilon }}}_e& =\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\partial v_x}{\partial x}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\partial v_y}{\partial y}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\partial v_z}{\partial z}}\right)}^{\mathrm{2}}\right.\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial v_x}{\partial y}}+{\displaystyle \frac{\partial v_y}{\partial x}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial v_x}{\partial z}}+\underset{\mathrm{¯}}{\underset{\mathrm{¯}}{{\displaystyle \frac{\partial v_z}{\partial x}}}}\right)}^{\mathrm{2}}\\ & {\left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial v_y}{\partial z}}+\underset{\mathrm{¯}}{\underset{\mathrm{¯}}{{\displaystyle \frac{\partial v_z}{\partial y}}}}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}\end{array}$$

The underlined terms are omitted in the BP approximation as explained below. The four equations (Eqs. 1 and 8) with the four unknowns (v_x, v_y, v_z and p) form the FS equation system. Various approximations to this set of equations exist ranging from shallow ice approximation (SIA) to different types of one layer and multilayer approximations for non-horizontal shear components of the stress tensor (Hutter, 1983; Hindmarsh, 2004; Bueler and Brown, 2009; Ahlkrona et al., 2013). The FS problem has a saddle point character which is difficult to solve (Benzi et al., 2005). In order to reduce the computational demand and derive a better-posed problem, Blatter (1995) and Pattyn (2003) developed the so-called BP approximation scheme which is valid for most of the ice sheet domain and widely used in ice flow models. The main assumption of their model is that the vertical component of the momentum balance is approximated as hydrostatic (i.e. $\partial t_{\mathrm{zz}}/\partial z={\mathit{\varrho }}_{\mathrm{i}}g$). This assumption is used to eliminate the pressure variable p from the FS equation system and reduces the FS problem with the two horizontal velocity components as unknown field variables. The further assumption that horizontal gradients of the vertical velocity are small compared to the vertical gradient of the horizontal velocity (i.e. $\partial v_z/\partial x\ll \partial v_x/\partial z$ and $\partial v_z/\partial y\ll \partial v_y/\partial z$) leads to a closed problem for the two horizontal velocity components. This approximation is equivalent to dropping the terms marked with an underline in Eqs. (8) and (9); it is termed LMLa in Hindmarsh (2004).

The assumptions to the BP scheme imply that the so-called bridging stresses (also known as vertical resistive stress; van der Veen and Whillans, 1989), i.e. the resistance to varying stress gradients in the direction of the ice flow, is neglected. Bridging effects are generally small and occur near the ice divide and at the grounding line (Pattyn, 2000). Since the BP schemes retain stresses of 𝒪(1) and 𝒪(ϵ) but delete stress terms of 𝒪(ϵ²) (Blatter, 1995), it is only valid to a certain aspect ratio and topographic variability. Once high-velocity gradients develop over short distances, BP may not provide an accurate solution.

The simplifications made by BP on the FS system reduces the computational demand by deriving a better-posed problem. By re-arranging the BP equations this system leads to a closed problem for the two horizontal velocity components v_x and v_y. In doing so, the BP approximations decouple the vertical velocity v_z and the pressure p from the full system. Both unknowns are recovered diagnostically by integrating the incompressibility equation (Eq. 1) and the vertical component of the momentum balance equation (Eq. 8), respectively. However, in our approach we follow an alternative way to align with the BP stress regime, thereby directly allowing us to compare FS and BP model simulations (see Sects. 3.1 and 3.2).

3 Model description and verification

3.1 Ice flow model

To solve the set of equations presented above (Eqs. 1 and 2 and boundary conditions), we use the COMmercial finite element SOLver (COMSOL) © Multiphysics (version 5.5, https://www.comsol.com/, last access: 1 May 2021). COMSOL provides an FS interface which is additionally manipulated in order to align with the BP stress regime. This is accomplished by discarding the corresponding terms (underlined terms in Eqs. 8 and 9) in the weak formulation but maintaining identical FS numerical details. Through this approach the BP stress regime does not further modify the default FS module; i.e. the BP version implemented in this study is not forming a multilayer model with the two horizontal velocity components as field variables. Therefore, our implemented BP approximation is hereafter termed BP-like. In doing so, both flow modes (FS and BP-like) are running through the same (FS) solving algorithm. Hereby we ensure that differences in the computed results due to numerical issues are largely eliminated. However, the disadvantage is that we do not benefit from the reduced complexity of the BP approximation. As a consequence of the avoided decoupling, $\partial v_z/\partial z$ in Eq. (9) is not replaced by $-\left(\partial v_x/\partial x+\partial v_y/\partial y\right)$ as commonly done in HO approximations (e.g. Pattyn, 2003).

3.2 Numerics

The mass balance and momentum equations constitute a non-linear equation system. Independent of the employed flow regime, a damped Newton method (which is COMSOL's default solver) is applied in a fully coupled manner to solve for the velocity vector v and the pressure p. The unstabilized Stokes equation (Eq. 8) is subject to the Babuska–Brezzi condition, which states that the basis functions for p have to be of lower order than for v. This issue could be circumvented with suitable stabilization techniques. Here, we test the numerical robustness of our results by employing two different discretization schemes. (1) We use P2+P1 Lagrange elements which consist of quadratic basis functions for the velocity and linear for the pressure. (2) The second scheme uses linear elements for both the velocity components and the pressure field (P1+P1 Lagrange elements), and numerical stabilization is achieved with streamline diffusion (Galerkin least-squares (GLS); Hauke and Hughes, 1994). Due to the equal-order interpolation, the latter is computationally more efficient but less accurate and sensitive to the GLS stabilization parameters (Helanow and Ahlkrona, 2018). In the following, the schemes will be referred as P2P1 and P1P1GLS, respectively.

To further test the numerical robustness we compare two different implementations (called “strong” and “weak”) for the Dirichlet condition (Eq. 7) at the ice base. The strong imposition applies the constraints pointwise at each node and, therefore, modifies the underlying variational formulation. The weak method is rather a component of the variational formulation by adding exterior facet integrals. For the latter, we use the so-called Lagrange multiplier method (Babuska, 1973; Verfürth, 1986; Urquiza et al., 2014) in which the Dirichlet condition is multiplied with the Lagrange multiplier λ.

The flow regime equations are forming a saddle point problem, for which either a direct solver or special preconditioning is needed. The verification tests ISMIP-HOM (see Sect. 3.3) are solved with the direct solver MUMPS (https://graal.ens-lyon.fr/MUMPS/; last access: 1 January 2022). For the NEGIS region experiments (Sect. 4), the degrees of freedom (DOFs) increase rapidly with increasing spatial resolution, and a direct solver would be uncommon for such large numbers of DOFs. For the larger problems we rely on the iterative GMRES solver (Saad, 2003) which is accelerated with appropriate preconditioners. Simulations that make use of the strong imposition rely on a domain decomposition solver with an overlapping additive Schwarz method (ASM; Widlund and Toselli, 2005) which requires less computation time and working memory compared to MUMPS. Unfortunately, the ASM preconditioner shows a very high computational demand for simulations that employ the weak imposition. Since the involved Lagrange multiplier induces a zero on the diagonal of the system matrix, several common preconditioners are not applicable. Therefore, we employ the Vanka algorithm (Vanka, 1986; John and Matthies, 2001) which is specifically designed for large indefinite problems with a saddle point character. Based on our simulations, we found Vanka to be very memory efficient and computationally fast for large problems, although it requires more Newton and GMRES iterations compared to ASM (see Sect. 4 and Tables 1 and 2).

Regarding the desired consistency between the simulation results, we must clarify the different usages of the linear solvers: the coarser models show a poor performance with the ASM and Vanka preconditioner. At the same time, the higher resolutions run out of memory with the MUMPS solver. Depending on modelling domain and employed discretization, MUMPS runs out of memory at around 5 million DOFs. For the largest DOF model that could be solved with the direct solver MUMPS on our cluster system, we also performed a simulation with the iterative linear solvers ASM and Vanka. A comparison of both simulations reveals differences well below 10⁻⁷ m a⁻¹ in the surface velocity, $v_{\mathrm{s}}=|\mathit{v}{|}_{\mathrm{s}}$; differences between MUMPS and Vanka are below 10⁻³ m a⁻¹. Consequently, we assume that all applied linear solvers provide comparable results.

3.3 Verification by benchmarks

The ISMIP-HOM benchmark targets a comparison for three-dimensional ice flow models arising from basal topography undulations, from basal sliding and by changing the length scale. From the ISMIP-HOM set-up (Pattyn et al., 2008) we select the three-dimensional and non-transient experiments Experiment A (“ice flow over a bumpy bed”) and Experiment C (“ice-stream flow I”). We rely on this set-up to verify our flow implementations, in particular the BP-like implementation. Experiment A is a parallel-sided slab with a series of sinusoidal oscillations in the basal topography with a frequency of $\mathit{\omega }=\mathrm{2}\mathit{\pi }/L$, with L being the domain length. In this set-up no basal sliding is prescribed (i.e. v_b=0). Experiment C is a parallel-sided slab with a series of sinusoidal oscillations on the basal drag parameter β² (Eq. 6) with the frequency ω. At the lateral margins, periodic boundary conditions are applied. The length scale is subsequently reduced as $L=\mathrm{160},\mathrm{80},\mathrm{40},\mathrm{20},\mathrm{10}$ and 5 km so that full-Stokes effects can develop due to the higher aspect ratio $H/L$, with H being the mean ice thickness. The domain is discretized with $n_x=n_y=L/\mathrm{41}$ grid points in the horizontal plane. In the vertical direction, 15 layers refined to the base are chosen. We run the benchmarks with a length scale of 5 km with the linear (P1P1GLS) and second-order (P2P1) element functions. Additionally, Experiment C is run with the strong and weak imposition of the basal friction law. All other length scales are run with P1P1GLS and for Experiment C with the strong imposition of the friction law. The ISMIP-HOM experiments use an ice density of ${\mathit{\varrho }}_{\mathrm{i}}=\mathrm{910}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. For a detailed description of the experimental set-up, please refer to Pattyn et al. (2008).

Our results for FS and BP-like show good agreements with the corresponding models for FS and non-FS (nFS; nFS includes BP) for the A and C test cases conducted (see Figs. S2 and S3). However, the original exercise reveals distinct differences between FS and all other approximations on short length scales (Fig. 1). In particular, for Experiment C the “double-peak” in the surface velocity is present in the FS solution but absent in the BP-like regime, which is consistent with the original ISMIP-HOM results. Independent of the element functions used and friction law imposition, we found a good agreement with the original ISMIP-HOM contributions and within our own model versions. Generally, the P2P1 model reveals slightly higher surface velocities than P1P1GLS. The different friction law implementations show very minor differences. From the resulting behaviour, the implementations of the FS and BP-like stress regimes are verified.

Figure 1 Results of the ISMIP-HOM experiments A (a) and C (b) for the length scale L=5 km. Surface velocity component v_x (m a⁻¹) at $y=L/\mathrm{4}$. Values computed in this study for FS and BP-like are compared to the BP model “rhi2” and the FS model “oga1” from the original ISMIP-HOM benchmark (Pattyn et al., 2008). Please note that the BP-like P1P1GLS and “rhi2” BP solutions in Experiment A overlap on each other; in Experiment C the “rhi2” BP solution is overlaid by BP-like P2P1 strong and BP-like P2P1 weak.

4 Description of NEGIS region experiments

We select two different regions of the NEGIS (Fig. 2) to investigate the effect of FS and BP on its ice dynamics. The first region focuses on an ice-stream-like flow regime with abrupt changes in slow and fast flow of the southern branch (SB) of NEGIS (it is termed “ice stream”; see red line in Fig. 2). The second region also includes the western branch (WB) of NEGIS and the outlet regions of 79NG and ZI (it is termed “outlet”; see cyan line in Fig. 2). For the ice-stream experiment a region upstream of the grounding line and downstream of the ice divide is selected to avoid grounding line treatment and to focus on an ice-stream-like flow regime. The outlet region captures the grounding zone of 79NG and the marine terminus of ZI. As we rely on the BedMachine v3 dataset (BM; Morlighem et al., 2017), a small area at 79NG becomes afloat. In both set-ups, the geometry of the NEGIS regions is laterally defined along flow lines that were deduced from the MEaSUREs surface velocity dataset (Joughin et al., 2016, 2018).

Bed and surface topographies are taken from BM and the GIMP dataset (Howat et al., 2014), respectively. As consistency is the primary target of the model set-up to study the response of both stress regimes to different initial and boundary conditions, we do not perform an initial relaxation of the free surfaces. A relaxation run based on, for example, the BP-like scheme and subsequent simulation runs with relaxed geometry with BP-like and FS would likely experience a “shock” (Goelzer et al., 2018) for the FS model because the geometry is consistent to the BP-like physics (and in the same way also for an FS relaxation run). Performing relaxation simulations for BP-like and FS individually would result in different initial geometries and also different initial conditions (e.g. velocity). This would be an undesired effect for our consistency approach which would make a clear comparison difficult. Here, we followed the strategy to have physical parameters and the geometry input equal in the FS and BP-like to allow for a better comparison.

Figure 2 Overview of the NEGIS region and modelling domain. (a) Observed surface velocities (Joughin et al., 2016, 2018). (b) Bed topography (Morlighem et al., 2017). The red and cyan lines in (a) and (b) delineate the ice-stream and outlet modelling domains, respectively. White lines in (a) indicate flux gate locations at 79NG and ZI (see Fig. 9 and related text below). WB and SB indicate the western and southern branches of the NEGIS. Symmetry, free-slip and stress in (b) indicate the boundary condition types used as lateral margins. Background image is a RADARSAT Mosaic (Joughin, 2015; Joughin et al., 2016).

In the NEGIS experiments we consider a Budd-like friction law (Budd et al., 1979) with a linear (m=1) and non-linear relationship (m=3). The basal drag parameter β² in the basal boundary condition (Eq. 6) is assumed to be

$$\begin{array}{}\text{(10)}& {\mathit{\beta }}^{\mathrm{2}}=k^{\mathrm{2}}{N}^{\mathrm{1}/m}|{\mathit{v}}_{\mathrm{b}}{|}^{\mathrm{1}/m-\mathrm{1}},\end{array}$$

where k is the basal friction coefficient found by an inversion technique (see below and Appendix A). The effective normal pressure is the difference between the ice pressure P_i and the basal water pressure P_w:

$$\begin{array}{}\text{(11)}& N=P_{\mathrm{i}}-P_{\mathrm{w}}.\end{array}$$

The basal water pressure is computed in marine parts, i.e. where the ice base is below the sea level (h_b<z_sl), according to Huybrechts (1992):

$$\begin{array}{}\text{(12)}& {\displaystyle }{P}_{\mathrm{i}}=-(\mathbf{t}\cdot \mathit{n})\cdot \mathit{n},\text{(13)}& {\displaystyle }{P}_{\mathrm{w}}=-min({\mathit{\varrho }}_{\mathrm{w}}g{h}_{\mathrm{b}},\mathrm{0}),\end{array}$$

where ${\mathit{\varrho }}_{\mathrm{w}}=\mathrm{1023}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ is the density of the ocean water and h_b the bed elevation. Here, P_i is the normal stress instead of the approximated ice overburden pressure. The assumptions made about the water pressure imply that the base is perfectly connected to the ocean at any location in the domain that is below sea level. Overall, this might be incorrect, but it forms an appropriate assumption in the absence of an additional hydrological model.

Solving a subset of an ice stream poses unknown boundary conditions in the interior of the ice sheet. A simple approach would be prescribing the measured surface velocities as a depth-averaged velocity profile. However, we choose boundary conditions that are free to adjust during the solution process (Fig. 2b). Inflow boundaries are chosen to have a symmetry boundary condition, which implies $\mathit{v}\cdot \mathit{n}=\mathrm{0}$ and vanishing shear stresses. A free-slip condition is chosen at lateral along-flow boundaries and a normal stress condition at the outflow boundaries. Dependent on the set-up, the outflow boundaries are land- or marine-terminating fronts of the glacier or located within the ice. For those outflow boundaries we chose stress boundary condition $\mathbf{t}\cdot \mathit{n}=-F_{\mathrm{0}}\mathit{n}$. For outflow boundaries located within the ice (e.g. at 79NG) we choose $F_{\mathrm{0}}={\mathit{\varrho }}_{\mathrm{i}}g(z_{\mathrm{s}}-z)+\mathit{\delta }$; for land- or marine-terminating fronts (e.g. ZI) we set $F_{\mathrm{0}}=min\left({\mathit{\varrho }}_{\mathrm{i}}g\right(z_{\mathrm{s}}-z),\mathrm{0})+\mathit{\delta }$, where z_s is the surface elevation and z the vertical coordinate. Here, δ is an additional stress that is tuned to match the observed velocities at those boundaries. The additional stress is set to zero except at 79NG and ZI, where δ is 0.5 and 0.3 MPa at 79NG and ZI, respectively.

Dependent on the set-up, we conduct experiments ranging from a horizontal mesh resolution of l=12.8 to 0.1 km. Please note that l=0.15 km is the nominal resolution of the BM dataset. The horizontal grid sizes span a typical range of resolutions used in ice sheet projections (e.g. Aschwanden et al., 2019; Goelzer et al., 2020; Seroussi et al., 2020) and beyond. The three-dimensional meshes comprise between 1960 and 32 864 970 prism elements, respectively, resulting in 5280 to 72 416 740 DOFs for the mesh cases. For a number of simulations, the enhancement factor in Eq. (5) is varied as E=0.1, 1, 3 and 6 to investigate whether simulation differences stem from basal sliding (small E) or internal velocity (large E). In the following E=1 will be referred as the “classical case”. The upper and lower values of E are extreme but not too far away from currently used values (Ma et al., 2010; de Boer et al., 2015). To test the sensitivity on the assumed friction law we run simulations with m=1 and m=3 (Eq. 10). The numerical robustness is tested by the aforementioned numerical discretization schemes and the implementations of the basal boundary condition. We do not run all possible combinations in all set-ups as it is computationally too expensive. Overall, computational times varied between a few seconds and several hours. A broad overview of computational costs is given in Tables 1 and 2. If not mentioned otherwise all set-ups employ 11 vertical layers (refined to the base), which is sufficient compared to higher resolutions (Fig. S1).

Figure 3 Overview of input parameters and exemplary simulation results of the ice-stream NEGIS region. The upper row shows slope of bed topography (a), basal friction coefficient for E=1 inferred with ISSM (b) and simulated surface velocity with FS, E=1, m=1 and with FS P1P1GLS strong (c) for 6400 m horizontal resolution. The lower rows show the same fields but for 100 m horizontal resolution. Background image is a RADARSAT Mosaic (Joughin, 2015; Joughin et al., 2016).

For the parameters A(T,W) (Eq. 5) and k² (Eq. 10) we utilize external products and capacities.

1. To rely on a realistic thermal state but avoiding an intensive thermal spin-up, we make use of an enthalpy field from a palaeo-climatic spin-up simulated by the ice sheet model SICOPOLIS with SIA (Rückamp et al., 2019). The field is bilinearly interpolated to the COMSOL meshes. The temperature-dependent part of the rate factor in Eq. (5) for cold ice is computed following Cuffey and Paterson (2010, Table 3.4 therein), and the water-content-dependent part for temperate ice is computed following Lliboutry and Duval (1985).

2. We make use of the inversion capability in the Ice-sheet and Sea-level System Model (ISSM; Larour et al., 2012) to obtain a spatially varying basal friction coefficient k² (Figs. A1 and A2, see Appendix A). The ISSM inversion is performed with l=0.8 km with BP and for the individual E values. Similarly, the SICOPOLIS temperature field is interpolated to the three-dimensional ISSM mesh. Subsequently, the inferred fields for k² are bilinear interpolated to the COMSOL meshes.

Table 1 Overview of mesh characteristics and consumed computational resources of the ice-stream region. Number of non-linear (nlin) iterations (damped Newton's method), number of linear (lin) iterations, computational time and physical RAM usage are listed for E=0.1, m=1, FS and P1P1GLS strong; values do not vary much among other E and stress regime settings. Solver and cluster settings (OpenMP (OMP) shared-memory threads and distributed-memory tasks (MPI)) are kept constant between each resolution and stress regime. Computational time includes the solver procedure (e.g. assembling the stiffness matrix, load vector) and outside routines (e.g. construction of mesh, saving results). Generally, COMSOL is run in hybrid mode, meaning the process is distributed to a user-defined number of OMP threads and MPI tasks. Without testing the performance thoroughly we achieve a good performance by OMP parallelism per physical node and MPI communication between physical nodes (the code is simply run in OMP mode if MPI is equal to one). The listed RAM usage refers to one MPI task. Please note that the cluster settings (MPI and OMP) have been chosen in a way that the job returns a solution within the cluster “wall time” but do not indicate any performance scaling.

^a Simulations are run on the AWI cluster Cray CS400. Each production node is equipped with 64 GB RAM and 2×18 core CPUs (Intel Xeon Broadwell CPU E5-2697).
^b Simulations are run at the North German Supercomputing Alliance (HLRN) cluster Lise. Each production node has 362 GB RAM and 2×48 core CPUs (Intel Cascade Lake Platinum 9242 (CLX-AP)).
^c Domain decomposition solver with an overlapping multiplicative Schwarz method. Number of sub-domains are set equal to the number of MPI tasks.

Table 2 Overview of mesh characteristics and consumed computational resources of the outlet region similar to Table 1. Number of non-linear (nlin) iterations (damped Newton's method), number of linear (lin) iterations, computational time and physical RAM usage are listed for E=1, m=3, FS and P1P1GLS weak.

All simulations are run on the AWI cluster Cray CS400. Each production node is equipped with 64 GB RAM and 2x18 core CPUs (Intel Xeon Broadwell CPU E5-2697).

5 Results of the NEGIS region experiments

5.1 Ice-stream region

An overview of input parameters and FS simulation results of the ice-stream region for l=6400 and 100 m is provided in Fig. 3. The slope of the bed topography (Fig. 3a, d) and the friction coefficient (Fig. 3b, e) of the finer resolution reveals more small-scale features and stronger amplitudes, particularly a rapidly varying bed. Despite both compared resolutions relying on the same inferred basal friction coefficient and input geometry, the interpolation to the computational grids produces striking differences. As a consequence the simulated surface velocity fields differ quantitatively (Fig. 3c, f). Both resolutions reproduce the NEGIS to some degree, but the finer resolution features an extended fast flow region (≳350 m a⁻¹) and more pronounced shear margins (i.e. sharper velocity gradients). As the ISSM simulations inferring the friction coefficient are conducted with a 0.8 km horizontal grid resolution and with the BP ice flow approximation, the COMSOL results among all applied resolutions and flow regimes are not expected to match the observed velocity field. The purpose of the transferred friction coefficient from the ISSM inversion to COMSOL is to roughly reproduce the observed NEGIS flow pattern. A better fit to observations could certainly be achieved by performing the inversion for each setting separately; however, our experiments are designed to focus on stress regime inter-comparisons and not on an assessment of observations.

In the following, the FS solution is considered as the exact solution. Therefore, our analysis builds on presenting differences of the BP-like solution to FS. The differences of surface velocities, $\mathrm{\Delta }{v}_{\mathrm{s}}^{\text{FS}-\text{BP-like}}=v_{\mathrm{s}}^{\text{FS}}-v_{\mathrm{s}}^{\text{BP-like}}$, between FS and BP-like are displayed for the 100 m grid resolution in Fig. 4. Spatial differences reveal a different pattern among the E values. For stiffer ice (E=0.1) the BP-like solutions exhibit higher flow speeds which are almost confined to the fast flow region (Fig. 4a); this area partly coincides with a smooth bed and surface topography (see Fig. 3). Overall, discrepancies between FS and BP-like increase with higher velocities. The largest impact of using a BP-like solution results from the surface flow velocities, which are up to 43 m a⁻¹ faster (relative difference is about 19 %) compared to FS. For the “classical case” (E=1) the pattern is similar but with reduced magnitudes. In this set-up, the BP-like velocities are up to 19 m a⁻¹ higher (relative difference is about 6 %). With softer ice (E>1), the pattern in surface velocity differences changes, and, eventually, for very soft ice (E=6) differences are more pronounced in regions with higher topographic variability (compare Fig. 3) and exhibit a rippled pattern. Velocity differences are scattered around the FS solution by ±80 m a⁻¹; the spread in differences is more pronounced for softer ice (see scatter plots in Fig. 4).

Figure 4 Difference of surface velocity $\mathrm{\Delta }{v}_{\mathrm{s}}^{\text{FS}-\text{BP-like}}=v_{\mathrm{s}}^{\text{FS}}-v_{\mathrm{s}}^{\text{BP-like}}$ for a resolution of 100 m. (a–d) Differences for E=0.1, 1, 3 and 6, respectively. Background image is a RADARSAT Mosaic (Joughin, 2015; Joughin et al., 2016).

Based on the spatial differences, we compute a spatially averaged surface velocity for all parameter settings to obtain an integrative overview (Fig. 5). The analysis generally reveals that both the absolute spatially averaged surface velocity (Fig. 5a–d) and relative spatially averaged surface velocity differences (Fig. 5e–h) start to diverge somewhere below 1.6 and 0.8 km horizontal grid resolution (which is in the range of the mean ice thickness of 1330 m). Away from this threshold towards finer resolutions, solution differences increase between the stress regimes. For the finer resolutions BP-like produces higher mean surface velocities compared to FS. At the finest resolution (below the nominal resolution of the BM dataset), the mean velocity seems to reach a converged state for the BP-like stress regime (less pronounced for E=0.1), but FS reveals a small dip towards lower velocities. The tendency of the diverging relative differences at the finest resolution is therefore likely a result from the still changing FS solution under grid refinements. The obtained tendency below 150 m is intriguing as the differences solely arise from grid refinements and not due to the input geometry.

Relative maximum differences between FS and BP-like are obtained for the stiff ice case (E=0.1). A discussion on why stiff ice is more sensitive will come in Sect. 5.5. The spatially averaged ice flow velocity difference is up to 7 m a⁻¹ for E=0.1; the relative error between FS and BP-like is up to ∼ 5.8 %. The resulting difference for E=1 shows the same trend but is of lower magnitude (up to ∼ 1.5 %). With softer ice (E>1) the differences decrease. Interestingly, both stress regimes reveal an intermediate peak in mean surface velocity at l=1600 m over all E values.

Figure 5 Integrative overview of spatially averaged surface velocities from FS and BP-like stress regimes. The upper row (a–d) shows surface velocity mean for E=0.1, 1, 3 and 6, respectively. The lower row (e–h) shows the relative difference of BP-like to FS for E=0.1, 1, 3 and 6, respectively. The lighter background indicates the local maximum and minimum relative differences. Please note the logarithmic scale of the x axis and the different scale of the y axis in the upper row.

5.2 Outlet region

As for the ice-stream region, an overview of input parameters and FS simulation results of the outlet NEGIS region are compared for l=6400 m and 150 m with P1P1GLS weak and m=3 (Fig. 6). Similarly, the slope of the bed topography (Fig. 6a, d) and the friction coefficient (Fig. 6b, e) of the finer resolution reveals more small-scale features and stronger amplitudes. In both resolutions, the simulated surface velocity fields agree fairly well and reproduce the NEGIS (Fig. 6c, f). However in shear margins and gradients in the along-flow direction are more pronounced at higher resolution.

Figure 6 Overview of input parameters and simulation results of the outlet NEGIS region. The upper row shows (a) slope of bed topography, (b) basal friction coefficient for E=1 and (c) simulated surface velocity with FS P1P1GLS weak and m=3 for 6400 m horizontal resolution. The lower rows show the same fields but for 150 m horizontal resolution. Background image is a RADARSAT Mosaic (Joughin, 2015; Joughin et al., 2016).

Differences in surface velocities between FS and BP-like are displayed for the 150 m grid resolution, E=1 and P1P1GLS in Fig. 7. Overall, spatial differences reveal a similar pattern among the employed friction exponents and friction implementations. The largest differences occur in the vicinity of the 79NG grounding line (up to −570 m a⁻¹, see Table 3). For a power friction law (m=3) the differences are larger compared to the linear friction law (m=1). The P1P1GLS weak set-up produces slightly larger differences than P1P1 strong.

Figure 7 Difference of surface velocity $\mathrm{\Delta }{v}_{\mathrm{s}}^{\text{FS}-\text{BP-like}}=v_{\mathrm{s}}^{\text{FS}}-v_{\mathrm{s}}^{\text{BP-like}}$ in the outlet domain for a resolution of 150 m and the P1P1GLS discretization. (a, b) Differences for m=1 and the strong and weak friction implementations, respectively. (c, d) Differences for m=3 and the strong and weak friction implementations, respectively. Please note the different colour bar compared to Fig. 4. Background image is a RADARSAT Mosaic (Joughin, 2015; Joughin et al., 2016).

Table 3 Maximum absolute velocity differences $(max(\left|\mathrm{\Delta }{v}_{\mathrm{s}}^{\text{FS}-\text{BP-like}}\right|)$, unit is m a⁻¹) and their corresponding relative velocity differences ($max\left(\right|\mathrm{\Delta }{v}_{\mathrm{s}}^{\text{FS}-\text{BP-like}}\left|\right)/v_{\mathrm{s}}^{\mathrm{FS}}$, in brackets, unit is  %) for different areas in the outlet domain. Results are listed for l=150 m, P1P1GLS and different set-ups for the basal boundary condition.

We compute relative spatially averaged surface velocity differences for all parameter settings to obtain an integrative overview (Fig. 8). The analysis discloses the clear trend that BP-like simulations lead to higher velocities than FS on average. Overall, the differences increase with higher resolution independent of the employed numerical characteristics. Interestingly, differences between FS and BP-like are more pronounced when using a power friction law compared to a linear friction law (up to 17 % compared to 8 %, respectively). Furthermore, the weak implementation produces larger differences compared to the strong scheme. The P2P1 discretization reveals differences that are in between P1P1GLS strong and P1P1GLS weak results. Differences still tend to increase further below the highest employed resolution. Notably, the area-averaged differences are higher in the outlet region than in the ice-stream region (see Fig. 5f).

Figure 8 Relative surface velocity differences of area-averaged surface velocities between the FS and BP-like stress regimes in the outlet domain. (a) Linear friction law (m=1). (b) Non-linear friction law (m=3). The blue and red circles show two simulations from a “FS without (wo) bridging stresses” and LTSML-like regimes simulated with P1P1GLS weak (see “Discussion”). For these two stress regimes, differences are calculated relative to FS.

5.3 Impact on ice discharge

To estimate the impact on ice discharge, we calculated the ice mass flux along two flux gates of 79NG and ZI (see white lines in Fig. 2a). The flux gate of 79NG corresponds to the grounding line of the BM dataset, while the flux gate of ZI was shifted approx. 3 km upstream from the calving front. The ice mass flux is calculated as $Q={\int }_{\mathrm{\Gamma }}\overline{v}H{\mathit{\varrho }}_{\mathrm{i}}\mathrm{d}s$, where $\overline{\mathit{v}}$ indicates the depth-averaged ice velocity, and the usage of the upper-case Q indicates integration over the flux gate area Γ with the line segments ds. On average the total simulated flux for FS across 79NG's grounding line is 14.7 and 14.2 Gt a⁻¹ at the flux gate at ZI's terminus. Figure 9 displays relative ice mass flux differences between FS and BP-like. Overall, the BP-like model results in higher discharge than FS, whereas the differences show a different magnitude: at 79NG, BP-like ice discharge is up to 40 % (m=1) and 50 % (m=3) higher than FS. In contrast, a moderate difference of about 5 % in both friction laws is found for ZI.

Figure 9 Relative ice mass flux differences between FS and the BP-like stress regimes ($\mathrm{\Delta }{Q}^{\text{FS}-\text{BP-like}}=(Q^{\text{FS}}-Q^{\text{BP-like}})/Q^{\text{FS}})$. The location of the flux gates at 79NG and ZI is indicated in Fig. 2. On average the total simulated flux for FS at 79NG is 14.7 and at ZI 14.2 Gt a⁻¹. The triangles in (b) indicate simulations using P1P1GLS weak with a resolution of l=200 m but a re-gridded geometry from a simulation with l=1600 m (see Sect. 6).

5.4 Free surface evolution

In this section, we investigate how the velocity differences are represented in the solution of the ice surface position. The evolution of the free surfaces is given by the kinematic boundary condition

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial z_{\mathrm{s}}}{\partial t}}=v_x{n}_x+v_y{n}_y+v_z{n}_z+{\dot{a}}_{\mathrm{s}},\end{array}$$

where ${\dot{a}}_{\mathrm{s}}(x,y)$ is the accumulation–ablation function at the surface. Here, we assume no melting or accumulation, i.e. ${\dot{a}}_{\mathrm{s}}(x,y)=\mathrm{0}$, and the so-called emergence velocity (v⋅n) determines the surface elevation change, $\partial z_{\mathrm{s}}/\partial t$. Since we do not run a forward model the computed free surface evolution is representative for the first time step.

Figure 10 shows the emergence velocity for FS and the difference between FS and BP-like. The FS emergence velocity reveals a somewhat noisy pattern. Due to the hyperbolic nature of Eq. (14) numerical stabilization is necessary (Donea and Huerta, 2003; Riviere, 2008) for solving the free surface equation which certainly will smooth out sharp gradients and oscillations. However, the FS emergence velocity shows values in the range of $\pm \mathrm{50}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$ (higher at some isolated grid nodes), which is high but not very extreme. Considering values for surface slopes and velocities of the observational datasets (BM and MEaSUREs) of up to 0.1 and 1000 m a⁻¹, respectively, leads to values of about 100 m a⁻¹ for the horizontal part of the emergence velocity (i.e. v_xn_x+v_yn_y). Assuming an initial time step of 0.05 years the estimated maximum height change corresponds roughly to ±2.5 m. Assuming ongoing elevation changes $\partial z_{\mathrm{s}}/\partial t$ in that area of about 1 to 2 m a⁻¹ (derived from laser altimetry, pers. comm. Veit Helm, AWI, and older data published in Helm et al., 2014), the estimated maximum ratio between emergence velocity and height changes $(\mathit{v}\cdot \mathit{n}/(\partial z_{\mathrm{s}}/\partial t\left)\right)$ is about 25. A transient relaxation run would smooth the surface slopes. However, the high values are not critical as we are interested in the differences between FS and BP-like.

The emergence velocity differences between FS and BP-like reveals a higher emergence velocity in FS, particularly at the grounding line at 79NG. According to Eq. (14), BP-like would therefore compute a surface elevation that is lower than the FS surface within the first time step. Although we are not running a time-dependent simulation and, therefore, not capturing grounding line migration, we evaluate the buoyancy balance from the stationary set-up. The ice is assumed to be grounded if the normal stress at the ice base is larger than the sea water pressure, i.e. P_i>P_w (e.g. Durand et al., 2009). BP-like reveals a grounding line that lies further upstream (≈ 0.5 km) than the one calculated from FS. However both are still located within the hinge zone of 79NG.

Figure 10 Evaluation of the kinematic boundary condition (emergence velocity v⋅n) at the surface for l=150 m. Surface mass balance is assumed to be zero. The upper row shows (a) FS emergence velocity and (b) emergence velocity differences between FS and BP-like $(\mathit{v}\cdot \mathit{n}{|}_{\mathrm{FS}}-\mathit{v}\cdot \mathit{n}{|}_{\mathrm{BP}-\mathrm{like}})$ for m=1 and P1P1GLS strong. Panels (c) and (d) show the equivalent fields for m=3 and P1P1GLS weak. Background image is a RADARSAT Mosaic (Joughin, 2015; Joughin et al., 2016).

5.5 Relevant physical processes

The reason for the difference in the simulated stress regime is a complex composition of ice dynamic processes and how they are resolved in the ice flow model, such as the basal boundary condition and the internal deformation. In addition, basal drag is coupled to internal deformation via the effective viscosity, making it difficult to separate the two entangled processes. Nevertheless, we strive to give a flavour of the origins of the model differences by analysing key processes of ice dynamics below.

Figures 11 and 12 show the difference in characteristic quantities between FS and BP-like for the ice-stream and outlet domain, respectively. First of all we would expect that FS would lead to larger ice velocities than BP due to the effective deformation rate (Eq. 9). FS captures the entire stress tensor, which is equivalent to making ice softer on average for the same flow parameters (assuming that strain-rate components shared by FS and BP-like have the same magnitude). Indeed, the depth-averaged viscosity, $\overline{\mathit{\eta }}$, is basically lower for FS compared to BP-like (Figs. 11g–i and 12c, f). For E=0.1 the depth-averaged viscosity is lower in some larger areas, while for E=6 the pattern changes towards several smaller patches; the magnitude of $\overline{\mathit{\eta }}$ among the E values changes by one order.

However, since FS generally reveals lower velocities than BP-like, the softer ice in FS appears to be compensated by another process. A fundamental control on ice flow is the basal drag indicating the glacier's behaviour for fast flow. The spatial differences of basal drag inferred from the two models is generally larger for FS compared to BP (Figs. 11d–f and 12b, e). The spatial patterns generally coincide with the differences in surface velocity (Figs. 4a and 7). In the ice-stream region, the difference is largest for E=0.1 (up to 140 kPa between FS and BP-like). In this set-up, it seems that the basal drag is an essential control of the ice flow. With softer ice, the difference in the basal drag decreases. In the outlet domain, the basal drag difference is more pronounced for the weak implementation of the basal boundary condition with m=3 compared to the strong implementation with m=1. At the 79NG and ZI area, the maximum difference in basal drag is 316 kPa (99 %) and 121 kPa (73 %) for the weak implementation with m=3; for the strong set-up with m=1 it is 227 kPa (72 %) and 96 kPa (55 %), respectively.

Figure 11 Comparison of characteristic quantities for a resolution of l=100 m, m=1 and P1P1GLS strong at the ice-stream domain. Difference (FS minus BP-like solution) of basal velocities v_b (a–c) and basal drag τ_b (d–f). Relative difference of depth-averaged viscosity $\overline{\mathit{\eta }}$ (g–i) between FS and BP-like. The domain-averaged mean viscosity for FS and E=0.1, 1 and 6 are 10.8, 3.6 and 1.1×10¹⁴ Pa s, respectively; for BP-like they are 12.6, 4.1 and 1.3×10¹⁴ Pa s, respectively. The equivalent fields for E=3 reveal patterns and magnitudes that range from E=1 to E=6. Background image is a RADARSAT Mosaic (Joughin, 2015; Joughin et al., 2016).

Figure 12 Comparison of characteristic quantities for a resolution of l=150 m and E=1 at the outlet domain. The upper row (a–c) shows absolute difference of basal velocities v_b (a), basal drag τ_b (b) and relative difference of depth-averaged viscosity $\overline{\mathit{\eta }}$ (c) for m=1 and P1P1GLS strong. The domain-averaged mean viscosities for FS are 3.00×10¹⁴ and for BP-like 3.40×10¹⁴ Pa s. The lower row (d–f) shows the equivalent fields but for m=3 and P1P1GLS weak. The domain-averaged mean viscosity for FS is 3.11×10¹⁴ and for BP-like 3.39×10¹⁴ Pa s. Background image is a RADARSAT Mosaic (Joughin, 2015; Joughin et al., 2016).

Inspired by the findings of Gudmundsson (2003) and Hindmarsh (2004), deduced from idealized set-ups, we explore whether the criteria when a BP solution is invalid (BP invalid in regions of high slip ratio, high aspect ratio and high topographic variations) hold in our realistic set-up. In Figs. 13 and 14 we present two-dimensional histograms of the relative velocity difference between FS and BP-like for discrete bins of FS slip ratio $v_{\mathrm{b}}/v_{\mathrm{s}}$, the local aspect ratio, $\mathit{ϵ}=H/l$, and the bed slope, $\left|\mathrm{\nabla }{h}_{\mathrm{b}}\right|$ (see detailed line plots in Figs. S4–S7). Simulated differences emerge with increasing aspect ratio, ϵ. With increasing aspect ratio neglecting terms of 𝒪(ϵ²), i.e. $\partial v_z/\partial x$ and $\partial v_z/\partial y$, the BP-like (and also BP) model becomes problematic and the solution inaccurate. Similarly, the difference between FS and BP is most pronounced in regions with a high slip ratio (i.e. near-plug-flow conditions prevail), while with a higher resolution a larger number of grid nodes are involved in producing this error. At a first glance, the distribution of the relative error with respect to the bed slope reveals larger errors at smaller slopes. However, the distribution is certainly biased by large errors occurring in the fast-flowing and smooth trunk of the NEGIS region (see Fig. 2). By neglecting the bed slope bin from 0 to 0.05, the higher bed undulations lead to larger errors at higher resolutions; however, the bed slope seems to play a secondary role compared to, for example, the slip ratio as fewer grid nodes are involved. The analysis underlines that the source of stress regime deviations occurs in regions of high slip ratio, high aspect ratio and high topographic variability. A similar conclusion for mountain glaciers was drawn in Le Meur et al. (2004) by comparing FS and SIA.

Figure 13 Two-dimensional histograms for E=1, m=1, P1P1GLS strong and four different grid resolutions for the ice-stream region. Counts of grid nodes per relative velocity error bin and per (a) bed slope bin, (b) local aspect ratio bin and (c) FS slip ratio bin. The sizes of the triangles are normalized by the total number of horizontal grid nodes of each mesh.

Figure 14 Two-dimensional histograms for E=1, m=3, P1P1GLS weak and four different grid resolutions for the outlet region. Counts of grid nodes per relative velocity error bin and per (a) bed slope bin, (b) local aspect ratio bin and (c) FS slip ratio bin. The sizes of the triangles are normalized by the total number of horizontal grid nodes of each mesh.

6 Discussion

At the ice-stream region, the quantitative comparison of both stress regimes reveals that, particularly in extremes cases (E=0.1 and E=6), the BP-like solution produces higher surface velocities than FS (E=0.1) and a large spread around the FS solution (E=6). However, these extreme cases at such a high resolution are to date rarely used in the community. For E=1 (“classical case”) and E=3, we found distinct velocity differences, but those are very low (about 1.5 % in spatially averaged surface velocities). Although these velocity differences are small, they tend to increase further at the highest grid resolution used. We cannot rule out the possibility that at a higher resolution the disagreements would be more substantial. However, since FS is computationally expensive, it is satisfying to note that the differences between the approximated model and the exact solution are small in resolutions that are currently being run within the ice sheet modelling community. Therefore, it might be favourable to use a BP stress regime instead of FS for most ice-stream modelling studies.

However, once the outlet regions of NEGIS are included (i.e. 79NG and ZI), the analysis shows that differences between FS and BP-like are large in the 79NG grounding zone. This finding is consistent with the work by Morlighem et al. (2010) at PIG. The overestimation of BP-like surface velocities to FS in our study agrees with their necessity of a lower inferred basal drag for FS to reproduce observed velocities compared to BP. They attribute this behaviour to the developing bridging stresses in FS (t_zz about 2 % larger in FS compared to BP) due to the rising bed towards the grounding line, which in turn causes the reduction in basal drag. However, in our simulations a clear connection to bridging stresses could not be identified. Figure 15a and c show the relative differences of t_zz between FS and BP-like. Although the differences are in a similar range (±2 %) as in Morlighem et al. (2010), there is no clear trend that could explain the lower velocities in FS. However, at PIG the bed rises towards the grounding line, while at 79NG the bed mainly slopes down towards the grounding line. This geometric setting might control how the stresses develop. In our case, the differences mainly stem from the basal drag (Figs. 11d–f and 12b, e) – that is for most areas higher in FS compared to BP-like. The regions where BP-like simulates a lower basal drag than FS overlap with regions where the assumption that horizontal gradients of the vertical velocity are small compared to the vertical gradient of the horizontal velocity is invalid (Fig. 15b, d). These terms are of similar order, and therefore, dropping stress terms of 𝒪(ϵ²), i.e. $\partial v_z/\partial x$ and $\partial v_z/\partial y$, leads to a lower basal drag in the BP-like case. Consequently, we assume that basal sliding is overestimated in BP-like simulations.

Figure 15 Comparison of (a, c) relative difference of t_zz at the ice base between FS and BP-like and (b, d) ratio of vertical gradient of the horizontal velocity to horizontal gradients of the vertical velocity $\left(\right|(\partial v_x/\partial z)/(\partial v_z/\partial x)|+|(\partial v_y/\partial z)/(\partial v_z/\partial y)\left|\right)$ in absolute terms from the FS solution. Panels (a) and (b) show results for l=150 m, E=1, m=1 and P1P1GLS strong. Panels (c) and (d) show results for l=150 m, E=1, m=3 and P1P1GLS weak. Background image is a RADARSAT Mosaic (Joughin, 2015; Joughin et al., 2016).

In order to highlight that bridging stresses play a minor role in our set-up, we have run two other stress regime versions (executed with P1P1GLS weak and up to a resolution of l=400 m). The first version keeps the bridging stresses in the third equation of the momentum balance (Eq. 8; this is similar to LTSML (van der Veen and Whillans, 1989; Hindmarsh, 2004; Kleiner and Humbert, 2014)); the second version is an FS model without (“wo”) bridging stresses. The area-averaged results reveal that LTSML is close to BP-like, while “FS wo bridging stresses” is rather similar to FS (Fig 8b).

In our diagnostic simulations, we compute lower ice mass fluxes with FS compared to BP-like. This supports the tendency that the sea-level rise estimate of FS models is lower than the one of simpler models (Seddik et al., 2012, 2017; Favier et al., 2014). However, we observe a different response at 79NG and ZI. At 79NG, ice discharge differences are up to 50 %, while at ZI, ice discharge differences are about 5 %. Also, at 79NG, differences seem to increase below the highest resolution employed here, while at ZI, the anticipated discharge difference increase appears to be moderate. A thorough analysis of why these two regions behave differently revealed no significant differences between ice dynamic processes (e.g. basal drag shows a similar decrease). However, a potential source might be the geometric setting; e.g. the data coverage of available in situ retrieved data leads to different topographic variability in both regions. To test whether the response behaviour is an effect of topographic variability, we re-run a simulation with P1P1GLS weak, m=3 and l=200 m but using the geometry from a re-gridded l=1600 m simulation. Re-gridding from a coarser geometry acts as smoothing and leads to a less variable geometry. The simulated discharge reveals that 79NG and ZI respond somewhat similarly (see triangles in Fig. 9b). However, whether the different response at 79NG and ZI with the high-resolution bedrock stems from in situ data coverage or the fact that ZI has generally a smoother geometry than 79NG must be examined in further studies.

We tested the numerical robustness of the underlying finite element method by employing different discretization schemes and analysing how the non-penetration condition at the ice base is enforced. In general, the different schemes show a similar velocity difference pattern between FS and BP-like. For m=1 and m=3 the area-averaged velocity differences for l=150 m are within −6.1 % to −7.6 % and −13.1 % to −16.2 %, respectively. Therefore, the choice of the employed friction law plays a larger role than the choice of numerical characteristics tested in this study.

However, some tests regarding the numerical robustness are missing. Because we keep the number of layers constant over the entire modelling domains, large variations in the aspect ratio occur when refining the mesh. The element aspect ratio is critical for the numerical error, and the classical BP system would suffer less under bad element aspect ratios (Larour et al., 2012, Sect 5.1.3 therein). In order to assess the impact of element aspect ratios, we carried out a study in which we varied the number of vertical layers from 3 to 11. For this test we re-run the simulations up to a resolution of l=400 m with the P1P1GLS strong discretization scheme and m=3 applied to the outlet domain. The area-averaged surface velocity differences between FS and BP-like are very small (Fig. S8). This demonstrates that the comparison is not suffering under bad element aspect ratios and that the numerical errors behave similarly in FS and BP-like.

The basal boundary conditions are retrieved by inverting for the basal friction coefficient, with BP employing a resolution of l=800 m. That means that the boundary conditions are not consistent with FS or with other resolutions. In order to study whether this inconsistency influences the results, we repeat a set of experiments with two different fields for the friction coefficient. To this end, we run two ISSM friction inversions with (1) E=1, l=200 m, m=3 and BP and (2) E=1, l=800 m, m=3 and FS and feed the inferred friction fields to COMSOL. The area-averaged differences between FS and BP-like of the three set-ups show very minor differences (Fig. S9).

By implementing FS and BP-like solvers using the COMSOL finite element package, we have developed a useful tool for exploring differences between FS and BP-like without complications due to differing numerics. The developed model is also very flexible in order to test various ice flow components (e.g. by simply including bridging stresses to convert the BP-like model to a LTSML-like model). In our model set-up, we solve for the same physical problem as in the “classical” BP scheme but without taking advantage of the reduced number of equations to be solved. Regarding the model physics, our results compare well with other BP implementations, but our BP-like model is more expensive. Therefore, the developed BP-like model is not intended for long time integrations or large ensembles.

One limitation of our study is that the simulations are not prognostic; i.e. we have not investigated how the solution differences propagate to or interact with other components of an ice sheet model, e.g. by coupling to the ice thickness evolution. The calculation of ice discharge and the emergence velocity demonstrates that the stress regimes show distinct discrepancies in a key process, but whether these differences are accumulated in the glaciers mass balance over transient runs, eventually leading to substantially different estimates of ice mass changes, or fade out must be investigated in future studies.

A process not considered here is the thermo-mechanical coupling. For instance, within a thermo-mechanical model, the feedback of temperature on the ice viscosity, in turn, changes the ice flow. Additionally, there are several feasible processes of how the stress regime differences can modify the thermal regime: (1) the different representation of the basal drag presumably generates a different amount of frictional heating, (2) the different vertical deformation produces a different internal heat source rate, and (3) the internal advection scheme leads to different heat transport. Though the role of these feedbacks is not explored here, there are indications that small initial differences become much larger over long time integrations (Zhang et al., 2015).

7 Conclusions

We compared two approaches to represent ice flow dynamics, namely the FS and the BP-like formulation, for two different regions of the NEGIS. Both stress regimes are implemented within the same dynamic ice flow model to enable a consistent comparison in which numerical differences are largely eliminated. Our comparison experiments between FS and BP-like stress regimes reveal that BP-like tends to overestimate surface velocities and the resulting ice discharge. The models disagreements start at a horizontal resolution of around 1000 m and continuously increase with finer horizontal resolution (down to 100 m).

In the ice-stream region, when considering a rheology with much softer ice (E=6), an enhanced spread in surface velocity differences between FS and BP-like was modelled. In contrast, stiffer (E=0.1) ice leads to enhanced ice flow velocities in BP-like compared to FS. The classical rheology case (E=1) reveals very moderate differences between FS and BP-like. However, once the outlet regions of NEGIS are included, i.e. 79NG and ZI, ice flow treatment with FS is essential, and important surface velocity differences occur (up to 570 m a⁻¹ at the grounding line of 79NG). Numerical characteristics play a minor role compared to the choice of the friction type: model differences are particularly important when a power-law friction is applied (as opposed to a linear approach). The basal drag and topographic undulations are identified as the general causes for the model disagreements.

Whether FS will reduce uncertainties of future sea-level projections must be examined further in subsequent transient studies. Although computer power is increasing, transient simulations in FS with a horizontal resolution of about 100 m in dynamic, relevant regions are still challenging. Our diagnostic simulations and previous studies indicate that FS and BP differ significantly in regions with a grounding zone, and results from non-FS models should therefore be viewed with caution.

Appendix A: Inversion with ISSM

The simulations presented make use of a basal friction coefficient, k², that is inferred by an inversion method. For the inversion of the basal friction coefficient, we operate the Ice-Sheet and Sea-level System Model (ISSM; Morlighem et al., 2010; Larour et al., 2012), an open source finite element flow model appropriate for continental-scale and outlet glacier applications. The set-up is as described above; i.e. (1) domain outline and geometry dataset is similar as in the COMSOL simulations (see Fig. 2), (2) ice rheology is taken from SICOPOLIS spin-up (Rückamp et al., 2019) and interpolated onto the computational mesh, and (3) the friction law follows the form in Eq. (6). But for the ease of computational time, model calculations are performed on a structured finite element grid with a horizontal resolution of 1 km and the BP stress regime. The inversion is conducted separately for each enhancement factor E=0.1, 1, 3 and 6 and the sliding exponents m=1 and m=3.

Within the inverse problem a cost function (J) that measures the misfit between observed, ${\mathit{v}}^{\mathrm{obs}}=(v_x^{\mathrm{obs}},v_y^{\mathrm{obs}})$, and simulated, $\mathit{v}=(v_y,v_y)$, velocities is minimized. We use the observed velocities from the MEaSUREs project (Joughin et al., 2016, 2018) as targets within the inversion. The cost function is composed of two terms which fit the velocities in fast- and slow-moving areas. A third term is a Tikhonov regularization to avoid oscillations due to over fitting. The cost function is defined as follows:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}J(\mathit{v},k)& ={\mathit{\gamma }}_{\mathrm{1}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }(v_x-v_y^{\mathrm{obs}}{)}^{\mathrm{2}}+v_y-v_y^{\mathrm{obs}}{)}^{\mathrm{2}}\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{s}}\\ & +{\mathit{\gamma }}_{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{{\mathrm{\Gamma }}_{\mathrm{s}}}{\int }\ln {\left({\displaystyle \frac{\sqrt{v_x^{\mathrm{2}}+v_y^{\mathrm{2}}}+\mathit{\epsilon }}{\sqrt{{v_x^{\mathrm{obs}}}^{\mathrm{2}}+{v_y^{\mathrm{obs}}}^{\mathrm{2}}}+\mathit{\epsilon }}}\right)}^{\mathrm{2}}\\ \mathrm{d}{\mathrm{\Gamma }}_{\mathrm{s}}& +{\mathit{\gamma }}_{\mathrm{t}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }\mathrm{\nabla }k\cdot \mathrm{\nabla }k\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{b}},\end{array}\end{array}$$

where ε is a minimum velocity used to avoid singularities, and Γ_s and Γ_b are the ice surface and ice base, respectively. An L-curve analysis was performed to pick the Tikhonov parameter γ_t (not shown). We obtained a good agreement with the observed velocities by choosing with ${\mathit{\gamma }}_{\mathrm{1}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}$, ${\mathit{\gamma }}_{\mathrm{2}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{4}}$ and ${\mathit{\gamma }}_{\mathrm{t}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ (see Figs. A1 and A2).

Figure A1 Overview of the ISSM inversion performed with the different values for E=0.1, 1, 3 and 6 at the ice-stream region. (a–d) Simulated surface velocities. (e–h) Velocity difference between observed and simulated velocities. (i–l) Inferred basal friction coefficient k².

Figure A2 Overview of the ISSM inversion performed with the different values for m=1 and m=3 at the outlet region. The upper row shows simulated surface velocities (a), velocity difference between observed and simulated velocities (b), and the inferred basal friction coefficient k² (c) for m=1. The lower row shows the equivalent fields for m=3.