The Utrecht Finite Volume Ice-Sheet Model (UFEMISM) version 2.0 – Part 1: Description and idealised experiments

Abstract

Projecting the anthropogenic mass loss of the Greenland and Antarctic ice sheets requires models that can accurately describe the physics of flowing ice and its interactions with the atmosphere, the ocean, and the solid Earth. As the uncertainty in many of these processes can only be explored by running large numbers of simulations to sample the phase space of possible physical parameters, the computational efficiency and user-friendliness of such a model are just as relevant to its applicability as is its physical accuracy. Here, we present and verify version 2.0 of the Utrecht Finite Volume Ice-Sheet Model (UFEMISM). UFEMISM is a state-of-the-art finite-volume model that applies an adaptive grid in both space and time. Since the first version published 2 years ago, v2.0 has added more accurate approximations to the Stokes flow, more sliding laws, different schemes for calculating the ice thickness rates of change, a more numerically stable time-stepping scheme, more flexible and powerful mesh generation code, and a more generally applicable discretisation scheme. The parallelisation scheme has changed from a shared-memory architecture to distributed memory, enabling the user to utilise more computational resources. The version control system (git) includes automated unit tests and benchmark experiments to aid with model development, as well as automated installation of the required libraries, improving both user comfort and reproducibility of results. The input/output (I/O) now follows the NetCDF-4 standard, including automated remapping between regular grids and irregular meshes, reducing user workload for pre- and post-processing. These additions and improvements make UFEMISM v2.0 a powerful, flexible ice-sheet model that can be used for long palaeoglaciological applications, as well as large ensemble simulations for future projections of ice-sheet retreat, and that is ready to be used for coupling within Earth system models.

1 Introduction

One of the most worrisome and, at the same time, most uncertain possible long-term consequences of anthropogenic climate change is mass loss of the Greenland and Antarctic ice sheets, leading to global sea-level rise (Oppenheimer et al., 2019; van de Wal et al., 2019; Fox-Kemper et al., 2021). Projections of the Antarctic contribution to sea-level rise in 2100 under RCP8.5, which were studied in the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6), range from $-$2.5 to $+$17 $\mathrm{cm}$ (Seroussi et al., 2020), with a possible high-end value of $+$59 $\mathrm{cm}$ (van de Wal et al., 2022) and consequently much more on longer timescales. Part of this large uncertainty stems from poorly constrained physical properties and processes in the Antarctic ice-sheet system, including subglacial conditions (e.g. Kazmierczak et al., 2022; Berends et al., 2023a), basal sliding (Sun et al., 2020), interactions between the ice shelf and the ocean in the sub-shelf cavity (e.g. Burgard et al., 2022; Berends et al., 2023b), calving (e.g. Crawford et al., 2021) and ice-dynamical processes (e.g. Rückamp et al., 2022). However, even in idealised experiments where all these quantities are known perfectly, different ice-sheet models can predict rates of sea-level rise that differ by a factor of 3 (Cornford et al., 2020). This represents the uncertainty arising from the numerical models themselves, which disagree on how to translate the physical equations to computer code. These model differences include the way the momentum balance (typically represented by the Stokes equations) is approximated, the choice of grid, the numerical treatment of discontinuities in basal friction and melt rates at the grounding line, and the way the model is initialised. Sampling both this model-intrinsic uncertainty and the uncertainty in the physical properties and processes of the actual ice sheet requires ice-sheet models that have the computational power and the flexibility to perform large numbers of simulations at an adequate resolution to capture these processes. To meet this challenge, many research groups working on ice-sheet modelling have recently directed their efforts at creating new, more powerful ice-sheet models (e.g. Pattyn, 2017; Hoffman et al., 2018; Quiquet et al., 2018; Lipscomb et al., 2019; Robinson et al., 2020; Berends et al., 2022).

Here, we present version 2.0 of the Utrecht Finite Volume Ice-Sheet Model (UFEMISM). The main distinguishing feature of the model is its dynamic adaptive mesh. This approach was pioneered by Durand et al. (2009) and Gladstone et al. (2010) and has since been applied in BISICLES (Cornford et al., 2013), ISSM (dos Santos et al., 2019), and (in glacier-scale applications) Elmer/ice (Todd et al., 2018). This structure allows the model to resolve the grounding line at high ($<$ 5 $\mathrm{km}$) resolutions during multi-millennial simulations. Since the publication of v1.0, many new features have been added to UFEMISM, and many existing features have been improved in terms of power, flexibility, and user-friendliness. In Sect. , we provide the physical equations for ice flow that are solved by the model. This includes several approximations to the Stokes equations for the momentum balance (Sect. ), several sliding laws (Sect. ), a new numerical scheme for treating basal friction at the grounding line (Sect. ), different temporal discretisation schemes to calculate the ice geometry rates of change (Sect. ), and a new adaptive time-stepping scheme (Sect. ). In Sect. , we describe several improvements that were made to the model code. This includes a change from a shared-memory to a distributed-memory implementation (Sect. ) and a thoroughly reworked I/O module that now follows the NetCDF-4 standard (Unidata, 2023) and is much more flexible and user-friendly (Sect. ). It also includes a version control system that includes automated unit tests and benchmark experiments to aid in developing robust code and automated installation of external libraries to improve user-friendliness and the reproducibility of results (Sect. ). In Sect. , we present the results of a number of idealised-geometry experiments to verify the new model physics and numerics.

This paper, part 1, focuses on the basic mathematics and physics of the model and their verification in idealised benchmark experiments. Part 2, which was submitted for review and publication separately (Bernales et al., 2025), focuses on model additions required for the application of UFEMISM to realistic ice sheets such as those in Greenland and Antarctica. It includes descriptions of the routines for inverting subglacial bed roughness and ocean temperatures in shelf cavities, different sub-shelf melt parameterisations, initialisation approaches, and future projections of mass loss.

2 Model description

2.1 General

UFEMISM is a large-scale ice-sheet model. It solves different approximations of the Stokes equations to calculate the flow velocities of the ice. These are combined with the ice accumulation and loss rates at the surface and the basal and lateral boundaries of the ice sheet to find the thinning and thickening rates of the ice, which are integrated through time to find the evolution of the ice sheet. Note that hereafter, we refer to UFEMISM version 1.0 as “v1.0” and to version 2.0 as “v2.0”.

The main distinguishing feature of UFEMISM compared to many other ice-sheet models is the use of a dynamic adaptive grid. The two-dimensional plane on which the model operates is discretised as an irregular triangular mesh, an example of which is shown in Fig. .

Figure 1

A demo mesh generated by UFEMISM for the Antarctic ice sheet, using a 10 $\mathrm{km}$ grounding-line resolution and a resolution of up to 200 $\mathrm{km}$ for the ice-sheet interior and the open ocean.

[Figure omitted. See PDF]

Earlier research on ice-sheet modelling has shown that the accuracy of a numerical model is particularly sensitive to the spatial resolution of the grid around the grounding line (Durand et al., 2009; Gladstone et al., 2012; Pattyn et al., 2012). There, the discontinuous basal friction, which is non-zero underneath the grounded ice but zero underneath the floating ice, causes strong gradients in englacial stresses and therefore the ice geometry. Although different solutions have been presented in earlier literature to produce accurate results at coarser resolutions (see Sect. ), the resolutions required at the grounding line are still much higher than those needed in the slow-moving interior of the ice sheet, in order to achieve the same level of accuracy in the ice thickness evolution. As the demand for both the temporal coverage and the number of ice-sheet simulations increases, computational efficiency becomes a more important property of ice-sheet models. Using a uniform high resolution over the entire ice sheet when it is only needed in the relatively small area around the grounding line is therefore undesirable. UFEMISM solves this problem by using a mesh that has a high resolution only where needed and a low resolution where possible. This is the “adaptive” part of the mesh. However, as the ice-sheet geometry changes over the course of a model simulation, the location of the grounding line changes as well. This means that, after a while, the grounding line might no longer be located within the high-resolution area of the mesh. A possible solution would be to use a mesh with a high resolution over a wider area, enveloping the expected future migration of the grounding line. While this is a feasible approach for century-scale simulations, doing this for the multi-millennial applications for which UFEMISM is also intended would mean creating a mesh with a very large high-resolution area, thus offsetting the benefits of the adaptive mesh. Instead, UFEMISM periodically checks the mesh fitness to the modelled ice-sheet geometry and, if needed, automatically creates a new mesh that conforms to the new ice-sheet geometry (with a high-resolution area around the new grounding-line position), remapping the model data from the old mesh to the new one. This is the “dynamic” part of the mesh. Berends et al. (2021) showed that this approach results in good computational performance, with no significant loss of accuracy.

While the general principles of the dynamic adaptive mesh have not changed significantly in v2.0 with respect to v1.0, the way these principles are implemented has changed in several ways. The new mesh generation code, the scheme used to discretise the partial differential equations of the model on the mesh, and the scheme used to remap data from one mesh to another are presented in Appendices , , and , respectively.

2.2 Momentum balance

UFEMISM v2.0 includes solvers for several different approximations to the Stokes equations, which neglect increasingly more terms in these equations. Of these approximations, the Blatter–Pattyn approximation (BPA; Blatter, 1995; Pattyn, 2003), which is described in Sect. , neglects the fewest terms. The depth-integrated viscosity approximation (DIVA; Goldberg, 2011; Sect. ), the shallow-shelf approximation (SSA; Morland, 1987; Sect. ), the shallow-ice approximation (SIA; Morland and Johnson, 1980; Sect. ), and the hybrid SIA–SSA (Bueler and Brown, 2009; Sect. ) can all be derived by neglecting more and more terms. For a comprehensive description of the Stokes equations and a derivation of the different approximations, we recommend reading Greve and Blatter (2009).

2.2.1 Blatter–Pattyn approximation

The BPA arises from the Stokes equations by assuming hydrostatic equilibrium and neglecting the stresses arising from horizontal variations in the vertical velocity (i.e. $\frac{\partial w}{\partial x}\ll \frac{\partial u}{\partial z}$, $\frac{\partial w}{\partial y}\ll \frac{\partial v}{\partial z}$; Pattyn, 2003; see Table  for the definitions of the physical quantities and constants). This means that the pressure $p$ and the vertical velocity $w$ disappear as degrees of freedom from the momentum balance so that only the horizontal velocities $uv$ remain to be solved for. The BPA produces ice velocities that are generally very close to those from the Stokes equations (Pattyn et al., 2008) so that it is generally able to describe the large-scale evolution of continental ice-sheet–shelf systems such as the Antarctic ice sheet. Deviations from the full-Stokes solution are more noticeable in e.g. thermomechanically coupled problems of ice streams (Schoof and Mantelli, 2021), advection problems of tracers (Jouvet et al., 2020), and flow at ridges and domes (Seddik et al., 2011). While less computationally expensive to solve than the Stokes equations, the BPA is still significantly slower than DIVA or the hybrid SIA–SSA owing to the fact that, where those approximations either parameterise or neglect vertical variations in horizontal velocities or strain rates, the BPA solves for such variations explicitly. This requires the model to discretise the vertical dimension as well, whereas DIVA and the hybrid SIA–SSA operate in the two-dimensional plane, yielding a system of linear equations that is larger by a factor equal to the number of vertical layers in the model.

Table 1

Model symbols, units, and default values where applicable.

Symbol	Description	Units	Value
$A$	Glen's flow law factor	${\mathrm{Pa}}^{-\mathrm{n}}{\mathrm{yr}}^{-\mathrm{1}}$
$b$	Bedrock elevation	$\mathrm{m}$
$\mathit{\beta }$	Basal friction coefficient	$\mathrm{Pa}{\mathrm{m}}^{-\mathrm{1}}\mathrm{yr}$
${\dot{\mathit{\epsilon }}}_{\mathrm{e}}$	Effective strain rate	${\mathrm{yr}}^{-\mathrm{1}}$
$g$	Gravitational acceleration	$\mathrm{m}{\mathrm{s}}^{-\mathrm{2}}$	9.81
$H$	Ice thickness	$\mathrm{m}$
$\mathit{\eta }$	Effective viscosity	$\mathrm{Pa}\mathrm{yr}$
$n$	Glen's flow law exponent		3
$\mathit{\rho }$	Density of ice	$\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$	910
$s$	Surface elevation	$\mathrm{m}$
$\mathit{u}$	Horizontal ice velocity vector	$\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$
$u$	Horizontal ice velocity in $x$ direction	$\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$
$v$	Horizontal ice velocity in $y$ direction	$\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$
$w$	Vertical ice velocity	$\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$
$\mathit{\zeta }$	Scaled vertical coordinate		0–1

The set of partial differential equations that must be solved in order to find the three-dimensional horizontal ice velocities $uv$ reads $$\begin{array}{}\text{1a}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial }{\partial x}}\left[\mathrm{2}\mathit{\eta }\left(\mathrm{2}{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right)\right]& +{\displaystyle \frac{\partial }{\partial y}}\left[\mathit{\eta }\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right]\\ & +{\displaystyle \frac{\partial }{\partial z}}\left[\mathit{\eta }{\displaystyle \frac{\partial u}{\partial z}}\right]=\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial x}},\end{array}\text{1b}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial }{\partial y}}\left[\mathrm{2}\mathit{\eta }\left(\mathrm{2}{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial u}{\partial x}}\right)\right]& +{\displaystyle \frac{\partial }{\partial x}}\left[\mathit{\eta }\left({\displaystyle \frac{\partial v}{\partial x}}+{\displaystyle \frac{\partial u}{\partial y}}\right)\right]\\ & +{\displaystyle \frac{\partial }{\partial z}}\left[\mathit{\eta }{\displaystyle \frac{\partial v}{\partial z}}\right]=\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial y}}.\end{array}\end{array}$$

The effective viscosity $\mathit{\eta }$ is related to the effective strain rate ${\dot{\mathit{\epsilon }}}_{\mathrm{e}}$ by Glen's flow law (Glen, 1955): 2$$\mathit{\eta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\mathrm{1}/n}{\dot{\mathit{\epsilon }}}_{\mathrm{e}}^{\frac{\mathrm{1}-n}{n}}.$$

The flow factor can be set to a uniform fixed value (as is done in the idealised experiments presented here) or can be calculated from the ice temperature, following the Arrhenius-type relation provided in Berends et al. (2021).

The effective strain rate ${\dot{\mathit{\epsilon }}}_{\mathrm{e}}$ is given by 3$$\begin{array}{rl}{\dot{\mathit{\epsilon }}}_{\mathrm{e}}=& [{\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v}{\partial y}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)}^{\mathrm{2}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u}{\partial z}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial v}{\partial z}}\right)}^{\mathrm{2}}{]}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

At the ice surface, the zero-stress boundary condition reads $$\begin{array}{}\text{4a}& {\displaystyle }& {\displaystyle }\mathrm{2}{\displaystyle \frac{\partial s}{\partial x}}\left(\mathrm{2}{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right)+{\displaystyle \frac{\partial s}{\partial y}}\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)-{\displaystyle \frac{\partial u}{\partial z}}=\mathrm{0},\text{4b}& {\displaystyle }& {\displaystyle }\mathrm{2}{\displaystyle \frac{\partial s}{\partial y}}\left(\mathrm{2}{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial u}{\partial x}}\right)+{\displaystyle \frac{\partial s}{\partial x}}\left({\displaystyle \frac{\partial v}{\partial x}}+{\displaystyle \frac{\partial u}{\partial y}}\right)-{\displaystyle \frac{\partial v}{\partial z}}=\mathrm{0}.\end{array}$$

