ISMIP-HOM benchmark experiments using Underworld

. Numerical models have become an indispensable tool for understanding and predicting the flow of ice sheets and glaciers. Here we present the full-Stokes software package Underworld to the glaciological community. The code is already well established in simulating complex geodynamic systems. Advantages for glaciology are that it provides a full-Stokes solution for elasto-visco-plastic materials and includes mechanical anisotropy. Underworld uses a material point method to track the full history information of Lagrangian material points, of stratigraphic layers and of free surfaces. We show that Underworld successfully reproduces the results of other full-Stokes models for the benchmark experiments of the ISMIP-HOM project. Furthermore , we test FE meshes with different geometries, and highlight the need to be able to adapt the finite-element grid to discontinuous interfaces between materials with strongly different properties, such as the ice-bedrock boundary.

Underworld includes mechanical anisotropy (Moresi et al., 2006;Sharples et al., 2016). It employs the Material Point Method (MPM) (Sulsky et al., 1994;Moresi et al., 2003) where Lagrangian material points are combined with a Finite Element mesh. First and foremost, these material points allow tracking of the strain history and rheological or physical changes on distinct Lagrangian points. Further, tracking of the material points allow us to understand the deformation of individual volumes or layers within the ice sheet and the evolution of the surface. Particles can also be used to record the crystallographic preferred orientation (CPO), and thus the local mechanical anisotropy of the material. This way, the mechanical anisotropy can evolve as a result of the local deformation. The combination of both, anisotropic rheology and particle tracking has potentialis advantageous for the modeling of large-scale folds of stratigraphic layers observed in ice sheets (Wolovick, 2014;NEEM community members, 2013;Bons et al., 2016;Cavitte et al., 2016;Leysinger-Vieli et al., 2018), in particular when the folding is a result of the anisotropic rheology of ice (Bons et al., 2016).
Finally, Underworld can be coupled with other models to investigate surface effects, such as geological sedimentation and erosion, and processes that affect the base of the model, such as mantle deformation and heat flux processes (Salles et al., 2016;Bahadori et al., 2022). These) that have their equivalents at the surface of ice sheets in the form of snow precipitation, ablation and both surface and basal melting (e.g. Jacobson and Raymond, 1998;Smith-Johnsen et al., 2020). melting. For all these reasons Underworld, which is already well established for the simulation of complex tectonic processes (for instance Sandiford et al., 2020;Carluccio et al., 2019;Capitanio et al., 2019;Korchinski et al., 2018), surface processes ( Bahadori et al., 2022REF) and long-term ground water motion (Mather et al., 2022), seems also well suited to simulate ice-sheet and glacier flow.
Any numerical model needs to be validated or benchmarked. The Ice Sheet Model Intercomparison Project for Higher-Order Models (ISMIP-HOM;  and supplement or https://frank.pattyn.web.ulb.be/ismip/welcome.html; http://homepages.ulb.ac.be/~fpattyn/ismip/) provides tests for the comparison of computational ice-sheet flow models for different purposes. "Higher-order" here refers to models that go beyond the shallow-ice approximation (SIA) up to full Stokes solutions (as Underworld does).
ISMIP-HOM includes both 2D and 3D experiments. The flow law is Glen's law with a stress exponent n=3 and, in one experiment, Newtonian flow. In this paper we publish the results for the full suitesuites of experiments of the benchmarkstandard. We focus on three issues: (i) the viability of the results as compared to solutions provided by other models, (ii) the computation time, and (iii) the influence of the geometry of the underlying Finite Element grid. The tests are performed using the 2.10 release of the software package Underworld. Finally, we provide one example of how mechanical anisotropy and tracking of the stratigraphy can be incorporated in Underworld to illustrate the potential of Underworld to simulate mechanically complex systems and the resulting structures within a glacier or ice sheet.

Governing equations
The solution in Underworld is based on the Stokes equation of slow flow of a Newtonian incompressible fluid: ∂ v i ∂ x i =0. (2) Here σij is the absolute stress tensor, P the pressure, g the gravitational acceleration and v the velocity (see Tables 1 and 2 for symbols used). Simulations are based on Glen's flow law for viscous flow (Glen, 1955), according to which the strain rate (ε ij) is proportional to the deviatoric stress (τ ij ) to the power n, the stress exponent. This flow law can be written as: where A is the temperature-dependent rate factor, τII the second invariant of the deviatoric stress tensor τ ij (Nye, 1953).
Based on Newtonian flow, where τ ij =2ηε ij , we define an effective viscosity ηice after Eq. 3 as:  x, y axes parallel and vertical to the tilted surface, referred to as 'horizontal' / 'vertical' m

Characteristics of Underworld with regard to specific challenges in the modeling of ice flow
Underworld is designed to solve some of the special problems relevant to modeling geodynamic processes. Identical problems arise in the modeling of ice. Some of these challenges are: a the modeling of discontinuities in the material properties at layer boundaries, for instance at the ice-rock and ice-air interfacesinterface; Underworld addresses these issues with the so-called material point method (MPM) (Sulsky et al., 1994;Moresi et al., 2003), which is closely related to the venerable particle-in-cell (PIC) method. MPM uses a Eulerian Finite Element mesh in order to calculate the incremental development of the velocity field and other field variables, such as, for example, temperature and pressure. In the MPM method, Lagrangian material points ('particles') carry the density, viscosity, thermal conductivity, and other relevant material parameters. They thereby record the history at their current location at every time step and some historical properties like the stress at previous time step for simulations of viscoelastic deformations (Farrington et al., 2014). The underlying mesh provides solutions for the incremental movement of the material points. The method is advantageous has large advantages in the modeling of the emergence of structures (e. g. folding, see Mühlhaus et al., 2002) or where very strong deformation is involved, as in the deformation across shear zones or near the base of an ice sheet. In MPM, the mesh does not carry any history information other than deformation of the boundary and therefore can be re-meshed at any time as required and without loss of accuracy.
There is an unavoidable smoothing which comes from the coarseness of the computational mesh relative to the material point density (Moresi et al., 2003). While material boundaries are represented by a continuous interpolant on the grid, they are necessarily discrete in case of particles. This can lead to fluctuations in the solution close to sharp rheological or mechanical boundaries (Yang et al., 2020), for instance at the interface between ice and underlying rock. In the ice itself, a change in mechanical parameters is usually more gradual and is controlled, for example, by the temperature gradient.
Another complication in the numerical modeling of ice flow is the highly anisotropic behavior of ice, created by the nearorthotropic properties of the ice crystal. The possibility to model anisotropic flow is built into Underworld (Mühlhaus et al., 2004). Like any other local material property, the orientation of the anisotropy can be stored on the particle level.
Underworld offers a variety of possible solvers, including the well-known MUMPS, LU and multi-grid methods. We have carried out brief comparative precision tests with these solvers and could not find any difference in precision to the standard solver which is based on the multi-grid method. Throughout this study we will generally use MUMPS for 2D models and multi-grid for 3D models. The exception are tests of the computation time, for which we contrast a variety of solvers.

Description of experiments
The ISMIP-HOM benchmark experiments have been described in detail in Pattyn and Payne (2006) and . We perform all the experiments as described in these publications. For simplicity we also apply the same alphabetical numbering scheme and refer to the experiments as Experiment 'A' to Experiment 'F'. As we strive to focus on the essentials in the following descriptions, we refer to the originalindicated publications for further technical details if needed. Experiments 'A', 'B', 'C' and 'F' are three-dimensional. Experiments 'D' and 'E' consider only two spatial dimensions, and Experiment 'B' has an additional version in 2D. Experiments 'A' through 'D' are performed for a variety of horizontal system dimensions L, with L = 5, 10, 20, 40 80 or 160 km. Experiment 'A' through 'E' use a flow law based on n = 3, experiment 'F' applies Newtonian flow where n = 1 (Table 1).
We tested the influence of the mesh geometry on the results and the CPU time consumption as a function of the total degree of freedom using Experiment 'D' and the 2D-version of Experiment 'B'.

Experiments A and B
A and B consider a slab of ice with a mean ice thickness H = 1000 m, lying on a sloping bed with a mean slope α = 0.5°. The bedrock topography consists of a series of sinusoidal bumps (Experiment A) or ripples (Experiment B) with an amplitude of 500 m (Fig. 1). The minimum thickness of the rock layer is 500 m, and the total height of the model is 2000 m. The flow of ice is governed by Eq.
(3). The bedrock viscosity is constant, and ice is frozen to the bedrock. Relevant material parameters are compiled in Table 1.
The surface elevation is described by the formula: bedrock topography for Experiment A is described using zs by: and for Experiment B by:

Experiments C and D
Experiments C and D are similar to Experiments A and B, although the topography of the bedrock is flat. Instead, the coefficient of basal friction β 2 varies in a sinusoidal manner. The ice thickness H is constant at 1000 m. The slope angle α of the ice surface and of the underlying bedrock is 0.1°. Ice flow is governed by Eq. 3 and the material parameters are summarized in Table 1. According to the benchmark specifications, Experiment C is run exclusively in 3D, and Experiment D exclusively in 2D.
The basal friction coefficient relates the basal drag τ b to the basal velocity v b by With ω=2 π /L the coefficient of basal friction β 2 for Experiment C is defined as: (10) Figure 2 shows the basal drag and the basal velocity calculated from Eq. (10) (Experiment D) and Eq. (8). Here, the velocity is calculated for a constant basal shear stress τ b =ρ ice g H sin (α ), according to the shallow-ice approximation (SIA) (Hutter, 1983). Notice the singularity in the velocity field ( Fig. 2b), which develops because β 2 =0 at x=¾ L.

Experiment E
Experiment E is a two-dimensional diagnostic experiment along the central flowline of a glacier in the European Alps (Haut Glacier d'Arolla). The basic experiment and geometry are described in Blatter et al. (1998) and by Pattyn et al. (2002). The general geometry of the glacier profile as used in the experiment is shown in Fig. 3. The experiment is run with two different basal conditions: 1) the ice is everywhere frozen to the ground (β 2 = ∞), or 2) a zone of zero traction (β 2 = 0) between x = 2200m and x = 2500m exists. Compare Eq. (8) for the meaning of the basal friction coefficient β 2 .

Experiment F
Experiment F is a prognostic experiment in which a free ice surface relaxes until a steady state is reached for a zero surface mass balance. The slab of ice is resting on a bed with a mean slope α = 3°. The bedrock plane parallels the surface, but is perturbed by a Gaussian bump. The initial bedrock (B 0 ) and surface (S 0 ) elevation are described by: Experiment 'F' applies a Newtonian (n = 1) flow law, given in Table 1, so that the effective viscosity ηeff = (2 An=1) -1 . The experiment is run with two different values for the slip ratio c, so that c = 1 and c = 0 is applied. c is used to describe the basal friction coefficient β 2 (see Eq. (8)) by (13)

The FE mesh
The model domain is discretized by quadrilateral Q1/P0 elements, where velocity is continuously linear and pressure is discontinuously constant. The pressure grid is offset from the velocity grid. A direct comparison of the pressure at fixed locations to the results of other models is therefore prone to some interpolation error. Periodicity of the in-and outflow boundaries is applied in all experiments except experiment 'E'.
The ice-rock interface is defined by either of the 2 following methods, depending on the type of experiment: i) by particles, or ii) directly by the mesh geometry.

