Hostname: page-component-cd9895bd7-p9bg8 Total loading time: 0 Render date: 2024-12-23T02:03:28.746Z Has data issue: false hasContentIssue false
  • English
  • Français

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

Published online by Cambridge University Press:  04 June 2024

Constantijn J. Berends Open the ORCID record for Constantijn J. Berends [Opens in a new window] ,
Roderik S. W. van de Wal Open the ORCID record for Roderik S. W. van de Wal [Opens in a new window]  and
Paul A. Zegeling
Constantijn J. Berends*
Affiliation:
Institute for Marine and Atmospheric research Utrecht (IMAU), Utrecht University, Utrecht, The Netherlands
Roderik S. W. van de Wal
Affiliation:
Institute for Marine and Atmospheric research Utrecht (IMAU), Utrecht University, Utrecht, The Netherlands Faculty of Geosciences, Department of Physical Geography, Utrecht University, Utrecht, The Netherlands
Paul A. Zegeling
Affiliation:
Faculty of Science, Department of Mathematics, Utrecht University, Utrecht, The Netherlands
*
Corresponding author: Constantijn J. Berends; Email: [email protected]
Rights & Permissions [Opens in a new window]

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.

Keywords

Glacial rheologyice dynamicsice-sheet modelling
Type
Article
Information
Journal of Glaciology , First View , pp. 1 - 15
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
Copyright © The Author(s), 2024. Published by Cambridge University Press on behalf of International Glaciological Society

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, Reference Fox-Kemper2021). 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, Reference Goelzer2020; Seroussi and others, Reference Seroussi2020; Aschwanden and others, Reference Aschwanden, Bartholomaus, Brinkerhoff and Truffer2021). 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, Reference Seroussi2013; Babaniyi and others, Reference Babaniyi, Nicholson, Villa and Petra2021), the subglacial topography (Perego and others, Reference Perego, Price and Stadler2014), the subglacial hydrology (Kazmierczak and others, Reference Kazmierczak, Sun, Coulon and Pattyn2022), the sub-shelf ocean conditions and melt rates (Favier and others, Reference Favier2019; Jourdain and others, Reference Jourdain2020; Berends and others, Reference Berends, Stap and van de Wal2023a), and the subglacial friction and sliding law (Sun and others, Reference Sun2020; Berends and others, Reference Berends, van de Wal, van den Akker and Lipscomb2023b). Another part stems from uncertainty in the future change in climate and mass balance (Goelzer and others, Reference Goelzer2020; Seroussi and others, Reference Seroussi2020). 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, Reference Pattyn2008; Bernales and others, Reference Bernales, Rogozhina, Greve and Thomas2017; Rückamp and others, Reference Rückamp, Kleiner and Humbert2022), and different numerical treatments of basal friction (Feldmann and others, Reference Feldmann, Albrecht, Khroulev, Pattyn and Levermann2014; Leguy and others, Reference Leguy, Lipscomb and Asay-Davis2021) and sub-shelf melt (Cornford and others, Reference Cornford2020; Leguy and others, Reference Leguy, Lipscomb and Asay-Davis2021; Berends and others, Reference Berends, Stap and van de Wal2023a) 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. ub = 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, Reference Pattyn2003; Lipscomb and others, Reference Lipscomb2019), as well as discretisation schemes for the momentum balance itself (e.g. Bueler and Brown, Reference Bueler and Brown2009; Berends and others, Reference Berends, Goelzer, Reerink, Stap and van de Wal2022), 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, Reference Pattyn2008). 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, Reference Morland and Johnson1980). 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, Reference Paterson1994) 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, Reference Bueler and Brown2009; CISM, Lipscomb and others, Reference Lipscomb2019; IMAU-ICE, Berends and others, Reference Berends, Goelzer, Reerink, Stap and van de Wal2022), 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, Reference Arakawa and Lamb1977), 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, Reference Martin2011; Yelmo, Robinson and others, Reference Robinson2019; IMAU-ICE, Berends and others, Reference Berends, Goelzer, Reerink, Stap and van de Wal2022). 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 ua 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 (nz-1)-by-n z, while $M_{\partial /\partial z}^{c\to a}$ is n z-by-(nz-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 (nz-1)-by-(nz-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 L2-norm of the difference between two successive solutions of u is less then 10−8 of the L2-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 (Reference Fornberg2006). 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−2. 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−16 Pa−3 yr−1. 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 L2-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 nz = 8 (red), nz = 16 (blue), nz = 32 (green), and nz = 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 = 214 = 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 101–102 (e.g. 81 layers for SICOPOLIS and 121 layers for PISM in their respective ISMIP6 Antarctica set-ups; Seroussi et al., Reference Seroussiin 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, Reference Pattyn2003) 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, Reference Pattyn2008) 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 (Reference Pattyn2008) 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, Reference Berends, Goelzer, Reerink, Stap and van de Wal2022) showed that the SIA, the hybrid SIA/shallow shelf approximation (Bueler and Brown, Reference Bueler and Brown2009), and the depth-integrated viscosity approximation (Goldberg, Reference Goldberg2011) 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 (Reference Pattyn2008). 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, Reference Pattyn2008), 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 (Reference Pattyn2008) 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 (Reference Pattyn2008) 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 (Reference Pattyn2008) 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 (Reference Pattyn2008) 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.

Supplementary material

The supplementary material for this article can be found at https://doi.org/10.1017/jog.2024.45.

Code and data availability