A similar dynamical boundary condition at the ice base includes a basal friction term: $$\begin{array}{}\text{5a}& {\displaystyle }\begin{array}{rl}\mathrm{2}{\displaystyle \frac{\partial b}{\partial x}}\left(\mathrm{2}{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right)& +{\displaystyle \frac{\partial b}{\partial y}}\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\\ & -{\displaystyle \frac{\partial u}{\partial z}}+{\displaystyle \frac{\mathit{\beta }}{\mathit{\eta }}}u=\mathrm{0},\end{array}\text{5b}& {\displaystyle }\begin{array}{rl}\mathrm{2}{\displaystyle \frac{\partial b}{\partial y}}\left(\mathrm{2}{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial u}{\partial x}}\right)& +{\displaystyle \frac{\partial b}{\partial x}}\left({\displaystyle \frac{\partial v}{\partial x}}+{\displaystyle \frac{\partial u}{\partial y}}\right)\\ & -{\displaystyle \frac{\partial v}{\partial z}}+{\displaystyle \frac{\mathit{\beta }}{\mathit{\eta }}}v=\mathrm{0}.\end{array}\end{array}$$

UFEMISM currently does not include a stress boundary condition at the ice front for any of the momentum balance approximations. Instead, it uses the “infinite-slab” approach, where the momentum balance is solved on the entire grid, including ice-free cells (which the solvers assume are covered by a very thin (0.1 $\mathrm{m}$ by default) layer of ice), and a simple Neumann boundary condition is applied at the domain boundary. While we have not (yet) tested this in UFEMISM, earlier experiments with IMAU-ICE (Berends et al., 2022) showed that the effect of this approach on the velocity solution is generally very small. It should also be noted that this approach is only used in the momentum balance; when solving for conservation of mass, the model does account for ice-free vertices so that a calving front is explicitly included.

In order to solve the BPA, the vertical dimension must be discretised as well. This is not straightforward, as the surface and base of the ice are generally not flat and evolve over time so that these surfaces will generally move in between grid points in the vertical. In UFEMISM, as in most other ice-sheet models that solve a three-dimensional version of momentum balance, this problem is solved by introducing a terrain-following coordinate transformation, which is described in Appendix .

2.2.2 Depth-integrated viscosity approximation

DIVA, which is the default option for the momentum balance approximation in v2.0, arises by neglecting the stresses that arise from vertical variations in the horizontal strain rates in the BPA (i.e. $\frac{\partial }{\partial z}(\frac{\partial u}{\partial x},\frac{\partial u}{\partial y},\frac{\partial v}{\partial x},\frac{\partial v}{\partial y})\approx \mathrm{0})$ and integrating the resulting equations vertically. This means that, while the BPA is solved in three dimensions, DIVA operates in the two-dimensional plane, greatly reducing the computational expense of solving it. In the Ice Sheet Model Intercomparison Project for Higher-Order Models (ISMIP-HOM; Pattyn et al., 2008) experiments, DIVA produces velocities that agree well with the Stokes solution down to horizontal scales of about 20 $\mathrm{km}$ for basal topographical features (Berends et al., 2022; Robinson et al., 2022; this study, Sect. ).

The partial differential equations of DIVA read $$\begin{array}{}\text{6a}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial }{\partial x}}\left[\mathrm{2}\mathit{\eta }H\left(\mathrm{2}{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right)\right]& +{\displaystyle \frac{\partial }{\partial y}}\left[\mathit{\eta }H\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right]\\ & -{\mathit{\beta }}_{\text{eff}}u=\mathit{\rho }gH{\displaystyle \frac{\partial s}{\partial x}},\end{array}\text{6b}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial }{\partial y}}\left[\mathrm{2}\mathit{\eta }H\left(\mathrm{2}{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial u}{\partial x}}\right)\right]& +{\displaystyle \frac{\partial }{\partial x}}\left[\mathit{\eta }H\left({\displaystyle \frac{\partial v}{\partial x}}+{\displaystyle \frac{\partial u}{\partial y}}\right)\right]\\ & -{\mathit{\beta }}_{\text{eff}}v=\mathit{\rho }gH{\displaystyle \frac{\partial s}{\partial y}}.\end{array}\end{array}$$

Here, ${\mathit{\beta }}_{\text{eff}}$ is a term describing both basal friction and vertical shear stress:

7$${\mathit{\beta }}_{\text{eff}}={\displaystyle \frac{\mathit{\beta }}{\mathrm{1}+\mathit{\beta }{F}_{\mathrm{2}}}}.$$

The integral term $F_{\mathrm{2}}$, which can be thought of as a (scaled) depth integral of the inverse viscosity, is defined as 8$$F_n=\underset{b}{\overset{h}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\eta }}}{\left({\displaystyle \frac{s-z}{H}}\right)}^n\mathrm{d}z.$$

Note that, in Eq. (), $n=\mathrm{2}$; Eq. () lists the general form because elsewhere in DIVA, $F_{\mathrm{1}}$ appears as well. A comprehensive derivation of these and other required equations, including a step-by-step approach for how to solve them numerically, can be found in Lipscomb et al. (2019).

While the mathematical derivation is too cumbersome to include here, it can be shown that, in the absence of horizontal strain (i.e. $\frac{\partial u}{\partial x}\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}\frac{\partial v}{\partial y}=\mathrm{0})$, DIVA is identical to the SIA. In a preliminary experiment, we used DIVA to perform the moving-margin experiment from EISMINT-1 (Huybrechts et al., 1996), which describes a roughly Greenland-sized, idealised, circular, polythermal (though not thermomechanically coupled) ice sheet lying on a flat bed, achieving a steady state through a simple, spatially variable mass balance. The resulting ice sheet, which is dominated by vertical shearing, was nearly identical to that produced by the SIA, with a difference of only a few metres in ice thickness, concentrated near the ice divide and the ice margin.

A more practical advantage of DIVA, which has previously been pointed out by Robinson et al. (2022), is that the system of linear equations that must be solved (Eq. ) is almost identical to that of the SSA (Eq. ). Ice-sheet models that already contain code to solve the SSA can therefore be altered to solve DIVA instead with relatively little effort, including only a few simple calculations to evaluate Eq. (), altering the friction term that enters into the system of linear equations.

2.2.3 Shallow-shelf approximation

The SSA arises by neglecting all vertical variations in the BPA (i.e. $\frac{\partial u}{\partial z}\frac{\partial v}{\partial z}\approx \mathrm{0})$, leaving only the membrane stresses, and then vertically integrating the result. This is generally accepted to be a valid approximation in areas of negligible basal shear stress, such as ice shelves, and in well-lubricated, fast-flowing ice streams.

The partial differential equations of the SSA read $$\begin{array}{}\text{9a}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial }{\partial x}}\left[\mathrm{2}\mathit{\eta }H\left(\mathrm{2}{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right)\right]& +{\displaystyle \frac{\partial }{\partial y}}\left[\mathit{\eta }H\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right]\\ & -\mathit{\beta }u=\mathit{\rho }gH{\displaystyle \frac{\partial s}{\partial x}},\end{array}\text{9b}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial }{\partial y}}\left[\mathrm{2}\mathit{\eta }H\left(\mathrm{2}{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial u}{\partial x}}\right)\right]& +{\displaystyle \frac{\partial }{\partial x}}\left[\mathit{\eta }H\left({\displaystyle \frac{\partial v}{\partial x}}+{\displaystyle \frac{\partial u}{\partial y}}\right)\right]\\ & -\mathit{\beta }v=\mathit{\rho }gH{\displaystyle \frac{\partial s}{\partial y}}.\end{array}\end{array}$$

Neglecting the same strain rates reduces the expression for the effective strain rate that is used in Glen's flow law to

10$$\begin{array}{rl}{\dot{\mathit{\epsilon }}}_{\mathrm{e}}=& [{\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v}{\partial y}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial v}{\partial y}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)}^{\mathrm{2}}{]}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

It should be noted that v1.0 further simplified Eq. () by neglecting the gradients in the effective viscosity (after expanding the derivative outside the square brackets using the product rule). While this made the numerical solver more stable (and also significantly faster, requiring fewer non-linear viscosity iterations to converge), it was later discovered that this could lead to significant errors in the velocity and the ice thickness evolution. V2.0 therefore solves the SSA (and DIVA) without this simplification, gaining physical accuracy at the cost of computational performance. Including these additional terms necessitated a change in the discretisation scheme so that in v2.0 the ice velocities are defined on the triangle centres, whereas in v1.0 they were defined on the edges. The new discretisation scheme is presented in more detail in Appendix .

2.2.4 Shallow-ice approximation

The SIA arises by neglecting the membrane stresses in the BPA (i.e. $\frac{\partial u}{\partial x}\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}\frac{\partial v}{\partial y}\approx \mathrm{0})$, leaving only the vertical shear strain rates. This is generally accepted to be a valid approximation for the thick, slow-moving ice in the interior of the Greenland and Antarctic ice sheets, where the flow is dominated by deformation due to vertical shearing, rather than by basal sliding. These assumptions simplify the Stokes equations to $$\begin{array}{}\text{11a}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial u}{\partial z}}\right)=\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial x}},\text{11b}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial v}{\partial z}}\right)=\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial y}}.\end{array}$$

Similarly, the effective strain rate that is used in Glen's flow law reduces to

12$${\dot{\mathit{\epsilon }}}_{\mathrm{e}}={\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u}{\partial z}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial v}{\partial z}}\right)}^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}.$$

Substituting Eq. () into Eq. (), and assuming a stress-free boundary condition at the ice surface and a no-slip boundary condition at the ice base, leads to the following analytical solution for the vertical profile of the horizontal ice velocity: 13$$\mathbf{u}\left(z\right)=-\mathrm{2}\left(\mathit{\rho }g{)}^n{\left|\mathrm{\nabla }s\right|}^{n-\mathrm{1}}\mathrm{\nabla }s\underset{b}{\overset{z}{\int }}A\right(T\left(\mathit{\zeta }\right)\left)\right(s-\mathit{\zeta }{)}^n\mathrm{d}\mathit{\zeta }.$$

Note that it is not possible to include a sliding law in UFEMISM v2.0 when using only the SIA; for this, the hybrid SIA–SSA must be used (see Sect. ).

2.2.5 Hybrid shallow-ice and shallow-shelf approximation

In this approach, the SIA and SSA are solved separately, following the approach proposed by Bueler and Brown (2009). Based on the observation that the flow regime in most areas of an ice sheet is generally dominated by either vertical shear (described by the SIA) or basal sliding (described by the SSA), the two solutions are then simply added together to find an approximation of the flow of the entire ice sheet. This approach produces accurate results in terms of large-scale ice flow (e.g. Bueler and Brown, 2009; Berends et al., 2022) but starts to deviate significantly from the Stokes solution earlier than DIVA as the length scale decreases (Berends et al., 2022; this study).

2.2.6 Vertical velocities

The assumption that glacial ice is incompressible is expressed mathematically as

14$${\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial w}{\partial z}}=\mathrm{0}.$$

The BPA, DIVA, and the (hybrid) SIA–SSA only solve for the horizontal velocities $uv$. From those, the horizontal divergence $\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}$ can be calculated. Integrating Eq. () in the vertical dimension then yields the vertical velocity $w$: 15$$w\left(z\right)=w(z=b)-\underset{b}{\overset{z}{\int }}\left({\displaystyle \frac{\partial u}{\partial x}}\left(\mathit{\zeta }\right)+{\displaystyle \frac{\partial v}{\partial y}}\left(\mathit{\zeta }\right)\right)\mathrm{d}\mathit{\zeta }.$$

Here too, the terrain-following coordinate transformation must be applied before evaluating the vertical integral. The way this is done in UFEMISM is described in Appendix .

2.3 Sliding laws

UFEMISM v2.0 includes a number of different sliding laws for the user to choose from, which relate the basal shear stress ${\mathit{\tau }}_{\mathrm{b}}$ to the basal velocity $u_{\mathrm{b}}$ through the basal friction coefficient $\mathit{\beta }=\frac{\left|{\mathit{\tau }}_{\mathrm{b}}\right|}{\left|u_{\mathrm{b}}\right|}$. All sliding laws are presented here as they are coded in the model, with the basal friction coefficient $\mathit{\beta }$ expressed as a function of the basal speed $u_{\mathrm{b}}=\left|u_{\mathrm{b}}\right|$. The first option is a Weertman-type sliding law (Weertman, 1957):

16$$\mathit{\beta }=C_{\mathrm{w}}{u}_{\mathrm{b}}^{\frac{\mathrm{1}}{m}-\mathrm{1}}.$$

Here, $C_{\mathrm{w}}$ represents the (spatially variable) subglacial bed roughness.

The second option is a Coulomb-type sliding law (Iverson et al., 1998): 17$$\mathit{\beta }=N\text{tan}\mathit{\phi }{u}_{\mathrm{b}}^{-\mathrm{1}}.$$

Here, $N$ is the effective pressure between the ice and the bedrock, which is equal to the ice overburden pressure minus the subglacial water pressure. Currently, the subglacial water pressure is defined simply as 96 % of the ice overburden pressure, following Winkelmann et al. (2011), optionally scaled with a bedrock elevation-dependent parameterisation developed for Antarctica by Martin et al. (2011); the addition of a more elaborate subglacial hydrology model to UFEMISM is planned for future work. The (spatially variable) till friction angle $\mathit{\phi }$ is a measure of the subglacial bed roughness.

The third option is a Budd-type sliding law, proposed by Bueler and van Pelt (2015): 18$$\mathit{\beta }=N\text{tan}\mathit{\phi }{\displaystyle \frac{u_{\mathrm{b}}^{q-\mathrm{1}}}{u_{\mathrm{0}}^q}}.$$

Here, $u_{\mathrm{0}}$ is a transition velocity, with a default (configurable) value of 100 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$. Note that this is a Budd-type sliding law (i.e. a power-law dependence on velocity, scaled with the effective pressure) for the current choice of exponent $q=\mathrm{0.3}$; for $q=\mathrm{1}$, this becomes a regularised Coulomb sliding law, with no dependence on velocity. This sliding law was the only option in UFEMISM v1.0 (Berends et al., 2021).

The fourth option is the hybrid sliding law proposed by Tsai et al. (2015), as formulated by Asay-Davis et al. (2016): 19$$\mathit{\beta }=\text{min}\left({\mathit{\alpha }}^{\mathrm{2}}N,{\mathit{\beta }}^{\mathrm{2}}{u}_{\mathrm{b}}^{\frac{\mathrm{1}}{m}}\right)u_{\mathrm{b}}^{-\mathrm{1}}.$$

Note that here the (spatially variable) subglacial bed roughness is described by two separate parameters: ${\mathit{\alpha }}^{\mathrm{2}}$ for the Coulomb-type part of the friction and ${\mathit{\beta }}^{\mathrm{2}}$ (which is not the square of the basal friction coefficient $\mathit{\beta }$ but a confusingly named separate entity, which we maintain for the sake of consistency with earlier literature) for the Weertman-type part.

The final option is the hybrid sliding law proposed by Schoof (2005), as formulated by Asay-Davis et al. (2016): 20$$\mathit{\beta }={\displaystyle \frac{{\mathit{\beta }}^{\mathrm{2}}{u}_{\mathrm{b}}^{\frac{\mathrm{1}}{m}}{\mathit{\alpha }}^{\mathrm{2}}N}{{\left[{\mathit{\beta }}^{\mathrm{2}m}{u}_{\mathrm{b}}+({\mathit{\alpha }}^{\mathrm{2}}N{)}^m\right]}^{\frac{\mathrm{1}}{m}}}}{u}_{\mathrm{b}}^{-\mathrm{1}}.$$