FE grids with a rectangular hull
In most experiments we assume that the bedrock is identical to the basal grid boundary. In experiments with a flat bedrock topography ('C', 'D' and 'F') this means that the resulting shape of the mesh is rectangular.
An exception is the 2D-version of Experiment 'B': in order to evaluate the impact of the mesh geometry and of a particle representation of the bedrock material, we define the sinusoidal bedrock topography using particles on a FE grid with a rectangular hull. Different materials are represented by particles with different rheological properties. Particles are assigned to either ice or bedrock depending whether their depth exceeds the local ice thickness H. As pointed out, the material point method can lead to spurious fluctuations in gradients at material boundaries, if particles of different materials are both located in one element.
In order to test and improve this behavior we define three different internal grid geometries and compare the smoothness of the resulting basal shear stress τb. These geometries are (i) a classic rectangular grid, (ii) a grid, where the grid resolution increases in the vicinity of the ice-bedrock interface, and (iii) a grid, where the mesh perfectly fits the rock surface ( Fig. 4).

Rectangular grid geometry
The benchmark assumes a constant height of the model for each experiment, while the width is varied in a series of runs during such an experiment. It is typically expected that the accuracy of the solution and the computation time are most optimal if the aspect ratio of cells is close to one.

Structurally conforming grid
We adapt the rectangular mesh to the underlying topography defined by the ice-rock boundary by vertically shifting its vertices, so that the model resolution is significantly increased close to the bedrock surface.
We assume here that the grid in y-direction originates at Y0 = 0 and ends in Y1 = 1. The rock surface is defined by a function s(x). Δy = s(p1) -p2 is the vertical distance between the rock surface and p2. Both, resolution as well as the geometry of individual cells, are adapted to the interface line as shown in Fig. 3b. m is a structural adaption parameter, controlling the intensity of the mesh deformation. It is worth pointing out that the resolution in x is not affected by this geometry.

