Articles | Volume 13, issue 3
Geosci. Model Dev., 13, 1155–1164, 2020
Geosci. Model Dev., 13, 1155–1164, 2020

Development and technical paper 12 Mar 2020

Development and technical paper | 12 Mar 2020

A new open-source viscoelastic solid earth deformation module implemented in Elmer (v8.4)

A new open-source viscoelastic solid earth deformation module implemented in Elmer (v8.4)
Thomas Zwinger1, Grace A. Nield2,3, Juha Ruokolainen1, and Matt A. King2 Thomas Zwinger et al.
  • 1CSC–IT Center for Science Ltd., Espoo, Finland
  • 2Surveying and Spatial Sciences, School of Technology, Environments and Design, University of Tasmania, Hobart, Australia
  • 3Department of Geography, Durham University, Durham, UK

Correspondence: Thomas Zwinger (


We present a new, open-source viscoelastic solid earth deformation model, Elmer/Earth. Using the multi-physics finite-element package Elmer, a model to compute viscoelastic material deformation has been implemented into the existing linear elasticity solver routine. Unlike approaches often implemented in engineering codes, our solver accounts for the restoring force of buoyancy within a system of layers with depth-varying density. It does this by directly integrating the solution of the system rather than by applying stress-jump conditions in the form of Winkler foundations on inter-layer boundaries, as is usually needed when solving the minimization problem given by the stress divergence in commercial codes. We benchmarked the new model with results from a commercial finite-element engineering package (ABAQUS, v2018) and another open-source code that uses viscoelastic normal mode theory, TABOO, using a flat-earth setup loaded by a cylindrical disc of 100 km in diameter and 100 m in height at the density of ice. Evaluating the differences in predicted surface deformation at the centre of the load and two distinctive distances (100 and 200 km), average deviations of 7 and 2.7 cm of Elmer/Earth results to ABAQUS and TABOO, respectively, were observed. In view of more than 100 cm maximum vertical deformation and the different numerical methods and parameters, these are very encouraging results. Elmer is set up as a highly scalable parallel code and distributed under the (L)GPL license, meaning that large-scale computations can be made without any licensing restrictions. Scaling figures presented in this paper show good parallel performance of the new model. Additionally, the high-fidelity ice-sheet code Elmer/Ice utilizes the same source base as Elmer and thereby the new model opens the way to undertaking high-resolution coupled ice-flow–solid-earth deformation simulations, which are required for robust projections of future sea-level rise and glacial isostatic adjustment.

1 Introduction

Reconstructing ice-sheet history and predicting ice-sheet response to changes in climate are imperative for accurately predicting future ice-mass loss and hence sea-level rise. An important component of ice-sheet evolution is the isostatic response of the solid earth that occurs as a result of changes in the mass of the ice sheet. Over glacial cycles the waxing and waning of ice sheets causes the underlying earth to deform as the ice loading at the surface grows and shrinks. This deformation occurs both instantaneously as an elastic response and over longer timescales as the viscous mantle flows back to previously glaciated regions in order to regain gravitational equilibrium. How fast or slowly the earth deforms depends on the underlying mantle viscosity, and, although typically thought to occur over several thousands of years (Whitehouse2018, and references therein), recent studies have shown regions undergoing much more rapid (decadal) rebound in response to present-day changes (Nield et al.2014; Barletta et al.2018).

This isostatic response of the bedrock can strongly influence ice-sheet dynamics. Deformation of the earth changes the elevation of the ice sheet which in turn affects the surface temperature and the rate of accumulation or ablation. Solid earth deformation also alters the gradient of the bedrock on which the ice sheet rests, particularly at the periphery, altering the internal forces as well as the driving stress and therefore the flow of the ice sheet (Le Meur and Huybrechts1996; Adhikari et al.2014). In marine-grounded ice sheets lying on a reverse slope bed (e.g. West Antarctica) these effects can be critical. As the grounding line retreats further along the reverse slope into deeper water, ice flux across the grounding line increases leading to increased loss (Schoof2007). However, bedrock uplift can have a stabilizing effect by reducing the slope of the reverse bed and thereby slowing the retreat of the grounding line (Gomez et al.2010, 2013).

Including the isostatic response of bedrock in an ice-sheet model is therefore crucial to obtaining accurate predictions of ice-sheet mass balance, and there are several methods which can be used. Computing the isostatic response with a self-gravitating viscoelastic spherical earth is the most accurate, but most computationally expensive, method. Several simple approximations are often made using models with a combinations of local lithosphere or elastic lithosphere with diffusive asthenosphere or relaxing asthenosphere (Le Meur and Huybrechts1996; Rutt et al.2009). Of these, Le Meur and Huybrechts (1996) found the best performing is the “ELRA” (e.g. Greve2001) model (elastic lithosphere with relaxing asthenosphere) which is widely used in ice-sheet modelling, mainly due to its simplicity and fast computations. However, Bueler et al. (2007) found significant differences in resulting bed elevation and ice-sheet thickness when using a model with ELRA compared to a spherical self-gravitating model due to the shortcomings of using a constant relaxation time for the mantle as opposed to mode-dependent relaxation times (Peltier1974).

A further improvement to an ice-sheet model can be made by coupling a model of solid earth deformation to the ice-sheet model. Studies have demonstrated that the feedback between the two systems can have large impacts on ice-sheet evolution (Gomez et al.2013; de Boer et al.2014). Using a coupled model Gomez et al. (2015) showed a reduced estimate of Antarctic ice-mass loss compared with a model without solid earth effects included. However, due to the large computational expense of these models, they remain at a relatively low resolution both spatially and temporally therefore omitting short wavelength and short timescale deformations. A recent study by Larour et al. (2019) showed that models need kilometre-scale resolution in the horizontal components to accurately predict ice-sheet evolution in the region of ice-sheet mass change, particularly for the short wavelength elastic component of solid earth deformation. This demonstrates the clear need for a full Stokes ice-sheet model capable of computing high resolution solid earth rebound.

Wu (2004) presented a recipe to adapt existing commercial finite-element codes to compute earth deformation as a response to ice loads, both for flat-earth as well as spherical self-gravitating setups. Finite elements have the advantage that they in general can use unstructured meshes in order to provide the resolution needed in regions where either physics or geometry demand it while keeping the model size limited. Many finite-element packages also include versatile solution methods that often also work in parallel computing environments – an essential feature to address continental-size problems at high resolution.

2 Mathematical and numerical model

The implementation of the viscoelastic rheology and additional force terms, to a large extent, follows the one suggested by Wu (2004). Adopting their notation, we start from the viscoelastic stress tensor, τ defined by the differential equation

(1) τ t = τ 0 t + μ ν τ - Π 1 ,

with the stress τ0 in the case of incompressibility given by

(2) τ 0 = Π 1 + 2 μ ϵ ,

where Π denotes the isotropic part of the Cauchy stress, i.e., the pressure. In the derivatives of Eqs. (1) and (2), t stands for time, 1 denotes the unit-tensor, μ the shear modulus and ν is the viscosity. The strain-tensor ϵ written in terms of the displacement d denotes as

(3) ϵ = sym ( d ) = 1 2 d + ( d ) T .

The linearized equation of motion for solid earth deformation (Wu2004) is given by

(4) τ - ( ρ 0 g 0 d ˙ ) - ρ 1 g 0 - g 0 ϕ 1 = 0 .

Where ρ0 and g0 are hydrostatic background density and gravity, respectively, and ρ1 is the perturbed density. The direction of g0 is in negative radial direction. According to Wu (2004,  Sect. 3) a flat-earth model is derived from Eq. 4 by assuming incompressibility and ignoring self-gravitational effects (i.e., redistribution of mass), making the third and fourth terms vanish. Further, sphericity is ignored, leading to changes aligned with the unit vector of a Cartesian system in vertical direction, ez. This leads to the equation of motion for a non-self-gravitating flat-earth model with layer-wise constant material. It reduces to a balance between the divergence of the stress (first term) and a restoring force due to the advection of pre-stress of the material (Wu2004)

(5) τ - ρ g ( e z d ) = 0 .

Here, ρ=ρ0 and g=||g0|| is the magnitude of the local acceleration by gravity, which points into the negative direction of ez.

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 (Wu2004) 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

(6) τ t τ i + 1 - τ i Δ t , Π t Π i + 1 - Π i Δ t .

Here, i is the current, and i+1 the implicit time step as well as Δt=ti+1-ti the time-step size between. The solution of the time-evolution problem then reads

(7) - 1 Π i + 1 + 2 μ Φ ϵ i + 1 = - Φ Π i + 2 μ Φ ϵ i - Φ τ i ,

with ϕ=1/(1+(μ/ν)Δt). The balance Eq. (5) of linear momentum is then solved for the new time step:

(8) τ i + 1 ( d ) - ρ g e z d i + 1 = 0 .

The weak formulation then results from the integral over the whole domain Ω (with its confining surface ∂Ω) using the test and weighting function vectors u,vH1:

(9) Ω τ ( u ) ( v ) d V - Ω ( τ ( u ) n ) v d A - Ω ρ g e z u v d V = 0 .

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, t=τ(u)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

(10) Ω τ ( u ) ϵ ( v ) d V - Ω ( τ ( u ) n ) v d A - Ω ρ g e z u v d V = 0 .

The system is completed by boundary conditions that are either provided by a value for any component of the stress vector, t=τ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 Benchmark tests

Benchmark tests are performed in order to validate the new implementation of Elmer/Earth in comparison to two other codes: ABAQUS and TABOO. We force the models with changing surface load, representing an idealized ice loading experiment. Specific geometry, earth structure and ice loading for the benchmarking case are described in Sect. 3.3. The two other codes are briefly introduced in the following sections.

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 (Wu2004). 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; Spada2003) 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 (Peltier1974). 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 Clark1976).

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


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 ν=1×1044Pa s in the lithosphere, thereby enforcing an approximately elastic behaviour on the timescale of the load. This can be justified by the Maxwell time tM=ν/μ being of the order of 1033s, 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×1018 and 1×1022Pa s, respectively, and the elastic and density parameters are depth-averaged values from the Preliminary Reference Earth Model (Dziewonski and Anderson1981; 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−3. 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.


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.


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 εd=1×10-7 .

4 Comparison of results

Comparing the results of the benchmarking exercise with two models that use different methods gives us confidence in the implementation of the new Elmer code. Figure 4 shows displacement with time at three locations – the centre of the disc (indicated by 0 km) and at 100 and 200 km distance from the centre of the disc.

Figure 4Comparison of results for deformation at the load centre (0 km), 100 and 200 km for Elmer/Earth, ABAQUS and TABOO.


The displacement curves for all three models over major parts of the simulation agree to within an order of 10 cm (see Fig. 5) in relation to a maximum deformation of 1.1 m by ABAQUS at the centre. The largest difference is observed at the centre of the disc where the Elmer/Earth model deforms slightly less than ABAQUS and almost insignificantly more as TABOO but reaches this deformation more quickly than the other codes (i.e. it has a faster relaxation time). As a consequence, Fig. 5 shows differences in vertical displacement between models (also between ABAQUS and TABOO) to be largest in the very beginning (when applying the load) and around the time of sudden unloading.

Figure 5Difference in deformation of Elmer/Earth relative to ABAQUS and TABOO at the load centre (0 km), 100 and 200 km.


The small differences between the results could be caused by several factors. Mesh differences between Elmer/Earth and ABAQUS are the likely cause of some small differences with ABAQUS having a regular grid mesh and Elmer having a finer mesh at the centre of the disc. There seems to be a correlation of the resolution in the centre with the displacement in both FEM-based models. It seems that the ABAQUS model setup does not provide enough horizontal mesh resolution at the centre, where the load is applied. This is confirmed by results obtained with mesh 2 (half mesh size) in Elmer/Earth, which produced displacements even larger than the one with the constant 10 km mesh from ABAQUS (see Sect. 5).

The deformation calculated by TABOO is less than Elmer/Earth and ABAQUS at each location. This may be due to the fundamental differences in the computation methods employed by the TABOO code, implementing normal mode methods rather than finite-element methods. Furthermore, TABOO computes deformation on a self-gravitating solid earth, whereas ABAQUS and Elmer do not include self-gravitation, which would result in some differences between these models. Nevertheless, the differences observed in the displacement curves are still within an acceptable tolerance.

5 Performance and accuracy of the Elmer/Earth deformation model

In order to obtain some insight into parallel performance as well as the dependency on the mesh resolution of Elmer/Earth, three meshes with different resolutions and mesh partitions (4, 16 and 32) have been created (see Table 2). Partitioning of the meshes has been performed by the mesh-conversion program ElmerGrid (part of the Elmer installation) using the METIS k-way partitioning scheme (Karypis and Kumar1998).

Table 2Parameters of the meshes and their partitions used for Elmer/Earth test runs.

Download Print Version | Download XLSX

Identical numerical parameters and methods, as described in Sect. 2, were applied throughout all runs.

5.1 Strong and weak scaling

Tests were performed on the Linux cluster raijin (Australian National Computational Infrastructure2017), 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 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.


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

6 Conclusions

We presented a newly implemented viscoelastic addition to the linear elasticity solver of the open-source finite-element package Elmer and its application to a flat-earth model. Robust projection of future ice-sheet change depends on coupled solid earth and ice dynamic processes at high spatial resolution, and Elmer/Earth provides a new open-source capability in conjunction with the existing ice-sheet model Elmer/Ice (Gagliardini et al.2013). Elmer/Earth, on its own, provides a new tool for modelling viscoelastic solid earth deformation due to surface loading changes.

For the time being, Elmer/Earth is a so-called flat-earth model (Wu2004). In its current state it ignores sphericity and self-gravitational effects and neglects accounting for the deformation induced by the redistribution of ocean water masses. This introduces certain limitations on its applicability (Wu and Johnston1998). Consequently, future applications of this particular model version should be confined to regional studies of ice sheets or highly localized loads, such as glaciers and ice caps.

We benchmarked Elmer/Earth with another FEM code, ABAQUS, as well as a spherical viscoelastic normal mode code, TABOO, and these comparisons show good agreement in the range of deviation in solution method as well as numerical approaches.

Scaling figures presented in Sect. 5 are what one would expect from other parallel performance tests of Elmer. A good performance tuning strategy will have to make sure that a good ratio between partition size (i.e., computation mainly bounded by memory access) and communication between the different MPI tasks is obtained. OpenMP multi-threading is in principle available for certain modules in Elmer but is not implemented for the linear elasticity solver; however, it might be a potential way to boost performance within a single node (Byckling et al.2017).

Code availability

Elmer (version 8.4) is available for download under GitHub (, last access: 4 March 2020). The revision (SHA-1 14c19b6) used in this study can be retrieved from (last access: 6 November 2019). TABOO is an open-source code available for download under GitHub (, last access: 4 March 2020). In this study we used version v1.1 (SHA-1 6163bec), which can be downloaded from (last access: 6 November 2019). ABAQUS is proprietary software and needs a purchased license. We used the ABAQUS 2018 release in this study. Information on how to obtain the software can be found under (last access: 6 November 2019).

Video supplement

The animation (, Zwinger2019) shows the temporal deformation of the benchmark case for the reference run (on mesh 1) as discussed in the article. Deformations in the video shown are exaggerated by a factor of 10 000.

Author contributions

TZ helped developing and implementing the model setup and performing the computations in Elmer. GAN contributed to the design of the benchmark setup and performed the computations in ABAQUS and TABOO. JR implemented the altered model equations in the source code of Elmer. MAK conceived the study and consulted in the model implementation and contributed to the design of the benchmark test. All authors contributed to the paper.

Competing interests

The authors declare that they have no conflict of interest.


Development of the viscoelastic model was supported under the Australian Research Council's Special Research Initiative for Antarctic Gateway Partnership (Project ID SR140300001) and Discovery Project DP170100224. Part of the work of Thomas Zwinger was enabled by a visiting scientist scholarship from UTAS. This research was undertaken with the assistance and resources from the National Computational Infrastructure (NCI Australia), an NCRIS-enabled capability supported by the Australian Government. We want to express our gratitude to Peter Råback (CSC) for solving a problem with post-processing of Elmer/Earth data and Fredrik Robertsén (CSC) for the discussion on scalability test results. We are grateful to Giorgio Spada for making TABOO open source. We want to thank the two reviewers and the editor for constructive suggestions to improve the quality of this paper.

Financial support

This research has been supported by the Australian Research Council (grant nos. SR140300001 and DP170100224).

Review statement

This paper was edited by Thomas Poulet and reviewed by PingPing Huang and Surendra Adhikari.


Adhikari, S., Ivins, E. R., Larour, E., Seroussi, H., Morlighem, M., and Nowicki, S.: Future Antarctic bed topography and its implications for ice sheet dynamics, Solid Earth, 5, 569–584,, 2014. a

Australian National Computational Infrastructure: Online description of raijin system, available at: (last access: 19 July 2019), 2017. a

Barletta, V., Bevis, M., Smith, B., Wilson, T., Brown, A., Bordoni, A., Willis, M., Khan, S., Rovira-Navarro, M., Dalziel, I., Smalley, R., Kendrick, E., Konfal, S., Caccamise, D., Aster, R., Nyblade, A., and Wiens, D.: Observed rapid bedrock uplift in Amundsen Sea Embayment promotes ice-sheet stability, Science, 360, 1335–1339,, 2018. a

Bueler, E., Lingle, C., and Brown, J.: Fast computation of a viscoelastic deformable Earth model for ice-sheet simulations, Ann. Glaciol., 46, 97–105,, 2007. a

Byckling, M., Kataja, J., Klemm, M., and Zwinger, T.: OpenMP SIMD Vectorization and Threading of the Elmer Finite Element Software, in: Proceedings 13th International Workshop on OpenMP, Springer Lecture Notes, 123–137,, 2017. a, b

de Boer, B., Stocchi, P., and van de Wal, R. S. W.: A fully coupled 3-D ice-sheet–sea-level model: algorithm and applications, Geosci. Model Dev., 7, 2141–2156,, 2014. a

Dziewonski, A. and Anderson, D.: Preliminary reference Earth model, Phys. Earth Planet. In., 25, 297–356,, 1981. a

Eisenstat, S., Elman, H., and Schultz, M.: Variational iterative methods for nonsymmetric systems of linear equations, SIAM J. Numer. Anal., 20, 345–357, 1983. a

Farrell, W. and Clark, J.: On Postglacial Sea Level, Geophys. J. Roy. Astr. S., 46, 647–667,, 1976. a

Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a new-generation ice sheet model, Geosci. Model Dev., 6, 1299–1318,, 2013. a

Geuzaine, C. and Remacle, J.-F.: Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, Int. J. Numer. Meth. Eng., 11, 1309–1331,, 2009. a

Gomez, N., Mitrovica, J., Huybers, P., and Clark, P.: Sea level as a stabilizing factor for marine-ice-sheet grounding lines, Nat. Geosci., 3, 850–853,, 2010. a

Gomez, N., Pollard, D., and Mitrovica, J.: A 3-D coupled ice sheet – sea level model applied to Antarctica through the last 40 ky, Earth Planet. Sc. Lett., 384, 88–99,, 2013. a, b

Gomez, N., Pollard, D., and Holland, D.: Sea-level feedback lowers projections of future Antarctic Ice-Sheet mass loss, Nat. Communi., 6, 8798,, 2015. a

Greve, R.: Glacial Isostasy: Models for the Response of the Earth to Varying Ice Loads, in: Continuum Mechanics and Applications in Geophysics and the Environment Springer, Berlin, Germany, 307–325,, 2001. a

Hibbitt, D., Karlsson, B., and Sorensen, P.: Getting Started with ABAQUS, Version (6.14), hibbitt, Karlsson & Sorensen, Inc. Edn., 2016. a

Karypis, G. and Kumar, V.: Multilevel k-way Partitioning Scheme for Irregular Graphs, J. Parallel Distr. Com., 48, 96–129,, 1998. a

Larour, E., Seroussi, H., Adhikari, S., Ivins, E., Caron, L., Morlighem, M., and Schlegel, N.: Slowdown in Antarctic mass loss from solid Earth and sea-level feedbacks, Science, 364, 7908,, 2019. a

Le Meur, E. and Huybrechts, P.: A comparison of different ways of dealing with isostasy: examples from modeling the Antarctic ice sheet during the last glacial cycle, Ann. Glaciol., 23, 309–317,, 1996. a, b, c

Nield, G., Barletta, V., Bordoni, A., King, M., Whitehouse, P., Clarke, P., Domack, E., Scambos, T., and Berthier, E.: Rapid bedrock uplift in the Antarctic Peninsula explained by viscoelastic response to recent ice unloading, Earth Planet. Sc. Lett., 397, 32–41,, 2014. a

Peltier, W.: The impulse response of a Maxwell earth, Rev. Geophys. Space Ge., 12, 649–669,, 1974. a, b

Råback, P., Malinen, M., Ruokolainen, J., Pursula, A., and Zwinger, T.: Elmer Models Manual, available at: (last access: 5 March 2020), 2019. a

Rutt, I., Hagdorn, M., Hulton, N., and Payne, A.: The Glimmer community ice sheet model, J. Geophys. Res., 114, F02004,, 2009. a

Schoof, C.: Ice Sheet Grounding Line Dynamics: Steady States, Stability, and Hysteresis, J. Geophys. Res., 112, F03S28,, 2007. a

Spada, G.: The theory behind TABOO, available at: (last access: 5 March 2020), 2003. a

Spada, G., Antonioli, A., Boschi, L., Brandi, V., Cianetti, S., Galvani, G., Giunchi, C., Perniola, B., Piana Agostinetti, N., Piersanti, A., and Stocchi, P.: TABOO, User Guide, available at: (last access: 5 March 2020), 2003. a

Steffen, H., Kaufmann, G., and Wu, P.: Three-dimensional finite-element modeling of the glacial isostatic adjustment in Fennoscandia, Earth Planet. Sc. Lett., 250, 358–375,, 2006. a, b

Whitehouse, P. L.: Glacial isostatic adjustment modelling: historical perspectives, recent advances, and future directions, Earth Surf. Dynam., 6, 401–429,, 2018. a

Wu, P.: Using commercial finite element packages for the study of earth deformations, sea levels and the state of stress, Geophys. J. Int., 158, 401–408,, 2004.  a, b, c, d, e, f, g, h, i

Wu, P. and Johnston, P.: Validity of Using Flat-Earth Finite Element Models in the Study of Postglacial Rebound, in: Dynamics of the Ice Age Earth, Trans Tech Publications Ltd, Switzerland, 191–202, 1998. a

Zwinger, T.: Deformation in flat-earth model Elmer/Earth,, 2019. a

Short summary
We present a newly developed flat-earth model, Elmer/Earth, for viscoelastic treatment of solid earth deformation under ice loads. Unlike many previous approaches with proprietary software, this model is based on the open-source FEM code Elmer, with the advantage for scientists to apply and alter the model without license constraints. The new-generation full-stress ice-sheet model Elmer/Ice shares the same code base, enabling future coupled ice-sheet–glacial-isostatic-adjustment simulations.