Improvements on the discretisation of boundary conditions to the momentum balance for glacial ice

Abstract

The flow of glacial ice is typically approximated as a nonNewtonian viscous fluid, with the momentum balance described by (an approximation to) the Stokes equations, and the nonlinear rheology described by a flow law. The most commonly used rheological law for glacial ice, Glen's flow law, yields infinite viscosity in the case of zero deformation, which can be the case at the ice surface. This poses a problem when solving the momentum balance numerically. We show that two commonly-used discretisation schemes for the boundary conditions at the ice surface and base, which yield proper numerical convergence when applied to simpler problems, produce poor numerical convergence and large errors, when used to solve the momentum balance with Glen's flow law. We show that a discretisation scheme based on the concept of ghost nodes, which substitutes the boundary conditions directly into the momentum balance equations, yields second-order numerical convergence and errors that can be up to four orders of magnitude smaller. These results are robust across different momentum balance approximations. We show that the improved boundary conditions are particularly useful for solving the 3-D higher-order Blatter-Pattyn Approximation (BPA). In general, this work underlines the importance of thoroughly verifying the numerical solvers used in ice-sheet models, before applying them to future projections of ice-sheet mass loss.

1. Introduction

Sea-level rise caused by large-scale mass loss of the Greenland and Antarctic ice sheets is one of the most worrying potential consequences of unmitigated anthropogenic climate change (Fox-Kemper and others, 2021). Numerical models of ice-sheet flow are the most commonly used tools to project this mass loss into the future, but these projections contain substantial uncertainties (Goelzer and others, 2020; Seroussi and others, 2020; Aschwanden and others, 2021). Part of these uncertainties arises from uncertainties in physical properties of the present-day ice-sheets, such as the temperature and viscosity of the ice (Seroussi and others, 2013; Babaniyi and others, 2021), the subglacial topography (Perego and others, 2014), the subglacial hydrology (Kazmierczak and others, 2022), the sub-shelf ocean conditions and melt rates (Favier and others, 2019; Jourdain and others, 2020; Berends and others, 2023a), and the subglacial friction and sliding law (Sun and others, 2020; Berends and others, 2023b). Another part stems from uncertainty in the future change in climate and mass balance (Goelzer and others, 2020; Seroussi and others, 2020). However, a substantial contribution to the uncertainty stems from the way the physics of flowing ice are described mathematically, and solved numerically. Different approximations to the momentum balance (Pattyn and others, 2008; Bernales and others, 2017; Rückamp and others, 2022), and different numerical treatments of basal friction (Feldmann and others, 2014; Leguy and others, 2021) and sub-shelf melt (Cornford and others, 2020; Leguy and others, 2021; Berends and others, 2023a) at the grounding line all contribute to the spread in modelled projections of ice-sheet mass loss. Whereas uncertainties in the physical properties of the ice can be reduced by improving the quantity and quality of observations, uncertainties in the numerical representation of physical processes must be reduced by the effort of the ice-sheet modelling community.

The majority of ice-sheet models approximate the flow of glacial ice as a nonNewtonian viscous fluid, numerically solving (a simplified approximation to) the Stokes equations to calculate the ice velocity. The (nonlinear) relation between the strain rate and the effective viscosity is described by the flow law defining the rheological properties. All the different approximations to the Stokes equations take the form of an elliptic partial differential equation (PDE), with one or more of the ice velocity components u, v, w as the unknowns to be solved for. As the Stokes equations neglect all momentum terms, there is no time dependence in the equations, so that the velocity of the ice at any point in time depends only on the ice geometry (disregarding the possible dependence of the viscosity on the ice temperature, impurities, damage, or other quantities). This means that boundary conditions only need to be prescribed at the geometrical boundary of the ice sheet. At the ice base, englacial temperatures that are well below the (pressure-corrected) melting point are often assumed to imply that the ice is frozen to the substrate, i.e. u_b = 0, which numerically takes the form of a Dirichlet boundary condition. However, at the ice surface and margin, and at the sliding parts of the ice base, a zero-stress boundary condition applies. Numerically, this boundary condition takes the form of a Neumann boundary condition, where not the value of the unknown itself is specified, but that of its gradient. In the case of sliding at the ice base, possibly the magnitude of the basal velocity is involved as well, leading to a Robin, or mixed boundary condition. While several commonly used ice-sheet models provide the physical equations describing the boundary conditions to their momentum balance approximation (e.g. Pattyn, 2003; Lipscomb and others, 2019), as well as discretisation schemes for the momentum balance itself (e.g. Bueler and Brown, 2009; Berends and others, 2022), few provide details on the discretisation scheme for the boundary conditions. As Neumann boundary conditions can be discretised in different ways, the available literature is often insufficient to reproduce the numerical solvers used in these ice-sheet models.

Here, we investigate three different discretisation schemes for the zero-stress Neumann boundary condition to the momentum balance. Starting with the simplest possible case, we numerically solve the shallow ice approximation (SIA) in the vertical column and compare the results to the analytical solution. We show that not all discretisation schemes for the Neumann boundary condition perform equally well. Particularly, the more obvious choice of using one-sided finite differences to approximate gradients at the domain boundary, results in large errors and poor numerical convergence when Glen's flow law is used to describe the ice rheology. Instead, a ghost-point scheme, which allows the Neumann boundary condition to be substituted directly into the momentum balance equation, produces better results, with numerical convergence of the expected order, and errors that can be up to four orders of magnitude smaller. We demonstrate that the unintuitively poor performance of the one-sided finite difference schemes is caused by a singularity in Glen's flow law, which predicts viscosities that diverge to infinity as the strain rates approach zero. When we instead use a simpler flow law, which predicts finite viscosity everywhere, the one-sided finite difference schemes perform as expected. In a next step, we perform a similar set of experiments with the Blatter-Pattyn approximation (BPA), applied to Experiments A and C of the Ice-Sheet Model Intercomparison Project for Higher-Order Models (ISMIP-HOM; Pattyn and others, 2008). Here we find the same results: the one-sided finite difference schemes perform as expected when applied to the simplified flow law, but fail when applied to Glen's flow law, whereas the ghost-point scheme performs well in both cases.

2. Shallow ice approximation

2.1 Physics

