2.1 Implementation in Elmer/Earth

Elmer/Earth is based on the open-source finite-element package Elmer (Råback et al., 2019). In order to build a flat-earth model as described in the previous section, Eq. (1) has been added to the existing linear elasticity solver of Elmer. In the case of incompressibility, the additional variable of pressure, Π, has been introduced to the solver. This avoids the singularity of the compressible formulation in the case of the Poisson ratio approaching 0.5.

Many commercial codes lack an implementation of the second term in Eq. (5), which implies a transformation of the stress to reduce the formulation to only the first term. As a consequence of this stress transformation, additional jump conditions in the form of Winkler foundations (Wu, 2004) have to be imposed on internal boundaries that mark a jump in either the gravity or the density. This can be inconvenient in building the model, as a detailed description of the setup may contain boundaries for more than 10 layers.

Here we take advantage of the accessibility of the source code of Elmer by including this term in the weak formulation that uses the viscoelastic stress. The second term in Eq. (5) thereby contributes to the stiffness matrix. Naturally, the formulation still needs a layered structure of the model, i.e., material parameters are kept constant for certain layers. This can be easily achieved as Elmer allows material parameters to be prescribed as well as body forces (in our case gravity), on the basis of elements or even integration points (in addition to nodal values). This means that we are able to impose discontinuities in parameters over elements anywhere in the discretized computing domain without placing Winkler foundation boundaries at layer interfaces. In other words, no boundary conditions have to be set at internal layer boundaries. By including this term in the weak formulation of the problem, the method then automatically applies the restoring force on element boundaries with jumps in material properties or gravity, without the need to place boundaries in the mesh.

Discretization of the time derivatives for stress and pressure (in the case of incompressible material) is implemented by the first-order implicit difference

Here, i is the current, and i+1 the implicit time step as well as $\mathrm{\Delta }t=t^{i+\mathrm{1}}-t^i$ the time-step size between. The solution of the time-evolution problem then reads

with $\mathit{\varphi }=\mathrm{1}/(\mathrm{1}+(\mathit{\mu }/\mathit{\nu }\left)\mathrm{\Delta }t\right)$. The balance Eq. (5) of linear momentum is then solved for the new time step:

The weak formulation then results from the integral over the whole domain Ω (with its confining surface ∂Ω) using the test and weighting function vectors $\mathit{u},\mathit{v}\in H^{\mathrm{1}}$:

Note that the divergence of the stress tensor has been partially integrated, leading – according to Green's theorem – to a term that integrates the stress vector, $\mathit{t}=\mathit{\tau }\left(\mathit{u}\right)\cdot \mathit{n}$, over ∂Ω with its surface normal n. Taking additionally into account that τ(u) is a symmetric tensor, only the symmetric part of sym(∇v)=ϵ(v) contributes to the first integral, leading to the symmetric stiffness matrix in the weak formulation

The system is completed by boundary conditions that are either provided by a value for any component of the stress vector, $\mathit{t}=\mathit{\tau }\cdot \mathit{n}$, in the second integral (Neumann condition) of Eq. (10) or by imposing a value for any component of the deformation vector, d (Dirichlet condition).

Equation (10) is solved using the standard Galerkin method with first-order basis functions in the case of the benchmark described in Sect. 3. Apart from this particular choice, Elmer provides a variety of possible basis functions left to the choice of the user. The iteration for the viscous contribution is computed on the Gaussian integration points. In the case of incompressibility, stabilization has to be applied by the residual free bubble method.

3.1 Reference model ABAQUS

We use the finite-element software package ABAQUS (Hibbitt et al., 2016; software version 2018) to construct a model to verify the results of the new viscoelastic solver implemented in Elmer. We choose this approach to replicate the geometry and equations implemented in the Elmer/Earth model as fully as possible. The model is a 3-D flat-earth model which computes the solid earth deformation in response to a changing surface load using the approach of Wu (2004). Buoyancy forces are accounted for by applying Winkler foundations to layer boundaries within the model where a density contrast occurs between two layers, and at the surface (Wu, 2004). The model has a large lateral extent to prevent boundary effects in the area of interest (Steffen et al., 2006) and has zero displacement imposed on its lateral and bottom boundaries. The model includes layers from the surface of the earth to the core–mantle boundary with parameters shown in Table 1.

3.2 Reference model TABOO

TABOO is an open-source post-glacial rebound calculator (Spada et al., 2003; Spada, 2003) that computes the deformation of the earth in response to a changing surface (glacial) load. The TABOO model assumes a spherically symmetric, incompressible earth with a Maxwell viscoelastic rheology (non-rotating, self-gravitational). TABOO implements the classical viscoelastic normal mode method commonly used in studies of glacial isostatic adjustment (Peltier, 1974). There are several inbuilt solid earth models available in TABOO with a specific earth structure and parameters and we use one of these for our synthetic benchmarking case study (Table 1, Sect. 3.3). Deformation is computed up to a user-specified spherical harmonic degree, and we chose 2048 (equivalent to approximately 10 km).