The Matlab code for performing the experiments and creating the figures presented here, is available in the Supplementary Material.

Acknowledgements

We would like to thank Frank Pattyn, Jorjo Bernales, and Tim van den Akker for their helpful comments during the execution of this project.

Author contributions

CJB performed the experiments and analysed the data. CJB wrote the draft of the manuscript; all authors contributed to the final version.

Financial support

This publication was supported by PROTECT. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 869304, PROTECT contribution number 102.

Competing interests

None.

Appendix A – Irregular vertical grid

Many numerical ice-sheet models use an irregular vertical grid, with the layers being more closely-spaced near the ice base, to more accurately capture the high strain rates there. We found that the discretisation issues presented in this work are present regardless of whether the vertical grid is regular or irregular, and so we chose to use a regular grid in the main work, as the equations are slightly simpler. Here, we demonstrate that using an irregular grid does not improve the results – in fact, for the two-point and three-point one-sided schemes, it arguable leads to even less accurate results.

For these experiments, we replace Eqn (6) with the following expression:

(A1)$$z_a^k = b + H\left({1-\displaystyle{{R^{( k-1) /n_z-1}-1} \over {R-1}}} \right).$$

Here, the grid ratio R approximates the ratio between the first and last grid spacings:

(A2)$$\displaystyle{{z_a^2 -z_a^1 } \over {z_a^{n_z} -z_a^{n_z-1} }}\approx R.$$

For example, for R = 0.1, the grid points at the ice base are spaced approximately 10 times closer than those at the ice surface. Shown in Fig. 13 are the relative velocity errors at the ice surface in the 1-D SIA experiment with Glen's flow law, for different values of R, for the two-point and three-point one-sided schemes, and the ghost-node scheme (we did not derive the expressions for the four-point and five-point one-sided schemes on an irregular grid). All solutions were calculated with n z = 1, 024 grid points.

Figure 13. Relative velocity error at the ice surface for different values of the irregular grid ratio R, for (A) the two-point one-sided scheme, (B) the three-point one-sided scheme, and (C) the ghost-node scheme.

For the two-point and three-point one-sided schemes, the velocity error actually increases for smaller values of R. This is because, for the same number of grid points, using a smaller spacing at the ice base implies using a larger spacing at the ice surface. Therefore, the errors in the discretisation of the surface boundary condition, which dominates the total error, gets larger with decreasing values of R. The ghost-node scheme shows a less clear behaviour, with increasing errors for both large and small values of R, and a tentative optimum around R = 1 (i.e. a regular grid). This suggests that there is no added value in using an irregular grid to solve the momentum balance. Of course, this might be very different for e.g. a thermodynamical model, which typically also has the strongest gradients near the ice base, and does not have any singularities in the solution, so that an irregular grid might well be of added value there.

Shown in Fig. 14 is the convergence of the velocity error with the grid resolution for the SIA experiment with Glen's flow law, for the two-point one-sided scheme, the three-point one-sided scheme, and the ghost-node scheme, for different values of the grid ratio R. This shows that the convergence differences between the different discretisation schemes remain qualitatively unchanged when using an irregular grid.

Figure 14. Convergence of the velocity error with the grid resolution for the SIA experiment with Glen's flow law, for (A) the two-point one-sided scheme, (B) the three-point one-sided scheme, and (C) the ghost-node scheme, for different values of the grid ratio R (red: R = 0.1, green: R = 1, blue: R = 10).

Appendix B – The regularisation term in Glen's flow law

It is common practice for numerical ice-sheet models to add a small regularisation term to the computation of the effective strain rate, so that it never becomes zero, and the effective viscosity therefore never becomes infinite. Here, we investigate the effect of changing the value of the regularisation term $\varepsilon _0^2$ in Eqn (19). Shown in Fig. 15 are the relative velocity errors at the ice surface in the 1-D SIA experiment, for different values of $\varepsilon _0^2$, for the five different discretisation schemes, all calculated with n z = 1,024 grid points. Using a smaller number (e.g. n z = 64) yields qualitatively the same results, with all errors becoming uniformly larger.

Figure 15. Relative error at the ice surface vs the regularisation term ${\boldsymbol \varepsilon }_0^2$, for the five different discretisation schemes.

The convergence rates for the different schemes as a function of the regularisation term are shown in Fig. 16.

Figure 16. Convergence rates with the vertical resolution vs the regularisation term ${\boldsymbol \varepsilon }_0^2$, for the five different discretisation schemes.

The ghost-node scheme and the three-point and four-point schemes yield relative errors that are independent of the regularisation term as long as $\varepsilon _0^2 < 10^{{-}18}$. In this range, they display first-order convergence. For larger values of $\varepsilon _0^2$, they initially seem to produce better results and higher-order convergence, before the errors start to increase again when $\varepsilon _0^2 > 10^{{-}10}$. This is because, for small values of $\varepsilon _0^2$, these schemes underestimate the velocities. When $\varepsilon _0^2 > 10^{{-}10}$, the lower effective viscosity values start producing unphysical viscous deformation in the upper part of the ice column, increasing the velocity. This physical error compensates the discretisation error, apparently reducing the total error. The two-point scheme displays increasing errors when $\varepsilon _0^2 < 10^{{-}30}$, accompanied by a breakdown in the convergence rate. This is likely due to round-off errors, as the effective viscosity increases at the last staggered gridpoint becomes many orders of magnitude larger than elsewhere in the ice column. Interestingly, the five-point scheme shows almost the opposite behaviour, with small errors and second-order convergence when $\varepsilon _0^2 < 10^{{-}30}$, and large errors and poor/no convergence for larger values. It is unclear exactly why this should be the case. Lastly, the ghost-node scheme robustly produces small velocity errors and second-order convergence for all values of $\varepsilon _0^2 < 10^{{-}10}$. For larger values, the physical errors introduced by the reduced effective viscosity become apparent, just as for the other schemes.