The SIA neglects all viscous stresses except the vertical shear stress (Morland and Johnson, 1980). This approximation holds when the aspect ratio of the ice-sheet geometry is small, i.e. when the characteristic length of horizontal features is either very small or very large with respect to the ice thickness, and when sliding velocities are small compared to flow velocities due to vertical shear. This is the case for large areas of the interior of the Greenland and Antarctic ice sheets.

Choosing an x-coordinate parallel to the ice flow, the SIA reads:

(1)$$\displaystyle{\partial \over {\partial z}}\left({\eta \displaystyle{{\partial u} \over {\partial z}}} \right) = \rho g\displaystyle{{\partial h} \over {\partial x}}.$$

The symbols used throughout this work are defined in Table 1. Glen's flow law (Paterson, 1994) relates the effective viscosity η to the effective strain rate $\dot{\varepsilon }$:

(2)$$\eta = \displaystyle{1 \over 2}A( {T^\ast } ) ^{{-}1/n}\dot{\varepsilon }^{( 1-n) /n}, \;$$

(3)$$\dot{\varepsilon } = \left[{\displaystyle{1 \over 4}{\left({\displaystyle{{\partial u} \over {\partial z}}} \right)}^2} \right]^{1/2}.$$

Table 1. Symbols, notation, and units used in this work

Here, A(T*) is a temperature-dependent flow rate factor. When choosing a no-slip boundary condition at the ice base (i.e. u(z = b) = 0), and a zero-stress boundary condition at the ice surface (i.e. (∂u/∂z)(z = h) = 0), this nonlinear partial differential equation has the following analytical solution:

(4)$$u( z ) = {-}2( {\rho g} ) ^n\left\vert {\displaystyle{{\partial h} \over {\partial x}}} \right\vert ^{n-1}\displaystyle{{\partial h} \over {\partial x}}\mathop \smallint \limits_b^z A( {T^\ast ( {{z}^{\prime}} ) } ) ( {h-{z}^{\prime}} ) ^nd{z}^{\prime}.$$

For the case of isothermal ice, where A(T*(z)) = A, this simplifies to the widely used:

(5)$$u( z ) = {-}2( {\rho g} ) ^n\left\vert {\displaystyle{{\partial h} \over {\partial x}}} \right\vert ^{n-1}\displaystyle{{\partial h} \over {\partial x}}\displaystyle{A \over {n + 1}}( {H^{n + 1}-{( {h-z} ) }^{n + 1}} ).$$

2.2 Numerical solution

This analytical solution can be approximated numerically. As is common practise in many ice-sheet models (e.g. PISM, Bueler and Brown, 2009; CISM, Lipscomb and others, 2019; IMAU-ICE, Berends and others, 2022), we use a staggered-grid approach to solve the momentum balance, where material properties such as the englacial temperature, flow rate factor, strain rate, and effective viscosity, are staggered with respect to fluxes such as the ice velocities. This is illustrated for the SIA, which only concerns the vertical dimension z, in Figure 1.

Figure 1. The vertically staggered grid used to solve the SIA. The ice velocity u is defined on the regular grid (solid circles), while the effective viscosity η is defined on the staggered grid (open circles). In some literature, the staggered grid would be indicated by ‘half-indexing’, e.g. η_k−(1/2), η_k+(1/2).

The regular grid, which we shall from here on call the a-grid (Arakawa and Lamb, 1977), has n _z nodes spaced at regular intervals of Δz, with the first (k = 1) and last (k = n _z) nodes coinciding respectively with the ice base and the ice surface, so that:

(6)$$z_a^k = b + H\displaystyle{{k-1} \over {n_z-1}}, \;$$

(7)$$\Delta z = \displaystyle{H \over {n_z-1}}.$$

Many ice-sheet models use an irregular vertical grid, with thinner, more closely-spaced layers near the ice base to more accurately capture the higher strain rates there (e.g. PISM, Martin and others, 2011; Yelmo, Robinson and others, 2019; IMAU-ICE, Berends and others, 2022). We include some additional experiments in Appendix A where we investigate the effect of such an irregular grid.

The staggered grid, which we shall from here on call the c-grid, has n _z − 1 nodes, which lie halfway between the a-grid nodes:

(8)$$z_c^k = \displaystyle{{z_a^k + z_a^{k + 1} } \over 2}.$$

Note that we use the subscript to indicate the grid, and the superscript to indicate the node index. We approximate the vertical profile of the horizontal ice velocity u(z) by discretising it on the a-grid, whereas the effective viscosity η is discretised on the c-grid:

(9)$$u_a^k \approx u( {z_a^k } ) , \;$$

(10)$$\eta _c^k \approx \eta ( {z_c^k } ).$$

Gathering all the n _z values of $u_a^k$ together yields the discretised velocity vector u _a:

(11)$${\boldsymbol u}_a = \big [ {u_a^k } \big ] = \big [ {u_a^1 , \;u_a^2 , \;\ldots , \;u_a^{k-1} , \;u_a^k , \;u_a^{k + 1} , \;\ldots , \;u_a^{n_z-1} , \;u_a^{n_z} } \big ] ^T.$$

Similarly, the discretised effective viscosity vector is expressed as:

(12)$$ \eqalign{ {\boldsymbol \eta }_c = [ {\eta_c^k } ] = [ {\eta_c^1 , \;\eta_c^2 , \;\ldots , \;\eta_c^{k-1} , \;\eta_c^k , \;\eta_c^{k + 1} , \;\ldots , \;\eta_c^{n_z-2} , \;\eta_c^{n_z-1} } ] ^T}.$$

To approximate the gradient operators ∂/∂z between the two grids, we use a two-point central differencing scheme:

(13)$$\displaystyle{{\partial f} \over {\partial z}}_{a\to c}^k = \displaystyle{{\,f_a^{k + 1} -f_a^k } \over {\Delta z}}, \;$$