Note that the terms on the right-hand side of Eq. () are again the ${\mathit{\beta }}^{\mathrm{2}}$ term from Eq. (). In the idealised-geometry experiments presented here, the bed roughness is spatially uniform. For applications to realistic ice sheets, UFEMISM v2.0 includes routines for inverting the bed roughness by nudging. These are presented in part 2 of this work (Bernales et al., 2025).

2.4 Sub-grid friction scaling

UFEMISM v1.0 used a grounding-line flux condition (Schoof, 2007; Pollard and DeConto, 2012) to improve grounding-line migration. While the flux condition generally seems to produce more accurate results in unbuttressed geometries (e.g. Pattyn et al., 2012), extending this solution to geometries where buttressing plays a significant role has proved problematic (Reese et al., 2018). Furthermore, while the implementation of this scheme in v1.0 performed well in idealised geometries, it frequently resulted in numerical instability in the more complex geometries encountered in e.g. the Antarctic ice sheet. Therefore, in UFEMISM v2.0 the flux condition has been replaced by a sub-grid friction scaling scheme, following the approach used in ISSM (Seroussi et al., 2014), PISM (Feldmann et al., 2014), Elmer/ice (Gagliardini et al., 2016), CISM (Leguy et al., 2021), and IMAU-ICE (Berends et al., 2022). Here, the area fraction of each mesh triangle (where the velocities are defined), which is covered by grounded ice, is calculated by bilinearly interpolating the thickness above floatation on the three vertices spanning the triangle. The basal friction coefficient $\mathit{\beta }$, which is calculated using the sliding law, is then multiplied by this grounded fraction before being used to solve the momentum balance. This approach is much more numerically stable, does not require any special treatment to include buttressing, and works well in both idealised and realistic geometries.

2.5 Conservation of energy

The way the heat equation inside the ice is approximated and discretised is unchanged from UFEMISM v1.0. The approximation, which is based on Greve (1997), includes terms describing horizontal and vertical advection, vertical diffusion, and internal strain heating, with the annual mean temperature of the ocean and atmosphere and the geothermal heat flux serving as boundary conditions. Horizontal diffusion and the possible formation of liquid water inside the ice column are neglected. The governing equations and their discretisation (which uses an explicit scheme for the horizontal advective terms and an implicit scheme for the vertical advective and diffusive terms) are provided and verified in the EISMINT-1 benchmark experiments (Huybrechts et al., 1996) by Berends et al. (2021). Unless otherwise specified by the user, ice temperature affects the ice flow factor through an Arrhenius-type relation, following Huybrechts (1992).

2.6 Conservation of mass

After the momentum balance has been solved to find the ice velocities, the condition of conservation of mass can be used to find the rates of change in the ice geometry. Conservation of ice mass for a shallow layer of incompressible ice in the two-dimensional plane is expressed mathematically as

21$${\displaystyle \frac{\partial H}{\partial t}}=-\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{\mathit{u}}H\right)+m.$$

Here, $m$ is the net mass balance, including terms at the ice base, the ice surface, and the lateral boundaries, while $\stackrel{\mathrm{‾}}{u}$ is the vertically averaged, horizontally vector-valued ice velocity. Since UFEMISM always assumes a uniform, constant ice density, vertical variations in the horizontal velocities are not needed to solve the continuity equation. Equation () is discretised spatially using the finite-volume scheme that gave UFEMISM its name, which is derived in Appendix , resulting in the following equation: 22$${{\displaystyle \frac{\partial H}{\partial t}}}^{\mathrm{i}}=-{\mathbf{M}}_{\text{divQ}}{H}^{\mathrm{i}}+m^i.$$

Here, the ice thickness vector ${\mathit{H}}^{\mathrm{i}}$ contains the values of $H$ on all the vertices $i$. ${\mathbf{M}}_{\text{divQ}}$ is a matrix whose coefficients depend on the mesh geometry and the ice velocities, which can be multiplied by the ice thickness vector ${\mathit{H}}^{\mathrm{i}}$ to find the ice flux divergence $\mathrm{\nabla }\cdot (uH{)}^{\mathrm{i}}={\mathbf{M}}_{\text{divQ}}{\mathit{H}}^{\mathrm{i}}$. UFEMISM v2.0 offers three different options to discretise Eq. () temporally: an explicit scheme, an implicit scheme, and a semi-implicit scheme. In all three cases, the thickness rate of change $\frac{\partial H}{\partial t}$ is discretised using a simple first-order scheme. In the explicit scheme, all terms on the right-hand side of Eq. () are defined at time $t$: 23$${\displaystyle \frac{H^{i,t+\mathrm{\Delta }t}-H^{i,t}}{\mathrm{\Delta }t}}=-{\mathbf{M}}_{\text{divQ}}{H}^{i,t}+m^{i,t}.$$

Rearranging the terms yields the following expression, which can be evaluated to find $H^{i,t+\mathrm{\Delta }t}$: 24$$H^{i,t+\mathrm{\Delta }t}=(I-{\mathbf{M}}_{\text{divQ}}\mathrm{\Delta }t)H^{i,t}+m^{i,t}\mathrm{\Delta }t.$$

In the implicit scheme, the ice thickness on the right-hand side of Eq. () is defined at time $t+\mathrm{\Delta }t$: 25$${\displaystyle \frac{H^{i,t+\mathrm{\Delta }t}-H^{i,t}}{\mathrm{\Delta }t}}=-{\mathbf{M}}_{\text{divQ}}{H}^{i,t+\mathrm{\Delta }t}+m^{i,t}.$$

Rearranging the terms yields the following matrix equation that can be solved for $H^{i,t+\mathrm{\Delta }t}$: 26$$(I+{\mathbf{M}}_{\text{divQ}}\mathrm{\Delta }t)H^{i,t+\mathrm{\Delta }t}=H^{i,t}+m^{i,t}\mathrm{\Delta }t.$$

Lastly, the semi-implicit scheme is derived by defining the ice thickness on the right-hand side of Eq. () as the weighted average of $H^{i,t}$ and $H^{i,t+\mathrm{\Delta }t}$: 27$$\begin{array}{rl}{\displaystyle \frac{H^{i,t+\mathrm{\Delta }t}-H^{i,t}}{\mathrm{\Delta }t}}=& -{\mathbf{M}}_{\text{divQ}}(f_{\mathrm{s}}{H}^{i,t+\mathrm{\Delta }t}\\ & +(\mathrm{1}-f_{\mathrm{s}})H^{i,t})+m^{i,t}.\end{array}$$

Here, using a coefficient $f_{\mathrm{s}}=\mathrm{0}$ implies a fully explicit scheme, $f_{\mathrm{s}}=\mathrm{1}$ implies a fully implicit scheme, $\mathrm{0}<f_{\mathrm{s}}<\mathrm{1}$ implies a semi-implicit scheme, and $f_{\mathrm{s}}>\mathrm{1}$ implies an over-implicit scheme. Rearranging the terms yields the following matrix equation that can be solved for $H^{i,t+\mathrm{\Delta }t}$: 28$$\begin{array}{rl}(I+f_{\mathrm{s}}{\mathbf{M}}_{\text{divQ}}\mathrm{\Delta }t)H^{i,t+\mathrm{\Delta }t}& =(I-\mathrm{\Delta }t(\mathrm{1}-f_{\mathrm{s}}\left){\mathbf{M}}_{\text{divQ}}\right)\\ & \times H^{i,t}+m^{i,t}\mathrm{\Delta }t.\end{array}$$

Note that the (semi-)implicit schemes are only implicit in terms of the ice thickness. The flux divergence is still computed based on the velocity solution at the time step, making even the implicit scheme technically semi-implicit. Recent work by Bueler (2023) has looked into the possibilities of using fully implicit solvers for coupled momentum–mass conservation equations, but these solvers have not (yet) been implemented in UFEMISM.

2.7 Time stepping

In v2.0, we use the predictor–corrector (PC) time-stepping scheme by Robinson et al. (2020). While the SIA and SSA both have well-defined critical time steps, no such condition has yet been derived for DIVA or the BPA. The PC scheme essentially operates by calculating two solutions of $H^{i,t+\mathrm{\Delta }t}$: one with an explicit solution of the ice velocity $u$ and one with a pseudo-implicit solution. The difference between the two solutions of $H^{i,t+\mathrm{\Delta }t}$ is a measure of the temporal discretisation error, which can be used to adapt the time step; if the error is found to be increasing, the time step is decreased and vice versa. Robinson et al. (2022) showed that this scheme is particularly suitable to DIVA (and, by extension, the BPA), where the error is less sensitive to larger time steps than in the hybrid SIA–SSA due to the weaker dependence of the velocity on the local surface slope.

A time step in the PC scheme consists of three parts: the predictor step, the update step, and the corrector step. First, in the predictor step, the “predicted” ice thickness is calculated, based on the current ice thickness and the current velocity solution:

29$$\begin{array}{rl}{H}_{\text{pred}}^{t+\mathrm{\Delta }t}& =H^t+\mathrm{\Delta }{t}^t\left[\right(\mathrm{1}+{\displaystyle \frac{{\mathit{\zeta }}_t}{\mathrm{2}}}\left){\displaystyle \frac{\partial H}{\partial t}}\right(H^t,u^t)\\ & -{\displaystyle \frac{{\mathit{\zeta }}_t}{\mathrm{2}}}{\displaystyle \frac{\partial H}{\partial t}}(H^{t-\mathrm{\Delta }t},{\mathit{u}}^{t-\mathrm{\Delta }t})].\end{array}$$

Here, ${\mathit{\zeta }}_t=\frac{\mathrm{\Delta }{t}^t}{\mathrm{\Delta }{t}^{t-\mathrm{\Delta }t}}$ is the ratio between the current and the previous time steps.

Then, in the update step, a new ice velocity solution is calculated for the predicted ice thickness: 30$${\mathit{u}}^{t+\mathrm{\Delta }t}=\mathit{u}\left(H_{\text{pred}}^{t+\mathrm{\Delta }t}\right).$$

Lastly, in the corrector step, the “corrected” ice thickness is calculated, based on the current ice thickness and the new velocity solution: 31$$\begin{array}{rl}{H}_{\text{corr}}^{t+\mathrm{\Delta }t}& =H^t+{\displaystyle \frac{\mathrm{\Delta }{t}^t}{\mathrm{2}}}\left[{\displaystyle \frac{\partial H}{\partial t}}\right(H^t,u^t)\\ & +{\displaystyle \frac{\partial H}{\partial t}}\left(H_{\text{pred}}^{t+\mathrm{\Delta }t},{\mathit{u}}^{t+\mathrm{\Delta }t}\right)].\end{array}$$

The discretisation error $\mathit{\tau }$ in the ice thickness is estimated based on the difference between the predicted and the corrected ice thicknesses: 32$${\mathit{\tau }}^{t+\mathrm{\Delta }t}={\displaystyle \frac{{\mathit{\zeta }}_t\left(H_{\text{corr}}^{t+\mathrm{\Delta }t}-H_{\text{pred}}^{t+\mathrm{\Delta }t}\right)}{(\mathrm{3}{\mathit{\zeta }}_t+\mathrm{3})\mathrm{\Delta }{t}^t}}.$$

The time step is then adapted based on the maximum discretisation error: 33$$\mathrm{\Delta }{t}^{t+\mathrm{\Delta }t}={\left({\displaystyle \frac{\mathit{ϵ}}{\text{max}\left|{\mathit{\tau }}^{t+\mathrm{\Delta }t}\right|}}\right)}^{(k_{\mathrm{I}}+k_{\mathrm{p}})}{\left({\displaystyle \frac{\mathit{ϵ}}{\text{max}\left|{\mathit{\tau }}^t\right|}}\right)}^{-k_{\mathrm{p}}}\mathrm{\Delta }{t}^t.$$

Here, $\mathit{ϵ}$ is the target truncation error in the ice thickness rate of change (configurable, default value 3 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$), and $k_{\mathrm{I}}=\mathrm{0.2}$ and $k_{\mathrm{p}}=\mathrm{0.1}$ are tuning parameters (values taken from Robinson et al., 2020).

It should be noted that Eqs. () and () involve different realisations of the ice thickness rate of change $\frac{\partial H}{\partial t}$. However, the equations for the different ice thickness schemes (Eqs. , , and ) yield $H^{t+\mathrm{\Delta }t}$ for a given value of $\mathrm{\Delta }t$. In the model, the current ice thickness $H^t$ is subtracted from that and the remainder divided by $\mathrm{\Delta }t$ to find $\frac{\partial H}{\partial t}$. Later on, $\frac{\partial H}{\partial t}$ is then adapted by the predictor–corrector scheme to yield the “final” value of $\frac{\partial H}{\partial t}$ that is used by the model.

3 Code

3.1 Parallelisation

A major change in v2.0 with respect to v1.0 is the switch from a shared-memory architecture, where all parts of the memory are accessible via a common bus to all computing cores, a distributed-memory architecture, which involves communication between memory-separated computing nodes. Memory access within shared-memory nodes outperforms message passing between separated-memory nodes, which implies that, all else being equal, v2.0 would be (slightly) slower than v1.0. However, the shared-memory architecture can only run on the number of cores within a single multi-core, shared-memory node (typically 32 or 64 on many high-performance scientific computing systems). The distributed-memory architecture is not limited in this way, allowing the user to scale up to far larger numbers of cores if necessary. With distributed-memory MPI, the code path and the communication paradigm stay the same whether running on a single-node or a multi-node configuration. However, inter-nodal communication is usually much slower than intra-nodal communication, which might cause an observable slowdown in the algorithm when moving to multiple nodes.

Solving the matrix equation representing the momentum balance is currently the most computationally demanding part of the model by far, often accounting for more than 80 % of the total computation time of a simulation (when using DIVA; the hybrid SIA–SSA and the BPA are more expensive to solve and would account for an even larger fraction). UFEMISM uses the PETSc library (Balay et al., 2021) for this. Most of the other operations that require data exchange between processes (e.g. remapping and calculating gradients) are represented by matrix operations, which are also handled by PETSc. In cases where the user requires a process to access data from another process, UFEMISM offers a set of standardised routines that, in turn, use the MPI API to facilitate this. For example, gathering a distributed array to a single process would require allocating memory on that process (possibly after performing an MPI reduction to determine how much memory is needed), calling one of the MPI_gather routines, and, finally (optionally), deallocating the memory for the distributed array. We combined these steps into a single subroutine to ease the workload of an aspiring UFEMISM developer. Currently, MPI_gather routines are only used for I/O and boundary communications when necessary for domain-wide computations that are not supported by PETSc.

We performed some simple simulations to assess the scalability of UFEMISM v2.0. These consist of the spin-up phase of the (modified, plan-view) Marine Ice Sheet Intercomparison Project (MISMIP) experiment, using DIVA with an 8 $\mathrm{km}$ resolution at the grounding line, for a period of 10 000 years. These simulations were run on the Snellius supercomputer on the AMD Rome 7H12 nodes (with 128 cores each). These preliminary results show that v2.0 does not yet scale well when using more than 32 cores, as shown in Fig. . Reducing the grounding-line resolution to 2 $\mathrm{km}$ to increase the size of the problem does not result in better scaling. This likely has to do with the way data communication between processes is handled by PETSc, which could be improved by paying more attention to the way the model domain is partitioned over the processes and the way PETSc decides which data should be communicated. These improvements are reserved for future work. Another contributing factor could be that the model setup used for the scaling test was too “small” (i.e. had too few vertices), causing the communication latencies between cores to dominate the total computation time. This is supported by the slowdown observed at 64 cores. Unfortunately, the time spent on communications is not (yet) measured separately in v2.0. However, it should be noted that v2.0 in its current form is already capable of performing multi-millennial simulations of the Antarctic ice sheet, solving DIVA with a grounding-line resolution of $<$ 5 $\mathrm{km}$ across selected basin-scale regions (e.g. the Amundsen Sea drainage basin), on a dual-core, consumer-grade laptop (Macbook Pro M2 2023), within 24 $\mathrm{h}$ of wall-clock time (Bernales et al., 2025). Large-scale practical applications of the model are therefore already feasible even without these future improvements.