3.3 Test model setup

In order to test and compare the newly built Elmer/Earth model, a simple benchmark case has been set up for each of the models presented in Sect. 3.1 and 3.2. The benchmark case consists of a simple one-dimensional earth structure with parameters varying in the radial direction only, loaded and unloaded with a disc of ice. The models in Elmer/Earth and ABAQUS both use a flat-earth approximation, whereas TABOO is a fully spherical model. The effects of sphericity are negligible for the size of load we use for our benchmarking case. None of the models solve the “sea-level equation” (Farrell and Clark, 1976).

For the flat-earth approximation, the three-dimensional model domain stretches 4000 km in each horizontal direction from the centre of the ice load. This distance is 80 times the diameter of the test load, which is more than sufficient to allow mantle deformation below the ice load (Steffen et al., 2006). With depth, the model extends from the earth’s surface at a radius of 6371 km to the core–mantle boundary with a total depth of 2891 km.

Geometry construction and meshing for Elmer/Earth simulations was achieved using the open-source software Gmsh (Geuzaine and Remacle, 2009). The lateral mesh resolution for the ABAQUS model is a constant 10 km, whereas it varies for Elmer/Earth from 10 km for the area over which the load is applied to 200 km, increasing linearly, at the lateral domain boundaries (see Fig. 1a). The vertical resolution increases with depth as shown in Fig. 1b. The TABOO model has a spectral resolution equivalent to 10 km.

Figure 1Top and side view of the reference run Elmer/Earth mesh (mesh1). The different layers corresponding to varying material parameters shown in panel (b) are given in Table 1. Annotated coordinates are in kilometres.

Download

The earth structure used for the benchmarking case is one that is included as part of the TABOO package and is summarized in Table 1. The solid earth model consists of an elastic lithosphere, a viscoelastic upper mantle divided into three layers, and a viscoelastic lower mantle. Elmer/Earth applies incompressibility throughout the whole column and an extremely high viscosity of $\mathit{\nu }=\mathrm{1}\times {\mathrm{10}}^{\mathrm{44}}$ Pa s in the lithosphere, thereby enforcing an approximately elastic behaviour on the timescale of the load. This can be justified by the Maxwell time $t_{\mathrm{M}}=\mathit{\nu }/\mathit{\mu }$ being of the order of 10³³ s, which indicates that viscous effects would only be significant at timescales several order of magnitudes larger than the timing of the load signal.

The viscosity of the upper and lower mantle is set to 1×10¹⁸ and 1×10²² Pa s, respectively, and the elastic and density parameters are depth-averaged values from the Preliminary Reference Earth Model (Dziewonski and Anderson, 1981; PREM). These parameters can easily be assigned to layers in both ABAQUS and Elmer. The relatively low value for the upper mantle helps to shorten the timescales for the benchmark test.

For the benchmark case we compute the deformation caused by an instantaneously imposed ice load at t=0. Starting from an equilibrium bedrock with zero deformation, an ice load is instantaneously applied at the centre of the domain at the very beginning of the simulation. It is a 100 km diameter disc of 100 m height with a prescribed constant density of 917 kg m⁻³. The load is maintained for 100 years after which it is instantaneously removed and the rebound computed for a further 100 years. The result on the vertical plane of symmetry from the reference run described in Sect. 5 is shown in Fig. 2.

Figure 2Cross section of the reference run with Elmer/Earth (mesh1) showing the vertical deformation at 99 years into the simulation at maximum deformation. Deformation is shown as colour texture as well as isoline (white in 0.1 m spacing). The boundaries between the lithosphere and upper and lower mantle (as given in Table 1) are annotated as grey lines. Annotated coordinates are in kilometres.

Download

The temporal evolution of the vertical displacement of the reference Elmer/Earth run (mesh1) over a line at the surface from the centre to the margin (0–200 km) is depicted in Fig. 3.

Figure 3Temporal evolution of vertical displacement of the reference Elmer/Earth run (mesh1) over a line at the surface from the centre to the margin.

Download

Table 1Properties of the different layers in the flat-earth model benchmark. Vertical distances are with respect to earth's centre. The ABAQUS reference model uses a material model with a constant Poisson ratio of 0.49 throughout the whole domain.

Download Print Version | Download XLSX

3.4 Numerical settings in Elmer/Earth

For all runs of Elmer/Earth presented in Sects. 4 and 5, the same numerical methods and parameters have been applied. A time-step size for the implicit backward differentiation formula (BDF) of the equivalent of 1 year has been chosen – in Sect. 5 we discuss the impact in accuracy by halving this time-step size. The resulting system matrix of the linear elasticity solver was first pre-conditioned using an ILU (incomplete lower–upper) factorization of first-order degree (ILU1, in Elmer terminology). To obtain a solution, its inverse was approximated using the GCR (generalized conjugate residual) Krylov subspace method (see, e.g., Eisenstat et al., 1983). A convergence criterion was applied for the relative norm of the solution vector between two iteration steps of ${\mathit{\epsilon }}_{\mathit{d}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{7}}$ .

