The Earth is subjected to 100,000 year cycles of glaciation and deglaciation. The deformations induced by glacial and oceanic loading and the continuous attempt at recovery of the isostatic equilibrium within the solid Earth, are referred to as Glacial Isostatic Adjustment (GIA). This process is ongoing still and yields a large contribution to present day surface deformation and sea level change in formerly glaciated areas. In order to accurately model GIA, the lateral viscosity variations within the interior of the Earth are accounted for (Kaufmann et al., 2000; Steffen et al., 2006; Wu and van der Wal, 2003). Additionally, an increased level of accuracy is obtained by adopting a combination of linear and non-linear viscoelasticity as demonstrated by Barnhoorn et al. (2011); Forno and Gasperini (2007); van der Wal et al. (2013);Wu and Wang (2008). 

During the Last Glacial Maximum, the British-Irish and Fennoscandian Ice Sheets covered large parts of Northwestern Europe. The interior of the Earth in this area is known to consist of material of very heterogeneous tectonic origin (Artemieva et al., 2006). Additionally, research in this area is promoted by the availability of the independent regional ice model Bradley2018 (Bradley, personal communication), an RSL observation database for the Rhine-Meuse Delta (Hijma and Cohen, 2019), and a collection of GPS derived uplift rates throughout Europe (Teferle, personal communication). 

At the Astrodynamics and Space Missions research group of Delft University of Technology, a 3D GIA FEM model has been developed to model GIA in Antarctica (Blank et al., 2017). This model follows the work by Wu (2004) and van der Wal et al. (2013), and is complemented with an iterative algorithm to solve the sea level equation in accordance with Kendall et al. (2005). This research aims to provide a single GIA model best suited for the prediction of GIA induced vertical surface deformation in Northwestern Europe, by adapting the existing model. In doing so, a better understanding of the interior of the Earth in Northwestern Europe can be achieved.

The response of the Earth is dictated by the composite rheology creep flow laws for olivine (Hirth and Kohlstedt, 2003). By varying the grain size as well as the water content of the mantle material, and by implementing a global temperature model of the Earth’s interior, four 3D composite rheology Earth models are obtained. The fifth Earth model considered is the radially symmetric VM5a viscosity profile developed by Peltier et al. (2015) in conjunction with the global ice model ICE-6G_C. The performance of all five Earth model configurations in combination with both the ICE-6G_C model and the Bradley2018 model is analysed in terms of relative sea level and uplift rates.

It is found that the Bradley2018 model is the preferred ice model for GIA modelling in Northwestern Europe. The ICE-6G_C model outperforms the Bradley2018 model at far-field RSL sites, which is attributed to its superior representation of global eustatic sea level rise. The 3D composite rheologies lead to improved fits to RSL observations for the majority of the investigated measurement sites compared to the 1D scenario. The dry 4 mm grain size rheology yields the best overall performance out of all rheological configurations considered. A preference towards wet rheology exists in regions of Sveco-Norwegian tectonic origin. The strongest rheology is preferred in the mid-west of Scotland. No definitive connection is found between the local tectonic origin and preferred rheology fromRSL simulations. It is believed that this analysis may benefit from the inclusion of laterally varying grain sizes and water content inferred from geophysical observations, as well as the extension of the variable space for the water content.

For both ice models, an improved fit to observed uplift rates can be obtained through the application of a 3D composite rheology. The GPS derived uplift rates can be reproduced best using the dry 10 mm grain size rheology in combination with the Bradley2018 ice history. This model is deemed to be best suited for simulation of GIA induced uplift rates in Northwestern Europe. The second-best performance in terms of uplift rate is found using the 1D Earth model. In Scandinavia the 4 mm dry rheology yields uplift rates equal to roughly half the observed uplift rates, while the uplift rates for the 10 mm and 4 mm wet rheologies are near-zero. In the far-field, where other surface deformation mechanisms may infer a larger deformation rate than GIA (Fokker et al., 2018), no model could reproduce the observed uplift rates. 

The presence of a high viscosity anomaly beneath Eastern Fennoscandia is captured by the 3D rheologies and results in a shift of the centre of positive and negative uplift rates. As the spatial distribution of minima and maxima in both uplift rates and RSL change rates is sensitive to the inclusion of a 3D rheology, this should be accounted for in future regional sea level change and surface deformation projections.