Figure 2

Strong scaling for UFEMISM v2.0 with the 10 000-year initialisation phase of the (modified, plan-view) MISMIP experiment, run with DIVA (see also Sect. ). The domain consists of approximately 13 000 triangles. The full model is the sum of the ice dynamic and non-ice-dynamic components. I/O was disabled for these scaling tests. With more than 32 cores, a slowdown instead of a speedup is visible.

[Figure omitted. See PDF]

3.2 I/O

All output files of v2.0 are now NetCDF-4 standard (Unidata, 2023), and all input files are NetCDF too (replacing the small number of text-based files that v1.0 inputted and outputted). UFEMISM's NetCDF input routines automatically interact with the routines for remapping data between Cartesian grids (typical of ice-sheet-specific data, e.g. BedMachine; Morlighem et al., 2019), long–lat grids (typical of global climate model output), and triangular meshes (e.g. output from other UFEMISM simulations). The user can provide input data in any of those formats, and UFEMISM will automatically detect the type of grid (by parsing the names of the dimensions of the NetCDF file), choose the appropriate remapping function for that grid, and remap the data to the model mesh. The sparse matrices representing the remapping operators (commonly known as “`weights”; see Appendix ) are stored in memory so that if more data are read from the same input file later on, the matrix is reused instead of needing to calculate it again. All of this is done automatically, requiring no user intervention. Currently, second-order conservative remapping is used by default; with a single keyword, this can be changed to e.g. bilinear or nearest-neighbour interpolation. Input files that do not cover the entire computational domain are extrapolated on a nearest-neighbour basis; routines for applying a user-defined missing value or doing a linear or Gaussian extrapolation instead exist and can be easily integrated here. Projection parameters specified in the header of the NetCDF file are not read. UFEMISM assumes that input grids use the same projection as the model itself (i.e. the ISMIP standard projections for Greenland and Antarctica). Converting between different projections therefore must be done by the user before providing files to the model.

UFEMISM produces output on both the model mesh and a Cartesian grid (with a user-defined resolution). The former is useful for detailed post-processing or visualisation, while the latter can be conveniently used for cursory inspection of model output using any NetCDF viewing software, as well as for coupling to other models where square-grid input is more convenient. The user can specify which data fields should be included in the output files within the model configuration files; the full list of the 100$+$ fields (both two- and three-dimensional) that the user can choose from can be found in the NetCDF output module. Adding a new field requires about 10 lines of new code. The standard output includes all the data required for a “perfect restart” (so that e.g. running one simulation of 200 years yields identical results to two subsequent simulations of 100 years), which is necessary for script-based coupling to other models. Additionally, UFEMISM generates a separate NetCDF file with time series of domain-integrated quantities (e.g. mass balance components, ice volume).

3.3 Version control

UFEMISM is maintained on GitHub (https://github.com/IMAU-paleo/UFEMISM2.0, last access: 16 June 2025). GitHub Actions (https://docs.github.com/en/actions, last access: 16 June 2025) has been set up to automatically perform all the unit tests that have been built in for the routines interfacing with OpenMPI and PETSc: the NetCDF I/O routines, mesh generation, remapping, and PDE discretisation. This enables the user to quickly diagnose any problems occurring in the model. A number of benchmark experiments have been set up similarly, which are automatically run when Git branches are merged. Figures for these experiments, following the style of the publications where these benchmark experiments were first presented (e.g. Pattyn et al., 2008, for the ISMIP-HOM experiments) are created automatically. The UFEMISM GitHub repository also features integration with the nix package manager (https://nixos.org/, last access: 16 June 2025). This should allow the user to install all the required libraries (OpenMPI, PETSc, NetCDF) with their transient dependencies, using the exact version numbers for each of them, with a single command.

4 Idealised-geometry experiments

4.1 ISMIP-HOM

The Ice-Sheet Model Intercomparison Project for Higher-Order Models (ISMIP-HOM; Pattyn et al., 2008) contains several experiments to benchmark the velocities produced by the momentum balance in an idealised geometry. These experiments describe a slab of ice on a sloping bed. In experiments A and B, no sliding is allowed, and periodic undulations are superimposed on the flat bed slope, either in both the along-slope and the cross-slope directions (experiment A) or in only the along-slope direction (experiment B). In experiments B and C, the bedrock remains a flat slope, but sliding is now allowed, with the basal friction coefficient varying periodically in both the along-slope and the cross-slope directions (experiment C) or in only the along-slope direction (experiment D). Six different versions of each experiment exist, differing in the wavelength of the bedrock undulations or the friction variations, with values ranging between 160 and 5 $\mathrm{km}$. Decreasing the wavelength increases the aspect ratio of the ice geometry, making the more simplified momentum balance approximations, such as the SIA and SSA, less accurate. The experimental setup is described in full by Pattyn et al. (2008).

Figure 3

Ice surface velocities calculated by UFEMISM with the hybrid SIA–SSA, DIVA, and the BPA in the six different versions of ISMIP-HOM's experiment A (periodic bedrock undulations in both directions), compared to the model ensemble by Pattyn et al. (2008), which is divided into the full-Stokes model ensemble (green) and the higher-order model ensemble (blue).

[Figure omitted. See PDF]

Figure 4

Ice surface velocities calculated by UFEMISM with the hybrid SIA–SSA, DIVA, and the BPA in the six different versions of ISMIP-HOM's experiment C (flat, sloping bed; periodic variations in friction in both directions), compared to the model ensemble by Pattyn et al. (2008), which is divided into the full-Stokes model ensemble (green) and the higher-order model ensemble (blue).

[Figure omitted. See PDF]

The velocities calculated by UFEMISM v2.0 for ISMIP-HOM experiments A and C using the hybrid SIA–SSA, DIVA, and the BPA are compared to the ISMIP-HOM ensemble by Pattyn et al. (2008) in Figs.  and .

In experiment C (Fig. ), which concerns sliding over a bed with spatially varying roughness, all three approximations result in velocities that agree well with the ensemble, with only the BPA solution lying (slightly) outside the ensemble range, differing from the full Stokes solution by up to 13 %. In experiment A (Fig. ), which concerns viscous, non-sliding flow over an undulating bed, the hybrid SIA–SSA starts to diverge from the ensemble with the increasing aspect ratio of the geometry at spatial scales of about 80 $\mathrm{km}$. UFEMISM's solution to the BPA lies within the ensemble for all spatial scales. DIVA produces a relative velocity error of about 17 % in the $L$ $=$ 40 $\mathrm{km}$ experiment, which increases to 25 % and 40 % for the $L$ $=$ 20 $\mathrm{km}$ and $L$ $=$ 10 $\mathrm{km}$ experiments, respectively. The choice of what level of error is acceptable is, to some extent, subjective. Considering the inter-model spread in ensembles of realistic experiments (e.g. ISMIP6-Antarctica; Seroussi et al., 2020) and the fact that ISMIP-HOM's experiment A has rather extreme subglacial topography, we believe it is justified to use DIVA in settings where subglacial topographical features have a typical length scale larger than 20 $\mathrm{km}$. Of course, when it becomes computationally feasible to use the BPA in large-scale realistic experiments, this is to be preferred.

Figure 5

(a) Transects of the ice sheet at the end of each of the three phases for the 10 $\mathrm{km}$ simulation. (b) Grounding-line position over time, using different grounding-line resolutions.

[Figure omitted. See PDF]

4.2 MISMIP

To demonstrate the effectiveness of our sub-grid basal friction scaling scheme at resolving grounding-line migration, we performed an experiment along the lines of the Marine Ice Sheet Intercomparison Project (MISMIP; Pattyn et al., 2012). The original experiment describes a flowline over a simple linear slope, which is subjected to a spatially uniform positive mass balance. Rather than transforming this one-dimensional flowline into a two-dimensional flowband, we have opted to extrude the one-dimensional geometry radially to create a circular, cone-shaped island. This results in the formation of a circular, dome-shaped ice sheet, which flows radially outward, feeding into an ice shelf that extends outward to infinity. While this means the resulting grounding-line position no longer matches the (semi-)analytical solution provided by Pattyn et al. (2012), it offers the advantage of checking the full two-dimensional stress balance (instead of only the $x$ component). The experimental protocol consists of stepwise decreasing/increasing the flow parameter $A$ in Glen's flow law, resulting in an advance/retreat of the grounding line. After being spun up to a steady state, a single advance–retreat cycle should, physically, result in the same grounding-line position as before. Any remaining difference in position, i.e. grounding-line hysteresis where there should be none, must therefore be a numerical path dependency, which the original MISMIP study showed could be significant (up to several hundred kilometres) in models that do not pay special attention to the way the discontinuous friction at the grounding line is handled (Pattyn et al., 2012).

We performed simulations with grounding-line resolutions of 10, 8, 5, and 4 $\mathrm{km}$ using DIVA. We start with a 10 000-year spin-up phase, with a uniform flow factor of $A$ $=$ 10⁻¹⁶ ${\mathrm{Pa}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$. We then decrease the flow factor to $A$ $=$ 10⁻¹⁷ ${\mathrm{Pa}}^{-\mathrm{3}}{\mathrm{yr}}^{-\mathrm{1}}$ for a period of 10 000 years, resulting in an advance of the grounding line by about 200 $\mathrm{km}$. Finally, we revert the flow factor back to its original value, causing the grounding line to retreat again. While the original experiment involves several more decreases in the flow factor before moving on to the stepwise increases, only a single decreased/increased step is sufficient to assess the level of unwanted numerical grounding-line hysteresis, which is what we aim to investigate here. The results of this experiment are shown in Fig. . Figure a shows transects of the ice sheet at the end of each of the three phases (spin-up, advance, retreat) for the 10 $\mathrm{km}$ simulation, while Fig. b shows the position of the grounding line over time for all three resolutions. The difference in the grounding-line position between the end of the spin-up phase at 10 $\mathrm{kyr}$ and the end of the retreat phase at 30 $\mathrm{kyr}$ is smaller than twice the grounding-line resolution in all simulations. Note that all these simulations were performed with the dynamic adaptive mesh, while in v1.0, a mesh update would result in a small but noticeable “jump” in the grounding-line position (Berends et al., 2021, their Fig. 10b; note that the study used the hybrid SIA–SSA instead of DIVA, the flux condition scheme instead of the sub-grid friction scaling scheme, and much coarser resolutions of 64–16 $\mathrm{km}$). Some improvements to the remapping scheme in v2.0 (see Appendix ) have greatly reduced this problem. Lastly, the sub-grid friction scaling scheme in v2.0 results in a more symmetrical, circular grounding line (not shown) than the flux condition scheme in v1.0. A well-known (but, as far as we are aware, never published) issue with flux condition schemes in square-grid models is the “octagonal” grounding line that can sometimes appear (in square-grid models; on unstructured grids, the grid dependency is often less obvious); a similar undesirable dependency on the grid geometry could sometimes be seen in v1.0.

4.3 MISMIP$+$

The third Marine Ice Sheet Model Intercomparison Project (MISMIP$+$; Asay-Davis et al., 2016) investigates the retreat of an ice stream feeding into a buttressed shelf. In the steady state, the ice stream flows down an 80 $\mathrm{km}$ wide, $\sim$ 500 $\mathrm{km}$ long fjord. The grounding line rests on a retrograde slope, which is kept stable by the strongly buttressed ice shelf. In the experiment, the ice sheet starts from a steady state and is subjected to a strong sub-shelf melt forcing. The resulting loss of buttressing causes the grounding line to retreat by about 50 $\mathrm{km}$ over the course of the 100-year simulation. The experimental setup is described by Asay-Davis et al. (2016), while the results of the intercomparison are presented by Cornford et al. (2020). The resulting grounding-line retreat was found to vary by a factor of approximately 3 between different models. A large part of this spread was attributed to (small) differences in initial conditions, as well as the choice of sliding law (the experimental protocol allows one to choose between three different sliding laws).

Figure 6

(a) Mid-stream transects of the ice sheet at the beginning (solid lines) and the end (dashed lines) of the 100-year retreat simulation at different resolutions (see legend). (b) Mid-stream grounding-line position over time at different resolutions (see legend) compared to the Cornford et al. (2020) model ensemble (mean shown by solid black line, spread shown by grey-shaded area).

[Figure omitted. See PDF]

We have performed the MISMIP$+$ experiment ice1r (100 years of increased-melt forcing) with UFEMISM v2.0 using the Schoof sliding law (Eq. , chosen here because, out of the three options in the MISMIP$+$ protocol, we find it results in the best numerical stability) and DIVA at grounding-line resolutions ranging from 5 $\mathrm{km}$ to 500 $\mathrm{m}$. Glen's flow law parameter $A$ has been tuned separately for each simulation to achieve a stable mid-stream grounding-line position at $x$ $=$ 450 $\mathrm{km}$. The results of these simulations are compared to the model ensemble results by Cornford et al. (2020) in Fig. . The UFEMISM results lie well within the Cornford et al. (2020) ensemble range. Note that these simulations were all performed with the dynamic adaptive mesh. In the 500 $\mathrm{m}$ simulation, the mesh was updated about once every model year on average, at no significant computational expense (as the computation time is dominated by solving the momentum balance). The solution does not seem to converge to a unique value with increasing resolution, which might be explained by the fact that the flow factor $A$ is tuned for each experiment individually. When we repeat the simulations with the same flow factor for every resolution (not shown), the resulting grounding-line retreat curves in Fig. b are more parallel but at the cost of an increased spread in initial positions (though still mostly within $\pm$ 5 $\mathrm{km}$ of the 450 $\mathrm{km}$ target).

5 Discussion and conclusions

We presented version 2.0 of UFEMISM and verified it in a number of benchmark experiments with idealised geometries. We have shown that the model is able to solve different approximations to the Stokes equations and integrate the resulting thinning rates through time to model the evolving ice geometry on a dynamic adaptive mesh. The results lie within the published model ensembles for all these experiments. These verified model capabilities provide the groundwork for the realistic applications presented in part 2 of this work (Bernales et al., 2025).

The numerical stability of the model has been greatly improved. This includes the new time-stepping scheme and the switch from a simple successive over-relaxation scheme to PETSc for solving the matrix equations. While these changes have generally improved the computational performance of the model, a direct comparison between v1.0 and v2.0 is complicated by the changes that have been made to the model physics and discretisation, such as the un-simplification of the SSA, the change from a grounding-line flux condition to a sub-grid friction scaling scheme, and the change from defining velocities on the edges to the triangle centres. Comparing the performance is further complicated by the absence of several new features in v1.0 that are required for realistic simulations of the Greenland or Antarctic ice sheets. For example, v1.0 lacks the modules for inverting the basal friction and the sub-shelf melt so that it cannot start from the same steady state as v2.0. Initialising the model with a spin-up instead, using simple parameterisations for the basal friction and melt, would lead to a very different, generally smoother initial ice geometry, which would artificially increase the stability of the model and inflate its performance.

The ISMIP-HOM experiments presented here, as well as the work by Rückamp et al. (2022), demonstrate the importance of considering the model's approximation to the Stokes equations when moving to high resolutions. At the high resolutions that UFEMISM can now achieve, topographical features can be resolved that would invalidate the underlying assumptions of DIVA. However, solving the BPA can easily require 50 times more computation time than solving DIVA, which would be unfeasible for many practical applications. Improving the model's performance when using large numbers of cores, as mentioned before, could be a way to solve this problem. Another approach could be to reduce the size of the physical problem by moving to regional ice-sheet modelling, limiting the model domain to a single drainage basin. In preparation for such an approach, the code of UFEMISM's routines for solving the ice thickness equation has been written in such a way as to easily allow the user to define regions where the ice thickness should not change.

The infinite-slab approach used by UFEMISM to simplify the momentum balance at the calving front, while not expected to greatly affect the solution, is outdated and should be replaced by a proper stress boundary condition in future work.

The current version of the model does not yet scale well, which is a major remaining point of improvement. We suspect part of this problem lies with the way PETSc is implemented in UFEMISM and, consequently, the way it handles inter-process communication. Although the (simple) mesh partitioning scheme that was created for v 1.0 (Berends et al., 2021) generally results in good load balancing, we suspect that currently, a lot of computation time is wasted by PETSc determining what data it should communicate (i.e. figuring out the non-zero structures of the different sub-matrices), when this information can already be determined a priori from the mesh connectivity. However, even with this sub-optimal performance, the model is already capable of performing high-resolution ($<$ 5 $\mathrm{km}$), multi-millennial simulations of the Antarctic ice sheet (Bernales et al., 2025), within a few hours on a consumer-grade (dual-core) laptop (although moving to even higher resolutions would currently still require the user to wait several days for the simulation to complete). Improving this part of the model's performance should be the focus of future work.

Appendix A Mesh generation

UFEMISM uses an extended version of Ruppert's algorithm (Ruppert, 1995) to iteratively refine a simple initial mesh until it meets the requirements of the ice-sheet geometry. In Ruppert's original algorithm, the mesh is inspected to find “bad” triangles, which are triangles whose smallest internal angle lies below a certain threshold value (typically 25°). These triangles are then “split”, meaning that a new vertex is added at that triangle's circumcentre, and the Delaunay triangulation is updated to include the new vertex. In UFEMISM, Ruppert's algorithm is extended to additionally mark as bad any triangle whose longest leg exceeds the maximum resolution for the area of the domain where that triangle lies. For example, if the grounding line passes through a triangle whose longest leg exceeds the user-defined maximum grounding-line resolution, that triangle is marked as bad, even if it meets Ruppert's original smallest-angle criterion.

While the general functionality of the mesh generation code has not fundamentally changed from v1.0, the way meshes are refined is quite different now. In v2.0, the mesh generation code is provided with data fields of bedrock elevation and ice thickness, which can be defined either on a square grid or on a mesh. This geometry is then “reduced” to obtain a list of [$x,y$] points that together span the grounding line (and similarly for the calving front). This is illustrated in Fig. .

Figure A1

The grounding line of the ice sheet can be represented by a series of short line segments. This grounding line was created from the BedMachine Antarctica dataset (Morlighem et al., 2019) at 40 $\mathrm{km}$ resolution so that the individual segments are at most $\mathrm{40}\sqrt{\mathrm{2}}\mathrm{km}$ long.

[Figure omitted. See PDF]

This line is provided as input to the mesh generation code, which simply checks which triangles cross with any section of the line and splits them if necessary. An advantage of this approach is that the code paths for generating a mesh based on an ice-sheet geometry that is provided on a square grid (e.g. BedMachine; Morlighem et al., 2019) and for a geometry provided on a mesh (e.g. during a mesh update in a UFEMISM simulation) are identical from the point where these geometries are reduced to lines.

Figure A2

Each row shows how the mesh refinement algorithm refines an existing mesh (first column) with a refinement forcing (second column) to produce a new mesh (third column). (a) Starting with the five-vertex, four-triangle “dummy” mesh, the line refinement algorithm is provided with a series of short line segments spanning a simple circle. (b) The mesh is further refined (to an even higher resolution) over a series of short line segments spanning a semicircle. (c) The mesh is further refined over two points. (d) A dummy mesh is refined over a series of line segments spanning the Antarctic grounding line, yielding a mesh that would be more suitable for the ice-sheet model.

[Figure omitted. See PDF]

In addition to the line-based mesh refinement code, v2.0 also contains point-based and polygon-based refinement routines. The point-based routine can be used to obtain a high resolution at a certain location of interest, for example an ice-core site. The polygon-based routine can be used to increase the mesh resolution over a certain ice-sheet section, e.g. the Pine Island Glacier drainage basin. The point-based and line-based refinement are illustrated in Fig. .

Through the config file, the user can set separate maximum resolutions for the entire domain, for grounded ice, floating ice, (a band of specified width around) the grounding line, the calving front, the grounded ice margin, and the coastline.

Appendix B Discretisation

The discretisation scheme used in v1.0, described in Berends et al. (2021), which was based on neighbour functions, has been replaced by a least-squares-based scheme based on Syrakos et al. (2017). The advantage of this new scheme is that it can be easily extended to work on different Arakawa grids (a benefit, since v2.0 makes a lot more use of staggering than v1.0 did due to the change in the definition of the velocities from mesh edges to mesh triangles) and higher orders of accuracy and that it can be coded much more elegantly.

Let $f:R^{\mathrm{2}}\to R$ be a function defined on the model domain, and let $f_{\mathrm{a}}{f}_{\mathrm{b}}{f}_{\mathrm{c}}$ be its discretised approximations on the mesh vertices (equivalent to the Arakawa A grid), triangles (B grid), and edges (C grid), respectively. For convenience, the discretised approximations to the partial derivatives of $f$ on the different grids are written as $f_{x,a}=(\frac{\partial f}{\partial x}{)}_{\mathrm{a}}{f}_{yy,c}=(\frac{{\partial }^{\mathrm{2}}f}{\partial x^{\mathrm{2}}}{)}_{\mathrm{c}}$. These partial derivatives can be expressed as linear combinations of $f_{\mathrm{a}}{f}_{\mathrm{b}}{f}_{\mathrm{c}}$, e.g.

B1$$f_{x,a}={\mathbf{M}}_{x,a,a}{f}_{\mathrm{a}}.$$

Here, ${\mathbf{M}}_{x,a,a}$ is an nV-by-nV matrix (with nV being the number of vertices in the mesh). In the notation convention used here, $M$ has three subscript indices. The first indicates the operation represented by $M$, $x$ for $\frac{\partial }{\partial x}$, $yy$ for $\frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}$, and so on, and $m$ for mapping $f$ between the different Arakawa grids. The second and third indices represent the source and destination Arakawa grids, respectively. For example, $M_{m,a,b}$ maps a data field from the vertices to the triangles.

B1 First-order regular grid

Syrakos et al. (2017) describe a (weighted) least-squares scheme for discretising the gradient operator on an unstructured grid. Let $f_{\mathrm{a}}^i{f}_{x,a}^i{f}_{y,a}^i$ be the values of the function $f$ and its first partial derivatives on vertex $i$. The value $f_{\mathrm{a}}^j$ of $f$ on vertex $j$, which neighbours vertex $i$, can then be expressed as a Taylor expansion of $f$ around $i$:

B2$$f_{\mathrm{a}}^j=f_{\mathrm{a}}^i+\mathrm{\Delta }{x}_j{f}_{x,a}^i+\mathrm{\Delta }{y}_j{f}_{y,a}^i+O\left(\mathrm{\Delta }{x}_j^{\mathrm{2}},\mathrm{\Delta }{y}_j^{\mathrm{2}}\right).$$

Here, $\mathrm{\Delta }{x}_j,\mathrm{\Delta }{y}_j$ is the displacement between vertices $j$ and $i$, as illustrated in Fig. .

Figure B1

Illustration showing part of a mesh. Vertex $i$ is indicated by the green dot, while its six neighbours are shown in orange. Vertex $j$ is one of its neighbours, with the displacement $\mathrm{\Delta }{x}^j,\mathrm{\Delta }{y}^j$ shown by the red arrows.

[Figure omitted. See PDF]

If $i$ has $n$ neighbours, this results in the following system of $n$ linear equations (defining $\mathrm{\Delta }{f}_{\mathrm{a}}^j\equiv f_{\mathrm{a}}^j-f_{\mathrm{a}}^i$, dropping the truncation error $O(\mathrm{\Delta }{x}_j^{\mathrm{2}},\mathrm{\Delta }{y}_j^{\mathrm{2}})$, and introducing the vertex weights $w_j$ for the weighted least-squares approximation): B3$$\begin{array}{rl}& \underset{W}{\underbrace{\left[\begin{array}{cccc}{w}_{\mathrm{1}}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& w_{\mathrm{2}}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\dots }& w_n\end{array}\right]}}\underset{\mathrm{b}}{\underbrace{\left[\begin{array}{c}\mathrm{\Delta }{f}_{\mathrm{a}}^{\mathrm{1}}\\ \mathrm{\Delta }{f}_{\mathrm{a}}^{\mathrm{2}}\\ \mathrm{⋮}\\ \mathrm{\Delta }{f}_{\mathrm{a}}^n\end{array}\right]}}\\ & =\underset{W}{\underbrace{\left[\begin{array}{cccc}{w}_{\mathrm{1}}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& w_{\mathrm{2}}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\dots }& w_n\end{array}\right]}}\underset{A}{\underbrace{\left[\begin{array}{cc}\mathrm{\Delta }{x}_{\mathrm{1}}& \mathrm{\Delta }{y}_{\mathrm{1}}\\ \mathrm{\Delta }{x}_{\mathrm{2}}& \mathrm{\Delta }{y}_{\mathrm{2}}\\ \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{\Delta }{x}_n& \mathrm{\Delta }{y}_n\end{array}\right]}}\underset{z}{\underbrace{\left[\begin{array}{c}{f}_{x,a}^i\\ f_{y,a}^i\end{array}\right]}}.\end{array}$$