Appendix C – Solving the Blatter-Pattyn approximation

C.1 Terrain-following coordinate transformation

In order to solve the BPA, it is desirable that the ice base and ice surface both coincide with a node. This not possible to achieve with a Cartesian coordinate system for any 3-D ice-sheet geometry. We solve this problem by introducing a terrain-following coordinate transformation:

(C1a)$$\hat{x}( {x, \;y, \;z} ) = x, \;$$
(C1b)$$\hat{y}( {x, \;y, \;z} ) = y, \;$$
(C1c)$$\zeta ( {x, \;y, \;z} ) = \displaystyle{{h( {x, \;y} ) -z} \over {H( {x, \;y} ) }}.$$

With this transformation, ζ = 0 at the ice surface, and ζ = 1 at the ice base. Applying the transformation to the gradient operators ∂/∂x, ∂/∂y, ∂/∂z yields:

(C2a)$$\displaystyle{\partial \over {\partial x}} = \displaystyle{\partial \over {\partial \hat{x}}} + \displaystyle{{\partial \zeta } \over {\partial x}}\displaystyle{\partial \over {\partial \zeta }}, \;$$
(C2b)$$\displaystyle{\partial \over {\partial y}} = \displaystyle{\partial \over {\partial \hat{y}}} + \displaystyle{{\partial \zeta } \over {\partial y}}\displaystyle{\partial \over {\partial \zeta }}, \;$$
(C2c)$$\displaystyle{\partial \over {\partial z}} = \displaystyle{{\partial \zeta } \over {\partial z}}\displaystyle{\partial \over {\partial \zeta }}.$$

Applying the chain rule to Eqns (C1a-c), the gradients ∂ζ/∂x, ∂ζ/∂x, ∂ζ/∂x of ζ are:

(C3a)$$\displaystyle{{\partial \zeta } \over {\partial x}} = \displaystyle{1 \over H}\left({\displaystyle{{\partial h} \over {\partial x}}-\zeta \displaystyle{{\partial H} \over {\partial x}}} \right), \;$$
(C3b)$$\displaystyle{{\partial \zeta } \over {\partial y}} = \displaystyle{1 \over H}\left({\displaystyle{{\partial h} \over {\partial y}}-\zeta \displaystyle{{\partial H} \over {\partial y}}} \right), \;$$
(C3c)$$\displaystyle{{\partial \zeta } \over {\partial z}} = \displaystyle{{-1} \over H}.$$

C.2 Discretising the BPA

We solve the BPA on a 3-D grid, which is constructed by vertically extruding a square horizontal grid. In the horizontal plane, we consider a regular a-grid, and a b-grid that is staggered in both the x- and y-directions relative to the a-grid. In the vertical dimension, we consider a regular k-grid, and a staggered ks-grid. Of the four possible combinations this offers, we use three: the ak-grid (horizontally regular, vertically regular), the bk-grid (horizontally staggered, vertically regular), and the bks-grid (horizontally staggered, vertically staggered). This is illustrated in Fig. 17. Note that the index k is oriented positive in ζ, so that k = 1 now lies at the ice surface, and k = n z at the ice base.

Figure 17. The different staggered 3D-grids. Note that the ‘real’ vertical dimension is displayed pointing upwards; because ζ = 0 at the ice surface, this means that vertical layer k + 1 lies below layer k.

The ice thickness H, the bedrock elevation b, and the surface elevation h are defined on the ak-grid; the horizontal ice velocity components u, v are defined on the bk-grid; and the strain rates ∂u/∂x, ∂u/∂y, ∂u/∂z, ∂v/∂x, ∂v/∂y, ∂v/∂z and the effective viscosity η are defined on both the ak-grid and the bks-grid.

The discretised approximation to the BPA contains one linear equation for every degree of freedom, meaning that there are 2n xn yn z linear equations. Note that both the a-grid and the b-grid have n x by n y nodes, which is necessary to implement the horizontal periodic boundary conditions of the ISMIP-HOM experiment. If a simple Dirichlet or Neumann boundary condition were to be used at the lateral domain borders, then the a-grid would have n x by n y nodes, whereas the b-grid would have (n x-1) by (n y-1) nodes, and there would be 2(n x − 1)(n y − 1)n z linear equations to solve for u, v.

We first define the discretised velocity vector υbkq, and the discretised effective viscosity vectors η ak and η bks:

(C4a)$${\boldsymbol \upsilon }_{bkq} = [ {\upsilon_{bkq}^{r_{bkq}( {i, j, k, q} ) \;} } ].$$
(C4b)$$\upsilon _{bkq}^{r_{bkq}( {i, j, k, q} ) \;} = \left\{{\matrix{ {u_{bk}^{r_{bk}( {i, j, k} ) } , \;\;\; q = 1, \;} \cr {v_{bk}^{r_{bk}( {i, j, k} ) } , \;\;\; q = 2, \;} \cr } } \right.$$
(C5)$${\boldsymbol \eta }_{ak} = [ {\eta_{ak}^{r_{ak}( {i, j, k} ) \;} } ] , \;$$
(C6)$${\boldsymbol \eta }_{bks} = [ {\eta_{bks}^{r_{bks}( {i, j, k_s} ) \;} } ].$$

Here, r bkq(i, j, k, q), r ak(i, j, k), and r bks(i, j, k s) are functions that produce a one-to-one mapping between the grid indices i ∈ [1, n x], j ∈ [1, n y], k ∈ [1, n z], k s ∈ [1, n z − 1], the velocity component index q ∈ [1, 2], and the matrix row index r. We choose the following mapping functions:

(C7a)$$r_{ak}( {i, \;j, \;k} ) = n_yn_z( {i-1} ) + n_z( {\,j-1} ) + k, \;$$
(C7b)$$r_{bk}( {i, \;j, \;k} ) = n_yn_z( {i-1} ) + n_z( {\,j-1} ) + k, \;$$
(C7c)$$r_{bks}( {i, \;j, \;k_s} ) = n_y( {n_z-1} ) ( {i-1} ) + ( {n_z-1} ) ( {\,j-1} ) + k_s.$$
(C7d)$$r_{bkq}( {i, \;j, \;k, \;q} ) = 2n_yn_z( {i-1} ) + 2n_z( {\,j-1} ) + 2( {k-1} ) + q.$$

Note though that the choice of these mapping functions does not affect the subsequent discretisation or solving scheme in any way, as long as the mapping is one-to-one, and one is careful to ensure the mapping is applied consistently everywhere.

Using these definitions of the vector forms, we can define the matrices representing the gradient operators in terrain-following coordinates. For example, the coefficients of the matrix $M_{\partial /\partial \hat{x}}^{ak\to bk}$, representing the gradient operator $\partial /\partial \hat{x}$ going from the ak-grid to the bk-grid, are given by:

(C8)$$\eqalign{& M_{\partial /\partial \hat{x}}^{ak\to bk} ( {r_{bk}( {i_{bk}, \;j_{bk}, \;k_{bk}} ) , \;r_{ak}( {i_{ak}, \;j_{ak}, \;k_{ak}} ) } ) \\ & \quad = \left\{\matrix{\displaystyle{{-1} \over {2\Delta x}}, \;\;\, {i_{bk} = i_{ak}, \;j_{bk} = j_{ak}, \;k_{bk} = k_{ak}, \;} \cr\cr \displaystyle{{-1} \over {2\Delta x}}, \;\, {i_{bk} = i_{ak} + 1, \;j_{bk} = j_{ak}, \;k_{bk} = k_{ak}, \;} \cr\cr \displaystyle{1 \over {2\Delta x}}, \;\;\, {i_{bk} = i_{ak}, \;j_{bk} = j_{ak} + 1, \;k_{bk} = k_{ak}, \;} \cr\cr \displaystyle{1 \over {2\Delta x}}, \;\;\, {i_{bk} = i_{ak} + 1, \;j_{bk} = j_{ak} + 1, \;k_{bk} = k_{ak}, \;} \cr\cr 0, \;\;\,\, {{\rm otherwise}.} \hfill} \right.}$$

This represents a simple two-sided differencing scheme. All other gradient operator matrices are similarly defined. Mapping operators between the different grids are likewise represented by matrices. For example, the coefficients of the matrix $M_{{\rm map}}^{ak\to bk}$, which represents the mapping operation from the ak-grid to the bk-grid, are given by:

(C9)$$\eqalign{& M_{{\rm map}}^{ak\to bk} ( {r_{bk}( {i_{bk}, \;j_{bk}, \;k_{bk}} ) , \;r_{ak}( {i_{ak}, \;j_{ak}, \;k_{ak}} ) } )\\ & \quad = \left\{\matrix{\displaystyle{1 \over 4}, \;\;\; {i_{bk} = i_{ak}, \;j_{bk} = j_{ak}, \;k_{bk} = k_{ak}, \;} \cr\cr \displaystyle{1 \over 4}, \;\;\; {i_{bk} = i_{ak} + 1, \;j_{bk} = j_{ak}, \;k_{bk} = k_{ak}, \;} \cr\cr \displaystyle{1 \over 4}, \;\;\; {i_{bk} = i_{ak}, \;j_{bk} = j_{ak} + 1, \;k_{bk} = k_{ak}, \;} \cr\cr \displaystyle{1 \over 4}, \;\; {i_{bk} = i_{ak} + 1, \;j_{bk} = j_{ak} + 1, \;k_{bk} = k_{ak}, \;} \cr\cr 0, \;\;\;\;\, {{\rm otherwise}.} \hfill} \right.}$$

All other mapping operator matrices are similarly defined. Lastly, we define mapping operators that map the velocity components u and v between the bkq-grid and the bk-grid. For example, the coefficients of the matrix $M_{{\rm map}}^{bku\to bk}$, which can be multiplied with the vector υbkq (which contains the values of both u and v on the bk-grid) to give u bk, are given by:

(C10)$$\eqalign{& M_{{\rm map}}^{bku\to bk} ( {r_{bk}( {i_{bk}, \;j_{bk}, \;k_{bk}} ) , \;r_{bkq}( {i_{bk}, \;j_{bk}, \;k_{bk}, \;q} ) } ) \cr & \quad = \left\{{\matrix{ {1, \;\;\;q = 1, \;} \cr \ \quad{0, \;\;\;{\rm otherwise}.} \cr } } \right.}$$

The different velocity components can be mapped between the bkq-grid and the bk-grid by the following matrix multiplications:

(C11a)$${\boldsymbol u}_{bk} = M_{{\rm map}}^{bku\to bk} {\boldsymbol \upsilon }_{bkq}, \;$$
(C11b)$${\boldsymbol v}_{bk} = M_{{\rm map}}^{bkv\to bk} {\boldsymbol \upsilon }_{bkq}, \;$$
(C11c)$${\boldsymbol \upsilon }_{bkq} = M_{{\rm map}}^{bk\to bku} {\boldsymbol u}_{bk} + M_{{\rm map}}^{bk\to bkv} {\boldsymbol v}_{bk}.$$

In order to solve the BPA, we need matrices representing the gradient operators in Cartesian coordinates ∂/∂x, ∂/∂y, ∂/∂z. These can be constructed by combining the matrices representing the gradient operators in terrain-following coordinates $\partial /\partial \hat{x}$, $\;\partial /\partial \hat{y}, \;\;\partial /\partial \zeta$, with the gradients of ζ, to represent the coordinate transformations given in Eqn (A3). For example, the matrix $M_{\partial /\partial x}^{ak\to bk}$, representing the gradient operator ∂/∂x going from the ak-grid to the bk-grid, is obtained by:

(C12)$$M_{\partial /\partial x}^{ak\to bk} = M_{\partial /\partial \hat{x}}^{ak\to bk} + D\left({{\displaystyle{{\partial {\boldsymbol \zeta }} \over {\partial {\boldsymbol x}}}}_{bk}} \right)M_{\partial /\partial \zeta }^{ak\to bk}.$$

All other matrices representing the gradient operators in Cartesian coordinates are obtained similarly. Using these definitions, the discretisation of the BPA is represented by the following matrix equation:

(C13)$$\eqalign{ A_{{\rm eq}1} = & M_{\partial / \partial x}^{ak\to bk} [ {2D( {{\boldsymbol \eta }_{ak}} ) ( {2M_{\partial /\partial x}^{bk\to ak} M_{{\rm map}}^{bku\to bk} + M_{\partial /\partial y}^{bk\to ak} M_{{\rm map}}^{bkv\to bk} } ) } ] \cr & + M_{\partial /\partial y}^{ak\to bk} [ {D( {{\boldsymbol \eta }_{ak}} ) ( {M_{\partial /\partial y}^{bk\to ak} M_{{\rm map}}^{bku\to bk} + M_{\partial /\partial x}^{bk\to ak} M_{{\rm map}}^{bkv\to bk} } ) } ] \cr &+ M_{\partial /\partial z}^{bks\to bk} [ {D( {{\boldsymbol \eta }_{bks}} ) M_{\partial /\partial z}^{bk\to bks} M_{{\rm map}}^{bku\to bk} } ].} $$
(C14)$$\eqalign{ A_{{\rm eq}2} = & M_{\partial /\partial y}^{ak\to bk} [ {2D( {{\boldsymbol \eta }_{ak}} ) ( {2M_{\partial /\partial y}^{bk\to ak} M_{{\rm map}}^{bkv\to bk} + M_{\partial /\partial x}^{bk\to ak} M_{{\rm map}}^{bku\to bk} } ) } ] \cr & + M_{\partial /\partial x}^{ak\to bk} [ {D( {{\boldsymbol \eta }_{ak}} ) ( {M_{\partial /\partial x}^{bk\to ak} M_{{\rm map}}^{bkv\to bk} + M_{\partial /\partial y}^{bk\to ak} M_{{\rm map}}^{bku\to bk} } ) } ] \cr& + M_{\partial /\partial z}^{bks\to bk} [ {D( {{\boldsymbol \eta }_{bks}} ) M_{\partial /\partial z}^{bk\to bks} M_{{\rm map}}^{bkv\to bk} } ].}$$
(C15)$$A = M_{{\rm map}}^{bk\to bku} A_{{\rm eq}1} + M_{{\rm map}}^{bk\to bkv} A_{{\rm eq}2}.$$
(C16)$${\boldsymbol b} = M_{{\rm map}}^{bk\to bku} \left({\rho g\displaystyle{{\partial h} \over {\partial x}}} \right)_{bk} + M_{{\rm map}}^{bk\to bkv} \left({\rho g\displaystyle{{\partial h} \over {\partial y}}} \right)_{bk}.$$
(C17)$$A\upsilon _{bkq} = {\boldsymbol b}. $$

In order to compute the stiffness matrix A, the effective viscosity η and the strain rates need to be computed on both ak-grid and the bks-grid. The horizontal stretch/shear strain rates ∂u/∂x, ∂u/∂y, ∂v/∂x, ∂v/∂y are computed on the ak-grid, and then mapped to the bks-grid:

(C18)$$\displaystyle{{\partial {\boldsymbol u}} \over {\partial {\boldsymbol x}}}_{ak} = M_{\partial /\partial x}^{bk\to ak} M_{{\rm map}}^{bku\to bk} {\boldsymbol \upsilon }_{bkq}, \;$$
(C19)$$\displaystyle{{\partial {\boldsymbol u}} \over {\partial {\boldsymbol x}}}_{bks} = M_{{\rm map}}^{ak\to bks} \displaystyle{{\partial {\boldsymbol u}} \over {\partial {\boldsymbol x}}}_{ak}.$$

The vertical shear strain rates ∂u/∂z, ∂v/∂z are computed on the bks-grid, and then mapped to the ak-grid:

(C20)$$\displaystyle{{\partial {\boldsymbol u}} \over {\partial {\boldsymbol z}}}_{bks} = M_{\partial /\partial z}^{bk\to bks} M_{{\rm map}}^{bku\to bk} {\boldsymbol \upsilon }_{bkq}, \;$$
(C21)$$\displaystyle{{\partial {\boldsymbol u}} \over {\partial {\boldsymbol z}}}_{ak} = M_{{\rm map}}^{bks\to ak} \displaystyle{{\partial {\boldsymbol u}} \over {\partial {\boldsymbol z}}}_{bks}.$$

The effective viscosity η is then calculated separately on both grids, using Eqn (31).

C.3 Boundary conditions to the BPA

As with the SIA, we explore three different ways to implement the zero-stress surface boundary conditions at the ice surface and base, which differ in the way they discretise the vertical shear strain rate ∂u/∂z: a two-point one-sided scheme, a three-point one-sided scheme, and a ghost-point scheme. Recall that the first equation of the zero-stress boundary condition at the ice surface reads:

(C22)$$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.$$

Transforming the vertical shear strain rate ∂u/∂z to terrain-following coordinates, the two-point one-sided scheme is discretised as follows (leaving out the discretisation of the horizontal stretch/shear strain rates for readability):

(C23)$$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 \zeta } \over {\partial z}}\displaystyle{{u_{bk}^{r_{bk}( {i, j, 2} ) } -u_{bk}^{r_{bk}( {i, j, 1} ) } } \over {\Delta \zeta }} = 0.$$

The three-point one-sided scheme reads:

(C24)$$ \eqalign{&\hskip-8pt 2\displaystyle{{\partial h} \over {\partial x}}\left[{2\displaystyle{{\partial u} \over {\partial x}} + \displaystyle{{\partial v} \over {\partial y}}} \right] \cr \quad &+ \displaystyle{{\partial h} \over {\partial y}}\left[{\displaystyle{{\partial u} \over {\partial y}} + \displaystyle{{\partial v} \over {\partial x}}} \right]-\displaystyle{{\partial \zeta } \over {\partial z}}\left[{\displaystyle{{-3} \over {2\Delta \zeta }}u_{bk}^{r_{bk}( {i, j, 1} ) } + \displaystyle{4 \over {2\Delta \zeta }}u_{bk}^{r_{bk}( {i, j, 2} ) } + \displaystyle{{-1} \over {2\Delta \zeta }}u_{bk}^{r_{bk}( {i, j, 3} ) } } \right] = 0.}$$

For the ghost-point scheme, we first expand the vertical shear stress term in the BPA, using the product rule:

(C25)$$\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 \eta } \over {\partial z}}\displaystyle{{\partial u} \over {\partial z}} + \eta \displaystyle{{\partial ^2u} \over {\partial z^2}} \cr & \quad = \rho g\displaystyle{{\partial h} \over {\partial x}}.}$$

