A generalised Fourier collocation for fast computation of glacial isostatic adjustment

Glacial isostatic adjustment (GIA) represents an important negative feedback on ice-sheet dynamics. The magnitude and time scale of GIA primarily depend on the upper mantle viscosity and the lithosphere thickness. These parameters have been found to vary strongly over the Antarctic continent, showing ranges of 10¹⁸ - 10²³ Pa s for the viscosity and 30 - 250 km for the lithospheric thickness. Recent studies show that coupling ice-sheet models to 3D GIA models capturing these spatial dependencies results in substantial differences in the evolution of the Antarctic Ice Sheet compared to the use of 1D GIA models, where the solid-Earth parameters are assumed to depend on the latitude but not on the longitude and the depth. However, 3D GIA models are computationally expensive and sometimes require an iterative coupling for the ice sheet and the solid-Earth solutions to converge. As a consequence, their use remains limited, potentially leading to errors in the simulated ice-sheet response and associated sea-level rise projections. Here, we propose to tackle this problem by generalising the Fourier collocation method for solving GIA proposed by Lingle and Clark (1985) and implemented by Bueler et al. (2007). The method allows for an explicit accounting of the effects of spatially heterogeneous viscosity and lithospheric thicknesses and is computationally very efficient. Thus, for a continental domain at relatively high spatial resolution (256 x 256 grid points) and a 1-year time step, the model runs with speeds of ca. 200 simulation years per second on a single CPU, while keeping the error low compared to 3D GIA models. As the time step is small enough, the need of an iterative coupling method is avoided, thus making the model easy to couple with ice-sheet models.