Using matrix notation, this equation reads $Wb=WAz$, which can be solved for $z$: B4$$z=(A^{\mathrm{T}}{W}^{\mathrm{T}}WA{)}^{-\mathrm{1}}{A}^{\mathrm{T}}{W}^{\mathrm{T}}Wb=Q{\mathit{\beta }}_{\mathrm{b}}.$$

Here, we have grouped the $A$ and $W$ terms into $Q=(A^{\mathrm{T}}{W}^{\mathrm{T}}WA{)}^{-\mathrm{1}}$ and ${\mathit{\beta }}_{\mathrm{b}}=A^{\mathrm{T}}{W}^{\mathrm{T}}Wb$. The symmetric 2-by-2 matrix $A^{\mathrm{T}}{W}^{\mathrm{T}}WA$, which needs to be inverted to find $Q$, is expressed as B5$$A^{\mathrm{T}}{W}^{\mathrm{T}}WA={\sum }_{c=\mathrm{1}}^n{w}_{\mathrm{c}}^{\mathrm{2}}.$$

Here, $c$ loops over all vertices that are connected to $i$ (the orange vertices in Fig. ). The second term, ${\mathit{\beta }}_{\mathrm{b}}$, is expressed as B6$${\mathit{\beta }}_{\mathrm{b}}={\sum }_{c=\mathrm{1}}^n{w}_{\mathrm{c}}^{\mathrm{2}}\left[\begin{array}{c}\mathrm{\Delta }{x}_{\mathrm{c}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\\ \mathrm{\Delta }{y}_{\mathrm{c}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\end{array}\right].$$

Once $Q$ has been calculated by inverting $A^{\mathrm{T}}{W}^{\mathrm{T}}WA$, the first partial derivative $f_{x,a}^i$ of $f$ on $i$ can be expressed as B7$$\begin{array}{rl}{f}_{x,a}^i& =Q(\mathrm{1},\mathrm{1}){\sum }_{c=\mathrm{1}}^n\left(w_{\mathrm{c}}^{\mathrm{2}}\mathrm{\Delta }{x}_{\mathrm{c}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\right)\\ & +Q(\mathrm{1},\mathrm{2}){\sum }_{c=\mathrm{1}}^n\left(w_{\mathrm{c}}^{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{c}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\right).\end{array}$$

Since we defined $\mathrm{\Delta }{f}_{\mathrm{a}}^j\equiv f_{\mathrm{a}}^j-f_{\mathrm{a}}^i$, this can be rewritten to read B8$$\begin{array}{rl}{f}_{x,a}^i& =-f_{\mathrm{a}}^i{\sum }_{c=\mathrm{1}}^n\left[w_{\mathrm{c}}^{\mathrm{2}}\left(Q\right(\mathrm{1},\mathrm{1})\mathrm{\Delta }{x}_{\mathrm{c}}+Q(\mathrm{1},\mathrm{2}\left)\mathrm{\Delta }{y}_{\mathrm{c}}\right)\right]\\ & +{\sum }_{c=\mathrm{1}}^n{f}_{\mathrm{a}}^c\left[w_{\mathrm{c}}^{\mathrm{2}}\left(Q(\mathrm{1},\mathrm{1})\mathrm{\Delta }{x}_{\mathrm{c}}+Q(\mathrm{1},\mathrm{2})\mathrm{\Delta }{y}_{\mathrm{c}}\right)\right].\end{array}$$

This means that the coefficients of the operator matrix ${\mathbf{M}}_{x,a,a}$ are given by B9$${\mathbf{M}}_{x,a,a}^{i,j}=\left\{\begin{array}{ll}\begin{array}{c}-{\sum }_{c=\mathrm{1}}^n\left[w_{\mathrm{c}}^{\mathrm{2}}\right(Q(\mathrm{1},\mathrm{1})\mathrm{\Delta }{x}_{\mathrm{c}}\\ +Q(\mathrm{1},\mathrm{2})\mathrm{\Delta }{y}_{\mathrm{c}}\left)\right]\end{array}& \text{if }i=j,\\ \begin{array}{c}{w}_j^{\mathrm{2}}\left(Q\right(\mathrm{1},\mathrm{1})\mathrm{\Delta }{x}_j\\ +Q(\mathrm{1},\mathrm{2})\mathrm{\Delta }{y}_j)\end{array}& \begin{array}{l}\text{if }j\text{ is}\\ \text{connected}\\ \text{to }i,\end{array}\\ \mathrm{0}& \text{otherwise.}\end{array}\right.$$

Similarly, the coefficients for $M_{y,a,a}$ are given by B10$$M_{y,a,a}^{i,j}=\left\{\begin{array}{ll}\begin{array}{c}-{\sum }_{c=\mathrm{1}}^n\left[w_{\mathrm{c}}^{\mathrm{2}}\right(Q(\mathrm{2},\mathrm{1})\mathrm{\Delta }{x}_{\mathrm{c}}\\ +Q(\mathrm{2},\mathrm{2})\mathrm{\Delta }{y}_{\mathrm{c}}\left)\right]\end{array}& \text{if }i=j,\\ \begin{array}{c}{w}_j^{\mathrm{2}}\left(Q\right(\mathrm{2},\mathrm{1})\mathrm{\Delta }{x}_j\\ +Q(\mathrm{2},\mathrm{2})\mathrm{\Delta }{y}_j)\end{array}& \begin{array}{l}\text{if }j\text{ is}\\ \text{connected}\\ \text{to }i,\end{array}\\ \mathrm{0}& \text{otherwise.}\end{array}\right.$$

The weights $w_j$ depend on the distance between $j$ and $i$: B11$$w_j={\displaystyle \frac{\mathrm{1}}{{\left|r_j-r_i\right|}^q}}.$$

Following Syrakos et al. (2017), we choose $q=\frac{\mathrm{3}}{\mathrm{2}}$.

B2 First-order staggered grid

The derivation in Sect.  holds for the case where both the function $f$ and its gradients $f_x{f}_y$ are defined on the same grid so that $f_i$ is known. However, if, for example, we want to calculate the first partial derivative of $f$ on the mesh triangles $f_{x,b}$ when $f$ itself is defined on the mesh vertices ($f_{\mathrm{a}}$), then this condition does not hold, and a slightly different derivation is needed.

Consider the Taylor series described by Eq. (). We once again write out the system of linear equations for $f$ on the collection of neighbouring points, but this time we do not introduce $\mathrm{\Delta }f$ so that we obtain the following expression:

B12$$\begin{array}{rl}& \underset{W}{\underbrace{\left[\begin{array}{cccc}{w}_{\mathrm{1}}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& w_{\mathrm{2}}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\dots }& w_n\end{array}\right]}}\underset{\mathrm{b}}{\underbrace{\left[\begin{array}{c}\mathrm{\Delta }{f}_{\mathrm{a}}^{\mathrm{1}}\\ \mathrm{\Delta }{f}_{\mathrm{a}}^{\mathrm{2}}\\ \mathrm{⋮}\\ \mathrm{\Delta }{f}_{\mathrm{a}}^n\end{array}\right]}}\\ & =\underset{W}{\underbrace{\left[\begin{array}{cccc}{w}_{\mathrm{1}}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& w_{\mathrm{2}}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\dots }& w_n\end{array}\right]}}\underset{A}{\underbrace{\left[\begin{array}{ccc}\mathrm{1}& \mathrm{\Delta }{x}_{\mathrm{1}}& \mathrm{\Delta }{y}_{\mathrm{1}}\\ \mathrm{1}& \mathrm{\Delta }{x}_{\mathrm{2}}& \mathrm{\Delta }{y}_{\mathrm{2}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{1}& \mathrm{\Delta }{x}_n& \mathrm{\Delta }{y}_n\end{array}\right]}}\\ & \times \underset{z}{\underbrace{\left[\begin{array}{c}{f}_{\mathrm{b}}^i\\ f_{x,b}^i\\ f_{y,b}^i\end{array}\right]}}.\end{array}$$

Following the same derivation as before, the symmetric 3-by-3 matrix $A^{\mathrm{T}}{W}^{\mathrm{T}}WA$ that needs to be inverted to find $Q$ is now given by B13$$A^{\mathrm{T}}{W}^{\mathrm{T}}WA={\sum }_{c=\mathrm{1}}^n{w}_{\mathrm{c}}^{\mathrm{2}}\left[\begin{array}{ccc}\mathrm{1}& \mathrm{\Delta }{x}_{\mathrm{c}}& \mathrm{\Delta }{y}_{\mathrm{c}}\\ \mathrm{\Delta }{x}_{\mathrm{c}}^{\mathrm{2}}& & \mathrm{\Delta }{y}_{\mathrm{c}}^{\mathrm{2}}\end{array}\right].$$

Similarly, ${\mathit{\beta }}_{\mathrm{b}}$ is now given by B14$${\mathit{\beta }}_{\mathrm{b}}={\sum }_{c=\mathrm{1}}^n{w}_{\mathrm{c}}^{\mathrm{2}}\left[\begin{array}{c}{f}_{\mathrm{a}}^c\\ \mathrm{\Delta }{x}_{\mathrm{c}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\\ \mathrm{\Delta }{y}_{\mathrm{c}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\end{array}\right].$$

This leads to the following expression for the coefficients of the matrices $M_{m,a,b}{M}_{x,a,b}{M}_{y,a,b}$: $$\begin{array}{}\text{B15}& {\displaystyle }{M}_{m,a,b}^{i,j}=\left\{\begin{array}{ll}\begin{array}{c}-{\sum }_{c=\mathrm{1}}^n\left[w_{\mathrm{c}}^{\mathrm{2}}\right(Q(\mathrm{1},\mathrm{1})\\ +Q(\mathrm{1},\mathrm{2})\mathrm{\Delta }{x}_{\mathrm{c}}\\ +Q(\mathrm{1},\mathrm{3})\mathrm{\Delta }{y}_{\mathrm{c}}\left)\right]\end{array}& \text{if }i=j,\\ \begin{array}{c}{w}_j^{\mathrm{2}}\left(Q\right(\mathrm{1},\mathrm{1})\\ +Q(\mathrm{1},\mathrm{2})\mathrm{\Delta }{x}_j\\ +Q(\mathrm{1},\mathrm{3})\mathrm{\Delta }{y}_j)\end{array}& \begin{array}{l}\text{if }j\text{ is}\\ \text{connected}\\ \text{to }i,\end{array}\\ \mathrm{0}& \text{otherwise,}\end{array}\right.\text{B16}& {\displaystyle }{M}_{x,a,b}^{i,j}=\left\{\begin{array}{ll}\begin{array}{c}-{\sum }_{c=\mathrm{1}}^n\left[w_{\mathrm{c}}^{\mathrm{2}}\right(Q(\mathrm{2},\mathrm{1})\\ +Q(\mathrm{2},\mathrm{2})\mathrm{\Delta }{x}_{\mathrm{c}}\\ +Q(\mathrm{2},\mathrm{3})\mathrm{\Delta }{y}_{\mathrm{c}}\left)\right]\end{array}& \text{if }i=j,\\ \begin{array}{c}{w}_j^{\mathrm{2}}\left(Q\right(\mathrm{2},\mathrm{1})\\ +Q(\mathrm{2},\mathrm{2})\mathrm{\Delta }{x}_j\\ +Q(\mathrm{2},\mathrm{3})\mathrm{\Delta }{y}_j)\end{array}& \begin{array}{l}\text{if }j\text{ is}\\ \text{connected}\\ \text{to }i,\end{array}\\ \mathrm{0}& \text{otherwise,}\end{array}\right.\end{array}$$ B17$$M_{y,a,b}^{i,j}=\left\{\begin{array}{ll}\begin{array}{c}-{\sum }_{c=\mathrm{1}}^n\left[w_{\mathrm{c}}^{\mathrm{2}}\right(Q(\mathrm{3},\mathrm{1})\\ +Q(\mathrm{3},\mathrm{2})\mathrm{\Delta }{x}_{\mathrm{c}}\\ +Q(\mathrm{3},\mathrm{3})\mathrm{\Delta }{y}_{\mathrm{c}}\left)\right]\end{array}& \text{if }i=j,\\ \begin{array}{c}{w}_j^{\mathrm{2}}\left(Q\right(\mathrm{3},\mathrm{1})\\ +Q(\mathrm{3},\mathrm{2})\mathrm{\Delta }{x}_j\\ +Q(\mathrm{3},\mathrm{3})\mathrm{\Delta }{y}_j)\end{array}& \begin{array}{l}\text{if }j\text{ is}\\ \text{connected}\\ \text{to }i,\end{array}\\ \mathrm{0}& \text{otherwise.}\end{array}\right.$$

B3 Second-order regular grid

Here, we extend the discretisation scheme by Syrakos et al. (2017) to include the second-order partial derivatives $f_{xx}{f}_{xy}{f}_{yy}$. First, we extend the Taylor expansion of $f$ around $i$ to include the following second-order terms:

B18$$\begin{array}{rl}{f}_{\mathrm{a}}^j& =f_{\mathrm{a}}^i+\mathrm{\Delta }{x}_j{f}_{x,a}^i+\mathrm{\Delta }{y}_j{f}_{y,a}^i+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }{x}_j^{\mathrm{2}}{f}_{xx,a}^i\\ & +\mathrm{\Delta }{x}_j\mathrm{\Delta }{y}_j{f}_{xy,a}^i+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }{y}_j^{\mathrm{2}}{f}_{yy,a}^i\\ & +O\left(\mathrm{\Delta }{x}_j^{\mathrm{3}},\mathrm{\Delta }{y}_j^{\mathrm{3}}\right).\end{array}$$

Writing out the system of linear equations for all neighbours of $i$ now yields the following expression: B19$$\begin{array}{rl}& \underset{W}{\underbrace{\left[\begin{array}{cccc}{w}_{\mathrm{1}}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& w_{\mathrm{2}}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\dots }& w_n\end{array}\right]}}\underset{\mathrm{b}}{\underbrace{\left[\begin{array}{c}\mathrm{\Delta }{f}_{\mathrm{a}}^{\mathrm{1}}\\ \mathrm{\Delta }{f}_{\mathrm{a}}^{\mathrm{2}}\\ \mathrm{⋮}\\ \mathrm{\Delta }{f}_{\mathrm{a}}^n\end{array}\right]}}\\ & =\underset{W}{\underbrace{\left[\begin{array}{cccc}{w}_{\mathrm{1}}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& w_{\mathrm{2}}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\dots }& w_n\end{array}\right]}}\\ & \times \underset{W}{\underbrace{\left[\begin{array}{ccccc}\mathrm{\Delta }{x}_{\mathrm{1}}& \mathrm{\Delta }{y}_{\mathrm{1}}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{x}_{\mathrm{1}}^{\mathrm{2}}& \mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{y}_{\mathrm{1}}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{1}}^{\mathrm{2}}\\ \mathrm{\Delta }{x}_{\mathrm{2}}& \mathrm{\Delta }{y}_{\mathrm{2}}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{x}_{\mathrm{2}}^{\mathrm{2}}& \mathrm{\Delta }{x}_{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{2}}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{2}}^{\mathrm{2}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{\Delta }{x}_n& \mathrm{\Delta }{y}_n& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{x}_n^{\mathrm{2}}& \mathrm{\Delta }{x}_n\mathrm{\Delta }{y}_n& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{y}_n^{\mathrm{2}}\end{array}\right]}}\\ & \times \underset{z}{\underbrace{\left[\begin{array}{c}{f}_{x,a}^i\\ f_{y,a}^i\\ f_{xx,a}^i\\ f_{xy,a}^i\\ f_{yy,a}^i\end{array}\right]}}.\end{array}$$

The symmetric 5-by-5 matrix $A^{\mathrm{T}}{W}^{\mathrm{T}}WA$ that needs to be inverted to find $Q$ is now given by B20$$\begin{array}{rl}& A^{\mathrm{T}}{W}^{\mathrm{T}}WA=\sum _{c=\mathrm{1}}^n{w}_{\mathrm{c}}^{\mathrm{2}}\\ & \left[\begin{array}{ccccc}\mathrm{\Delta }{x}_{\mathrm{c}}^{\mathrm{2}}& \mathrm{\Delta }{x}_{\mathrm{c}}\mathrm{\Delta }{y}_{\mathrm{c}}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{x}_{\mathrm{c}}^{\mathrm{3}}& \mathrm{\Delta }{x}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{c}}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{x}_{\mathrm{c}}\mathrm{\Delta }{y}_{\mathrm{c}}^{\mathrm{2}}\\ & \mathrm{\Delta }{y}_{\mathrm{c}}^{\mathrm{2}}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{x}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{c}}& \mathrm{\Delta }{x}_{\mathrm{c}}\mathrm{\Delta }{y}_{\mathrm{c}}^{\mathrm{2}}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{c}}^{\mathrm{3}}\\ & & \frac{\mathrm{1}}{\mathrm{4}}\mathrm{\Delta }{x}_{\mathrm{c}}^{\mathrm{4}}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{x}_{\mathrm{c}}^{\mathrm{3}}\mathrm{\Delta }{y}_{\mathrm{c}}& \frac{\mathrm{1}}{\mathrm{4}}\mathrm{\Delta }{x}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{c}}^{\mathrm{2}}\\ & & & \mathrm{\Delta }{x}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{c}}^{\mathrm{2}}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{x}_{\mathrm{c}}\mathrm{\Delta }{y}_{\mathrm{c}}^{\mathrm{3}}\\ & & & & \frac{\mathrm{1}}{\mathrm{4}}\mathrm{\Delta }{y}_{\mathrm{c}}^{\mathrm{4}}\end{array}\right].\end{array}$$