As with the SIA, we discretise ∂u/∂z and ∂2u/∂z 2 using standard three-point two-sided schemes:

(C26)$$\left({\displaystyle{{\partial u} \over {\partial z}}} \right)_{bk}^{r_{bk}( {i, j, k} ) } = \displaystyle{{\partial \zeta } \over {\partial z}}\left({\displaystyle{{u_{bk}^{r_{bk}( {i, j, k + 1} ) } -u_{bk}^{r_{bk}( {i, j, k-1} ) } } \over {2\Delta \zeta }}} \right).$$
(C27)$$\left({\displaystyle{{\partial^2u} \over {\partial z^2}}} \right)_{bk}^{r_{bk}( {i, j, k} ) } = \left({\displaystyle{{\partial \zeta } \over {\partial z}}} \right)^2\left({\displaystyle{{u_{bk}^{r_{bk}( {i, j, k + 1} ) } + u_{bk}^{r_{bk}( {i, j, k-1} ) } -2u_{bk}^{r_{bk}( {i, j, k} ) } } \over {\Delta \zeta^2}}} \right).$$

Substituting Eqn (C26) into Eqn (C22) yields the following expression for the value $u_{bk}^{r_{bk}( {i, j, k-1} ) }$ of u at the ghost node at k − 1 (note that, because of the terrain-following coordinate transformation, the ice surface now lies at the first node, rather than the last, as was the case for the SIA):