Grid fitted to ice-rock interface
The mesh defined by Eq. (10) does not guarantee that the ice-rock interface is aligned perfectly with finite element edges. This may still introduce stress perturbations in elements containing different materials with strong viscosity contrasts (Yang et al., 2020). Therefore it makes sense to apply another mesh structure whose mesh edges fit the ice-rock interface exactly.
We define the new vertical y-coordinate p ' 2 of a vertex point p = [p1, p2] from the regular mesh by The variable n denotes the nth node in the vertical y-direction, and n0 is a predefined node, which is relocated exactly to the rock surface.n t is the total number of vertical vertexes. The difference between Eq. (11) and (10) is that we fix n 0 in Eq. (11) while n 0 varies along x-direction for the case discussed in Eq. (10). As before, m is an adaption parameter which controls the intensity of the mesh deformation. The adaption of the grid geometry does not affect the position of nodes in x-direction.

FE grids with a non-rectangular hull
In all other experiments with an uneven bedrock topography (experiments 'A' and 'E' and the 3D-version of experiment 'B') we apply Eq. (11) with m = 0 to the lower system boundary. In case of experiment 'D', the interface ice-air is not flat, therefore particles represent ice and overlying air (Fig. 5).
8 220 225 Figure 5: Grid geometry and particle distribution used for experiment 'E', representing the glacier 'Haut Glacier d'Arolla'. Particles carry the rheological properties of ice and air. Blue: ice, red: air.