Similarly, ${\mathit{\beta }}_{\mathrm{b}}$ is now given by B21$${\mathit{\beta }}_{\mathrm{b}}={\sum }_{c=\mathrm{1}}^n{w}_{\mathrm{c}}^{\mathrm{2}}\left[\begin{array}{c}\mathrm{\Delta }{x}_{\mathrm{c}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\\ \mathrm{\Delta }{y}_{\mathrm{c}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\\ \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{x}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\\ \mathrm{\Delta }{x}_{\mathrm{c}}\mathrm{\Delta }{y}_{\mathrm{c}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\\ \frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\Delta }{f}_{\mathrm{a}}^c\end{array}\right].$$

Expressions for the coefficients of ${\mathbf{M}}_{x,a,a}{M}_{y,a,a}{M}_{xx,a,a}{M}_{xy,a,a}{M}_{yy,a,a}$ (which are now fourth-order accurate operators) can be derived similarly to before.

Appendix C Remapping

Because of the dynamic adaptive grid, data fields must often be remapped between square grids and (different) irregular triangular meshes. Extensive preliminary experiments have shown that only second-order conservative remapping results in accurate model results (e.g. ice thickness over time that matches the analytical solution in the Halfar dome experiment). Less accurate remapping schemes (nearest neighbour, bilinear, biquadratic, binning, Gaussian interpolation) all result in much more diffusion during each remapping operation and additionally violate conservation of mass and energy when remapping ice thickness and temperature, as these schemes are generally not conservative.

The mathematical theory behind conservative remapping is described by Jones (1999) and is relatively straightforward. However, Jones (1999) derived the equations in spherical coordinates, whereas UFEMISM uses Cartesian coordinates. Furthermore, UFEMISM uses a slightly different scheme, which conserves both local and global integrated values (the definition of “conservative” used by Jones), as well as extreme values (an important property, as we do not want to end up with negative ice thickness after remapping). We therefore provide a full derivation here.

C1 Theory

Let there exist two meshes that both cover the same domain $\mathrm{\Omega }$: a source mesh (indicated hereafter by the subscript $s$) and a destination mesh (subscript $d$). Suppose the source mesh is the one that existed before a mesh update and the destination mesh is the newly generated mesh. Let $f_{{\mathrm{s}}^a}$ be a discrete function defined on the vertices of the source mesh. The remapping problem then consists of finding a new discrete function $f_{{\mathrm{d}}^a}$, defined on the vertices of the destination mesh, such that $$\begin{array}{}\text{C1}& {\displaystyle }\int \underset{\mathrm{\Omega }}{\int }{f}_{{\mathrm{d}}^a}\mathrm{d}A=\int \underset{\mathrm{\Omega }}{\int }{f}_{{\mathrm{s}}^a}\mathrm{d}A,\text{C2}& {\displaystyle }\int \underset{A_{\mathrm{d}}^i}{\int }{f}_{{\mathrm{d}}^a}\mathrm{d}A=\int \underset{A_{\mathrm{d}}^i}{\int }{f}_{{\mathrm{s}}^a}\mathrm{d}A,\begin{array}{l}\text{where }{A}_{\mathrm{d}}^i\text{ values are the}\\ \text{Voronoi cells of the}\\ \text{vertices of the}\\ \text{destination mesh}\end{array}\text{C3}& {\displaystyle }\text{min}\left(f_{{\mathrm{s}}^a}\right)\le f_{{\mathrm{d}}^a}\le \text{max}\left(f_{{\mathrm{s}}^a}\right).\end{array}$$

Here, Eq. () implies conservation of the global integrated value, Eq. () implies conservation of local integrated values, and Eq. () implies conservation of extreme values.

Let $f(x,y)$ be a piecewise bilinear function, which is obtained from the discrete source function on the source mesh triangles $f_{{\mathrm{s}}^b}$ by bilinearly interpolating inside the triangles:

C4$$f(x,y)=f_{{\mathrm{s}}^b}+\left(x-x_{{\mathrm{s}}^b}\right){\left({\displaystyle \frac{\partial f}{\partial x}}\right)}_{{\mathrm{s}}^b}+\left(y-y_{{\mathrm{s}}^b}\right){\left({\displaystyle \frac{\partial f}{\partial y}}\right)}_{{\mathrm{s}}^b}.$$

Here, $x_{{\mathrm{s}}^b}{y}_{{\mathrm{s}}^b}$ are the coordinates of the geometric centre of source mesh triangle $s^b$. Note that $f_{{\mathrm{s}}^b}$ can be obtained from $f_{{\mathrm{s}}^a}$ using the operator matrices derived in Appendix : $$\begin{array}{}\text{C5}& {\displaystyle }& {\displaystyle }{f}_{{\mathrm{s}}^b}=M_{m,s^a,s^b}{f}_{{\mathrm{s}}^a},\text{C6}& {\displaystyle }& {\displaystyle }{\left({\displaystyle \frac{\partial f}{\partial x}}\right)}_{{\mathrm{s}}^b}=M_{x,s^a,s^b}{f}_{{\mathrm{s}}^a},\text{C7}& {\displaystyle }& {\displaystyle }{\left({\displaystyle \frac{\partial f}{\partial y}}\right)}_{{\mathrm{s}}^b}=M_{y,s^a,s^b}{f}_{{\mathrm{s}}^a}.\end{array}$$

The discrete function $f_{{\mathrm{d}}^a}$ on the vertices of the destination mesh is found by simply averaging $f(x,y)$ over the Voronoi cells $A_{{\mathrm{d}}^a}$ of the vertices of the destination mesh: C8$$f_{{\mathrm{d}}^a}={\displaystyle \frac{\mathrm{1}}{A_{{\mathrm{d}}^a}}}\int \underset{A_{{\mathrm{d}}^a}}{\int }f(x,y)\mathrm{d}A.$$

Note that, as Eq. () implies $\text{min}\left(f\right(x,y\left)\right)\le f_{{\mathrm{d}}^a}\le \text{max}\left(f\right(x,y\left)\right)$ and Eq. () implies $\text{min}\left(f_{{\mathrm{s}}^a}\right)\le f(x,y)\le \text{max}\left(f_{{\mathrm{s}}^a}\right)$, this implies $\text{min}\left(f_{{\mathrm{s}}^a}\right)\le f_{{\mathrm{d}}^a}\le \text{max}\left(f_{{\mathrm{s}}^a}\right)$, thus satisfying the conservation of extreme values required by Eq. (). Substituting Eq. () into Eq. () yields C9$$\begin{array}{rl}{f}_{{\mathrm{d}}_{\mathrm{a}}}& ={\displaystyle \frac{\mathrm{1}}{A_{{\mathrm{d}}^a}}}{\sum }_s^b[\int \underset{A_{s^b{d}^a}}{\int }(f_{{\mathrm{s}}^b}+\left(x-x_{{\mathrm{s}}^b}\right){\left({\displaystyle \frac{\partial f}{\partial x}}\right)}_{{\mathrm{s}}^b}\\ & +\left(y-y_{{\mathrm{s}}^b}\right){\left({\displaystyle \frac{\partial f}{\partial y}}\right)}_{{\mathrm{s}}^b}\left)\mathrm{d}A\right].\end{array}$$

Here, $A_{s^b{d}^a}$ indicates the area of overlap between the source mesh triangles $s^b$ and the destination mesh Voronoi cells $d^a$. This is illustrated in Fig. .

Figure C1

(a) The source mesh, with the triangle $j$ indicated. (b) The destination mesh, with the Voronoi cell of vertex $i$ indicated. (c) The area of overlap $A_{s^{bj}{d}^{ai}}$ between the source mesh triangle $j$ and the destination mesh vertex $i$ is indicated by the thick gold line. The perimeter of this area consists of sections of the perimeter of source mesh triangle $j$ and the Voronoi cell of destination mesh vertex $i$.

[Figure omitted. See PDF]

Equation () can be rearranged to read C10$$\begin{array}{rl}{f}_{{\mathrm{d}}_{\mathrm{a}}}& ={\displaystyle \frac{\mathrm{1}}{A_{{\mathrm{d}}^a}}}{\sum }_s^b[f_{{\mathrm{s}}^b}\int \underset{A_{s^b{d}^a}}{\int }\mathrm{d}A\\ & +{\left({\displaystyle \frac{\partial f}{\partial x}}\right)}_{{\mathrm{s}}^b}\int \underset{A_{s^b{d}^a}}{\int }\left(x-x_{{\mathrm{s}}^b}\right)\mathrm{d}A\\ & +{\left({\displaystyle \frac{\partial f}{\partial y}}\right)}_{{\mathrm{s}}^b}\int \underset{A_{s^b{d}^a}}{\int }\left(y-y_{{\mathrm{s}}^b}\right)\mathrm{d}A],\\ & {\displaystyle \frac{\mathrm{1}}{A_{{\mathrm{d}}^a}}}{\sum }_s^b\left[\right(f_{{\mathrm{s}}^b}-x_{{\mathrm{s}}^b}{\left({\displaystyle \frac{\partial f}{\partial x}}\right)}_{{\mathrm{s}}^b}\\ & -y_{{\mathrm{s}}^b}{\left({\displaystyle \frac{\partial f}{\partial y}}\right)}_{{\mathrm{s}}^b})\int \underset{A_{s^b{d}^a}}{\int }\mathrm{d}A\\ & +{\left({\displaystyle \frac{\partial f}{\partial x}}\right)}_{{\mathrm{s}}^b}\int \underset{A_{s^b{d}^a}}{\int }x\mathrm{d}A\\ & +{\left({\displaystyle \frac{\partial f}{\partial y}}\right)}_{{\mathrm{s}}^b}\int \underset{A_{s^b{d}^a}}{\int }y\mathrm{d}A].\end{array}$$

Since the area of overlap $A_{s^b{d}^a}$ between a triangle of the source mesh and a Voronoi cell of the destination mesh will generally be an irregularly shaped polygon, Eq. () is generally not easy to evaluate. However, the problem can be simplified by applying the divergence theorem, rewriting the three surface integrals in Eq. () into line integrals: $$\begin{array}{}\text{C11}& {\displaystyle }& {\displaystyle }\int \underset{A}{\int }\mathrm{d}A=\underset{\partial A}{\oint }x\mathrm{d}y,\text{C12}& {\displaystyle }& {\displaystyle }\int \underset{A}{\int }x\mathrm{d}A=-\underset{\partial A}{\oint }xy\mathrm{d}x,\text{C13}& {\displaystyle }& {\displaystyle }\int \underset{A}{\int }y\mathrm{d}A=\underset{\partial A}{\oint }xy\mathrm{d}y.\end{array}$$

Note that, as the perimeters of both the source mesh triangles and the destination mesh Voronoi cells are piecewise linear curves, the perimeter of the area of overlap $A_{s^b{d}^a}$ must therefore also be a piecewise linear curve. The expressions for the three line integrals along a straight line from $p=[x_p,y_p]$ to $q=\left[x_q,y_q\right]$ are given by $$\begin{array}{}\text{C14}& {\displaystyle }\underset{p}{\overset{q}{\int }}x\mathrm{d}y=x_p\mathrm{\Delta }y-y_p\mathrm{\Delta }x+{\displaystyle \frac{\mathrm{\Delta }}{x}}\mathrm{2}\mathrm{\Delta }y\left(y_q^{\mathrm{2}}-y_p^{\mathrm{2}}\right),\text{C15}& {\displaystyle }\begin{array}{rl}\underset{p}{\overset{q}{\int }}-xy\mathrm{d}x& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(x_p{\displaystyle \frac{\mathrm{\Delta }}{y}}\mathrm{\Delta }x-y_p\right)\left(x_q^{\mathrm{2}}-x_p^{\mathrm{2}}\right)\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\displaystyle \frac{\mathrm{\Delta }}{y}}\mathrm{\Delta }x\left(x_q^{\mathrm{3}}-x_p^{\mathrm{3}}\right),\end{array}\end{array}$$ C16$$\begin{array}{rl}\underset{p}{\overset{q}{\int }}xy\mathrm{d}y& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(x_p-y_p{\displaystyle \frac{\mathrm{\Delta }}{x}}\mathrm{\Delta }y\right)\left(y_q^{\mathrm{2}}-y_p^{\mathrm{2}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\displaystyle \frac{\mathrm{\Delta }}{x}}\mathrm{\Delta }y\left(y_q^{\mathrm{3}}-y_p^{\mathrm{3}}\right).\end{array}$$ Here, $\mathrm{\Delta }x=x_q-x_p,\mathrm{\Delta }y=y_q-q_p$. Substituting Eqs. ()–() into Eq. () yields C17$$\begin{array}{rl}{f}_{{\mathrm{d}}^a}=& {\displaystyle \frac{\mathrm{1}}{A_{{\mathrm{d}}^a}}}{\sum }_s^b[\left(f_{{\mathrm{s}}^b}-x_{{\mathrm{s}}^b}{\left({\displaystyle \frac{\partial f}{\partial x}}\right)}_{{\mathrm{s}}^b}-y_{{\mathrm{s}}^b}{\left({\displaystyle \frac{\partial f}{\partial y}}\right)}_{{\mathrm{s}}^b}\right)\\ & \times \underset{\partial A_{s^b{d}^a}}{\oint }x\mathrm{d}y-{\left({\displaystyle \frac{\partial f}{\partial x}}\right)}_{{\mathrm{s}}^b}\underset{\partial A_{s^b{d}^a}}{\oint }xy\mathrm{d}x+{\left({\displaystyle \frac{\partial f}{\partial y}}\right)}_{{\mathrm{s}}^b}\\ & \underset{\partial A_{s^b{d}^a}}{\oint }xy\mathrm{d}y].\end{array}$$

This implies that, in order to find the remapped value of $f$ on a destination vertex, we need to find all the source triangles overlapping with that vertex's Voronoi cell and calculate the three line integrals around the perimeter of the area of overlap between that source triangle and the destination Voronoi cell.

As can be seen from Eq. (), the remapped function $f_{{\mathrm{d}}^a}$ is a linear combination of the triangle source function values $f_{{\mathrm{s}}^b}$ and its gradients $\left(\frac{\partial f}{\partial x}{)}_{{\mathrm{s}}^b}\right(\frac{\partial f}{\partial y}{)}_{{\mathrm{s}}^b}$, which are in turn linear combinations of the vertex source function values $f_{{\mathrm{s}}^a}$. We can therefore rewrite Eq. () as a matrix equation. First, we define the three matrices ${\mathbf{B}}_{x\mathrm{d}y}$, ${\mathbf{B}}_{-xy\mathrm{d}x}$, and ${\mathbf{B}}_{xy\mathrm{d}y}$, which contain the line integrals around the areas of overlap between the source triangles $s^b$ and the destination Voronoi cells $d^a$: C18$${\mathbf{B}}_{x\mathrm{d}y}^{ij}=\underset{\partial A_{s^{bj}{d}^{ai}}}{\oint }x\mathrm{d}y,$$ $$\begin{array}{}\text{C19}& {\displaystyle }& {\displaystyle }{\mathbf{B}}_{-xy\mathrm{d}x}^{ij}=-\underset{\partial A_{s^{bj}{d}^{ai}}}{\oint }xy\mathrm{d}x,\text{C20}& {\displaystyle }& {\displaystyle }{\mathbf{B}}_{xy\mathrm{d}y}^{ij}=\underset{\partial A_{s^{bj}{d}^{ai}}}{\oint }xy\mathrm{d}y.\end{array}$$

Note that ${\mathbf{B}}_{x\mathrm{d}y}^{ij}$, ${\mathbf{B}}_{-xy\mathrm{d}x}^{ij}$, and ${\mathbf{B}}_{xy\mathrm{d}y}^{ij}$ are non-zero if and only if source triangle $j$ and destination Voronoi cell $i$ overlap.

These three matrices can be combined to yield the three remapping weight matrices ${\mathbf{W}}_{\mathrm{0}}$, ${\mathbf{W}}_{\mathrm{1},x}$, and ${\mathbf{W}}_{\mathrm{1},y}$: $$\begin{array}{}\text{C21}& {\displaystyle }& {\displaystyle }{\mathbf{W}}_{\mathrm{0}}^{ij}={\displaystyle \frac{{\mathbf{B}}_{x\mathrm{d}y}^{ij}}{A_{{\mathrm{d}}^{ai}}}},\text{C22}& {\displaystyle }& {\displaystyle }{\mathbf{W}}_{\mathrm{1},x}^{ij}={\displaystyle \frac{{\mathbf{B}}_{-xy\mathrm{d}x}^{ij}}{A_{{\mathrm{d}}^{ai}}}}-{\mathbf{W}}_{\mathrm{0}}^{ij}{x}_{s^{bj}},\text{C23}& {\displaystyle }& {\displaystyle }{\mathbf{W}}_{\mathrm{1},y}^{ij}={\displaystyle \frac{{\mathbf{B}}_{xy\mathrm{d}x}^{ij}}{A_{{\mathrm{d}}^{ai}}}}-{\mathbf{W}}_{\mathrm{0}}^{ij}{y}_{s^{bj}}.\end{array}$$

Substituting Eqs. ()–() into Eq. () yields C24$$f_{{\mathrm{d}}^a}={\mathbf{W}}_{\mathrm{0}}{f}_{{\mathrm{s}}^b}+{\mathbf{W}}_{\mathrm{1},x}{\left({\displaystyle \frac{\partial f}{\partial x}}\right)}_{{\mathrm{s}}^b}+{\mathbf{W}}_{\mathrm{1},y}{\left({\displaystyle \frac{\partial f}{\partial y}}\right)}_{{\mathrm{s}}^b}.$$

Substituting Eqs. ()–() into Eq. () yields C25$$\begin{array}{rl}{f}_{{\mathrm{d}}^a}& =\left({\mathbf{W}}_{\mathrm{0}}{M}_{m,s^a,s^b}+{\mathbf{W}}_{\mathrm{1},x}{M}_{x,s^a,s^b}+{\mathbf{W}}_{\mathrm{1},y}{M}_{y,s^a,s^b}\right)f_{{\mathrm{s}}^a}\\ & =M_{s^a,d^a}{f}_{{\mathrm{s}}^a}.\end{array}$$

Here, $M_{s^a,d^a}={\mathbf{W}}_{\mathrm{0}}{M}_{m,s^a,s^b}+{\mathbf{W}}_{\mathrm{1},x}{M}_{x,s^a,s^b}+{\mathbf{W}}_{\mathrm{1},y}{M}_{y,s^a,s^b}$ is an ${\text{nV}}_{\mathrm{d}}$-by-${\text{nV}}_s$ matrix that represents the second-order conservative remapping operation from the source mesh vertices to the destination mesh vertices.

C2 Implementation

In order to calculate the remapping matrix $M_{s^a,d^a}$, the three line integrals in Eqs. ()–() need to be calculated around the areas of overlap between all source mesh triangles and destination mesh Voronoi cells. While the line integrals themselves are simple enough (Eqs. –), determining which source triangles overlap with which destination Voronoi cells is not straightforward. Given the large numbers of vertices and triangles involved in high-resolution meshes (easily several tens of thousands of both), it is necessary to pay attention to computational efficiency.

The perimeter $\partial A_{s^{bj}{d}^{ai}}$ of the area of overlap $A_{s^{bj}{d}^{ai}}$ between source triangle $s^{bj}$ and destination Voronoi cell $d^{ai}$ consists of part of the perimeter $\partial A_{s^{bj}}$ of source triangle $s^{bj}$ and part of the perimeter $\partial A_{{\mathrm{d}}^{ai}}$ of destination Voronoi cell $d^{ai}$. This means that, in order to calculate the coefficients of the three matrices in Eqs. ()–(), it suffices to integrate once around every source triangle and around every destination Voronoi cell, carefully keeping track of the triangle or Voronoi cell of the opposite mesh with which it overlaps.

In UFEMISM, this is done using a collection of “line-tracing” subroutines. Given a line $[p,q]$, the model “traces” that line through a mesh and returns a list of all the Voronoi cells or triangles through which that line passes and the line integrals for all the individual line segments lying within them. Great care is taken to detect cases where the perimeters of source triangles and destination Voronoi cells coincide to prevent double-counting. By actively tracing the line and finding the index of the next triangle or Voronoi cell it crosses into using the connectivity lists of the triangle or cell it leaves, instead of performing a mesh-wide search operation every time, the computational expense is greatly reduced. Thus, calculating the remapping matrix only takes a fraction of the computation time required to create a new mesh.

Appendix D Terrain-following coordinate transformation

In order to solve the BPA, the heat equation, and conservation of mass, the vertical dimension must be discretised as well. In UFEMISM, this is done by introducing a terrain-following coordinate transformation: $$\begin{array}{}\text{D1a}& {\displaystyle }& {\displaystyle }\widehat{x}(x,y,z,t)=x,\text{D1b}& {\displaystyle }& {\displaystyle }\widehat{y}(x,y,z,t)=y,\text{D1c}& {\displaystyle }& {\displaystyle }\mathit{\zeta }(x,y,z,t)={\displaystyle \frac{s(x,y,t)-z}{H(x,y,t)}},\text{D1d}& {\displaystyle }& {\displaystyle }\widehat{t}(x,y,z,t)=t.\end{array}$$

Equation (c) implies that $\mathit{\zeta }=\mathrm{0}$ at the ice surface and $\mathit{\zeta }=\mathrm{1}$ at the ice base. Note that, in order to transform the heat equation, the time dimension is transformed as well. Applying this coordinate transformation results in the following expressions for the gradient operators: $$\begin{array}{}\text{D2a}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial x}}={\displaystyle \frac{\partial }{\partial \widehat{x}+\frac{\partial \mathit{\zeta }}{\partial x}\frac{\partial }{\partial \mathit{\zeta }}}},\text{D2b}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial y}}={\displaystyle \frac{\partial }{\partial \widehat{y}+\frac{\partial \mathit{\zeta }}{\partial y}\frac{\partial }{\partial \mathit{\zeta }}}},\text{D2c}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial z}}={\displaystyle \frac{\partial \mathit{\zeta }}{\partial z}}{\displaystyle \frac{\partial }{\partial \mathit{\zeta }}},\text{D2d}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}={\displaystyle \frac{\partial }{\partial \widehat{t}}}+{\displaystyle \frac{\partial \mathit{\zeta }}{\partial t}}{\displaystyle \frac{\partial }{\partial \mathit{\zeta }}}.\end{array}$$

