This paper presents a method that modifies commercial engineering-oriented finite element packages for the modelling of Glacial Isostatic Adjustment (GIA) on a self-gravitating, compressible and spherical Earth with 3-D structures. The approach, called the iterative finite element body and surface force (FEMIBSF) approach, solves the equilibrium equation for deformation using the ABAQUS finite element package and calculates potential perturbation consistently with finite element theory, avoiding the use of spherical harmonics. The key to this approach lies in computing the mean external body forces for each finite element within the Earth and pressure on Earth's surface and core-mantle boundary (CMB). These quantities, which drive the deformation and stress perturbation of GIA but are not included in the equation of motion of commercial finite element packages, are implemented therein. The method also demonstrates how to calculate degree-1 deformation directly in the spatial domain and Earth-load system for GIA models. To validate the FEMIBSF method, loading Love numbers (LLNs) for homogeneous and layered earth models are calculated and compared with three independent GIA methodologies: the normal-mode method, the iterative body force method and the spectral-finite element method. Results show that the FEMIBSF method can accurately reproduce the unstable modes for the homogeneous compressible model and agree reasonably well with the Love number results from other methods. It is found that the accuracy of the FEMIBSF method increases with higher resolution, but a non-conformal mesh should be avoided due to creating the so-called hanging nodes. The role of a potential force at the CMB is also studied and found to only affect the long-wavelength surface potential perturbation and deformation in the viscous time regime. In conclusion, the FEMIBSF method is ready for use in realistic GIA studies, with modelled vertical and horizontal displacement rates in a disc load case showing agreement with other two GIA methods within the uncertainty level of GNSS measurements.

It has been argued that viscous dissipation from mantle flow in response to surface loading during glacial cycles can result in short-term heating and thus trigger transient volcanism or changes in mantle properties, which may in turn affect mantle dynamics. Furthermore, heating near the Earth's surface can also affect the stability of ice sheets. We have studied the magnitude and spatial-temporal distribution of viscous heating induced in the mantle by the realistic ice model ICE-6G and gravitationally consistent ocean loads. Three types of mantle rheologies, including linear, non-linear and composite rheologies are considered to see if non-linear creep can induce larger viscous heating than linear rheology. We used the Coupled-Laplace-Finite-Element model of Glacial Isostatic Adjustment (GIA) to compute the strain, stress and shear heating during a glacial cycle. We also investigated the upper bound of temperature change and surface heat flux change due to viscous heating. We found that maximum viscous heating occurs near the end of deglaciation near the edge of the ice sheet with amplitude as high as 120 times larger than that of the chondritic radioactive heating. The maximum heat flux due to viscous heating can reach 30 mW m^-2, but the area with large heat flux is small and the timescale of heating is short. As a result, the upper bound of temperature change due to viscous heating is small. Even if 30 glacial cycles are included, the largest temperature change can be of the order of 0.3 °C. Thus, viscous heating induced by GIA cannot induce volcanism and cannot significantly affect mantle material properties, mantle dynamics nor ice-sheet stability. The Author(s) 2017. Published by Oxford University Press. All rights reserved.