Basal conditions
'Underworld' pre-implements no-slip and free-slip boundary conditions for flat system boundaries. However, experiments 'C' through 'E' require the usage of a friction law. We realize the basal drag requirement by a basal layer with Newtonian viscosity, whose viscosity is dependent on the basal friction coefficient β 2 .
Following the relation between the Newtonian viscosity of this basal layer (η b ) and β 2 is derived. For the sake of simplicity, we assume a flat horizontal surface, so that the usual notation can be used. In case of an uneven base x and y correspond to directions parallel or perpendicular to the lower system boundary.
Combining and integrating the relations τ =2 η bε (Newtonian flow) and ε x =∂ v x /∂ y (definition of the strain rate) leads to η b ( x )=h ⋅ β 2 (x ), with h the layer height. This relation is then used to define the local viscosity of the Newtonian layer.
At the top of this layer, the velocity condition defined in Eq. (8) is satisfied.

General performance tests
Below, we will first examine the CPU time consumption as a function of the grid resolution and the effect of the grid geometry on the smoothness of results.

CPU time consumption
The mesh resolution is one of the most important parameters that controls the precision of the solution. Since computation time increases with resolution, it is desirable to establish a relationship between these two quantities. Below we test and display the CPU time of the initial solution of the 2D-version of experiment 'B'. The computation time is not directly linked to the mesh resolution, but instead to the degrees of freedom (DOF) of the system. For 2D-Stokes problems with quadrilateral elements in 2D the degrees of freedom are 3N with N the total number of nodes . Figure 6 shows the relationship between the DOF and the computation time in a log-log plot for a series of 13 simulations with different resolutions and for different solvers (MUMPS, LU and multi-grid). The ratio of width to height of the grid cells remains constant in all experiments. On the hardware side, the simulations ran on a Symmetrical Multiprocessor (SMP) system with an Intel Xeon E5 processor and 32 GB RAM.
The interpolation between N and the computation time (s) expressed by a power regression is 0.00014 N 1.21 (multi-grid), 0.00037 N 1.06 (MUMPS) and 0.00013 N 1.15 (LU). In case of the multi-grid method outliers from the generally linear relation exist, which are related to the recursive refinement of the grid. Outliers do not exist for the direct solvers MUMPS and LU.
It is important to note, that -theoretically -the direct solvers should scale not better than N 2 . Therefore, the power regression cannotcan not be used in order to extrapolate these results to an arbitrary DOF. It can be seen in Fig. 6 that a slight upward curvature of the data with regard to the interpolation exists.