Applying the chain rule to Eq. (c) yields the following expressions for the gradients of $\mathit{\zeta }$: $$\begin{array}{}\text{D3a}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\zeta }}{\partial x}}={\displaystyle \frac{\mathrm{1}}{H}}\left({\displaystyle \frac{\partial s}{\partial x}}-\mathit{\zeta }{\displaystyle \frac{\partial H}{\partial x}}\right),\text{D3b}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\zeta }}{\partial y}}={\displaystyle \frac{\mathrm{1}}{H}}\left({\displaystyle \frac{\partial s}{\partial y}}-\mathit{\zeta }{\displaystyle \frac{\partial H}{\partial y}}\right),\text{D3c}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\zeta }}{\partial z}}={\displaystyle \frac{-\mathrm{1}}{H}},\text{D3d}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\zeta }}{\partial t}}={\displaystyle \frac{\mathrm{1}}{H}}\left({\displaystyle \frac{\partial s}{\partial t}}-\mathit{\zeta }{\displaystyle \frac{\partial H}{\partial t}}\right).\end{array}$$

The gradient operators in Eq. ()a–d can be represented by matrices, as derived in Appendix , by multiplying their untransformed equivalents by the gradients of $\mathit{\zeta }$, e.g.

D4$$M_{x,a,b}=M_{\widehat{x},a,b}+D\left({\displaystyle \frac{\partial \mathit{\zeta }}{\partial x}}\right)M_{\mathit{\zeta },a,b}.$$

