Finite-element modelling of glacial isostatic adjustment (GIA): Use of elastic foundations at material boundaries versus the geometrically non-linear formulation

Abstract

Abstract Commercial finite-element software packages like ABAQUS are widely used for geodynamic modelling, which usually requires considering isostatic effects, for example, to calculate displacements and stresses resulting from glacial isostatic adjustment (GIA). Since Wu (Geophysical Journal International, 158, 401–408, 2004) proposed that models created with commercial finite-element software need to implement elastic ( Winkler ) foundations at all external and internal material boundaries to account for restoring forces, this approach has been applied by many GIA and other geoscientific studies. However, there is no consensus about the necessity of implementing elastic foundations, which have the disadvantage that the stress output needs post-processing to obtain meaningful results. Here we demonstrate that the elastic-foundation approach was derived from an Eulerian formulation of the equation of motion for elastic and viscoelastic materials. Finite-element codes like ABAQUS, however, use a Lagrangian formulation, which renders the implementation of elastic foundations at all material boundaries unnecessary if the geometrically non-linear formulation (NLGEOM) is used. Results from viscoelastic half-space models show that for incompressible viscoelastic materials models with elastic foundations (but no NLGEOM) and models with NLGEOM (but no elastic foundations) yield vertical displacements underneath the load that differ by less than 1% from the analytical solution. Both models reach a state of isostatic equilibrium. In contrast, models without NLGEOM and elastic foundations do not reach isostatic equilibrium, i.e. the model surface continuously subsides under the load. Models with both NLGEOM and elastic foundations behave overly stiff and yield wrong displacement fields. Results from models based on Archimedes principle demonstrate that restoring forces are correctly calculated when using NLGEOM, which has the advantage that the stress output can be used without post-processing and that negative density contrasts (e.g. between salt and surrounding rocks) can be taken into account within the model domain, which is impossible when using elastic foundations.