Anisotropic metric-based mesh adaptation for ice flow modelling in Firedrake

Glaciological modelling is a computationally challenging task due to its high cost and complexity associated with large spatial-scale and long timescale simulations. In this paper we propose feature-based anisotropic mesh adaptation methods and demonstrate their effectiveness for time-dependent glaciological modelling on a Marine Ice Sheet Model Intercomparison Project (MISMIP+) experiment. Our methods use the Python-based Firedrake finite element library and the mmg2d remeshing software. We show that we are able to achieve solution accuracy comparable to uniform 0.5 km resolution mesh simulations using a sequence of adapted meshes with, on average, 30 times fewer vertices when adapting meshes based on the basal stress and 8–12 times fewer when adapting based on ice thickness and velocity. We further introduce a novel hybrid time-dependent fixed-point mesh adaptation algorithm that reaches mesh convergence approximately twice as fast compared with the existing global fixed-point algorithm. Since the fixed-point algorithms require that the problem is solved multiple times, the reported reduction in the number of vertices ultimately translates into a 3–6 times lower overall computational cost compared to uniform mesh simulations.

1 Introduction

Despite significant improvements in ice sheet models, the accurate prediction of the contribution of ice sheets to sea-level rise in the 21st century remains uncertain (Aschwanden et al., 2021). A major uncertainty pertains to correctly modelling the dynamics of the grounding line, which is where the ice transitions from being grounded to floating in the ocean (Durand and Pattyn, 2015). Early studies demonstrated that having appropriate refinement of the computational mesh near the grounding line is vital for obtaining reliable numerical results (Schoof, 2007; Vieli and Payne, 2005). This has been further confirmed by the Marine Ice Sheet Model Intercomparison Project (MISMIP+) experiments, highlighting the mesh resolution's importance regardless of the level of approximation in the force balance equations (Pattyn et al., 2012, 2013). While some studies have proposed methods to mitigate mesh dependency by introducing a smooth transition in the basal stress at the grounding line at the cost of model accuracy (Leguy et al., 2014), the majority have focused on reducing computational costs by numerical and meshing techniques while maintaining the model physics. Different strategies have been developed to achieve the latter, including models with explicit mesh movement to solve for the grounding line position (Vieli and Payne, 2005; Moreno-Parada et al., 2023) or fixed-grid models incorporating sub-grid parameterisation of the grounding line (Pattyn et al., 2006; Gladstone et al., 2012; Seroussi et al., 2014). However, these approaches rely on the hydrostatic criterion to define the grounding line position and thus cannot be readily applied to flow models that solve the complete force balance equations and treat the grounding line position as a contact problem (Gagliardini et al., 2016; de Diego et al., 2022).

On the other hand, mesh adaptation methods select finer mesh resolution in regions of interest (e.g. in the vicinity of the grounding line) and coarser elsewhere to preserve numerical resources. These methods have grown particularly popular in finite element and finite volume modelling due to their ability to utilise different types of meshes (structured, unstructured, or hybrid) and, through sophisticated error estimation and adaptation algorithms, to enhance computational efficiency and improve solution accuracy (Alauzet and Loseille, 2016). In the context of mesh adaptation, “solution accuracy” refers specifically to minimising discretisation error, i.e. the error that arises from approximating continuous mathematical equations by discrete numerical methods. Other sources of error are not addressed.

Different criteria have been used in the literature to define the regions of interest with desired mesh resolution. The distance to the grounding line, or alternatively distance from flotation, has been the most popular criterion to define areas warranting smaller mesh sizes (Durand et al., 2009; Goldberg et al., 2009; Gladstone et al., 2010; Cornford et al., 2013; Gudmundsson et al., 2012; Jouvet and Gräser, 2013; dos Santos et al., 2019). However, such criteria do not explicitly control the solution error and neglect the error contribution of domain regions away from the grounding line while refining other regions where this may not be necessary, hence making the simulation less efficient. Distance-from-flotation-based criteria have been shown to be inferior to refinement criteria based on an error estimator (dos Santos et al., 2019). Goldberg et al. (2009) use the jumps in strain rate at cell boundaries as a generic estimator of the numerical error to define areas that require finer mesh resolution. In their method, the total number of mesh nodes is constant and this generic estimator is not adapted to handle all flow regimes as it tends to increase the resolution at the shear margins, sometimes at the expense of the resolution at the grounding line. dos Santos et al. (2019) implement a true a posteriori error approximation (i.e. the actual numerical error is approximated after the solution has been obtained), the Zienkiewicz–Zhu error estimator, for the deviatoric stress tensor and ice thickness. They find that using the distance to grounding line and the error estimator criteria separately and in combination produces similar results on a benchmark steady-state problem but predict that the combined criterion would yield superior performance in simulations involving real ice sheets.

Existing ice sheet models utilising mesh adaptation also differ in the way they cope with grounding line movement in time-dependent problems. Some models use fixed meshes, which requires the a priori refinement of areas where the grounding line is susceptible to move through, while other models use mesh adaptation techniques that can involve moving meshes (Durand et al., 2009; Goldberg et al., 2009), mesh refinement and/or coarsening (Goldberg et al., 2009; Jouvet and Gräser, 2013; Cornford et al., 2013), or the entire remeshing of the domain (e.g. the Úa finite-element ice-flow model described in Gudmundsson et al., 2012). However, none of these studies consider the temporal distribution of spatial error associated with time-dependent mesh adaptations.

Mesh adaptation strategies have become popular in computational fluid dynamics to capture complex multiscale phenomena such as shock wave propagation (e.g. Frey and Alauzet, 2005; Alauzet et al., 2007). The general aim is to control the accuracy of the solution by adapting the size and shape (or anisotropy) of individual mesh elements. A core challenge is in finding efficient and reliable “estimators” of the numerical error which are used as the basis with which to define refinement criteria. Frey and Alauzet (2005) propose an estimator based on the interpolation error; the mesh size is then defined by an anisotropic metric map which equidistributes the error in the computational mesh. Such approaches are also known as feature-based approaches, since the constructed metric aims to capture the features of solution fields, or derived solution-dependent fields, from which meshes are to be adapted. While not based on a true a posteriori approximation error, it has been found effective in practice to control the numerical error and offers flexibility in combining metrics obtained for different variables. The generation of anisotropic meshes allows us to capture strongly directional processes and geometries even more efficiently. This is desirable in glacier shear margins, for example, where much finer resolution is required across the shear margin than along it. Such methods have been applied in finite element ice flow models to capture the dynamics of fast-flowing outlet glaciers using a metric defined from the observed surface speed while keeping overall computational costs low by coarsening other regions (Morlighem et al., 2010; Seddik et al., 2012; Gillet-Chaulet et al., 2012). These methods generate only a single mesh that is used throughout the simulation. However, they do not consider the non-linear and time-dependent coupling between the solution and the underlying mesh. This coupling suggests an iterative mesh adaptation procedure, which we will implement.

The purpose of this paper is twofold. Firstly, we aim to provide a self-contained description of anisotropic metric-based mesh adaptation methods suitable for ice sheet and glacier modelling. Secondly, we demonstrate and evaluate their ability to control solution accuracy while maintaining low computational cost in the context of grounding line dynamics modelling. We build on earlier applications of the method in glaciological modelling, and we implement a novel adaptation procedure appropriate for transient simulations which consists of a spatio-temporal error analysis and a generation of multiple meshes that control the error.

2 MISMIP+ experiment set-up

The numerical experiment from the third Marine Ice Sheet Model Intercomparison Project (MISMIP+; Asay-Davis et al., 2016; Cornford et al., 2020) features an idealised ice stream, which is an elongated region of fast-flowing ice within an ice sheet. While the ice sheet exhibits negligible basal sliding, the rapid ice stream motion is dominated by processes at the ice–bed interface. Viscous stresses can also be significant in some parts of the ice stream, such as near the grounding line, in shear margins, and where basal traction is low (Greve and Blatter, 2009; Stokes, 2018). The ice experiences a sudden change in flow regime at the grounding line, where the ice stream flows into the floating ice shelf, which is no longer in contact with the bed topography.

We follow the experiment design and prescribed parameter values as described in Asay-Davis et al. (2016), which places the ice stream in an elongated rectangular domain measuring 640 km in the x direction and 80 km in the y direction. The ice is flowing approximately parallel to the x axis over an idealised bedrock topography. The topography is described by a sixth-order polynomial in the x direction and an exponential in the y direction: prescribing an elongated central trough surrounded by steep walls. The bed topography slopes downwards throughout most of the domain but involves a retrograde slope at around x=450 km where the steady-state ice stream would ground. Beyond the grounding line is the ice shelf that is fed by the upstream flow of ice. The experiment prescribes a no-slip boundary condition at x=0, free-slip conditions at $y=-\mathrm{40}$ and y=40 km, and a fixed calving ice front at x=640 km where ice is removed from the domain. The prescribed topography and boundary conditions lead to solutions that are symmetrical about the middle of the domain, at y=0 km. Therefore, we choose to run the experiment in only half of the domain to preserve computational resources, as several participants in the intercomparison have also done (Cornford et al., 2020).

We focus on the Ice1 group of experiments from the MISMIP+ exercise as it produces more drastic glacier evolution and grounding line migration compared to other experiments in the exercise (Cornford et al., 2020). The Ice1 experiment runs for 200 years, where the first 100 years see a retreat of the glacier induced by ice shelf melting. The ice shelf melting is then removed in the last 100 years, when the glacier grows and re-advances. The drastic change in flow regime accompanied by the rapid migration of the grounding line in the Ice1 experiment presents an ideal test case for the application of mesh adaptation methods.

2.1 Solving equations of glacier flow