Here, $\mathbf{D}\left(f\right)$ represents a diagonal matrix with the elements of the vector $\mathit{f}$ on the diagonal; i.e. $D^{ij}={\partial }^{ij}{f}^i$ (with ${\partial }^{ij}$ being the Kronecker delta). By thus calculating the matrices for all the gradient operators, the stiffness matrix representing the momentum balance can be assembled.

The scaled vertical coordinate $\mathit{\zeta }$ is discretised using an irregular, log-linear grid: D5$${\mathit{\zeta }}^k=\mathrm{1}-{\displaystyle \frac{R^{\left(\frac{n-k}{n-\mathrm{1}}\right)}-\mathrm{1}}{R-\mathrm{1}}}k\in [\mathrm{1},n].$$

This implies that the ratio between the grid spacings at the ice surface and ice base is approximately equal to $R$, which is a configurable number with a default value of $R=\mathrm{10}$. This scheme results in improved accuracy of the solution near the ice base, where the strain rates (in the BPA) and the temperature gradients (in the heat equation) are the highest, without requiring additional vertical grid points. The number of vertical layers is configurable and is by default set to 12.

Appendix E Vertical ice velocities

Applying the terrain-following coordinate transformation from Appendix  to the expression for conservation of mass in Eq. () yields

E1$${\displaystyle \frac{\partial u}{\partial \widehat{x}}}+{\displaystyle \frac{\partial \mathit{\zeta }}{\partial x}}{\displaystyle \frac{\partial u}{\partial \mathit{\zeta }}}+{\displaystyle \frac{\partial v}{\partial \widehat{y}}}+{\displaystyle \frac{\partial \mathit{\zeta }}{\partial y}}{\displaystyle \frac{\partial v}{\partial \mathit{\zeta }}}+{\displaystyle \frac{\partial \mathit{\zeta }}{\partial z}}{\displaystyle \frac{\partial w}{\partial \mathit{\zeta }}}=\mathrm{0}.$$

The terms $\frac{\partial u}{\partial \widehat{x}}+\frac{\partial v}{\partial \widehat{y}}$ describe the divergence in the two-dimensional plane, in scaled coordinates: E2$${\displaystyle \frac{\partial u}{\partial \widehat{x}}}+{\displaystyle \frac{\partial v}{\partial \widehat{y}}}=\widehat{V}u.$$

Averaging this divergence over the Voronoi cell of a mesh vertex yields E3$$\stackrel{\mathrm{‾}}{\widehat{V}u}={\displaystyle \frac{\mathrm{1}}{A}}\int \underset{A}{\int }\left(\widehat{V}u\right)\mathrm{d}A.$$

By applying the divergence theorem, this integral can be transformed to a loop integral around the boundary of the Voronoi cell: E4$$\stackrel{\mathrm{‾}}{\widehat{V}u}={\displaystyle \frac{\mathrm{1}}{A}}\underset{\partial A}{\oint }(u\cdot \widehat{n})\mathrm{d}S.$$

Here, $\widehat{n}$ is the outward normal vector to the Voronoi cell boundary. Substituting this expression into Eq. () yields E5$${\displaystyle \frac{\partial w}{\partial \mathit{\zeta }}}={\displaystyle \frac{-\mathrm{1}}{\partial \mathit{\zeta }/\partial z}}\left[{\displaystyle \frac{\mathrm{1}}{A}}\underset{\partial A}{\oint }(u\cdot \widehat{n})\mathrm{d}S+{\displaystyle \frac{\partial \mathit{\zeta }}{\partial x}}{\displaystyle \frac{\partial u}{\partial \mathit{\zeta }}}+{\displaystyle \frac{\partial \mathit{\zeta }}{\partial y}}{\displaystyle \frac{\partial v}{\partial \mathit{\zeta }}}\right].$$

This expression can then be integrated over the transformed vertical dimension to find $w$: E6$$w\left(\mathit{\zeta }\right)=w(\mathit{\zeta }=\mathrm{1})-\underset{\mathrm{1}}{\overset{\mathit{\zeta }}{\int }}{\displaystyle \frac{\partial w}{\partial \mathit{\zeta }}}\mathrm{d}\mathit{\zeta }.$$

Note that the minus sign in Eq. () arises from the fact that $\mathit{\zeta }$ runs from 0 at the ice surface to 1 at the ice base, meaning that integrating upwards from the ice base means integrating in the negative $\mathit{\zeta }$ direction. The vertical velocity at the base is given by E7$$\begin{array}{rl}w(\mathit{\zeta }-\mathrm{1})& =w_{\mathrm{b}}=u_{\mathrm{b}}\left({\displaystyle \frac{\partial s}{\partial x}}-{\displaystyle \frac{\partial H}{\partial x}}\right)+v_{\mathrm{b}}\left({\displaystyle \frac{\partial s}{\partial y}}-{\displaystyle \frac{\partial H}{\partial y}}\right)\\ & +{\displaystyle \frac{\partial s}{\partial t}}-{\displaystyle \frac{\partial H}{\partial t}}.\end{array}$$

Appendix F Calculating the ice flux divergence operator

Conservation of ice mass for a shallow layer of ice in the two-dimensional plane is expressed mathematically as F1$${\displaystyle \frac{\partial H}{\partial t}}=-\mathrm{\nabla }\cdot \left(\mathbf{u}H\right)+m.$$

Here, $m$ is the net mass balance, including terms at the ice base, the ice surface, and the lateral boundaries. This equation is discretised spatially using the finite-volume scheme that gave UFEMISM its name. Averaging Eq. () over the Voronoi cell of vertex $i$ (the control volume of the finite-volume scheme) yields F2$${{\displaystyle \frac{\partial H}{\partial t}}}^i={\displaystyle \frac{-\mathrm{1}}{A^i}}\int \underset{A}{\overset{i}{\int }}\mathrm{\nabla }\cdot \left(\mathbf{u}H\right)\mathrm{d}A+m^i.$$

Using the divergence theorem, the double integral in Eq. () can be transformed: F3$${{\displaystyle \frac{\partial H}{\partial t}}}^i={\displaystyle \frac{-\mathrm{1}}{A^i}}\underset{\partial A^i}{\oint }\left(\mathbf{u}H\right)\cdot \widehat{n}\mathrm{d}S+m^i.$$

Here, $\widehat{n}$ is the outward unit normal to the boundary $\partial A^i$ of the Voronoi cell of vertex $i$. Let $(uH{)}^{ij}$ be average ice flux on the shared Voronoi cell boundary of vertices $i$ and $j$. Then the loop integral (Eq. ) can be transformed to the following sum: F4$${{\displaystyle \frac{\partial H}{\partial t}}}^i={\displaystyle \frac{-\mathrm{1}}{A^i}}\sum _{j=\mathrm{1}}^n\left[(\mathbf{u}H{)}^{ij}\cdot {\widehat{n}}^{ij}{L}^{ij}\right]+m^i.$$

Here, ${\widehat{n}}^{ij}$ is the unit normal vector pointing from vertex $i$ to vertex $j$, $L^{ij}$ is the length of their shared Voronoi cell boundary, and ${\sum }_{j=\mathrm{1}}^n$ sums over only the vertices $j$ that are connected to $i$. We then introduce an upwind scheme for the ice flux $\mathbf{u}H$: F5$$(\mathbf{u}H{)}^{ij}=\left\{\begin{array}{ll}{\mathbf{u}}^{ij}{\mathit{H}}^{\mathrm{i}}& \text{if, }{\mathbf{u}}^{ij}\cdot {\widehat{n}}^{ij}>\mathrm{0},\\ {\mathbf{u}}^{ij}{H}^j& \text{otherwise.}\end{array}\right.$$

This implies that, if the ice flows from vertex $i$ to vertex $j$, the ice thickness in vertex $i$ determines the flux and vice versa. This scheme offers better numerical stability than using the average ice thickness of $i$ and $j$ regardless of the flow direction.

Equations () and () imply that ${\frac{\partial H}{\partial t}}^i$ is a linear combination of the ice thicknesses ${\mathit{H}}^{\mathrm{i}}$. Equation () can therefore be represented by a matrix equation: F6$${{\displaystyle \frac{\partial H}{\partial t}}}^i=-{\mathbf{M}}_{\text{divQ}}{\mathit{H}}^{\mathrm{i}}+m^i.$$

Here, $M_{\mathrm{d}ivQ}$ is a matrix whose coefficients depend on the mesh geometry and the ice velocities, which can be multiplied by the ice thickness vector ${\mathit{H}}^{\mathrm{i}}$ to find the ice flux divergence $\mathrm{\nabla }\cdot \left(\mathbf{u}H\right)$. The coefficients of $M_{\mathrm{d}ivQ}$ are given by F7$${\mathbf{M}}_{\text{divQ}}^{ij}=\left\{\begin{array}{ll}\frac{\mathrm{1}}{A^i}{\sum }_{j=\mathrm{1}}^n\left[L^{ij}max\left({\mathbf{u}}^{ij}\cdot {\widehat{n}}^{ij},\mathrm{0}\right)\right]& \text{if }i=j,\\ \frac{L^{ij}}{A^i}min\left({\mathbf{u}}^{ij}\cdot {\widehat{n}}^{ij},\mathrm{0}\right)& \begin{array}{l}\text{if }i\text{ is}\\ \text{connected}\\ \text{to }j,\end{array}\\ \mathrm{0}& \text{otherwise.}\end{array}\right.$$

Code and data availability

The source code of UFEMISM v2.0, the scripts for compiling and running the model on a variety of computer systems, and the configuration files for all simulations presented here are freely available on GitHub: https://github.com/IMAU-paleo/UFEMISM2.0 (Berends et al., 2025). The exact version of the code that was used to produce the results presented here is archived at zenodo.org (10.5281/zenodo.8434469, Berends et al., 2025), though aspiring users are advised to check out the latest version from GitHub.

Author contributions

CJB wrote the new model code except for the distributed-memory code, which was developed by VA. CJB performed the experiments. VA and JAB set up the new version control system. CJB wrote the draft of the manuscript; all authors contributed to the final version.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We thank the reviewers, Thomas Zwinger and Thomas Albrecht, and the editor, Philippe Huybrechts, for their constructive criticism of the paper.

Financial support

Constantijn J. Berends 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. Jorge A. Bernales and Victor Azizi were supported by NWO under grant no. OCENW.KLEIN.515.

Review statement

This paper was edited by Philippe Huybrechts and reviewed by Thomas Zwinger and Torsten Albrecht.