Impact of the grid geometry
We tested the impact of the grid geometry based using the flow law parameters n=3 and A=10 -16 Pa -3 a -1 . In order to increase effects related to the geometry, the mesh resolution is intentionally relatively small with 128x64. Only 9 240 245 250 255 260 265 for the rectangular grid, a second, larger resolution is applied. The non-rectangular meshes use an adaption parameter m=0.25.
The results for the surface velocity are generally smooth and virtually identical in the three grid geometries. However, fluctuations of the shear stress can become considerable, especially close to the ice-bedrock interface ( Fig. 7), depending on the type of FE grid. These fluctuations are largest if a low-resolution rectangular mesh is used. Increasing the resolution of the rectangular mesh does not fully eliminate the fluctuations, but is a way to reduce the problem ( as highlighted incp Fig.  7c). A more conforming grid with an increased resolution around the rock surface (Eq. 8), reduces the fluctuations significantly, but does not completely eliminate them (Fig. 7b). The unambiguously smoothest results are achieved using a grid fully fitted to the rock surface after Eq. (10 (see Fig. 7a). This confirms that fluctuations become smaller when there is less mixing of materials with different properties within an element.
The absence of fluctuations in the stress field is a measure for the ability of the model to deal with discontinuous material boundaries. In the following simulations for Experiment B, we will therefore exclusively apply the fitted mesh. This is in line with the original benchmark paper by , who suggested to test the models with optimized settings.
The boundary between ice and the bedrock is usually the only non-planar material boundary in in glaciers and ice sheets. The other material boundary, between ice and air, can often be considered almost planar in the modeling of large ice sheets, while rheological changes within the ice itself can usually be treated as gradual, controlled for instance by the temperature field or the crystallographic fabric. It is an important conclusion from these results that a FE mesh fitted (even approximately) to the underlying bedrock topography can significantly improve the accuracy of ice flow simulations based on the material point method.

Impact of the grid resolution
Choosing a well suited mesh resolution is always a balance between precision requirements and computational limitations. The computation limits can be related either to the hardware or to general time constraints, which do not allow very long computation times. Also mesh geometry has an impact on the accuracy and may or may not require a larger resolution.
We performed the 2D experiments on a classical SMP workstation (Sec. 5.1). Given the existing time constraints and the number of simulations this allowed us to test resolutions of 64 x 32, 128 x 64 up to 256 x 128. Even at the relatively small resolution of 64 x 32 the solution for stress and velocity at the surface of the model are smooth and virtually identical to the results obtained with higher resolutions, independent of the mesh geometry (Sec. 5.2). This was different in case of of the interface between the ice and the underlying bedrock. While performance was well enough or very good with a resolution of 128 x 64 for grids adapted to the interface, a resolution of 256 x 128 was necessary to obtain usable results in case of a rectangular mesh. Two grid points above the interface the solution was smooth for all mesh geometries at 128 x 64 and didn't differ from higher resolutions.
Highest resolution 3D simulations ran mostly on a high performance computer cluster. Accuracy was systematically tested for Experiment B using a fitted grid, by comparing the results of the 2D and 3D experiments. Under perfect conditions 2D and 3D simulations which should yield identical results. We found that the necessary resolution is highly dependent on the system size, i. e. on the aspect ratio of the cells. For small system sizes of up to 20 km we found that a minimum resolution of 256 x 128 x 32 was necessary to produce similar results as the 2D model. Larger systems produced satisfactorily results already with a resolution of 128 x 64 x 8. Below we will compare the results of the individual benchmark experiments generated by Underworld to the results of alternative codes. Where applicable we run both, the three-dimensional and the two-dimensional version of the benchmark experiments. This applies to Experiments 'B' and 'C'. For the latter, Experiment 'D' is the associated 2D-setup. All results are compiled in the supplementarycomplimentary material to this paper. A note on the output-parameter Δp, which is the difference between the isotropic and the hydrostatic pressure: the curve progression is very similar to that of other Full Stokes model of the ISMIP-project, but show stronger fluctuations between adjacent evaluation points. This is due to the internal architecture of Underworld-experiments, where the mesh used for the calculation of pressure is a sub-mesh of the velocity-mesh with a staggered geometry. Fluctuations are a product of internal interpolation.

Experiments 'A' and 'B'
We fit the lower system boundary to the topography of the underlying bedrock in the 3D experiments by applying Eq. (11) with m = 0.2 to the lower system boundary. The 2D version of experiment 'B' is based on a rectangular grid, which is internally fitted to the rock surface as in Sec. 5.2 above. Other relevant parameters are given in Table 1.
Results of experiment 'A' are shown in Fig. S1 and S2  The surface velocity is controlled by the model width L and is in the range from a few m a -1 to more than 100 m a -1 . Figure 8 compares our results for Experiment 'B' to results compiled from 8 full-Stokes models that were previously published . A full comparison of the results is in the supplement to this paper.
In Experiment 'B', the shape of the horizontal surface velocity for L=5 km differs significantly from that for the other cases. The surface velocity is larger over the bump and thus anti-correlated with the ice thickness. Gagliardini et al. (2008) explain this as a mass conservation effect: horizontal flux cannot be balanced by vertical flux, because the vertical velocity would be to large for the given system size. This phenomenon does not occur in Experiment 'A', with the same system size, although the section is identical at the chosen location. This can be explained because ice can flow around the sides of the bump. To shed some light on the difference between 2D and 3D results, we ran 3D simulations with a resolution of 128 x 64 x 8 and of 256 x 128 x 32 for these system sizes and compared them to the results of the 2D simulation (256 x 128). Regarding the surface velocity, the higher-resolution 3D mesh yields indeed results closer to the 2D result than the lower-resolution simulation ( Figure S6 of the supplementary material). This is not the case for the basal shear stress, where the deviation from the 2D result is comparable for both resolutions.

Experiments 'C' and 'D'
Parameters for Experiments 'C' and 'D' are given in Table 1 for n=3. Results are summarized in Fig. S6 and S7 (Experiment 'C') and Fig. S8 and S9 (Experiment 'D') in the supplementary material.
In Experiment 'C' the surface velocity is strongly dependent on the system size L. In case of the smallest system size (L=5 km) the surface velocity is close to constant at ~16 m a -1 . With L=160 km it ranges from 8.8 m a -1 up to 122.5 m a -1 (Fig. S7, supplement). The shear stress results of the Underworld-simulations line up well with the comparison simulations. They show a sinusoidal curve with a maximum at x =0.25 and a minimum at x =0.75. With increasing model size, the maximum gets progressively flattened. The peak downwards at the singularity value x =0.75 is far less impacted by model size, which means it gets more dominant for bigger values of L.
In Experiment 'D', the maximum surface velocity increases with model width L and is in the range from 16 m a -1 to more than 235 m a -1 . Figure 9 showsFig.9shows the maximum and minimum surface velocity and the maximum shear stress τxy.
Results are generally in good agreement with the majority of model results they are compared to. However, the general variation between the results of the comparison models is generally larger than in the 2D version of Experiment 'B', the only other 2D experiment of the benchmark. The same statement applies to basal shear stress. Again, a full comparison of the results is included in the supplement to this paper.

Experiments 'E' and 'F'
Experiment 'E' is calculated for two different setups, where either the entire glacier is frozen to the underlying bedrock, or where a traction-free section close to the center exists. In both cases, the results of Underworld lies within the range of the comparison models (see Fig. S10 and S11 in the supplementary material). Both the calculated surface velocity and the basal shear stress are in the lower range of this range, which is true for the majority of models. In particular extreme values of the curve are a bit smoother than in some of the other Full Stokes models.
The prognostic Experiment 'F' calculates the surface velocity and the surface topography. Results are compiled in Fig. S12 and S13 in the supplementary material. Since the majority of potential comparison models are not capable to calculate a free surface, the results of Underworld can be compared to only two other software packages. Due to the special properties of the software, we model only the version of the experiment with a basal no-slip condition. Other values for the basal traction would either involve a very complex implementation or yield questionable results. We ran the model until no further change in topography or velocity occurred, and interpreted this as the stable state.
When we compare our results with those of the two reference models with respect to an analytical sample solution (Frank Pattyn, 2022, personal communication), qualitatively similar results are obtained. One of the reference models shows very good agreement with the theoretical surface elevation, but shows lower accuracy with respect to surface velocity. Underworld predicts the surface velocity the best of the three models, but tends to develop a slightly less extreme topography than the analytical solution.

Comparison experiments with anisotropic and isotropic ice
Ice is assumed to deform mostly by dislocation creep at natural strain rates (Budd and Jacka, 1989), whereby slip along the basal planes is much easier than slip along the other slip planes (Duval et al., 1983). This makes an ice single crystal mechanically effectively transversely isotropic, with the c-axis that is oriented perpendicular to the basal plane as the symmetry axis. The anisotropy of ice single crystals leads to a crystallographic preferred orientation (CPO) in a deforming aggregate of ice crystals (Alley, 1988;Llorens et al., 2017). The mechanical anisotropy of an aggregate may be of a lower symmetry than that of a single crystal but can be approximated to a first order as orthotropic (Gillet-Chaulet et al., 2006).
One of the features of Underworld, which makes it interesting for the modeling of ice, is its capability to model linear orthotropic viscosity, which includes transverse isotropy as a special case.. In order to demonstrate its effect we set up comparison experiments with anisotropic and isotropic ice rheology, based on the setup of Experiment 'B': ice flows over a sinusoidal surface, driven by a general 0.5° tilt of the model. Isotropic flow is governed by Eq. 4, using the parameters given in Table 1.
The orientation of the anisotropy is stored as the c-axis orientation on the level of the particles in Underworld. The aggregate anisotropy of one mesh element is calculated from the individual c-axis orientations of the cloud of particles in an element. Fabric evolution is simulated using the rotation of the c-axes in the flow field after e. g. Gillet-Chaulet et al. (2005) or Richards et. al. (2021). However, for simplicity, in our example here we assume that all basal planes are and remain horizontal. The anisotropy of ice, in terms of the ratio of resistance to slip parallel to crystallographic non-basal and basal planes, is about 60-80 in a single ice crystal (Duval et al., 1983). However, as we here simulate aggregates of crystals, we set the viscosity for shear non-parallel to the basal planes only ten times higher than the vicosity for shear parallel to the basal planes, which we assume equal to that used for the isotropic flow law. Figure 9 shows the shape of marker lines prior to deformation and after 750 a of flow for both, iso-and anisotropic ice models. Marker lines inherit their sinusoidal shape from the shape of the underlying topography and are then folded according to the localization of the highest shear rates. In isotropic ice (Fig. 9 b) the axial plane of the shear fold is mostly controlled by the underlying topography. In case of anisotropic ice (Fig. 9 c), the folding is more intense and the axial plane is sub-horizontal, showing that the anisotropy has a distinct effect on resulting fold structures. Figures 10 and 11 that show the velocity field and the strain rate field in both materials allow a better understanding of the deformation process. In case of isotropic ice, the flow field is controlled by the underlying topography. Hence the hill on the left acts as a bottleneck for the ice flow, and the hill sides funnel ice in and out of the bottleneck region. Consequently, the zone of high shear strain at the ice-rock interface extends lang the crest of the bedrock bump (Fig. 10 a) and the maximum velocity above it. Outside the bottle-neck region flow is comparatively evenly distributed.
In the case of anisotropic ice, the flow regime and thus the shear folding is quite different. Here, flow is strongly controlled by the inherent anisotropy and far less by the bedrock topography. Looking at the velocity field (Figure 11 a), it becomes clear that the flow field is internally subdivided into a fast flowing upper part and slow flowing 'dead ice' in the lower part. Decoupling of the shallow and deep ice develops because the anisotropy facilitates shear along the horizontal basal planes. The resulting horizontal high strain zone spans the entire system and produces the distinct shear fold.

Experiments 'E' and 'F'
Experiment 'E' is calculated for two different setups, where either the entire glacier is frozen to the underlying bedrock, or where a traction-free section close to the center exists. In both cases, the results of Underworld lies within the range of the comparison models (see Fig. S10 and S11 in the supplementary material). Both the calculated surface velocity and the basal shear stress are in the lower range of this range, which is true for the majority of models. In particular extreme values of the curve are a bit smoother than in some of the other Full Stokes models.
The prognostic Experiment 'F' calculates the surface velocity and the surface topography. Results are compiled in Fig. S12 and S13 in the supplementary material. Since the majority of potential comparison models are not capable to calculate a free surface, the results of Underworld can be compared to only two other software packages. Due to the special properties of the software, we model only the version of the experiment with a basal no-slip condition. Other values for the basal traction would either involve a very complex implementation or result in questionable results. We ran the model until no further change in topography or velocity occurred, and interpreted this as the stable state. In comparison to the two reference models, very similar results are achieved. The topography and the surface speed have a similar general trend. Differences between Underworld and the reference models are similar to those between the two reference models itself. In case of surface topography, however, Underworld tends to develop a slightly less extreme topography.

Outlook
One of the great advantages of the Underword2 software package from the standpoint of the modeler are its great flexibility and traceability due to its hybrid nature as both, particle and Finite Element model. Another contribution to its flexibility is the easy extensibility of the core code and the interactive development thanks to the existence of a Python API. The compatibility of the API with the mathematical-scientific packages NumPy and SciPy provides easy access to a wide range of numerical resources. Anisotropic flow and heat flow are already part of the core package. Given the current debates in the ice community about the role of the anisotropy and hence of the fabric evolution of ice, or the ideas for improvements of the flow laws for ice (e.g. Kennedy et al, 2013;Llorens et al., 2017;Richards et al., 2021;Martín et al., 2009), this type of flexibility is an important precondition for a future proof numerical model for glacier ice. One example was the implementation of the orthotropic rheology and the fabric evolution in deformed ice following Gödert (2003) and Gillet-Chaulet et al. (2005) by the authors of this text. Possible improvements to the fabric evolution models -for instance by the SpecCAF model (Richards et. al, 2021) -can be easily implemented and tested.
A few desirable features are currently still missing, but are on the road map for future versions of Underworld. The soon-tobe released Underworld3 for instance allows greater flexibility of the mesh geometry, including the triangulation of shapes and areas with arbitrary geometry.

Conclusions
f[e] The software 'Underworld is well suited for solving deformation in complex geodynamic systems with non-linear elasto-visco-plastic materials, for which it provides a full-Stokes solution.
The Underworld software is designed to solve deformation in complex geodynamic systems with non-linear elasto-viscoplastic materials, for which it provides a full-Stokes solution. It is therefore Underworld is well suited for the modeling of glacier and ice-sheet flow, as it includes heat flow and anisotropic rheology. The combination of Lagrangian mass points (particles) and a Eulerian Finite Element solution allows thematerial point method allows accurate tracking of individual points as well as of inner and outer surfaces, such as deforming stratigraphic layers, but also of the thermal-mechanical properties in deforming materials. In case of large rheological differences across interfaces, the possibility to fit the grid to the interface greatly improves the accuracy of stress field, compared to other grid types. In case of ice flow experiments, it makes sense to fit the grid to the bedrock-ice boundary, which is advantageous for the tracking of, for example, deforming stratigraphic layers.
We compared results of Underworld simulations with those of other modeling approaches for the set of benchmark experiments provided by . Our results match the full-Stokes solutions that are compiled in that study. This means that Underworld is a viable alternative to other full-Stokes models, in particular where the Material Point Method is advantageous, such as when accurate tracking of material volumes or stratigraphic layers is desired. A further advantage is that, owing to the built-in Python API, Underworld is very flexible and can be extended to be applied to even the more complex processes which are involved in the flow of ice sheets and glaciers.
In case of large rheological differences across interfaces, the possibility to fit the grid to the interface greatly improves the accuracy of stress field, compared to other grid types. In case of ice flow experiments, it makes sense to fit the grid to the rock-ice boundary. Code/Data availability · The program code that defines the experiments discussed in this article is included in the supplement to this article.
· Graphical representations of the results as well as data files are found in the supplementary material to this article.
· In addition to the supplement, the complete data set is also available here: https://gitlab.com/underworld2/ISMIP-HOM-Data. Data files show output for all grid points in the horizontal for the domain from 0 to the model width L (or between 0 and 1 for scaled coordinates). All variables are taken either at the surface ys or at the basis of the ice at yb. The abbreviations for the study presented here and the comparison models are named accordingly (see Table 1). Because experimental data have been compiled by different authors, the abbreviations for this study are 'tsa1' (2D experiments of 'B' and 'D') and 'jla1' (all other experiments).

Data files for experiment A
For Experiment A the following output has been compiled in each file in that order: the normalized x position, the normalized y position, the horizontal velocity in x direction at the surface, the horizontal velocity in y direction at the surface, the vertical velocity at the surface, the basal shear stress in x direction, the basal shear stress in y direction and the difference between the isotropic and the hydrostatic stress at the basis of the ice. xẑ

Data files for Experiment B
For Experiment B the following output has been compiled in each file in that order: the normalized x position, the horizontal velocity at the surface, the vertical velocity at the surface, the basal shear stress and the difference between the isotropic and the hydrostatic stress at the basis of the ice.
For 3D experiments, values are given for a profile at ẑ=0.25: x v x ( y s ) v y ( y s ) τ xy ( y b )

Δ p
For 2D experiments: x v x ( y s ) v y ( y s ) τ xy ( y b ) Δ p

Data files for Experiment C
For Experiment C the following output has been compiled in each file in that order: the normalized x position, the normalized y position, the horizontal velocity in x direction at the surface, the horizontal velocity in y direction at the surface, the vertical velocity at the surface, the horizontal velocity in x direction at the base, the horizontal velocity in y direction at the surface, the basal shear stress in x direction, the basal shear stress in y direction and the difference between the isotropic and the hydrostatic stress at the basis of the ice. xẑ

Data files for Experiment D
For Experiment D the following output has been compiled in each file and in the given order: the normalized x position, the horizontal velocity at the surface, the vertical velocity at the surface, the horizontal velocity at the basis, the basal shear stress and the difference between the isotropic and the hydrostatic stress at the basis of the ice. x

Data files for Experiment E
For Experiment B the following output has been compiled in each file and in the given order: The normalized x position, the horizontal velocity at the surface, the vertical velocity at the surface, the basal shear stress and the difference between the isotropic and the hydrostatic stress at the basis of the ice.
x v x ( y s ) v y ( y s ) τ xy ( y b ) Δ p

Data files for Experiment F
For Experiment A the following output has been compiled in each file: The normalized x position, the normalized z position, the y position at the surface, the velocity in x direction, the velocity in y direction and the vertical velocity. xẑ y s v x ( y s ) v y ( y s ) v z ( y s ) 2 Explanations to the diagrams

Experiment A
This is a 3D ice-stream experiment of flow over a bumpy bedrock surface. The presented results have been calculated on an adapted grid (compare text) with a mesh resolution of 128 x 64 x 8. Compare text for further explanations and parameters.
The surface velocity is given in m/a, shear stress is given in kPa. Data is plotted in Figures S1 and S2. Each model size is shown separately and compared to previous results from full-Stokes solutions compiled in  (see Table S2).

Experiment B
This is a 3D and 2D ice-stream experiment of flow over a rippled bedrock surface. The presented results have been calculated on an adapted grid (compare text) with a mesh resolution of 128 x 64 x 8 (jla1) and 256 x 64 x 32 (tsa2) for 3D experiments. The resolution in 2D is 256 x 128. Compare main text for further explanations and parameters. The horizontal surface velocity is given in m/a, shear stress is given in kPa. Data is plotted in Figures S3 to S6 for 3D and 2D settings. Each model size is shown separately and compared to previous results from full-Stokes solutions compiled in  in case of the 3D experiments (see Table S2). Diagrams with results of the 3D experiments compare them to the results of the 2D setup.

Experiment C
This is a 3D experiment. The bedrock topography is flat, but the basal friction coefficient is prescribed by a sinusoidal function. The mesh resolution is 128 x 64 x 8. For further explanations see main text. Velocities are given in m/a, shear stress and pressure is given in kPa. Data is plotted in Figures S5 and S6. Each model size is shown separately and compared to previous results from full-Stokes solutions compiled in  (see Table S2).

Experiment D
This is a 2D ice-stream experiment. The results presented here have been calculated on an rectangular grid with a mesh resolution of 375 x 75. Parameters and explanations are given in the main text.
The horizontal velocity at the surface and at the basis of the ice sheet is given in m/a, shear stress in x-direction is given in kPa. Data is plotted in Figures S7, S8 and S9. Each model size is shown separately and compared to previous results from full-Stokes solutions compiled in  (see Table S2).

Experiment E
Exp. E is an experiment along the central flowline of a temperate glacier in the European Alps and thus essentially a 2D experiment. The model input consists of the longitudinal surface and bedrock profiles of Haut Glacier d'Arolla, Switzerland (Blatter et al., 1998 andPattyn, 2002). In a first experiment the ice is frozen to the ground. A second variant considers a narrow zone of zero traction close to the center. The mesh geometry has been adapted to the bedrock topography, the resolution is 50 x 50. Compare main text for parameters and further explanations. The horizontal velocity at the surface and at the basis of the ice sheet is given in m/a, shear stress in is given in kPa.
Data is plotted in Figures S11 and S12. Each model size is shown separately and compared to previous results from full-Stokes solutions compiled in  (see Table S2).

Experiment F
Experiment F is a 3D ice-stream experiment over a central Gaussian bump. The free surface is allowed to relax until a steady state is reached. The mesh resolution is 240 x 48 x 240, and the mesh geometry has been adapted to the bedrock topography. Experiments are run with a slip ratio c=0, which means that ice is effectively frozen to the ground. This experiment uses Newtonian viscosity, in difference to previous experiments. Compare main text for parameters and further explanations. Figures S12 and S13 compare the results to two full-Stokes solutions compiled in  and an analytical solution, contributed by Frank Pattyn (2022, personal communication). Surface height is given in m, velocity in m/a.

Comparison models
The data plots below compare the results of this study with previous full-Stokes solutions compiled in the reference paper by . Pattyn (2022, personal communication).