(C28)$$\eqalign{u_{bk}^{r_{bk}( {i, j, k-1} ) } & = u_{bk}^{r_{bk}( {i, j, k + 1} ) } \cr & \quad -\displaystyle{{2\Delta \zeta } \over {( {\partial \zeta /\partial z} ) }}\left({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]} \right).}$$

Substituting Eqn (C28) into Eqn (C27) to eliminate the ghost node yields the following expression for ∂2u/∂z 2 at the ice surface:

(C29)$$\eqalign{\left({\displaystyle{{\partial^2u} \over {\partial z^2}}} \right)_{bk}^{r_{bk}( {i, j, k} ) } = \left({\displaystyle{{\partial \zeta } \over {\partial z}}} \right)^2\displaystyle{2 \over {\Delta \zeta ^2}}\left[{u_{bk}^{r_{bk}( {i, j, k + 1} ) } - u_{bk}^{r_{bk}( {i, j, k} ) } - \displaystyle{{2\Delta \zeta } \over {( {\partial \zeta /\partial z} ) }}\left({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]} \right)} \right].}$$

Substituting Eqns (A22 and A29) into the BPA yields:

(C30)$$\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 \eta } \over {\partial z}}\left[{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]} \right]\cr & \quad + \left({\displaystyle{{\partial \zeta } \over {\partial z}}} \right)^2\displaystyle{{2\eta } \over {\Delta \zeta ^2}}\left[{u_{bk}^{r_{bk}( {i, j, 2} ) } -u_{bk}^{r_{bk}( {i, j, 1} ) } -\displaystyle{{2\Delta \zeta } \over {( {\partial \zeta /\partial z} ) }}\left({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]} \right)} \right] = \rho g\displaystyle{{\partial h} \over {\partial x}}.} $$