The description of the ice stream dynamics requires a momentum conservation equation, which describes how the velocity field u(x,t) evolves under the influence of forces. In this work we apply the shallow stream approximation, which yields the following depth-averaged momentum conservation equation:

$$\begin{array}{}\text{(1)}& \mathrm{\nabla }\cdot \left(H\mathbf{S}\right)+{\mathit{\tau }}_{\mathrm{b}}={\mathit{\rho }}_{\mathrm{I}}gH\mathrm{\nabla }s,\end{array}$$

where $H=H(\mathit{x},t)$ is the ice thickness, $\mathbf{S}=\mathbf{S}(\mathit{x},t)$ is the membrane stress tensor, ${\mathit{\tau }}_{\mathrm{b}}={\mathit{\tau }}_{\mathrm{b}}(\mathit{x},t)$ is basal friction, ρ_I is ice density, g is gravitational acceleration, and $s=s(\mathit{x},t)$ is surface elevation. The membrane stress and strain rates are related by Glen's flow power law. Furthermore, we require a description of how ice thickness H evolves in time, which is described by the following depth-averaged mass conservation equation:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial H}{\partial t}}+\mathrm{\nabla }\cdot \left(H\mathit{u}\right)=\dot{b},\end{array}$$

where $\dot{b}=\dot{b}(\mathit{x},t)$ is the climatic-basal mass balance rate. For a detailed discussion and derivation of the equations, we refer the reader to Greve and Blatter (2009).

To solve equations of glacier flow we use icepack, a Python library built on Firedrake, which includes relevant highly customisable glacier flow models (Shapero et al., 2021). Like most current ice sheet models, icepack implements the first-order explicit Euler approximation to the coupled model to solve equations of glacier flow. That is to say, at each time step the scheme uses ice geometry at the previous time step to solve the stress balance equations for ice velocity at the current time step, which is then used to solve the mass transport equation to evolve ice geometry in time (for details, see Shapero et al., 2021).

Icepack formulates the shallow stream model from a principle of least action, where the action functional consists of terms for viscosity, friction, gravity, and terminus. This is a generalisation of the shallow shelf approximation of the Stokes equations as it involves a bed friction term τ_b which is non-zero for grounded ice and zero for floating ice. Therefore, in order to produce a smooth transition and prevent shocks across the grounding line, we gradually reduce friction in its vicinity (Leguy et al., 2014). To do that, we use a modified form of the Schoof sliding law derived by Shapero et al. (2021). Basal shear stress, τ_b, is then given by

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{\mathrm{b}}=\left\{\begin{array}{ll}-\frac{{\mathit{\alpha }}^{\mathrm{2}}N|\mathit{u}{|}^{\mathrm{1}/\mathrm{3}}}{{\left[{\left({\mathit{\alpha }}^{\mathrm{2}}{\mathit{\beta }}^{-\mathrm{2}}N\right)}^{\mathrm{4}}+|\mathit{u}{|}^{\mathrm{4}/\mathrm{3}}\right]}^{\mathrm{1}/\mathrm{4}}}\frac{\mathit{u}}{\left|\mathit{u}\right|}& \mathrm{if}N>\mathrm{0},\\ \mathrm{0}& \mathrm{if}N\le \mathrm{0},\end{array}\right.\end{array}$$

where N is the effective pressure at the ice base, u is the horizontal ice velocity, α²=0.5, and β²=100 MPa m${}^{-\mathrm{1}/\mathrm{3}}$ a${}^{\mathrm{1}/\mathrm{3}}$ is the friction coefficient. The effective pressure N is defined as the difference in ice and water overburden pressure; thus τ_b is zero for floating ice (N≤0). Finally, icepack implements a Lax–Wendroff scheme for the mass transport model (Shapero et al., 2021). The computational cost for solving the stress balance equations is significantly higher than the mass transport equation.

The weak forms of the shallow stream and mass conservation equations are then discretised by Firedrake using finite element methods. To construct an efficient solver for the stress balance equation, we leverage their formulation as the derivative of a convex action functional, meaning that the Hessian is symmetric and positive-definite (Shapero et al., 2021). We therefore solve linear systems using a direct solver based on Cholesky factorisation (Amestoy et al., 2001, 2006). Backtracking line search is used to solve non-linear systems.

2.2 Uniform mesh refinement

In order to study the sensitivity of the solution to mesh resolution and validate the mesh adaptation results, we first run the MISMIP+ Ice1 experiment on a series of uniform structured meshes, with mesh resolutions of 4, 2, 1, 0.5, 0.25, and 0.125 km.

The Ice1 experiment starts from an initial state computed by running the model forward in time for a long time with constant climatic-basal mass balance rate $\dot{b}=\mathrm{0.3}$ m a⁻¹ until the sought quasi-steady state is reached. Following Shapero et al. (2021), we efficiently spin up the model by adopting a hierarchical uniform mesh refinement strategy. We begin with a structured uniform coarse mesh of 4 km resolution on which we compute a reasonable approximation of the solution. The mesh is then uniformly refined, and the simulation continues. This process is repeated five times, resulting in uniform meshes of 4, 2, 1, 0.5, 0.25, and ultimately 0.125 km step sizes in both x and y directions. The total spin-up time is 15 000 years. With each mesh refinement, the difference in ice volume and grounding line position at steady state diminishes. We conclude that the simulation has converged at a uniform mesh resolution of 0.25 km due to very small differences between the 0.5, 0.25, and 0.125 km resolution results at steady state. This matches convergence studies performed by MISMIP+ participants, who concluded that the 0.5 km mesh resolution is adequate for the experiments (Cornford et al., 2020, supplementary data). The initial quasi-steady-state ice thickness H, ice velocity u, and basal shear stress τ_b are shown in Fig. 1, alongside bed topography contours.

Figure 1 Initial steady-state (a) ice thickness h, (b) ice velocity u, and (c) basal shear stress τ_b of the MISMIP+ experiment with the sliding law of Eq. (3).

Once the steady-state solution had been computed, we initialise the Ice1 experiment by interpolating the steady-state ice thickness and velocity fields onto uniform-resolution meshes. In order to better isolate the differences in results due to different spatial discretisations, for each simulation we use the same time step size $\mathrm{\Delta }t=\mathrm{1}/\mathrm{24}$ a which satisfies the Courant–Friedrichs–Lewy (CFL) condition on the finest mesh. Results of the Ice1 experiment on each of the uniform meshes are shown in Fig. 2. Finer-resolution simulations exhibit faster ice volume loss and retreat of the grounding line, while the coarsest resolution (4 km) finds a new quasi-steady state. There is again a small difference between the 0.125, 0.25, and 0.5 km results in most of the temporal evolution. The largest differences appear in the position of the grounding line along the shear margin (the lateral boundary of the ice stream), where the grounding line position on a 0.5 km mesh fluctuates, and towards the end of the simulation in the grounding line position at the ice ridge near the domain boundary (Fig. 2d). Numbers of degrees of freedom and computation times associated with each simulation are given in Table 1, which ranges from 2.7 min for the coarsest-resolution simulation to 83.2 h for the finest-resolution simulation. All simulations were run in serial for the purpose of CPU time measurements, but they can be run in parallel.

Figure 2 Results of the uniform mesh convergence study for the Ice1 experiment for different mesh resolutions: (a) volume above flotation (V_f) over time, (b) grounding line position along the midchannel line (x_gl(y=0)) over time, and (c) grounding line position in the domain at t=100 a and (d) t=200 a.

Table 1 Number of vertices (N_v) and CPU time associated with Ice1 experiment simulations shown in Fig. 2 run on a single CPU on uniform meshes of varying resolution.

2.3 Mesh adaptation considerations

The convergence study of Fig. 2 confirmed that the MISMIP+ experiment is indeed highly sensitive to mesh resolution. The study, however, does not inform us of what parts of the domain are most sensitive to resolution, since we uniformly refined the mesh over the entire domain. If a certain problem is equally sensitive to resolution throughout the domain, then using a uniform mesh may be the most suitable choice. However, as discussed in detail in the Introduction, earlier studies do show that ice sheet simulations benefit from selective refinement, particularly near the grounding line (Schoof, 2007; Vieli and Payne, 2005). As general rules of thumb, we would expect that finer resolution is required in regions where the solution changes rapidly, near moving interfaces or free boundaries, around complex geometric features, and so on.

As the name suggests, feature-based mesh adaptation relies on identifying instructive problem-specific features, or sensors, to guide mesh adaptation. The choice of a sensor field is non-trivial, as any spatially varying scalar field qualifies. This includes the solutions themselves, i.e. the scalar-valued ice thickness and individual components of the vector-valued ice velocity, as well as any solution-dependent field, such as the components of the vector-valued basal shear stress. Multiple sensor fields may also be combined, both in space and in time.

In the next section we show how a metric is defined from a sensor field and how different choices of sensor fields lead to different adapted meshes. For a simpler example, we refer the reader to Appendix A, where we apply mesh adaptation methods to a familiar problem of a Poisson equation.

3 Anisotropic metric-based mesh adaptation

Consider a partial differential equation (PDE) problem written in the generic operator form F(u)=0 on a bounded polygonal domain Ω with an exact solution u∈V. Since the domain is bounded, we can prescribe a spatial discretisation in terms of a mesh, ℋ, which consists of a finite number of non-overlapping elements and a corresponding number of vertices N_v. The numerical solution u_h∈V_h, where V_h⊂V is a discrete subspace of V, is obtained here via finite element methods, with an associated approximation error u−u_h. The cost of computing the solution u_h on ℋ increases as N_v increases. However, while increasing N_v through the uniform refinement of a mesh would generally be expected to result in reduced errors, uniform refinement is expected to give a sub-optimal reduction in error as N_v increases. Mesh adaptivity seeks to address this: to more optimally and robustly link a decreasing error with an increasing N_v.

The general mesh adaptation problem considered in this paper can be formulated in an a priori way as follows:

$$\begin{array}{}\text{(4)}& \mathrm{Find}{\mathcal{H}}_{\mathrm{opt}}\mathrm{with}{N}_{\mathrm{v}}\mathrm{vertices}\mathrm{solving}\underset{\mathcal{H}}{min}\mathbf{E}\left(\mathcal{H}\right),\end{array}$$

where ℋ_opt is the optimal mesh that minimises some measure of the approximation error E(ℋ). The goal of mesh adaptation is therefore to achieve a minimal approximation error for a given computational cost, which can be directly related to the number of vertices, N_v, of ℋ. In the anisotropic mesh adaptation approaches presented here, the optimal mesh is found approximately and iteratively by changing the size, shape, and orientation of its elements.

3.1 Continuous metric framework: 2D overview

In order to formulate the mesh adaptation problem (Eq. 4) such that it is well-posed and that it can be analysed using familiar methods from variational calculus, Loseille and Alauzet (2011a) propose the continuous metric framework. The framework establishes a strong duality between discrete mesh elements and continuous mathematical objects stemming from Riemannian geometry. The key idea in Riemannian geometry is to define a local way of measuring distances in a space, which, unlike in Euclidean geometry, may vary between points x∈Ω. The notion of measuring distances arises from the definition of a Riemannian metric M(x), which is identified with a positive-definite matrix in ℝ^d×d at each point x of the d-dimensional domain. A spatially varying metric leads to local notions of geometric quantities such as length, volume, and angle. This is particularly useful in the context of anisotropic mesh adaptation, as it allows us to prescribe different sizes, orientations, and directional stretchings of individual mesh elements in the domain. We refer the reader to Loseille and Alauzet (2011a, b) for a detailed description and analysis of the framework. However, in order to be self-contained, we illustrate how Riemannian metrics can be used to drive mesh adaptation on a two-dimensional domain Ω⊂ℝ².

It follows from the spectral theorem that M admits a diagonalisation $\mathbf{M}=\mathbf{R}\mathrm{diag}({\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}}){\mathbf{R}}^{\mathsf{T}},$ where R=R(x) is a matrix of orthogonal and normalised eigenvectors of M, and λ₁ and λ₂ are the corresponding eigenvalues such that λ₁≥λ₂. By defining metric density ρ and anisotropy quotients r₁ and r₂ as