(14)$$\displaystyle{{\partial f} \over {\partial z}}_{c\to a}^k = \left\{{\matrix{ {\displaystyle{{\,f_c^k -f_c^{k-1} } \over {\Delta z}}, \;\;\;1 < k < n_z} \cr {0, \;\;\qquad\enspace\quad {\rm otherwise}.} \cr } } \right..$$

Note that the gradient of ∂f ^k/∂z _c→a is not defined on the first and last nodes.

These operators can be represented by matrices, which can be multiplied with the discretised vectors u_a and η_c to calculate their gradients:

(15)$$M_{\partial /\partial z}^{a\to c} = \left[{\matrix{ {\displaystyle{{-1} \over {\Delta z}}} & {\displaystyle{1 \over {\Delta z}}} & 0 & \cdots & 0 & 0 & 0 \cr\cr 0 & {\displaystyle{{-1} \over {\Delta z}}} & {\displaystyle{1 \over {\Delta z}}} & \cdots & 0 & 0 & 0 \cr\cr \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \cr\cr 0 & 0 & 0 & \cdots & {\displaystyle{{-1} \over {\Delta z}}} & {\displaystyle{1 \over {\Delta z}}} & 0 \cr\cr 0 & 0 & 0 & \cdots & 0 & {\displaystyle{{-1} \over {\Delta z}}} & {\displaystyle{1 \over {\Delta z}}} \cr\cr } } \right], \;$$

(16)$$M_{\partial /\partial z}^{c\to a} = \left[{\matrix{ 0 & 0 & 0 & \cdots & 0 & 0 & 0 \cr\cr {\displaystyle{{-1} \over {\Delta z}}} & {\displaystyle{1 \over {\Delta z}}} & 0 & \cdots & 0 & 0 & 0 \cr\cr 0 & {\displaystyle{{-1} \over {\Delta z}}} & {\displaystyle{1 \over {\Delta z}}} & \cdots & 0 & 0 & 0 \cr\cr \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \cr\cr 0 & 0 & 0 & \cdots & {\displaystyle{{-1} \over {\Delta z}}} & {\displaystyle{1 \over {\Delta z}}} & 0 \cr\cr 0 & 0 & 0 & \cdots & 0 & {\displaystyle{{-1} \over {\Delta z}}} & {\displaystyle{1 \over {\Delta z}}} \cr\cr 0 & 0 & 0 & \cdots & 0 & 0 & 0 \cr } } \right].$$

Note that the size of $M_{\partial /\partial z}^{a\to c}$ is (n_z-1)-by-n _z, while $M_{\partial /\partial z}^{c\to a}$ is n _z-by-(n_z-1). The discretised approximation to the SIA can now be expressed as a matrix equation:

(17)$$M_{\partial /\partial z}^{c\to a} D( {{\boldsymbol \eta }_c} ) M_{\partial /\partial z}^{a\to c} {\boldsymbol u}_a = A{\boldsymbol u}_a = {\boldsymbol b}\,. $$

Here, D(η_c) is an (n_z-1)-by-(n_z-1) diagonal matrix with the elements of η_c on the diagonal. The load vector b is an n _z-by-1 vector with every element having a value of ρg(∂h/∂x). The stiffness matrix A is equal to the product $M_{\partial /\partial z}^{a\to c} D( {{\boldsymbol \eta }_c} ) M_{\partial /\partial z}^{c\to a}$. We can write out the coefficients of the M-matrices to derive a single linear equation from this system:

(18a)$$M_{\partial /\partial z}^{c\to a} {\boldsymbol u}_a = \displaystyle{1 \over {\Delta z}}( {u_a^{i + 1} -u_a^i } ) , \;$$

(18b)$$D( {{\boldsymbol \eta }_c} ) M_{\partial /\partial z}^{a\to c} {\boldsymbol u}_a = \displaystyle{{\eta _c^i } \over {\Delta z}}( {u_a^{i + 1} -u_a^i } ) , \;$$

(18c)$$M_{\partial /\partial z}^{c\to a} D( {{\boldsymbol \eta }_c} ) M_{\partial /\partial z}^{a\to c} {\boldsymbol u}_a = \displaystyle{1 \over {{\rm \Delta }z^2}}[ {\eta_c^i ( {u_a^{i + 1} -u_a^i } ) -\eta_c^{i-1} ( {u_a^i -u_a^{i-1} } ) } ].$$

The nonlinear dependence of η on the strain rate ∂u/∂z is solved for by Picard iteration, iteratively solving for u based on the current estimate of η, and then recalculating η for the new solution of u. The iteration is stopped when the L²-norm of the difference between two successive solutions of u is less then 10⁻⁸ of the L²-norm of u itself.

We adopt the common practise of adding a small regularisation term $\varepsilon _0^2 = 10^{{-}20}$ to the effective strain rate in Glen's flow law:

(19)$$\dot{\varepsilon } = \left[{\displaystyle{1 \over 4}{\left({\displaystyle{{\partial u} \over {\partial z}}} \right)}^2 + {\rm \;}\varepsilon_0^2 } \right]^{{1 \over 2}}.$$

This prevents divide-by-zero errors when (∂u/∂z) = 0. In the two-point one-sided scheme, this is always the case at the last staggered node below the surface, $\partial u^{n_z-1}/\partial z_c$. Otherwise, this only happens in the very first viscosity iteration, as we choose an initial guess of u _a = 0; after that, (∂u/∂z _c) > 0 everywhere. Alternatively, enforcing $\dot{\varepsilon } \ge \varepsilon _0$ (thus changing only the values in the first iteration, or at the last staggered node) does not affect the results.

2.3 Boundary conditions

The first and last row of the stiffness matrix A and the load vector b must describe the boundary conditions at the base and surface of the ice, respectively. The no-slip condition at the base, i.e. u(z = b) = 0, is a Dirichlet boundary condition, which is represented simply by setting the diagonal element of the first row of A to unity, all other elements in the row to zero, and the first element of b to zero. The zero-stress boundary condition at the surface, i.e. (∂u/∂z)(z = h) = 0, is a Neumann boundary condition that is less trivial to implement.

The first option we explore, approximates the vertical shear strain rate at the ice surface, $\partial u^{n_z}/\partial z_a$, with a two-point one-sided differencing scheme:

(20)$${\rm \;}\displaystyle{{\partial u} \over {\partial z}}_a^{n_z} = \displaystyle{{u_a^{n_z-1} -u_a^{n_z} } \over {\Delta z}} = 0.$$

Since in this case, the ${\cal O}( {\Delta z^2} )$ terms in the Taylor expansion of ∂u/∂z around $z_a^{n_z}$ do not cancel out, this scheme is only first-order accurate (i.e. the truncation error is of order ${\cal O}( {\Delta z} )$).

We also explore three-point, four-point, and five-point one-sided differencing schemes to approximate the gradient:

(21)$$\displaystyle{{\partial u} \over {\partial z}}_a^{n_z} = \displaystyle{{-1} \over {2\Delta z}}u_a^{n_z-2} + \displaystyle{4 \over {2\Delta z}}u_a^{n_z-1} -\displaystyle{3 \over {2\Delta z}}u_a^{n_z} = 0, \;$$

(22)$$\eqalign{\displaystyle{{\partial u} \over {\partial z}}_a^{n_z} = \displaystyle{2 \over {6\Delta z}}u_a^{n_z-3} -\displaystyle{9 \over {6\Delta z}}u_a^{n_z-2} + \displaystyle{{18} \over {6\Delta z}}u_a^{n_z-1} -\displaystyle{{11} \over {6{\rm \Delta }z}}u_a^{n_z} = 0, \;}$$

(23)$$\eqalign{\displaystyle{{\partial u} \over {\partial z}}_a^{n_z} & = \,\displaystyle{{-3} \over {12\Delta z}}u_a^{n_z-4} + \displaystyle{{16} \over {12\Delta z}}u_a^{n_z-3} -\displaystyle{{36} \over {12\Delta z}}u_a^{n_z-2} \cr & \quad + \displaystyle{{48} \over {12\Delta z}}u_a^{n_z-1} -\displaystyle{{25} \over {12{\rm \Delta }z}}u_a^{n_z} = 0.}$$

It can be shown that these scheme are, respectively, second-order, third-order, and fourth-order accurate.

The last option we explore uses the concept of a ‘ghost node’. An in-depth explanation of the concept is provided by Fornberg (2006). We temporarily construct an additional node (the ‘ghost node’) on the a-grid outside the domain boundary, such that:

(24)$$z_a^{n_z + 1} = z_a^{n_z} + \Delta z = h + \Delta z.$$

We then formulate the discretisation of the SIA at the ice surface, without considering the boundary conditions just yet. Expanding Eqn (1) using the product rule yields:

(25)$$\displaystyle{{\partial \eta } \over {\partial z}}\displaystyle{{\partial u} \over {\partial z}} + \eta \displaystyle{{\partial ^2u} \over {\partial z^2}} = \rho g\displaystyle{{\partial h} \over {\partial x}}.$$

We discretise the first and second derivatives of u at the ice surface by using standard three-point two-sided differencing schemes:

(26)$$\displaystyle{{\partial u} \over {\partial z}}_a^{n_z} = \displaystyle{{u_a^{n_z + 1} -u_a^{n_z-1} } \over {2\Delta z}}, \;$$

(27)$$\displaystyle{{\partial ^2u} \over {\partial z^2}}_a^{n_z} = \displaystyle{{u_a^{n_z + 1} + u_a^{n_z-1} -2u_a^{n_z} } \over {\Delta z^2}}.$$

Substituting the zero-stress boundary condition, $( \partial u^{n_z}/\partial z_a) = 0$, into Eqn (26) implies that:

(28)$$u_a^{n_z + 1} = u_a^{n_z-1}.$$

Substituting this into Eqn (27) yields:

(29)$$\displaystyle{{\partial ^2u} \over {\partial z^2}}_a^{n_z} = \displaystyle{{2u_a^{n_z-1} -2u_a^{n_z} } \over {\Delta z^2}}.$$

Finally, substituting this into Eqn (25), the discretised expression for the SIA now reads:

(30)$$\eta _c^{n_z-1} \displaystyle{{2u_a^{n_z-1} -2u_a^{n_z} } \over {\Delta z^2}} = \rho g\displaystyle{{\partial h} \over {\partial x}}.$$

Note that $u_a^{n_z + 1}$ has disappeared from the expression. The ghost node is merely a tool to derive the expression, and does not appear in the final numerical solution.

Equation (30) can alternatively be derived by substituting the boundary condition from Eqn (28) directly into the discretised form of the SIA. Begin with the general form in Eqn (18c), and let i = n _z:

(31)$$\displaystyle{1 \over {{\rm \Delta }z^2}}[ {\eta_c^{n_z} ( {u_a^{n_z + 1} -u_a^{n_z} } ) -\eta_c^{n_z-1} ( {u_a^{n_z} -u_a^{n_z-1} } ) } ] = \rho g\displaystyle{{\partial h} \over {\partial x}}.$$

From Eqn (28), it also follows that $\eta _c^{n_z} = \eta _c^{n_z-1}$:

(32)$$\displaystyle{1 \over {{\rm \Delta }z^2}}[ {\eta_c^{n_z-1} ( {u_a^{n_z-1} -u_a^{n_z} } ) -\eta_c^{n_z-1} ( {u_a^{n_z} -u_a^{n_z-1} } ) } ] = \rho g\displaystyle{{\partial h} \over {\partial x}}.$$

This can be rearranged to read:

(33)$$\eta _c^{n_z-1} \displaystyle{{2u_a^{n_z-1} -2u_a^{n_z} } \over {\Delta z^2}} = \rho g\displaystyle{{\partial h} \over {\partial x}}.$$

Observe that this expression is identical to Eqn (30).

Lastly, it is noteworthy that we could alternatively keep both the general PDE in Eqn (31) and the boundary condition in Eqn (28) as separate linear equations, and include the ghost node as an additional degree of freedom. This means the matrix equation that must be solved has one additional row and column, so that we no longer need to manually eliminate the ghost node as a degree of freedom. It is trivially easy to write a program that solves this ‘extended’ matrix equation, and to find that this yields an identical solution. This illustrates the fundamental difference between the ghost-node scheme, and the different one-sided schemes. In the one-sided schemes, the linear equation representing the boundary condition replaces the linear equation representing the PDE, while in the ghost-node scheme, both linear equations are kept, and an additional degree of freedom (the ghost node) is introduced to allow both of them to be solved (although it is possible, as we have done here, to manually solve for the extra degree of freedom beforehand).

2.4 Experiments

We compare the numerical solution to the SIA with Glen's flow law to the analytical solution in the test case of an isothermal slab of ice with a thickness of H = 2000 m, lying on an inclined plane, sloping downward in the x-direction with a slope of (∂h/∂x) = −10⁻². We perform experiments for three different flow laws: Glen's flow law (Section 2.4.1), a linear, nondiverging flow law (Section 2.4.2), and an over-regularised variant of Glen's flow law (Section 2.4.3). We analyse the results of these experiments in Section 2.4.4.

2.4.1 Glen's flow law

For Glen's flow law, we use a uniform flow factor of A = 10⁻¹⁶ Pa⁻³ yr⁻¹. The analytical solution to the SIA for these parameters is shown in Figure 2.

Figure 2. The analytical solution to the SIA with Glen's flow law.

We use the discretisation schemes described in Section 2.2 to solve the SIA numerically. By doing so at increasing numbers of nodes, we can investigate how quickly the error with respect to the analytical solution decreases, and determine the rate of convergence. We do this for all three different options for discretising the Neumann boundary condition at the ice surface, described in Section 2.3. In order to analyse the convergence behaviour of the different discretisation schemes, we define the relative error in the velocity solution at the ice surface:

(34)$${\rm err\;}u = \left\vert {\displaystyle{{u^{k = 1}-u_{{\rm analytical}}( {z = h} ) } \over {u_{{\rm analytical}}( {z = h} ) }}} \right\vert.$$

Shown in Figure 3 are the velocity solutions resulting from the two-point one-sided scheme, for different numbers of grid points. As can be seen, the significant errors in the solution are not localised at the ice surface, but instead affect the solution throughout the ice column. This implies that calculating the error at a different point in the vertical column, or even as the L²-norm over the entire column, yields qualitatively the same results.

Figure 3. Velocity solutions resulting from the two-point one-sided scheme for the combination of the SIA with Glen's flow law, for n_z = 8 (red), n_z = 16 (blue), n_z = 32 (green), and n_z = 64 (yellow), compared to the analytical solution (solid black line).

2.4.2 Linear flow law

We perform a set of experiments where we solve the SIA with a different flow law. In this case, we define the effective viscosity as a simple function of the vertical coordinate z:

(35)$$\eta = a + z.$$

This expression yields a finite viscosity everywhere, with small values at the ice base and larger values at the ice surface. Note that the units in this expression are not consistent. It is not meant to represent a realistic flow law, but merely serves to create a mathematical problem that is qualitatively similar to the SIA, but without the diverging term in the differential equation.

The analytical solution to the SIA, combined with this linear, nondiverging flow law, and with the respective no-slip and zero-stress boundary conditions at the base and surface of the ice, reads:

(36)$$\eqalign{u( z ) = {-}\rho g\displaystyle{{\partial h} \over {\partial x}}( {a + H} ) \log ( {a + z} ) + \rho g\displaystyle{{\partial h} \over {\partial x}}( {z + ( {a + H} ) \log ( a ) } ) .}$$

The values for the different physical parameters of the experiment are provided in Table 2.

Table 2. Physical parameters used in the SIA experiment with the artificial flow law

The analytical solution to the SIA with this flow law, for these parameters, is shown in Figure 4.

Figure 4. The analytical solution to the SIA with the artificial flow law.

2.4.3 Over-regularised Glen's flow law

Here, we perform a set of experiments using Glen's flow law, but with a very large value of the regularisation term $\varepsilon _0^2 = 10^{{-}1}$ (see Eqn (19)). As with the linear flow law, this results in a finite viscosity everywhere, but additionally includes the nonlinearity of Glen's flow law. In this case, the analytical solution for the SIA with Glen's flow law is invalid. Instead we compare to a numerical solution with a very high resolution (found by using the ghost-node scheme with n _z = 2¹⁴ = 16, 384 grid points, which is shown in Fig. 5.

Figure 5. The high-resolution solution to the SIA with the over-regularised Glen's flow law.

2.4.4 Results

The relative errors in the ice velocity at the surface are shown as a function of the number of nodes, for the five respective discretisation schemes (the two-point, three-point, four-point, and five-point one-sided schemes, and the ghost-node scheme), in Figure 6. We have used values of the number of grid points n _z ranging from 16 to 1024. Typically, large-scale ice-sheet models applied to e.g. the Antarctic ice-sheet use numbers in the range of 10¹–10² (e.g. 81 layers for SICOPOLIS and 121 layers for PISM in their respective ISMIP6 Antarctica set-ups; Seroussi et al., in review).

Figure 6. Convergence of the relative error in the surface velocity with the number of nodes, for the combination of the SIA with (A) Glen's flow law, (B) the linear flow law, and (C) the over-regularised variant of Glen's flow law. The different graphs show the results for the five different ways of discretising the zero-stress boundary condition at the ice surface: the two-point (blue), the three-point (red), the four-point (green), the five-point (yellow) one-sided schemes, and the ghost-node scheme (purple). Log-linear curves (solid lines) are fitted to each set of experiments to calculate the order of convergence R, which is displayed in the legend.

For the case of Glen's flow law (Fig. 6a), the two-point, three-point, and four-point one-sided schemes display first-order convergence of the error with number of nodes. The five-point one-sided scheme and the ghost-node scheme display second-order convergence. The first-order convergence of the two-point one-sided scheme, and the second-order convergence of the ghost-point scheme, are to be expected from the order of accuracy of these respective schemes. However, the three- one-sided scheme is second-order accurate. The fact that this scheme still yields first-order convergence is therefore unexpected. The four-point and five-point one-sided schemes are, respectively, third-order and fourth-order accurate. However, as the discretisation of the SIA in the ice column (Eqn (18c)) is only second-order accurate, the total error will be dominated by the second-order term. For the five-point scheme, this is reflected in the results (which show second-order convergence), whereas the four-point scheme, like the three-point scheme, displays only first-order convergence.

For the case of the linear flow law (Fig. 6b), the two-point one-sided scheme still displays first-order convergence, while the five-point scheme and the ghost-node scheme still display second-order convergence. However, the three-point and four-point one-sided schemes now display second-order convergence.

For the case of the over-regularised variant of Glen's flow law (Fig. 6c), the results are qualitatively the same as for the linear flow law, with all schemes except the two-point scheme displaying second-order convergence. The only difference between this experiment and Glen's flow law, is that here there is no more singularity in the effective viscosity at the ice surface. This illustrates that it is likely the presence of this singularity that poses a problem to the different numerical solvers. The experiments with Glen's flow law also included a regularisation term, albeit with a much smaller value of $\varepsilon _0^2 = 10^{{-}20}$, which should prevent the effective viscosity from becoming infinite. In Appendix A, we present a number of simulations where we vary the value of $\varepsilon _0^2$ and study the resulting velocity errors and convergence rates of the different discretisation schemes. In summary, values of $\varepsilon _0^2$ that are small enough not to cause any significant viscous deformation in the upper part of the ice column, still yield high enough effective viscosities to cause issues with the different one-sided discretisation schemes. The ghost-node scheme, on the other hand, robustly produces small errors and second-order convergence for all values of $\varepsilon _0^2 = 10^{{-}10}$; for larger values, (unphysical) viscous deformation in the upper part of the ice column becomes noticeable.

3. Blatter-Pattyn approximation

3.1 Physics and numerical solution

The Blatter-Pattyn approximation (BPA; Pattyn, 2003) is a so-called higher-order approximation to the Stokes equations. It includes all viscous stresses except for those involving the vertical component of the ice velocity, which is still assumed to be negligibly small. Making no assumptions about the direction of the ice flow, the BPA reads:

(37a)$$\eqalign{\displaystyle{\partial \over {\partial x}}\left[{2\eta \left({2\displaystyle{{\partial u} \over {\partial x}} + \displaystyle{{\partial v} \over {\partial y}}} \right)} \right] + \displaystyle{\partial \over {\partial y}}\left[{\eta \left({\displaystyle{{\partial u} \over {\partial y}} + \displaystyle{{\partial v} \over {\partial x}}} \right)} \right] + \displaystyle{\partial \over {\partial z}}\left[{\eta \displaystyle{{\partial u} \over {\partial z}}} \right] = \rho g\displaystyle{{\partial h} \over {\partial x}},}$$

(37b)$$\eqalign{ \displaystyle{\partial \over {\partial y}}\left[{2\eta \left({2\displaystyle{{\partial v} \over {\partial y}} + \displaystyle{{\partial u} \over {\partial x}}} \right)} \right] + \displaystyle{\partial \over {\partial x}}\left[{\eta \left({\displaystyle{{\partial v} \over {\partial x}} + \displaystyle{{\partial u} \over {\partial y}}} \right)} \right] + \displaystyle{\partial \over {\partial z}}\left[{\eta \displaystyle{{\partial v} \over {\partial z}}} \right] = \rho g\displaystyle{{\partial h} \over {\partial y}}.}$$

Relative to the SIA, several additional terms now appear in the effective strain rate $\dot{\varepsilon }$, which were previously neglected:

(38)$$ \eqalign{ \dot{\varepsilon } = & \bigg[{{\left({\displaystyle{{\partial u} \over {\partial x}}} \right)}^2 + {\left({\displaystyle{{\partial v} \over {\partial y}}} \right)}^2} + \displaystyle{{\partial u} \over {\partial x}}\displaystyle{{\partial v} \over {\partial y}} + \displaystyle{1 \over 4}{\left({\displaystyle{{\partial u} \over {\partial y}} + \displaystyle{{\partial v} \over {\partial x}}} \right)}^2 \cr &\quad+ \displaystyle{1 \over 4}{\left({\displaystyle{{\partial u} \over {\partial z}}} \right)}^2 + \displaystyle{1 \over 4}{\left({\displaystyle{{\partial v} \over {\partial z}}} \right)}^2\bigg]^{1/2}.}$$

The zero-stress boundary condition at the ice surface now reads:

(39a)$$2\displaystyle{{\partial h} \over {\partial x}}\left[{2\displaystyle{{\partial u} \over {\partial x}} + \displaystyle{{\partial v} \over {\partial y}}} \right] + \displaystyle{{\partial h} \over {\partial y}}\left[{\displaystyle{{\partial u} \over {\partial y}} + \displaystyle{{\partial v} \over {\partial x}}} \right]-\displaystyle{{\partial u} \over {\partial z}} = 0, \;$$

(39b)$$2\displaystyle{{\partial h} \over {\partial y}}\left[{2\displaystyle{{\partial v} \over {\partial y}} + \displaystyle{{\partial u} \over {\partial x}}} \right] + \displaystyle{{\partial h} \over {\partial x}}\left[{\displaystyle{{\partial v} \over {\partial x}} + \displaystyle{{\partial u} \over {\partial y}}} \right]-\displaystyle{{\partial v} \over {\partial z}} = 0.$$

The boundary condition for sliding ice at the base reads:

(40a)$$2\displaystyle{{\partial b} \over {\partial x}}\left[{2\displaystyle{{\partial u} \over {\partial x}} + \displaystyle{{\partial v} \over {\partial y}}} \right] + \displaystyle{{\partial b} \over {\partial y}}\left[{\displaystyle{{\partial u} \over {\partial y}} + \displaystyle{{\partial v} \over {\partial x}}} \right]-\displaystyle{{\partial u} \over {\partial z}} + \displaystyle{\beta \over \eta } = 0, \;$$

(40b)$$2\displaystyle{{\partial b} \over {\partial y}}\left[{2\displaystyle{{\partial v} \over {\partial y}} + \displaystyle{{\partial u} \over {\partial x}}} \right] + \displaystyle{{\partial b} \over {\partial x}}\left[{\displaystyle{{\partial v} \over {\partial x}} + \displaystyle{{\partial u} \over {\partial y}}} \right]-\displaystyle{{\partial v} \over {\partial z}} + \displaystyle{\beta \over \eta } = 0.$$

For the case of ice that is frozen to the bedrock, preventing any basal sliding, the boundary condition is simply u = v = 0. This set of coupled, nonlinear partial differential equations has no known analytical solutions. The discretisation of the BPA and its boundary conditions follows the same general approach as that of the SIA, and is described in detail in Appendix A.

3.4 Experiments

3.4.1 Glen's flow law

We perform ISMIP-HOM experiment A and C (Pattyn and others, 2008) to investigate the different discretisations. These experiments concern an infinite slab of isothermal ice lying on an inclined plane. In experiment A, the ice is frozen to the base, and undulations in both horizontal directions are superimposed on the sloping bed, leading to nonnegligible horizontal stretch/shear strain rates. In experiment C, the bed is a simple flat, sloping plane, but basal sliding is allowed, with periodic spatial variations in the basal friction coefficient. Pattyn and others (2008) showed that the BPA yields results that are nearly identical to those of the Stokes equations in these settings, while later work (e.g. Berends and others, 2022) showed that the SIA, the hybrid SIA/shallow shelf approximation (Bueler and Brown, 2009), and the depth-integrated viscosity approximation (Goldberg, 2011) deviate significantly from the Stokes solution.

The geometry of experiment A is given by the following equations:

(41)$$h( {x, \;y} ) = 2000-x\tan \theta , \;$$

(42)$$H( {x, \;y} ) = 1000 + 500\sin \left({\displaystyle{{2\pi x} \over L}} \right)\sin \left({\displaystyle{{2\pi y} \over L}} \right), \;$$

(43)$$b( {x, \;y} ) = h( {x, \;y} ) -H( {x, \;y} ).$$

At the lateral domain boundaries, periodic boundary conditions apply, and a no-slip boundary condition is prescribed at the ice base. The parameters of the experiment are listed in Table 3.

Table 3. Parameters of the ISMIP-HOM A experiment with Glen's flow law

For experiment C, Eqn (42) is replaced by a uniform ice thickness of 1000 m, so that Eqn (43) yields a flat, sloping bed. The surface/bedrock slope θ is assigned a lower value of 0.1 degrees. The following expression describes the basal friction coefficient β:

(44)$$\beta ( {x, \;y} ) = 1000 + 1000\sin \left({\displaystyle{{2\pi x} \over L}} \right)\sin \left({\displaystyle{{2\pi x} \over L}} \right).$$

No analytical solutions exist for these experiments. Instead, model verification is done by comparing against the ensemble results by Pattyn and others (2008). We perform these experiments with three different discretisation schemes for the zero-stress boundary at the ice surface: the two-point scheme, the three-point scheme, and the ghost-node scheme. Implementations of the four-point and five-point schemes consistently failed to converge during the iteration over the nonlinear effective viscosity. We use a square grid of 40 by 40 nodes in the horizontal plane (identical to the models in the ensemble from Pattyn and others, 2008), and n _z ∈ [8, 32] nodes in the vertical column. Our discretisation of the BPA, and of the three different schemes to discretise the surface and basal boundary conditions, are derived in Appendix A. The results of our simulations are compared to the Pattyn and others (2008) ensemble in Figures 7 and 8.

Figure 7. Results of ISMIP-HOM Experiment A with the BPA, using Glen's flow law, with a horizontal length scale of L = 160 km. The ensemble results of Pattyn and others (2008) for the participating higher-order models (green) and full-Stokes models (blue) are shown by the shaded areas, with the solid line indicating the ensemble mean. The coloured lines show the results of our model at different numbers of vertical layers (see legend in Panel B), for the two-point one-sided scheme (panel A), the three-point one-sided scheme (panel B), and the ghost-point scheme (panel C). The small subpanels zoom in on the thin, slow-moving ice around x = L/4.

Figure 8. Results of ISMIP-HOM Experiment C with the BPA, using Glen's flow law, with a horizontal length scale of L = 160 km. The ensemble results of Pattyn and others (2008) for the participating higher-order models (green) and full-Stokes models (blue) are shown by the shaded areas, with the solid line indicating the ensemble mean. The coloured lines show the results of our model at different numbers of vertical layers (see legend in Panel B), for the two-point one-sided scheme (panel A), the three-point one-sided scheme (panel B), and the ghost-point scheme (panel C). The small subpanels zoom in on the thin, slow-moving ice around x = L/4.

In experiment C (Fig. 8), there is no visible difference between the three-point one-sided scheme and the ghost-node scheme, with both giving very accurate solutions even at coarse vertical resolutions. We suspect this is because the ice flow in experiment C is dominated by sliding rather than by vertical shearing, so that the horizontal velocities are nearly uniform in the vertical, implying that the errors in the velocity solution are dominated by the horizontal shearing terms.

We have additionally performed experiment A with a length scale of L = 5 km. The results of these simulations are shown in Figure 9.

Figure 9. Results of ISMIP-HOM Experiment A with the BPA, using Glen's flow law, with a horizontal length scale of L = 5 km. The ensemble results of Pattyn and others (2008) for the participating higher-order models (green) and full-Stokes models (blue) are shown by the shaded areas, with the solid line indicating the ensemble mean. The coloured lines show the results of our model at different numbers of vertical layers (see legend in Panel B), for the two-point one-sided scheme (panel A), the three-point one-sided scheme (panel B), and the ghost-point scheme (panel C). The small subpanels zoom in on the thin, slow-moving ice around x = L/4.

Since no analytical solutions exist, we instead compare to high-resolution numerical solutions to perform a convergence analysis. For the high-resolution solution, we use the same horizontal resolution (as we are only interested in the convergence with the vertical resolution), but a vertical grid of n _z = 128 layers. These are calculated with all three schemes separately. Based on the convergence analysis, the errors in the high-resolution solution should be at least an order of magnitude smaller than those in any of the other numerical solutions, small enough to not significantly affect the error estimates or the convergence analysis. The convergence analyses for all three experiments are shown in Figure 10.

Figure 10. Convergence of the relative error in the surface velocity with the number of nodes, for the combination of the BPA with Glen's flow law, for (A) ISMIP-HOM experiment A with a horizontal length scale of L = 160 km, (B) experiment C with a length scale of L = 160 km, and (C) experiment A with a length scale of L = 5 km. The different graphs show the results for the three different ways of discretising the zero-stress boundary condition at the ice surface: the two-point one-sided scheme (blue), the three-point one-sided scheme (red), and the ghost-node scheme (green). Log-linear curves (dashed lines) are fitted to each set of experiments to calculate the order of convergence R, which is displayed in the legend.

The results of these experiments are less straightforward to interpret than those of the 1-D SIA experiments. In the 1-D case, the strain rate is zero at the ice surface, and the effective viscosity diverges to infinity. In the 3-D case of the ISMIP-HOM experiments, the lateral variations in the geometry result in additional horizontal strain rates, so that the effective strain rate is never zero (although it can still become very small), and the effective viscosity is never infinite (although it can still become very large). This effect is more pronounced in experiment C, where the vertical shear strain rates are actually smaller than the horizontal stretch/shear strain rates. Furthermore, in experiment C there is an additional Neumann boundary condition at the ice base, which is discretised using the same scheme as that at the surface. As the effective viscosity at the ice base is much smaller, the convergence behaviour of the three-point scheme is likely different too. Additionally, due to the coordinate transformation, the vertical resolution also has a small effect on the accuracy of the discretisation of the horizontal strain rates, further complicating the analysis. And lastly, because of the bedrock undulations in experiment A, the constant number of vertical layers implies a different vertical resolution in different parts of the domain, which is further complicated by the fact that the horizontal and vertical strain rates are also laterally varying.

The difference in performance between the three-point, one-sided scheme and the ghost-node scheme is less clear than in the 1-D SIA experiments. For both flow laws, the ghost-node scheme produces smaller errors than the three-point scheme, with a more pronounced difference in experiment A (where the effective strain rates at the surface can still be very small) than in experiment C (where the horizontal strain rates are more pronounced).

The results for experiment A with a horizontal length scale of L = 5 km are qualitatively similar to those with L = 160 km. The convergence rate of the ghost-node scheme is still slightly higher than that of the three-point scheme, and the errors are approximately an order of magnitude smaller. Note that, because of the smaller horizontal length scale, the horizontal stretch and shear strain rates are much larger than before, so that the effective viscosity is much smaller. This likely explains why there is now less difference between the three-point scheme and the ghost-node scheme.

3.4.2 Linear flow law

As with the SIA, we performed an additional set of experiments where we replace Glen's flow law with an expression that does not depend on the strain rates. In order to maintain a uniform, high value of the effective viscosity at the ice surface, we here define the following expression:

$$\eta ( {x, \;y, \;z} ) = a + z + 1500-h( {x, \;y} ).$$

We have only performed simulations with the artificial flow law for experiment A with a length scale of L = 160 km. The parameters for this new experiment are listed in Table 4.

Table 4. Parameters for ISMIP-HOM experiment A with the linear flow law

The results of these experiments are shown in Figure 11.

Figure 11. Results of ISMIP-HOM Experiment A with the BPA, using the linear flow law, with a horizontal length scale of L = 160 km. The coloured lines show the results of our model at different numbers of vertical layers (see legend in Panel A), for the two-point one-sided scheme (panel A), the three-point one-sided scheme (panel B), and the ghost-point scheme (panel C). The small subpanels zoom in on the thin, slow-moving ice around x = L/4.

As before, we perform a convergence analysis by comparing to numerical solutions with n _z = 128 vertical layers. The results of this analysis are shown in Figure 12.

Figure 12. Convergence of the relative error in the surface velocity with the number of nodes, for the combination of the BPA with the linear flow law, for ISMIP-HOM experiment A. The three panels show the results for the three different ways of discretising the zero-stress boundary condition at the ice surface: (A) the two-point one-sided scheme, (B) the three-point one-sided scheme, and (C) the ghost-node scheme. Log-linear curves (dashed lines) are fitted to each set of experiments to calculate the order of convergence R, which is displayed in the legend.

Although the linear flow law removes some of the complexity from the experiment, it is still not straightforward to interpret. As before, the accuracy of the discretisation of the horizontal strain rates is affected by the vertical resolution because of the coordinate transformation, and both the horizontal and vertical strain rates, as well as the vertical resolution, are all laterally varying.

The two-point scheme still shows approximately first-order convergence. Both the three-point scheme and the ghost-node scheme show approximately second-order convergence, with the ghost-node scheme now producing slightly larger errors than the three-point scheme.

4. Conclusions and discussion

We have presented results of experiments solving the SIA and the BPA for ice sheets with idealised geometries. We have investigated different schemes for discretising the zero-stress boundary condition at the ice surface, and combined these schemes with different flow laws: Glen's flow law, which diverges to infinite viscosity when the strain rates approach zero, a linear, nondiverging flow law that predicts finite viscosities everywhere, and an over-regularised variant of Glen's flow law that is nonlinear, but no longer contains a singularity. We find that the two different one-sided finite difference schemes that we tested produce the expected convergence behaviour when combined with the different nondiverging flow laws. However, excepting the five-point scheme, they all produce much larger errors when Glen's flow law is used instead, while being reduced to linear convergence. A ghost-node scheme, which retains both the linear equations for the momentum balance and for the boundary conditions, performs well with all flow laws. These results hold for both the SIA and the BPA, with the caveat that the results of the BPA are less straightforward to interpret for a number of reasons (see Section 3.4.1). A solid mathematical explanation for why the different one-sided schemes perform poorly when the flow law contains a singularity, while the ghost-point scheme does not, remains elusive.

While the five-point one-sided scheme performed well with Glen's flow law in the 1-D SIA experiment, an implementation of this scheme (as well as one with the four-point scheme) in the BPA solver failed to converge in the nonlinear viscosity iteration. We cannot explain why this would be the case. However, as the ghost-node scheme produces similar good results in the SIA experiment, and also works well in the BPA, we do not find it important right now to investigate this issue further. It is worth mentioning here that the number of Picard iterations required to solve for the nonlinear effective viscosity is not significantly different for the different discretisation schemes, nor does there appear to be a significant difference in computation time, in any of the experiments presented here.

We have investigated the ghost-node scheme for the boundary conditions at the ice surface and the ice base. In realistic applications, similar boundary conditions need to be applied at the floating ice front. However, due to the water pressure on the submerged part of the front, the ice front is not stress-free, so that the strain rates will not tend to zero, and the effective viscosity will not diverge to infinity. While it would be interesting to investigate this further in future work, based on our findings here we do not expect the discretisation scheme there to have as much of an effect on the solution as it does at the ice surface.

The available literature on existing ice-sheet models (including both reviewed publications and unreviewed model documentation) rarely provides the discretisation scheme used for the boundary conditions to the momentum balance. However, the one-sided discretisation schemes are much more common in literature on numerical mathematics than the ghost-node scheme, which we have not seen mentioned in any literature on ice-sheet modelling. Our results show that the one-sided schemes, even though they might be expected a priori to produce acceptably small errors at modest numbers of vertical layers, do not always produce such results when combined with a diverging flow law, such as Glen's flow law. Depending on the vertical resolution used by a particular model, this can lead to significant biases in the velocity solution, which could then lead to biases in e.g. estimates of future sea-level contributions. This is especially important because the results of higher-order ice-sheet models, which must include these boundary conditions, are often used as benchmarks for simpler, vertically-integrated ice-sheet models, which do not. We hope these findings will motivate ice-sheet modellers to take extra care in verifying the numerical solvers used in their models, and to publish the results of these efforts.