Automated forward and adjoint modelling of viscoelastic deformation of the solid Earth

Robust models of viscoelastic Earth deformation under evolving surface loads underscore many problems in geodynamics and are particularly critical for paleoclimate and sea-level studies through their role in Glacial Isostatic Adjustment (GIA). A long-standing challenge in GIA research is to perform computationally efficient inversions for ice-loading histories and mantle structure using a physically realistic Earth model that incorporates three-dimensional viscosity variations and/or complex rheologies. For example, recent geodetic observations from melting ice sheets appear inconsistent with long-term sea-level records and have been used to argue for transient rheologies, generating debate in the literature and leaving large uncertainties in projections of future sea-level change. Here, we extend the applicability of G-ADOPT (a Firedrake-based finite element framework for geoscientific adjoint optimisation) to these problems. Our implementation solves the equations governing viscoelastic surface loading while naturally accommodating elastic compressibility, lateral viscosity variations, and non-Maxwell rheologies (including transience). We benchmark the approach against a suite of analytical and numerical test cases, demonstrating both accuracy and computational efficiency. Crucially, G-ADOPT enables automatic derivation of adjoint sensitivity kernels, allowing gradient-based optimisation strategies that are essential for high-dimensional inverse problems. Using synthetic Earth-like experiments, we illustrate its capability to reconstruct ice histories and recover mantle viscosity variations, providing a roadmap towards data assimilation and uncertainty quantification in GIA modelling and sea-level projections.