5.1 Strong and weak scaling

Tests were performed on the Linux cluster raijin (Australian National Computational Infrastructure, 2017), utilizing compute nodes, each equipped with two Intel Xeon Sandy Bridge (E5-2670, 2.6 GHz) processors summing up to 16 cores per compute node. The code was compiled using the Intel compiler suite (version 2019.2.187) with Open MP (OMP) enabled, mainly to activate utilization of OMP-SIMD instructions within the code (Byckling et al., 2017). CPU-specific optimization was enabled by compiler flags -O2 -march=sandybridge. Basic linear algebra libraries (Lapack, BLAS, ScaLapack) were linked in from the Intel MKL library. Message passing was enabled by linking to the Intel MPI library (version 5.1.0.097) provided on the system.

We want to emphasize that we only studied a limited set of problem sizes or computing resource configurations, and only single runs (no statistics) were performed. Results presented in the following thus have to be interpreted in view of the limitations. All runs performed are summarized in Table 3.

Table 3 Timings of different scalability test runs. All timings are given in seconds.

Download Print Version | Download XLSX

A comparison of a simulation performed with 16 cores (single compute node) with mesh2 (half size) and with 32 cores (two compute nodes) on mesh1 (reference) reveals a drop to 64 % of an ideal, linear weak scaling (increasing core numbers while maintaining the load/core) performance. This can be explained by adding additional latency to that part of the MPI communication that in the 32-core run has to be routed over the inter-nodal connection (Infiniband), whereas the 16-core run solely uses faster communication provided within a single compute node. Reassuringly, a similar value, namely 61 %, was obtained between runs on the double-size mesh (mesh3) with 64 cores on four compute nodes in relation to the reference problem (mesh1) run on 32 cores on two compute nodes. Studying the log files of the runs, it also becomes clear that the chosen GCR algorithm takes longer to converge with respect to the same convergence criteria if increasing the amount of mesh partitions. Another comparison with slightly less strict convergence criteria of the linear solution iteration algorithm led to a value of 84 %.

On the other hand, if looking at strong scalability (i.e., increasing core numbers while reducing load/core), doubling computational resources from 16 cores (single compute node) to 32 cores (inter-nodal) for the fixed-size smaller problem (mesh2) revealed a speedup of 1.66, which is below the ideal value of 2 (half wall-clock time, by doubling of cores). For the larger reference problem (mesh1), we achieve a speedup of 2.2 if increasing from 16 (single node) to a 32 core utilizing two compute nodes of the reference run. We have not investigated the particular cause of this super-linear scaling further but can speculate on it: reducing the memory or core needed improves the possibility of fitting more data into the cache and thereby enabling faster memory access (i.e., avoiding cache misses) and hence – despite the added latency from inter-nodal communication – allowing for a general acceleration.

Despite applying the same solution method, it is not really possible to compare the performance of Elmer/Earth to ABAQUS, since the latter was run on a different platform using a regular mesh of 10 km constant horizontal mesh size. Computational performance was not the main motivation behind using ABAQUS for the benchmarking exercise; rather we wanted to use a model that could best replicate the geometry and equations used. Nevertheless, it is interesting to note that the run time of ABAQUS was in the range of 6 h using 32 cores on a high-end workstation, hence about twice the time of Elmer/Earth reference run on the same amount of cores of a larger Linux cluster. These run times should not be used in a direct comparison for computational performance, since ABAQUS was run on a mesh significantly larger (600 k nodes) than the one of Elmer/Earth. However, TABOO is using a completely different model approach, such that any comparison would be obsolete.

5.2 Accuracy with respect to mesh and time-step size

We further studied the accuracy and consistency of Elmer/Earth results with respect to spatial and temporal discretization sizes. To that end, we ran the same numerical setup on all three meshes given in Table 2.

Figure 6Vertical deformation at the centre (0 km) of Elmer/Earth simulations using different spatial and temporal resolutions.

Download

Results are depicted in Fig. 6 and reveal that too low a spatial resolution (i.e., mesh2) – in that particular case in the horizontal as well as vertical direction – yields deformations that are too large. That might simply be because of too little resolution of the induced viscous deformation in under-resolved layers. The more finely resolved meshes (mesh1 and mesh3) show very little deviations in results, thus indicating consistency of the model beyond a resolution of about 5 km mesh size at the centre of the geometry and the vertical structure depicted in Fig. 1b. On the other hand, increasing temporal accuracy by reducing the time-step size from 1 year to half a year did not reveal any significant difference in the result for similar setups to the reference run (mesh1).