The boundary condition for nonfrozen ice at the ice base reads:

(C31)$$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 }u = 0.$$

Discretising the vertical shear strain rate using the two-point one-sided scheme yields:

(C32)$$\eqalign{ 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 \zeta } \over {\partial z}}\displaystyle{{u_{bk}^{r_{bk}( {i, j, n_z} ) } -u_{bk}^{r_{bk}( {i, j, n_z-1} ) } } \over {\Delta \zeta }} + \displaystyle{\beta \over \eta }u_{bk}^{r_{bk}( {i, j, n_z} ) } = 0.} $$

Similarly, the three-point one-sided scheme yields:

(C33)$$\eqalign{ 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 \zeta } \over {\partial z}}\left[{\displaystyle{3 \over {2\Delta \zeta }}u_{bk}^{r_{bk}( {i, j, n_z} ) } -\displaystyle{4 \over {2\Delta \zeta }}u_{bk}^{r_{bk}( {i, j, n_z-1} ) } + \displaystyle{1 \over {2\Delta \zeta }}u_{bk}^{r_{bk}( {i, j, n_z-2} ) } } \right] + \displaystyle{\beta \over \eta }u_{bk}^{r_{bk}( {i, j, n_z} ) } = 0.} $$

Lastly, the ghost-point scheme yields:

(C34)$$\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 \eta } \over {\partial z}}\left[{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]} \right] + \left({\displaystyle{{\partial \zeta } \over {\partial z}}} \right)^2\displaystyle{{2\eta } \over {\Delta \zeta ^2}} \cr & \quad \left[{u_{bk}^{r_{bk}( {i, j, n_z-1} ) } -u_{bk}^{r_{bk}( {i, j, n_z} ) } + \displaystyle{{2\Delta \zeta } \over {( {\partial \zeta /\partial z} ) }}\left({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]} \right) + \displaystyle{\beta \over \eta }u_{bk}^{r_{bk}( {i, j, n_z} ) } } \right] = \rho g\displaystyle{{\partial h} \over {\partial x}}.} $$

References

