Nordic countries experience post-glacial rebound, a movement where geographical contours slowly change elevation with respect to the mean sea level. The glacial isostatic adjustment (GIA) model aims to explain the phenomena, which combined with seismic data allows geoscientists to reconstruct elastic coefficients and viscosities of the Earth's lithosphere and upper mantle. The use of standard commercial codes are not adequate for GIA simulations and result in significant errors in the displacement field.
This negative outcome suggests the development of GIA codes that include advection of pre-stress in the model. The problem set up consists on a solid 2D elastic layer under a flat Earth approximation, described by three different models suggested by current studies in geophysics.For space discretization the mixed finite element method (mFEM) is used and efficient preconditioners are built for the resulting algebraic system in saddle point form. A three level GMRES iterative solution strategy is proposed, based on Schur Complement preconditioners coupled with Multigrid techniques. The implementation is presented as a ready-to-use toolbox that easily deals with problem parameters, geometries, compressible and fully incompressible materials and provides higher accuracy for the displacement field compared with the previously existent codes. It also can be easily extended to 3D geometries and allows the implementation of a viscoelastic mantle. The code is written in C++ using the deal.II library designed for FEM, permitting the use of readily-made software packages, such as Trilinos that are straightforwardly parallelizable.