$$\begin{array}{}\text{(5)}& \mathit{\rho }=\sqrt{\det \mathbf{M}}=\sqrt{{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}},r_{\mathrm{1}}=\sqrt{{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}}}=r_{\mathrm{2}}^{-\mathrm{1}},\end{array}$$

where anisotropy quotients measure the “stretching” strength in the directions of the eigenvectors in R, the diagonal decomposition can be written in a form that is more instructive in the context of mesh adaptation:

$$\begin{array}{}\text{(6)}& \mathbf{M}=\mathit{\rho }\mathbf{R}\left(\begin{array}{cc}{r}_{\mathrm{1}}^{-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& r_{\mathrm{2}}^{-\mathrm{1}}\end{array}\right){\mathbf{R}}^{\mathsf{T}}.\end{array}$$

In the continuous metric framework, the metric M prescribes the optimal size of the mesh element (through density ρ), its optimal shape (through quotients r₁, r₂), and its orientation (through eigenvectors in R). Thus, the metric-based approach introduces anisotropy in the resulting mesh.

Analogously to how the metric tensor M models a single mesh element, the metric tensor field $\mathcal{M}:\mathit{x}\in \mathrm{\Omega }\mapsto \mathbf{M}\left(\mathit{x}\right)$ models the entire mesh ℋ. To that end, the metric complexity

$$\begin{array}{}\text{(7)}& \mathcal{C}\left(\mathcal{M}\right)=\underset{\mathrm{\Omega }}{\int }\mathit{\rho }\left(\mathit{x}\right)\mathrm{d}\mathit{x}\end{array}$$

represents the integrated amount of resolution of ℳ; i.e. it provides an estimate for the number of vertices, N_v, in the mesh corresponding to that metric.

As introduced in George et al. (1991), the key idea of anisotropic metric-based mesh adaptation is to generate a mesh which has approximately unit-size element edge lengths when viewed in the prescribed Riemannian space. When viewed in the Euclidean space, the mesh is then appropriately adapted and anisotropic. Note that the unity criterion is generally impossible to satisfy exactly (consider discretising a rectangular domain with non-overlapping equilateral triangles in Euclidean space), so the criterion must be relaxed (for details, see Loseille and Alauzet, 2011a).

Following from the diagonalisation of M, the metric is commonly interpreted geometrically by the deformation of a unit circle with a metric map into an ellipse (see Fig. 1.1 in Loseille and Alauzet, 2011a). However, such geometric interpretation is unwieldy for visualising entire metric fields ℳ, and metric density and anisotropy quotient fields may be more appropriate, as shown in Fig. 3. The figure demonstrates how a high metric density ρ prescribes a small adapted mesh element size, and vice versa, while the anisotropy quotients r₁ and r₂ control the stretching of mesh elements in the directions of the eigenvectors in R. For example, the mesh adapted from the basal shear stress metric prescribes the finest resolution along the grounding line, where elements' shape and orientation follow the geometry of the grounding line. Furthermore, mesh elements along the shear margin are most strongly stretched along the x direction when adapted from the ice thickness metric.

Figure 3 Visualising components of metric fields constructed from the Hessian of the initial (left column panels) ice thickness, (middle column panels) ice velocity, and (right column panels) basal stress fields of MISMIP+ experiments (see Fig. 1): (top row panels) metric density ρ and (middle row panels) anisotropy quotient r₁. The corresponding adapted meshes are shown in the bottom row, including a zoom-in on different regions of the domain. Adjacent adapted edge lengths were not allowed to differ by more than 30 %. The zoomed-in subpanels maintain equal scaling of the x and y axes, while the domain is laterally shortened for legibility in other subpanels.

In cases when the exact solution u is not known, the main challenge lies in finding a reliable estimate of the approximation error u−u_h, which is then used in constructing the metric field ℳ. The metric then drives mesh adaptation in a way that reduces the estimated error. In what follows, we consider computing an error estimate based on interpolation error theory.

3.2 Error estimate based on the interpolation error

Let us consider the mesh adaptation problem stated in an a priori way in Eq. (4). It can be shown that the approximation error $\Vert u-u_h\Vert$ for elliptic problems is bounded by the interpolation error in approximating a continuous field u in V_h:

$$\begin{array}{}\text{(8)}& \Vert u-u_h\Vert \le c\Vert u-{\mathrm{\Pi }}_{h}u\Vert ,\end{array}$$

where Π_hu is the linear interpolant of u onto V_h and c is a mesh-independent constant (Céa's lemma, e.g. Ciarlet, 1991, 2002). In practice, the approach has been found to hold for hyperbolic problems as well (Frey and Alauzet, 2005). Since the solution u and its linear interpolant Π_hu are generally not known, the ideas presented in this work are instead applied a posteriori to the numerical solution u_h. In the context of general adaptation methods, we are interested in how the error scales with polynomial degree (p adaptation, e.g. Cuzzone et al., 2018; Kirby and Shapero, 2024) and mesh size (h adaptation), which we may readily change. In the former, the polynomial degree can be locally adjusted while maintaining a fixed mesh topology, which eliminates the need for mesh-to-mesh interpolation (see Sect. 3.6). However, higher-order approximations typically lead to increased computational cost and memory consumption. In contrast, modifying the local mesh element size provides greater control over computational efficiency by dynamically adjusting the number of degrees of freedom. In this paper we only consider h-adaptation approaches, but both are active and impactful fields of research.

In particular, it can be shown that the local interpolation error $\Vert u-{\mathrm{\Pi }}_{h}u\Vert$ in a single element can be related to the Hessian field of u (see e.g. Ciarlet, 1991). The result has been readily adopted within metric-based mesh adaptation research due to the (normally) low computational cost involved in computing the Hessian, the directional information contained within it, and the problem-independent nature of Hessian-based interpolation errors (Pain et al., 2001; Frey and Alauzet, 2005; Piggott et al., 2009b; Davies et al., 2011). It has also become a fundamental part of the continuous metric framework, with Loseille and Alauzet (2009) presenting the continuous interpolation error estimation involving a continuous linear interpolant π_Mu. The result was then extended in Alauzet and Olivier (2010), who presented a well-posed continuous formulation of a problem (Eq. 4) that minimises the interpolation error in the L^p norm:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathrm{Find}{\mathcal{M}}_{L^p}& \mathrm{which}\mathrm{minimises}{\mathbf{E}}_{L^p}\left(\mathcal{M}\right)\\ & ={\left(\underset{\mathrm{\Omega }}{\int }{\left(u-{\mathit{\pi }}_{\mathbf{M}}u\right)}^p\mathrm{d}\mathit{x}\right)}^{\frac{\mathrm{1}}{p}},\end{array}\end{array}$$

under the constraint that 𝒞(ℳ)=N_v.

For two-dimensional steady-state problems, Loseille and Alauzet (2011b) show that the optimisation problem in Eq. (9) has a unique solution which optimally controls the interpolation error in the L^p norm. The resulting optimal metric is given by

$$\begin{array}{}\text{(10)}& {\mathbf{M}}_{L^p}={\mathcal{C}}_s{\left(\underset{\mathrm{\Omega }}{\int }(\det \Vert \mathbf{H}\Vert {)}^{\frac{p}{\mathrm{2}p+\mathrm{2}}}\mathrm{d}\mathit{x}\right)}^{-\mathrm{1}}(\det \Vert \mathbf{H}\Vert {)}^{\frac{-\mathrm{1}}{\mathrm{2}p+\mathrm{2}}}\Vert \mathbf{H}\Vert ,\end{array}$$

where 𝒞_s is the target metric complexity and $\Vert \mathbf{H}\Vert$ is the absolute value of the Hessian matrix H=H(x) of u, obtained by taking the absolute value of its eigenvalues in the spectral decomposition. In the numerical demonstrations that follow, we construct the approximation of H using the robust mixed-L² projection recovery technique that was heavily influenced by the work of McRae et al. (2018). The approach ensures a continuous representation of second-order derivatives.

The analogous analysis for time-dependent problems was presented in Alauzet and Olivier (2010), which now involves not one but N_a∈ℕ meshes $\mathit{\{}{\mathcal{H}}_i{\mathit{\}}}_{i=\mathrm{1}}^{N_a}$, where each mesh is assigned to a different part of the simulation time interval. They show that the optimal local metric for a subinterval i is given by

$$\begin{array}{}\text{(11)}& {\mathbf{M}}_{L^p}^i=G_{\mathrm{st}}{\left(n{\left(\det \Vert {\mathbf{H}}_i\Vert \right)}^{\mathrm{2}}\right)}^{-\frac{\mathrm{1}}{\mathrm{2}p+\mathrm{2}}}\Vert {\mathbf{H}}_i\Vert ,\end{array}$$

where n is the number of time steps in each subinterval (here constant) and G_st is the global space–time normalisation constant

$$\begin{array}{}\text{(12)}& G_{\mathrm{st}}={\mathcal{C}}_{\mathrm{st}}{\left(n^{\frac{p}{p+\mathrm{1}}}\sum _{j=\mathrm{1}}^{N_a}\underset{\mathrm{\Omega }}{\int }{\left(\det \Vert {\mathbf{H}}_j\Vert \right)}^{\frac{p}{\mathrm{2}p+\mathrm{2}}}\mathrm{d}\mathit{x}\right)}^{-\mathrm{1}},\end{array}$$

where 𝒞_st is the specified space–time metric complexity. The space–time complexity provides an estimate for the average number of mesh vertices in the entire mesh sequence $\mathit{\{}{\mathcal{H}}_i{\mathit{\}}}_{i=\mathrm{1}}^{N_a}$. The number of vertices, however, may vary between meshes, as they are distributed in both space and time among individual meshes in the sequence according to Eq. (11) such that the optimal spatio-temporal error in the L^p norm is obtained.

3.3 Metric intersection

For many applications, adequately capturing the mesh resolution requirements needs information from multiple sources in the domain. For example, it may be desirable to consider all solutions of a multi-variable PDE of interest or a vector-valued field whose components may each instruct a different mesh resolution. Furthermore, a time-dependent field may require significantly different mesh resolution in time in order to adequately resolve its evolution. To this end, individual metrics may be constructed from each such sensor field at multiple time steps before they are combined to form a single metric that is representative of all those individual time steps. The combination of metrics is a non-trivial task, and the resulting metric field may vary significantly for different combination approaches and imply different computational cost. Throughout the paper we consider the most common method of combining metrics: metric intersection. An example of an intersected metric is shown in Fig. 4, where ice thickness and ice velocity metrics of Fig. 3 have been intersected. An inspection of the metric components reveals that the intersected metric indeed contains information from both of the individual metrics.

Figure 4 (a) Metric density ρ and (b) anisotropy quotient r₁ of a metric constructed by intersecting ice thickness and ice velocity metrics shown in Fig. 3, as well as (c) the corresponding adapted mesh.

Metric intersection produces a new metric tensor that captures the most restrictive (i.e. finest) resolution requirements from both original metrics. As described in Alauzet et al. (2007), the process involves transforming the two metrics to a common eigenbasis, taking the maximum eigenvalues from each pair of corresponding eigenvalues and transforming the metric into the original eigenbasis; i.e. given two metrics M₁ and M₂ in two-dimensional space, the intersection M_1∩2 is defined as

$$\begin{array}{}\text{(13)}& {\mathbf{M}}_{\mathrm{1}\cap \mathrm{2}}={\left({\mathbf{P}}^{-\mathrm{1}}\right)}^{\mathsf{T}}\left(\begin{array}{cc}\max \left({\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}\right)& \\ & \max \left({\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}\right)\end{array}\right){\mathbf{P}}^{-\mathrm{1}},\end{array}$$

where P is the matrix of normalised eigenvectors $\mathit{\{}{\mathit{e}}_i{\mathit{\}}}_{i=\mathrm{1},\mathrm{2}}$ of $\mathbf{N}={\mathbf{M}}_{\mathrm{1}}^{-\mathrm{1}}{\mathbf{M}}_{\mathrm{2}}$ and ${\mathit{\lambda }}_i^j={\mathit{e}}_i^{\mathsf{T}}{\mathbf{M}}_j{\mathit{e}}_i$. Related to the geometric interpretation of metrics with an ellipsoid (see Sect. 3), the metric intersection is visualised as producing a metric whose associated ellipsoid is the largest one that can fit within the ellipsoids of the original metrics (see Fig. 2 in Alauzet et al., 2007).

3.4 Metric post-processing

After the metric had been computed, it may be desired, or even required, to post-process the metric field before the adapted mesh is generated. This includes prescribing the minimum and maximum element edge lengths, as well as maximum tolerated element anisotropy, by appropriately truncating eigenvalues λ₁ and λ₂ (see Appendix A). Such constraints are particularly useful in order to incorporate existing physical knowledge into the mesh adaptation process, which may not be contained in the chosen sensor adaptation fields. Preventing large differences in element sizes and anisotropy between original and adapted meshes may also improve robustness.

Certain choices of sensor fields may lead to metric fields that exhibit a very sharp gradient, as is the case for the metric defined from the basal shear stress, shown in Fig. 3. Borouchaki et al. (1998) show that such large rates of variation may lead to poor quality of the adapted meshes. To that end, they introduce a metric gradation routine in order to control its variation and therefore smooth out the metric field.

We refer the reader to Dapogny et al. (2014b) for details on the metric gradation routine and edge length truncation implemented in the remeshing library mmg2d, which is used for all mesh adaptation results in this paper (see Sect. 3.7).

3.5 Time-dependent mesh adaptation algorithms

Mesh adaptation has been successfully integrated with finite element methods since the 1980s, when the literature almost exclusively focused on applying mesh adaptation to steady simulations (e.g. Peraire et al., 1987). The mesh adaptation procedure involved an iterative process which would typically begin with a coarse uniform mesh. The mesh is then iteratively adapted until the convergence of the steady-state solution is reached, according to a mesh convergence study. In this context, “mesh convergence” refers to the process where successive refinements of the mesh result in progressively smaller changes in the solution, indicating that the solution is approaching a stable state that is independent of further mesh refinement. For details on mesh convergence, see e.g. Loseille and Alauzet (2011a).

Early attempts at mesh adaptation for time-dependent simulations simply employed the above-described algorithm multiple times throughout the simulation. This approach is referred to as the classical mesh adaptation algorithm in Alauzet et al. (2007), and we refer to it as such throughout the paper. Namely, the mesh is adapted N_a∈ℕ times throughout the simulation, at times t_i, $i=\mathrm{0},\mathrm{\dots },N_a-\mathrm{1}$, based on the current solution state at that time. In this way, the simulation involves N_a different meshes rather than a single fixed mesh, which allows us to prescribe fine resolution in particular regions of the domain only at times when that is needed. After the mesh had been adapted, the solution and any other spatially varying fields in the PDE are then transferred onto it and the simulation resumes until the end simulation time T.

While suitable for steady simulations, the classical mesh adaptation algorithm exhibits several shortcomings specific to time-dependent simulations. Firstly, the approach does not control the temporal distribution of error since the metric is only normalised in space according to Eq. (10), as only the current solution is used to guide the mesh adaptation. Additionally, by not considering the future evolution of the solution, the adapted mesh may not adequately resolve the solution at subsequent simulation times. In such a case, the mesh is said to lag with respect to the solution. In the context of glaciological modelling, a common consequence of a lagging mesh is the migration of the grounding line out of the fine-resolution region of the mesh and into the coarse region, where the grounding line dynamics are no longer accurately captured (e.g. IGE-Elmer/Ice in Sun et al., 2020). While increasing the mesh adaptation frequency would alleviate the lag, doing so may introduce large errors due to frequent solution transfers between meshes (see Sect. 3.6).

To efficiently avoid the mesh lagging behind the solution, the adaptation algorithm requires a prediction of where fine resolution is likely to be required before the next mesh adaptation. An accurate prediction will ensure that the simulated phenomena remain well-resolved when reducing the total number of mesh adaptations N_a. Several prediction strategies have been proposed over the years, such as metric advection for advection-dominated problems (Wilson, 2010; Smith et al., 2016), and more generally applicable methods, such as the introduction of safety regions around the fine-resolution region of the mesh (Löhner and Baum, 1992). The latter has become the most common mesh adaptation criterion in marine ice sheet modelling, where a fine-resolution region was introduced within a set distance from the grounding line. However, such criteria do not explicitly control the solution error and neglect the error contribution of more distant regions while potentially refining regions where that is not necessary.

Here we adopt the global fixed-point mesh adaptation algorithm described in Alauzet and Olivier (2011) in order to remedy the stated challenges specific to time-dependent problems. The algorithm is presented as pseudo-code in Algorithm 1. First, the simulation time interval (0, T] is partitioned into N_a equally long subintervals:

$$\begin{array}{}\text{(14)}& (t_{j-\mathrm{1}},t_j]=\left(\right(j-\mathrm{1})T/N_a,jT/N_a],\end{array}$$

where j=1, …, N_a and $T/N_a$ is the fixed subinterval length. To each subinterval we assign a mesh ${\mathcal{H}}_j^{\left(k\right)}$, where k is the fixed-point iteration number. At each iteration k, the PDE problem is solved over the entire interval (0, T] on the sequence of meshes $\mathit{\{}{\mathcal{H}}_j^{\left(k\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$. From the obtained solutions we compute Hessian metrics for given discrete time steps and combine them in time according to Eq. (13). The computed metric fields are then space–time normalised according to Eq. (11). In such a way we obtain the metric field ${\mathcal{M}}_j^{\left(k\right)}$ associated with each subinterval (t_j−1, t_j]. Finally, the meshes are adapted to yield the next iteration's adapted mesh sequence $\mathit{\{}{\mathcal{H}}_j^{(k+\mathrm{1})}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$, and the simulation is restarted. The algorithm terminates when a set maximum number of iterations or prescribed convergence criteria have been reached.

Algorithm 1 Global fixed-point mesh adaptation algorithm.

for iteration k=1, …, k_max do

for time subinterval j=1, …, N_a do

Interpolate the initial condition/PDE solution onto ${\mathcal{H}}_j^{\left(k\right)}$

Solve PDE on subinterval (t_j−1, t_j]

Construct subinterval metric ${\mathcal{M}}_j^{\left(k\right)}$

end for

Normalise subinterval metrics $\mathit{\{}{\mathcal{M}}_j^{\left(k\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$ in space and time

Generate adapted meshes $\mathit{\{}{\mathcal{H}}_j^{(k+\mathrm{1})}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$

if converged then

Terminate

end if

end for

The global fixed-point adaptation algorithm offers several advantages over the classical mesh adaptation algorithm. Firstly, the algorithm remedies the lagging mesh problem, since the problem is first solved over the entire time interval. The constructed metric fields therefore contain a prediction of the solution evolution over the entire subinterval. Secondly, metric fields are normalised according to space–time normalisation in Eq. (11), which allows for the error to be controlled in both space and time. On the other hand, since meshes are adapted at the end of each iteration, each subinterval may start with a worse approximation of the initial state compared to the classical mesh adaptation algorithm, which adapts the subinterval mesh before the PDE is solved over the corresponding subinterval. This may imply significant additional computational cost if the initial mesh sequence does not accurately capture the modelled phenomena, which is often the case for coarse uniform meshes. As a result, mesh convergence may be slow and many fixed-point iterations may be required.

Here we propose a combination of the classical and global fixed-point mesh adaptation algorithms, which yields a better initial approximation of the optimal mesh sequence at a minimal additional cost. We achieve this by incorporating the classical algorithm within the first iteration of the global fixed-point algorithm. The latter iterations then proceed as in the global fixed-point iteration, without incorporating the classical mesh adaptation algorithm since the mesh sequence is already approximated well. The proposed hybrid algorithm is presented in Algorithm 2. Since the classical adaptation algorithm generates the adapted mesh sequence $\mathit{\{}{\mathcal{H}}_j^{\left(\mathrm{1}\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$ before the PDE is solved on each subinterval, the only added computational cost compared to the global fixed-point algorithm is that of adapting the meshes. This added cost is normally small compared to the cost of solving the PDE. Solution fields obtained from the first iteration of the hybrid algorithm are therefore more accurate than those from the global fixed-point algorithm.

Algorithm 2 Hybrid fixed-point mesh adaptation algorithm.

for iteration k=1, …, k_max do

for time subinterval j=1, …, N_a do

if k=1 then

Construct ℳ from $u(\mathit{x},t_{j-\mathrm{1}})$

Normalise ℳ in space

Generate adapted mesh ${\mathcal{H}}_j^{\left(k\right)}$ from ℳ

end if

Interpolate the initial condition/PDE solution onto ${\mathcal{H}}_j^{\left(k\right)}$

Solve PDE on subinterval (t_j−1, t_j]

Construct subinterval metric ${\mathcal{M}}_j^{\left(k\right)}$

end for

Normalise subinterval metrics $\mathit{\{}{\mathcal{M}}_j^{\left(k\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$ in space and time

Generate adapted meshes $\mathit{\{}{\mathcal{H}}_j^{(k+\mathrm{1})}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$

if converged then

Terminate

end if

end for

3.6 Interpolation between meshes

The final step of any mesh adaptation procedure is that of interpolating all spatially varying fields to the newly generated mesh. This is a crucial step in the mesh adaptation process, as the choice of transfer method may significantly affect the accuracy of the solution. Since the optimal choice of the transfer method largely depends on the specific problem at hand, we discuss it here with respect to glaciological applications. The process also depends on the type of mesh and mesh adaptivity used. For example, Cornford et al. (2013) use the more restrictive block-structured meshes in order to alleviate issues related to numerical diffusion and mass conservation that we discuss below for unstructured meshes.

The mesh adaptation toolkit Animate, described in Sect. 3.7, currently offers three methods for interpolating fields between different meshes: linear interpolation, Galerkin projection, and the bounded variant of the Galerkin projection (see Farrell, 2010, for an overview of these methods). The most popularly used interpolation method in the majority of applications is that of linear interpolation. However, linear interpolation does not preserve the integral of interpolated fields, and its application to ice flow problems would yield unphysical results due to conservation laws being violated. The use of Galerkin projection is therefore more appropriate since it is a conservative operator despite implying a larger computational cost. In particular, Galerkin projection may introduce new function extrema, which may be problematic in times of extreme glacier retreat when ice thickness is close to zero. In such cases, the interpolated ice thickness might become negative-valued. In order to prevent that, we use the bounded variant of the Galerkin projection for degree-1 basis (P₁) fields, which prevents the addition of new function extrema (Farrell et al., 2009; Farrell and Maddison, 2011). However, by bounding the interpolated solution, the operation introduces additional numerical diffusion which may diminish the accuracy of the overall simulation.

3.7 Software

Firedrake (Ham et al., 2023) is a Python-based finite element library for the numerical solution of PDEs utilised throughout this paper. It provides a high-level interface for the definition of finite element variational forms using the Unified Form Language (UFL; Alnæs et al., 2014) which are then compiled into efficient low-level code for the solution of the PDEs. Firedrake utilises PETSc for solving linear and non-linear systems, as well as for its underlying mesh concept (Balay et al., 1997, 2023; Lange et al., 2016).

Mesh generation in Firedrake can be performed using built-in classes; reading in meshes generated by external generators (such as gmsh; Geuzaine and Remacle, 2009); or using Netgen (Schöberl, 1997), an external generator that has been integrated into Firedrake. Recent developments involve the Animate library for metric-based anisotropic mesh adaptation (Wallwork et al., 2024a). Animate implements the presented Riemannian metric framework, Hessian recovery methods, and metric operations and interfacing with PETSc's Riemannian metric functionality (Wallwork et al., 2022). A separate library, Goalie, has been developed for time-dependent mesh adaptation and supports both feature- and goal-oriented adaptation approaches (Wallwork et al., 2024b). Both libraries have undergone development synchronously with the development of the work presented in this paper.

The process of locally adapting meshes with respect to metric fields is performed by the external library Mmg2d, which is part of the Mmg platform (Balarac et al., 2022; Dobrzynski and Frey, 2008) that has been integrated into PETSc. Mmg2d requires mesh topology and a metric field defined at each mesh node as input, and it returns the mesh where elements' size, shape, and orientation have been adapted to the given metric field. While the work presented in this paper focuses on 2D mesh adaptation, Animate and Goalie also support 3D mesh adaptation through the implementation of the Mmg3d and ParMmg mesh generation toolkits (Dapogny et al., 2014a).

4 Results

In this section we demonstrate the ability of a Hessian-based anisotropic mesh adaptation method to control the solution accuracy in the MISMIP+ Ice1 experiment. We consider constructing metric fields from the Hessian of ice thickness H and ice velocity u, their intersection, and the Hessian of basal stress τ_b. The basal stress field is computed from Eq. (3) and the hydrostatic approximation for the normal basal stress, and as such it depends on both H and u.

In the absence of a closed-form solution, we rely on the 0.25 km resolution results to define reference solutions u_ref(x,t) and H_ref(x,t). We validate the results of each simulation by computing the relative solution error in L² norm based on the reference solutions:

$$\begin{array}{}\text{(15)}& {\displaystyle }{\stackrel{\mathrm{̃}}{e}}_H={\displaystyle \frac{\Vert H_{\mathrm{int}}-H_{\mathrm{ref}}{\Vert }_{L^{\mathrm{2}}}}{\Vert H_{\mathrm{ref}}{\Vert }_{L^{\mathrm{2}}}}},\text{(16)}& {\displaystyle }{\stackrel{\mathrm{̃}}{e}}_{\mathit{u}}={\displaystyle \frac{\Vert \Vert {\mathit{u}}_{\mathrm{int}}\Vert -\Vert {\mathit{u}}_{\mathrm{ref}}\Vert {\Vert }_{L^{\mathrm{2}}}}{\Vert \Vert {\mathit{u}}_{\mathrm{ref}}\Vert {\Vert }_{L^{\mathrm{2}}}}},\end{array}$$

where H_int and u_int are solution fields interpolated onto the reference 0.25 km mesh and $\Vert \mathit{u}\Vert$ is the velocity magnitude. We further consider the maximum deviation of the volume above flotation,

$$\begin{array}{}\text{(17)}& \Vert \mathrm{\Delta }{V}_{\mathrm{f}}{\Vert }_{\mathrm{\infty }}=\underset{i}{max}\Vert V_{\mathrm{f}}\left(t_i\right)-V_{\mathrm{f},\mathrm{ref}}\left(t_i\right)\Vert ,\end{array}$$

and midchannel grounding line position,

$$\begin{array}{}\text{(18)}& \Vert \mathrm{\Delta }{x}_{\mathrm{gl}}{\Vert }_{\mathrm{\infty }}=\underset{i}{max}\Vert x_{\mathrm{gl}}\left(y=\mathrm{0},t_i\right)-x_{\mathrm{gl},\mathrm{ref}}\left(y=\mathrm{0},t_i\right)\Vert ,\end{array}$$

over all simulation years t_i, where V_f,ref and $x_{\mathrm{gl},\mathrm{ref}}(y=\mathrm{0})$ are the ice volume above flotation and grounding line position along the midchannel line, respectively, for the reference 0.25 km uniform mesh simulation (see Fig. 2).

4.1 Mesh adaptation strategy

In our mesh adaptation experiments we do not incorporate the gained insight into resolution sensitivity from Sect. 2.2. As such, we do not prescribe bounds on element edge lengths and element anisotropy. We use Mmg2d's default maximum metric gradation factor of 1.3, which prevents adjacent element edge lengths from differing by more than 30 %.

Crucially, when solving time-dependent mesh adaptation problems, we must decide how to split the time domain into N_a distinct subintervals. As discussed in Sect. 3.5, the choice should consider the balance between accurately resolving the ice stream in time and the added interpolation error and computational cost associated with each mesh adaptation. To that end, we choose to split the simulation time interval (0, 200 a] into 20 equally long subintervals, meaning that meshes are adapted every 10 years of simulation time. A constant time step $\mathrm{\Delta }t=\mathrm{1}/\mathrm{24}$ a is again used to ensure that the CFL condition is satisfied.

Similarly, we choose to consider solutions at every 1 year of simulation time in the construction of the metric fields. This provides us with 10 metric fields per subinterval, which is, again, a compromise between accurately representing the evolution of the ice stream and avoiding further computational cost associated with metric computations. The 10 metric fields are combined in time using the metric intersection of Eq. (13) to yield a single subinterval metric. This is repeated for each of the 20 subintervals. Afterwards, the 20 subinterval metrics are normalised in space and time according to Eq. (11).

It must be noted that there are certainly other viable choices of partitioning the time interval which we do not explore here. For example, considering that the thinning of the ice shelf and grounding line retreat is more extreme in the first half of the experiment (Fig. 2), a reasonable choice might involve more frequent adaptations in the first half of the time interval than in the second half. However, we do not incorporate such prior knowledge in the set-up of the experiments, and rely on space–time metric normalisation to distribute available degrees of freedom accordingly. Similarly, we always prescribe an equal target complexity to each subinterval metric in the first iteration of the hybrid fixed-point algorithm, which cannot utilise space–time normalisation. This will result in the numbers of vertices, N_v, of each mesh of the mesh sequence $\mathit{\{}{\mathcal{H}}_j^{\left(\mathrm{1}\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$ to be approximately equal.

For interpolating fields between meshes, we use the bounded variant of the Galerkin projection. As discussed in Sect. 3.6, its use is most appropriate in the experiment considered here as it is a conservative operation that does not introduce new extrema in the solution. Since the bounded variant of the Galerkin projection is only implemented for P₁ fields, we use P₁ functions to represent both ice velocity and ice thickness over elements. The leading-order error for this discretisation typically scales quadratically with the decreasing mesh size in the L² norm. However, the reduced regularity due to a discontinuity in the derivative of τ_b in practice results in a degradation of the convergence rate.

4.2 Comparing time-dependent mesh adaptation algorithms

We first wish to compare the mesh convergence rate of the global and hybrid fixed-point adaptation algorithms described in Sect. 3.5. To that end, we run experiments for the two approaches. The first iteration of the global fixed-point algorithm uses an initial mesh sequence $\mathit{\{}{\mathcal{H}}_j^{\left(\mathrm{1}\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$ of N_a=20 uniform 4 km resolution meshes, while the hybrid algorithm begins with meshes generated using the classical mesh adaptation algorithm. Meshes are adapted based on the Hessian of the basal stress τ_b, but the same conclusions follow from other choices of the adaptation sensor field, as discussed below. We repeat the procedure for two choices of target complexities 𝒞: 1600, which is comparable to the number of vertices of the uniform 4 km resolution mesh (see Table 1), and 3200.

As shown in Fig. 5, the global and hybrid mesh adaptation algorithms converge to qualitatively similar adapted meshes. However, the hybrid algorithm converges to a reasonable approximation of the optimal mesh twice as fast, in only three iterations, while the global algorithm converges in seven iterations, as shown in Fig. 6.

Figure 5 Comparing final adapted meshes obtained from the global fixed-point and hybrid mesh adaptation algorithms for target complexity C=1600. Final adapted meshes ℋ₉ corresponding to the (90, 100 a] subinterval are shown in (a) and (b), and final subinterval meshes ℋ₁₉ corresponding to the (190, 200 a] subinterval are shown in (c) and (d), with a zoom-in around the grounding line on the right. Meshes adapted using the global algorithm are shown in (a) and (c), while meshes from the hybrid algorithm are shown in (b) and (d).

Figure 6 Comparing solution convergence of global and hybrid mesh adaptation algorithms, with basal stress as the adaptation sensor field and target complexities 𝒞 of 1600 and 3200. The resulting maximum deviations from the reference solutions of the (a) volume above flotation, $\Vert \mathrm{\Delta }{V}_{\mathrm{f}}{\Vert }_{\mathrm{\infty }}$, and (b) midchannel grounding line position, $\Vert \mathrm{\Delta }{x}_{\mathrm{gl}}(y=\mathrm{0}){\Vert }_{\mathrm{\infty }}$, over all simulation years are shown.

In the first iteration of the global fixed-point algorithm, the initial mesh sequence of uniform 4 km resolution is not able to adequately resolve the ice stream dynamics, as shown previously in Sect. 2.2. The poor solution accuracy in the first iteration in turn leads to sub-optimal metric fields and the adapted mesh sequence used in the next iteration. This result is independent of the choice of sensor field or target complexity used in constructing the metric, since it is the initial uniform mesh sequence that is far from optimal. In comparison, deviations in the initial iteration of the hybrid fixed-point algorithm are up to an order of magnitude smaller since the ice stream is not modelled on the inadequate coarse uniform mesh sequence.

Furthermore, the two algorithms converge to adapted meshes with consistent vertex counts, as shown in Figs. 5 and 7. However, the meshes adapted in the first iteration of the hybrid algorithm are only normalised in space (see Algorithm 2), meaning that the number of vertices of each mesh in the mesh sequence $\mathit{\{}{\mathcal{H}}_j^{\left(\mathrm{1}\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$ is therefore approximately equal. As in the global fixed-point algorithm, the metric fields used to generate the next mesh sequence $\mathit{\{}{\mathcal{H}}_j^{\left(\mathrm{2}\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$ are constructed from solutions sampled from the entire time domain, so the metric fields are able to be normalised in time as well as in space. The adapted mesh sequences $\mathit{\{}{\mathcal{H}}_j^{\left(k\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$, k≥2 now contain a higher number of vertices in the first half of the simulation and lower number in the second half. This is expected since the grounding line migration is more pronounced in the first half of the simulation and requires a more widely refined mesh in order for the grounding line to remain within the finely refined mesh region as it migrates in time.

Figure 7 Temporal distribution of number of vertices, N_v, for several iterations of the hybrid algorithm and the final iteration of the global fixed-point algorithm for 𝒞=1600.

Overall, we have found that a reasonable mesh convergence may be expected in only a few iterations of the hybrid fixed-point algorithm. Namely, most of the effort is done in the first two iterations of the algorithm. Given a reasonable time interval partition and sufficient target complexity, the classical mesh adaptation algorithm embedded within the first iteration generates an adequate adapted mesh sequence $\mathit{\{}{\mathcal{H}}_j^{\left(\mathrm{1}\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$. The mesh sequence $\mathit{\{}{\mathcal{H}}_j^{\left(\mathrm{2}\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$ is then generated at the end of the first iteration from metric fields normalised in both space and time, which allocates additional resolution to subintervals where that is necessary. Solutions obtained in the second iteration are therefore even more accurate, as the error is optimally controlled in space and distributed in time. The mesh sequence generated at the end of the second iteration, $\mathit{\{}{\mathcal{H}}_j^{\left(\mathrm{3}\right)}{\mathit{\}}}_{j=\mathrm{1}}^{N_a}$, can therefore be assumed to be a reasonable approximation of the optimal mesh sequence.

4.3 Comparing adaptation sensor fields

Finally, we study the mesh convergence and solution accuracy in simulations utilising the hybrid fixed-point adaptation algorithm for different adaptation sensor fields. We repeat the mesh adaptation process for three global fixed-point iterations, which implies a total of four adaptations of each mesh. As we demonstrated in Sect. 4.2, this is enough to achieve reasonable mesh convergence. For each choice of sensor field, we prescribe the target metric complexity 𝒞 of 800 and repeat the simulation four more times such that the target complexity is doubled each time. The final simulation therefore specified the target complexity of 16𝒞=12 800. Adapted meshes ℋ₉ are shown in Fig. 8 for simulations specifying target complexity of 8𝒞, with their elements' aspect ratios shown in Fig. 9. The spatial error distribution at t=100 a on ℋ₉ for 8𝒞 is shown in Fig. 10, while the time-averaged errors are shown in Fig. 11.

Figure 8 Adapted meshes ℋ₉ associated with the (90, 100 a] subinterval, with a zoom-in around the grounding line on the right. Metrics were constructed from (a) H, (b) u, (c) intersected metric of H and u, and (d) τ_b. Shown grounding line contours are computed at the beginning and end of the subinterval from the reference 0.25 km uniform-resolution mesh simulation solutions.

Figure 8 reveals instructive areas of refinement and coarsening for each constructed metric field. When adapting based on ice thickness, the fine resolution is concentrated in two regions: along the shear margin and on the ice shelf near the grounding line. This means that the grounding line will easily migrate out of the fine-resolution region during glacier retreat, as shown by the grounding line contours of the reference solution in Fig. 8a. On the other hand, the grounding line will remain in the fine region during glacier advance as it migrates further downstream. In comparison, meshes adapted based on ice velocity are mostly refined in the vicinity of the grounding line and in the inner part of the ice shelf, where ice quickly speeds up as it is no longer influenced by friction at the bed. However, the region around the grounding line is not as finely resolved as in Fig. 8a; instead, the available resolution is more or less equally distributed over a much larger region. As expected, meshes adapted from the intersected metric of the two fields exhibit features of both. However, this also includes a less refined region upstream of the grounding line compared to that in Fig. 8b, which again leads to the migration of the grounding line into the coarser region of the mesh. Meanwhile, because the basal stress rapidly diminishes across the grounding line, the metric field computed from its Hessian will prescribe the finest mesh resolution in its vicinity, as shown in Fig. 8d. This ensures that the grounding line remains in the finely resolved region of the mesh throughout the subinterval. However, since basal stress is null in the ice shelf, the region of the mesh corresponding to it is coarsened. Variation in mesh cell sizes in the ice shelf is only due to the metric gradation routine.

A measure of cell aspect ratios is shown in Fig. 9, where a unit aspect ratio corresponds to a perfectly isotropic cell (i.e. an equilateral triangle). We observe that meshes whose metrics involve ice thickness contain two bands of highly anisotropic cells in the shear margin with a predominantly x-directional orientation separated by a band of isotropic cells. Since the variation in ice thickness is much sharper across the shear margin than along it, anisotropic elements help to efficiently capture such directional processes along the shear margin by avoiding unnecessary over-refinement in the x direction. Similarly, in order to capture the rapid variation in the basal stress across the grounding line, meshes adapted based on the basal stress have elements stretched in the direction of the grounding line along its entirety. On the other hand, meshes adapted from ice velocity contain mesh elements with lower anisotropy along the ice stream shear margin but strongly anisotropic elements in the x direction along the grounding line flank as the ice drastically speeds up across it.

The above observations are directly reflected in the spatial distribution of the solution errors, shown at t=100 a in Fig. 10. While the error in H along the shear margin is lowest on the mesh adapted from H, the same meshes lead to the highest errors in the vicinity of the grounding line and in the ice shelf due to the grounding line migrating out of the fine-resolution region (see Fig. 8a). A more widely refined mesh generated from u is able to most accurately model the sharp velocity changes in the ice shelf but does not leave enough resolution in the direct vicinity of the grounding line. Overall, the lowest error in both H and u is obtained on the mesh generated from τ_b, which contains enough resolution in the direct vicinity of the grounding line to be able to accurately represent it throughout the 10-year subinterval (see Fig. 8d). However, the poorest solution accuracy is obtained along the grounding line flank along the shear margin and in the ice ridge for all choices of adaptation sensor fields. This was also observed in Fig. 2, where the 0.5 km uniform-resolution results most deviate from the reference 0.25 km resolution solutions in these areas.

Figure 9 Aspect ratio of elements of adapted meshes ℋ₉ shown in Fig. 8, corresponding to the time subinterval (90, 100 a]. The meshes are adapted from the Hessian of (a) H, (b) u, (c) intersected metric of H and u, and (d) τ_b.

Figure 10 Spatial distribution of errors in H and u relative to 0.25 km uniform-resolution solutions shown at t=100 a for meshes ℋ₉ generated from the Hessian-based metrics of (c) H, (d) u, and (e) τ_b and 𝒞=6400.

Similar conclusions follow by investigating adapted meshes at other subintervals and error distribution at other times, which are then reflected in the time-averaged ${\stackrel{\mathrm{̃}}{e}}_H$ and ${\stackrel{\mathrm{̃}}{e}}_{\mathit{u}}$ shown in Fig. 11. We observe that ${\stackrel{\mathrm{̃}}{e}}_H$ decays at approximately a linear rate for all experiments, while the convergence of ${\stackrel{\mathrm{̃}}{e}}_{\mathit{u}}$ tends towards a quadratic rate for adaptive mesh experiments. At low target complexities, mesh sequences generated from the basal stress τ_b yield the most accurate solutions since they most accurately represent the evolving grounding line. However, the increase in accuracy gained by further refining such meshes is effectively capped, since the metric computed from τ_b does not prescribe refinement in the ice shelf. Nonetheless, we obtain a solution accuracy comparable to that of a uniform 0.5 km resolution mesh sequence which has on average nearly 30 times fewer vertices. On the other hand, errors computed on mesh sequences adapted from H, u, and their intersection are all similarly high for target complexities of 800 and 1600. Differences begin to be more pronounced for a target complexity of 3200, when the metric from u prescribes finer resolution in the vicinity of the grounding line compared to the metric from H. The errors begin to be comparable to that of the uniform 0.5 km uniform mesh simulation at the target complexity of 6400, with 8 to 12 fewer vertices. Unlike the meshes adapted from τ_b, the solution error for mesh sequences adapted from H and u keep decreasing for higher target complexities, as shown for the final target complexity of 12 800.

Figure 11 Accuracy of solutions computed on adapted meshes relative to 0.25 km uniform-resolution solutions averaged over the simulation time interval of modelled (a) ice thickness H and (b) ice velocity u for different sensor fields used in constructing the metric and for increasing specified target complexity: 𝒞, 2𝒞, 4𝒞, 8𝒞, and 16𝒞, where 𝒞=800. Solution accuracy of simulations involving uniformly refined meshes are shown for comparison.

4.4 CPU time efficiency

Even though we have been successful in outperforming uniform-resolution meshes at a fraction of the total number of vertices, this does not necessarily translate into a successful reduction in total computational time. Time-dependent mesh adaptation routines presented in this paper imply additional costs related to metric computation, mesh generation, and function interpolations between meshes. None of these are present in single, fixed-mesh simulations. Moreover, the procedure is iterative, which involves repeatedly solving the PDE problem and adapting the mesh sequence until the mesh sequence has converged. Over all iterations of the hybrid mesh adaptation algorithm, the PDE problem was solved 4 times, the mesh sequence was adapted 4 times, and ice velocity and thickness fields were each transferred 80 times between subsequent meshes. The final decision behind which adaptation method (if any) to employ should therefore consider the total computation time necessary to achieve a desired solution accuracy. In Fig. 12 we therefore plot the average ${\stackrel{\mathrm{̃}}{e}}_H$ and ${\stackrel{\mathrm{̃}}{e}}_{\mathit{u}}$ versus total CPU time for different choices of adaptation sensor field and for target complexities of 800, 3200, and 12 800. Thanks to PETSc's and Firedrake's logging infrastructure, we are able to easily extract the time involved in each of the above-mentioned components, which we show in Fig. 13. As in the uniform-resolution simulations, all simulations were run in serial for the purposes of CPU time measurements.

Figure 12 Performance and efficiency comparison of the hybrid mesh adaptation algorithm for different sensor fields and target complexities of 800, 3200, and 12 800. Experiments were run on a single CPU, and the algorithm terminated after three iterations. Average (a) ${\stackrel{\mathrm{̃}}{e}}_{\mathrm{H}}$ and (b) ${\stackrel{\mathrm{̃}}{e}}_{\mathit{u}}$ versus CPU time for different adaptation sensor fields: H, u, intersection of H and u, and τ_b, compared to uniform mesh refinement.

Figure 13 Distribution of CPU time across different computational tasks: PDE solve, mesh adaptation (including metric construction), and interpolation.

As shown in Fig. 13, solving the PDE problem accounted for approximately only 50 % of the computational time for 𝒞=800 and 75 % for 𝒞=12 800. This indicates that the efficiency of mesh adaptation methods is decreased for meshes with a small number of vertices due to the relatively high costs of mesh adaptation and solution interpolation. However, as the solver or mesh complexity increases, the benefits of mesh adaptation become more pronounced. The cost of mesh adaptation is particularly high for vector-valued sensor fields since the Hessian metric of each component must be computed separately and then intersected to form a single metric. The intersected metric of H and u therefore implies the highest computational cost, since it involves three Hessian field computations (one for the scalar field H and one for each component of the vector field u) and two intersections for each metric (one for intersecting components of u and one for intersecting H and combined components of u). As a result, the superior performance of the intersected metric field (H, u) observed in Fig. 11 is not reflected in the overall computational efficiency shown in Fig. 13. Meanwhile, adapting the meshes based on τ_b again clearly outperforms other sensor fields at low target complexities. Such adapted meshes achieve a solution accuracy comparable to the 0.5 km uniform-resolution mesh with a 6-fold reduction in total computational time. All other sensor fields considered in this section achieve the same accuracy with approximately a 3-fold reduction in total computational time.

5 Conclusions

In this study we have demonstrated the effectiveness of anisotropic metric-based mesh adaptation for ice flow modelling. All libraries used here are Python-based and built on top of Firedrake, which, in combination, provide an intuitive and accessible approach to finite element ice flow modelling and mesh adaptation.

Using the set-up of the Ice1 experiment of the Marine Ice Sheet Model Intercomparison Project (MISMIP+), we have shown that feature-based anisotropic mesh adaptation based on the Hessian of the solution-dependent field allows us to efficiently optimise the number of vertices in space and time. Our results confirm the importance of adequately modelling grounding line dynamics and the necessity of fine mesh resolution in its vicinity. Adapting the mesh sequence based on the basal stress is therefore the most efficient option, since the finest mesh resolution is prescribed along the grounding line at a particular expense of a coarsened ice shelf, where the basal stress is zero. Meshes adapted from ice velocity fields are the most efficient at capturing the ice shelf dynamics, while ice thickness prescribed the most resolution along the shear margin and in the region of the ice shelf near the grounding line. All choices of the adaptation sensor fields were able to achieve accuracy in ice thickness within a few per mille and the ice velocity magnitude within a few percent of the reference 0.25 km uniform-resolution simulation results. Specifically, we show that adapted mesh sequences are able to reach the same solution accuracy as the uniform 0.5 km resolution mesh but with 8–30 times fewer degrees of freedom. However, this translates into a 3–6 times lower computational cost, reflecting the cost of our iterative mesh adaptation procedure.

We emphasise the importance of the interpolation stage associated with each mesh adaptation. The bounded variant of the projection operator used in this paper introduces a potentially significant amount of numerical diffusion, which impairs the overall simulation accuracy. Future Animate development will focus on developing a post-processing routine in order to achieve a minimally diffusive projection, as described in Farrell et al. (2009).

While results here have been obtained at a significantly lower computational cost using mesh adaptivity, this study reveals promising possibilities for further time reduction. As shown in Sect. 4.4, the main contributor to the computational cost of the mesh adaptation methods employed here is their inherent iterative procedure. To that end, we have demonstrated one possible way of reducing the total number of iterations: by combining the classical mesh adaptation algorithm with the global fixed-point algorithm in order to obtain a more optimal initial mesh sequence. The presented hybrid algorithm required approximately 50 % fewer iterations in order to reach mesh convergence than the global fixed-point algorithm alone while still controlling spatial error and its temporal distribution. As noted in Sect. 4.1, further computational gain may be obtained from a strategic time interval partition, which may reduce the total number of mesh adaptations and solution transfers. Based on these findings, our future work is focused on developing new mesh adaptation schemes inspired by the “on-the-fly” nature of the classical mesh adaptation algorithm that will further reduce total simulation time. We expect that this will allow for even higher-fidelity large- and global-scale glaciological modelling, where each iteration of the fixed-point mesh adaptation algorithm implies significant computational cost. However, more research is needed to investigate the effectiveness of these methods in larger-scale and less idealised models.

Appendix A: Poisson equation

Following a similar example as in Piggott et al. (2009a), we demonstrate the anisotropic metric-based mesh adaptation presented in this paper on a relatively simple example of a Poisson equation. We consider a unit square domain $\mathrm{\Omega }=[\mathrm{0},\mathrm{1}{]}^{\mathrm{2}}$ and homogeneous Dirichlet boundary conditions. Using the method of manufactured solutions (Roache, 2001), we select the solution $u(x,y)=(\mathrm{1}-e^{-x/\mathit{ϵ}})(x-\mathrm{1})\sin \left(\mathit{\pi }y\right)$, where ϵ=0.01, which is the solution to the following steady Poisson problem:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\mathrm{\nabla }}^{\mathrm{2}}u=& \left(\left(\mathrm{2}/\mathit{ϵ}-\mathrm{1}/{\mathit{ϵ}}^{\mathrm{2}}\right)e^{-x/\mathit{ϵ}}-{\mathit{\pi }}^{\mathrm{2}}(x-\mathrm{1})\left(\mathrm{1}-e^{-x/\mathit{ϵ}}\right)\right)\\ & \sin \left(\mathit{\pi }y\right).\end{array}\end{array}$$

As shown in Fig. A1, the solution exhibits a sharp gradient in the x direction near the x=0 boundary since the term $\mathrm{1}-e^{-x/\mathit{ϵ}}$ decays rapidly for ϵ=0.01.

Figure A1 Contour plot of the manufactured solution $u(x,y)=(\mathrm{1}-e^{-x/\mathit{ϵ}})(x-\mathrm{1})\sin \left(\mathit{\pi }y\right)$ to the Poisson problem (Eq. A1).

The optimal adaptive mesh is found approximately and by iteration. The process begins with a (coarse) uniform mesh on which we compute an approximation of the solution u_h. The mesh is then adapted, and the solution u_h is again computed on this new, adapted mesh. The process repeats until a set maximum number of iterations, or some other criterion, is reached. We expect a reasonable metric definition to lead to a more optimal mesh, which in turn leads to increased accuracy, particularly in the first few iterations when we transition from a uniform mesh to adaptive meshes.

The first step of a feature-based mesh adaptation process is the construction of a metric field from the Hessian of a chosen feature, i.e. sensor field, as described in Sect. 3–3.4. In the relatively simple example of a Poisson problem (Eq. A1), the only choice of a sensor field is the solution itself. We can approximately control several parameters during metric construction, such as the maximum tolerated metric anisotropy, a_max. By choosing the value of a_max of 1 and 16 we obtain a pair of isotropic and anisotropic metric fields, respectively, whose components – metric density ρ and anisotropy quotients r₁ and $r_{\mathrm{2}}=r_{\mathrm{1}}^{-\mathrm{1}}$ – are shown in Fig. A2. Metric density is highest along the x=0 boundary and the y=0.5 midline in both metrics, but it is over an order of magnitude larger in the isotropic metric. As expected, anisotropy quotients of the isotropic metric are uniformly one, while the anisotropy quotients of the anisotropic metric vary spatially: r₁ (r₂) is $\mathrm{1}/\mathrm{16}$ (16) along the x=0 boundary and 16 ($\mathrm{1}/\mathrm{16}$) along the y=0.5 midline.

Figure A2 (a, b) Metric density and (c, d) anisotropy quotient r₁ of metrics defined from the Hessian of the solution of Eq. (A1), shown in Fig. A1, with maximum tolerated element anisotropy, a_max, of (a, c) 1 and (b, d) 16.

As described in Sect. 3.1, metric density ρ controls the size of the mesh elements during mesh adaptation, while anisotropy quotients r₁ and $r_{\mathrm{2}}=r_{\mathrm{1}}^{-\mathrm{1}}$ control their shape. The influence of metric components is directly seen on meshes adapted from these metrics, shown in Fig. A3. In both cases the region near the x=0 boundary is finely resolved, but the isotropic mesh requires nearly 4 times as many vertices as the anisotropic mesh to do so. This is directly reflected in the respective metric densities, where the density of the isotropic metric was over a magnitude larger than that of the anisotropic metric. Since the variation along the y direction is not as rapid as along the x direction, elements of the anisotropic mesh are stretched along the y direction, thus avoiding unnecessary over-refinement along the x=0 boundary. We further observe refinement along the y=0.5 midline, with the metric density and resulting number of elements again being higher in the isotropic case.

Figure A3 (a) Isotropic and (b) anisotropic meshes adapted based on the Hessian of the computed solution u_h of Eq. (A1). A close-up of the domain region $[\mathrm{0},\mathrm{0.15}]\times [\mathrm{0.3},\mathrm{0.7}]$ shows fine resolution along the x=0 boundary in both meshes. The isotropic mesh contains a much larger number of vertices N_v=1625 compared to the anisotropic mesh with only N_v=462, mainly due to the large number of isotropic elements along the x=0 boundary.

Since the analytical solution is available, we can easily compute the accuracy of the computed solution u_h by comparing it to the analytical solution u in L² norm. The error in the solutions computed on the isotropic and anisotropic meshes in Fig. A3 is roughly the same ($\mathrm{4.37}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{7.07}\times {\mathrm{10}}^{-\mathrm{5}}$, respectively), even though the anisotropic mesh has roughly 4 times fewer vertices. We can conclude that there was indeed no significant benefit from finely resolving the region near the x=0 boundary along the y direction.

We repeat the process again by considering intermediate values of a_max and by varying the total number of vertices of generated meshes. For comparison, we also solve the problem on uniform meshes. For each solution we compute the error and plot it in Fig. A4 against the number of vertices of the mesh on which it was computed. We observe that adaptive simulations achieved the same error as uniform simulations with over an order-of-magnitude fewer mesh vertices. In particular, there is a clear reduction in error for the same number of vertices when considering increasingly anisotropic meshes. However, we do not observe a clear improvement of increasing a_max from 8 to 16. Had the solution u had an even steeper gradient, we would expect to see benefit from further anisotropy.

Figure A4 Error in the numerical solution u_h to the Poisson problem given in Eq. (A1) against the number of vertices N_v of the mesh on which it was computed. Maximum tolerated metric anisotropy values $a_{\max }\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{4},\mathrm{8},\mathrm{16}\mathit{\}}$ were prescribed during metric construction.

Uniform refinement achieves the expected theoretical second-order convergence rate ($\mathcal{O}\left(N_{\mathrm{v}}^{-\mathrm{2}}\right)$), while adaptive refinement achieves slower convergence rates, between $\mathcal{O}\left(N_{\mathrm{v}}^{-\mathrm{1.35}}\right)$ and $\mathcal{O}\left(N_{\mathrm{v}}^{-\mathrm{1.48}}\right)$, although it should be noted that there is no longer a simple relationship between N_v and (minimum or average) element size for the adapted meshes. The errors in uniform and adaptive (isotropic and anisotropic) meshes would converge for N_v that is high enough, when even the uniform mesh is fine enough near the x=0 boundary but also other regions where such fine resolution is unnecessary. Anisotropic mesh adaptation prevents such unnecessary over-refinement by distributing the available resolution where it is necessary, according to the prescribed metric. Such an approach ensures that available computing resources are efficiently utilised.