Arakawa, A and Lamb, VR (1977) Computational design of the basic dynamical processes of the UCLA general circulation model. Methods in Computational Physics: Advances in Research and Applications 17, 173–265.Google Scholar
Aschwanden, A, Bartholomaus, TC, Brinkerhoff, DJ and Truffer, M (2021) Brief communication: a roadmap towards credible projections of ice sheet contribution to sea level. The Cryosphere 15, 5705–5715.CrossRefGoogle Scholar
Babaniyi, O, Nicholson, R, Villa, U and Petra, N (2021) Inferring the basal sliding coefficient field for the Stokes ice sheet model under rheological uncertainty. The Cryosphere 15, 1731–1750.CrossRefGoogle Scholar
Berends, CJ, Goelzer, H, Reerink, TJ, Stap, LB and van de Wal, RSW (2022) Benchmarking the vertically integrated ice-sheet model IMAU-ICE (version 2.0). Geoscientific Model Development 15, 5667–5688.CrossRefGoogle Scholar
Berends, CJ, Stap, LB and van de Wal, RSW (2023a) Strong impact of sub-shelf melt parameterisation on ice-sheet retreat in idealised and realistic Antarctic topography. Journal of Glaciology, 69, 1–15. doi: 10.1017/jog.2023.33CrossRefGoogle Scholar
Berends, CJ, van de Wal, RSW, van den Akker, T and Lipscomb, WH (2023b) Compensating errors in inversions for subglacial bed roughness: same steady state, different dynamic response. The Cryosphere 17, 1585–1600.CrossRefGoogle Scholar
Bernales, JA, Rogozhina, I, Greve, R and Thomas, M (2017) Comparison of hybrid schemes for the combination of shallow approximations in numerical simulations of the Antarctic Ice Sheet. The Cryosphere 11, 247–265.CrossRefGoogle Scholar
Bueler, E and Brown, J (2009) Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model. Journal of Geophysical Research 114, F03008, doi: 10.1029/2008JF001179CrossRefGoogle Scholar
Cornford, SL and 21 others (2020) Results of the third marine ice sheet model intercomparison project (MISMIP+). The Cryosphere 14, 2283–2301.CrossRefGoogle Scholar
Favier, L and 7 others (2019) Assessment of sub-shelf melting parameterisations using the ocean-ice-sheet coupled model NEMO(v3.6)-Elmer/Ice(v8.3). Geoscientific Model Development 12, 2255–2283.CrossRefGoogle Scholar
Feldmann, J, Albrecht, T, Khroulev, C, Pattyn, F and Levermann, A (2014) Resolution-dependent performance of grounding line motion in a shallow model compared with a full-Stokes model according to the MISMIP3d intercomparison. Journal of Glaciology 60, 353–360.CrossRefGoogle Scholar
Fornberg, B (2006) A pseudospectral fictitious point method for high order initial-boundary value problems. SIAM Journal on Scientific Computing 28, 1716–1729.CrossRefGoogle Scholar
Fox-Kemper, B and 17 others (2021) Ocean, Cryosphere and Sea Level Change. Climate Change 2021: The Physical Science Basis. Contribution of Working Group 1 to the Sixth Assessment Report of the Intergovernmental Panel on Climate change.Google Scholar
Goelzer, H and 41 others (2020) The future sea-level contribution of the Greenland ice sheet: a multi-model ensemble study of ISMIP6. The Cryosphere 14, 3071–3096.CrossRefGoogle Scholar
Goldberg, DN (2011) A variationally derived, depth-integrated approximation to a higher-order glaciological flow model. Journal of Glaciology 57, 157–170.CrossRefGoogle Scholar
Jourdain, NC and 6 others (2020) A protocol for calculating basal melt rates in the ISMIP6 Antarctic ice sheet projections. The Cryosphere 14, 3111–3134.CrossRefGoogle Scholar
Kazmierczak, E, Sun, S, Coulon, V and Pattyn, F (2022) Subglacial hydrology modulates basal sliding response of the Antarctic ice sheet to climate forcing. The Cryosphere 16, 4537–4552.CrossRefGoogle Scholar
Leguy, GR, Lipscomb, WH and Asay-Davis, XS (2021) Marine ice sheet experiments with the Community Ice Sheet Model. The Cryosphere 15, 3229–3253.CrossRefGoogle Scholar
Lipscomb, WH and 14 others (2019) Description and evaluation of the community ice sheet model (CISM) v2.1. Geoscientific Model Development 12, 387–424.CrossRefGoogle Scholar
Martin, MA and 6 others (2011) The Potsdam Parallel Ice Sheet Model (PISM-PIK) – Part 2: Dynamic equilibrium simulation of the Antarctic ice sheet. The Cryosphere 5, 727–740.CrossRefGoogle Scholar
Morland, LW and Johnson, IR (1980) Steady motion of ice sheets. Journal of Glaciology 25, 229–246.CrossRefGoogle Scholar
Paterson, WSB (1994) The Physics of Glaciers, 3rd Edn. Oxford, UK: Pergamon.Google Scholar
Pattyn, F (2003) A new three-dimensional higher-order thermomechanical ice sheet model: basic sensitivity, ice stream development, and ice flow across subglacial lakes. Journal of Geophysical Research 108, 2382. doi: 10.1029/2002JB002329CrossRefGoogle Scholar
Pattyn, F and 20 others (2008) Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP-HOM). The Cryosphere 2, 95–108.CrossRefGoogle Scholar
Perego, M, Price, S and Stadler, G (2014) Optimal initial conditions for coupling ice sheet models to Earth system models. Journal of Geophysical Research: Earth Surface 119, 1894–1917.CrossRefGoogle Scholar
Robinson, A and 5 others (2019) Description and validation of the ice-sheet model Yelmo (version 1.0). Geoscientific Model Development 13, 2805–2823.CrossRefGoogle Scholar
Rückamp, M, Kleiner, T and Humbert, A (2022) Comparison of ice dynamics using full-Stokes and Blatter-Pattyn approximation: application to the Northeast Greenland Ice Stream. The Cryosphere 16, 1675–1696.CrossRefGoogle Scholar
Seroussi, H and 5 others (2013) Dependence of century-scale projections of the Greenland ice sheet on its thermal regime. Journal of Glaciology 59, 1024–1034.CrossRefGoogle Scholar
Seroussi, H and 46 others (2020) ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century. The Cryosphere 14, 3033–3070.CrossRefGoogle Scholar
Seroussi, H and 52 others (in review) Evolution of the Antarctic Ice Sheet over the next three centuries from an ISMIP6 model ensemble. Earth's Future, in review.Google Scholar
Sun, S and 28 others (2020) Antarctic ice sheet response to sudden and sustained ice-shelf collapse (ABUMIP). Journal of Glaciology 66, 891–904.CrossRefGoogle Scholar
Figure 0

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

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

Figure 2

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

Figure 3

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

Figure 4

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

Figure 5

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

Figure 6

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

Figure 7

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.

Figure 8

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

Figure 9

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 10

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.

Figure 11

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.

Figure 12

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.

Figure 13

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

Figure 14

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.

Figure 15

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.

Figure 16

Figure 13. Relative velocity error at the ice surface for different values of the irregular grid ratio R, for (A) the two-point one-sided scheme, (B) the three-point one-sided scheme, and (C) the ghost-node scheme.

Figure 17

Figure 14. Convergence of the velocity error with the grid resolution for the SIA experiment with Glen's flow law, for (A) the two-point one-sided scheme, (B) the three-point one-sided scheme, and (C) the ghost-node scheme, for different values of the grid ratio R (red: R = 0.1, green: R = 1, blue: R = 10).

Figure 18

Figure 15. Relative error at the ice surface vs the regularisation term ${\boldsymbol \varepsilon }_0^2$, for the five different discretisation schemes.

Figure 19

Figure 16. Convergence rates with the vertical resolution vs the regularisation term ${\boldsymbol \varepsilon }_0^2$, for the five different discretisation schemes.

Figure 20

Figure 17. The different staggered 3D-grids. Note that the ‘real’ vertical dimension is displayed pointing upwards; because ζ = 0 at the ice surface, this means that vertical layer k + 1 lies below layer k.

Supplementary material: File

Berends et al. supplementary material

Berends et al. supplementary material
Download Berends et al. supplementary material(File)
File 10 MB