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 power-law and 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.

1 Introduction

Modelling the solid Earth's viscoelastic response under evolving surface loads is a foundational problem in geodynamics that has a rich history spanning more than a century of research (e.g. Gilbert, 1886; Haskell, 1935; Wu and Peltier, 1982; Mitrovica, 1996). Surface loads of geological interest can be of both natural and anthropogenic origin and operate across a wide range of spatial and temporal scales – from metres to thousands of kilometres, and from minutes to millions of years. Pertinent examples include coseismic deformation (e.g. Pollitz, 1996), ocean and body tides (e.g. Dehant et al., 1999; Lau et al., 2017; Matviichuk et al., 2020), reservoir impoundment (e.g. Chao et al., 2008), redistribution of surface and groundwater (e.g. Crittenden, 1963; Crowley et al., 2006; Austermann et al., 2020), volcanic loading (e.g. Moore, 1970; Rundle, 1982), and the formation and melting of ice sheets (e.g. Haskell, 1935; Cathles, 2016). Understanding these processes has required the development of increasingly complex models of deformation within Earth's interior and at its surface.

One of the most societally relevant applications of viscoelastic modelling is glacial isostatic adjustment (GIA), a process that links the evolution of ice sheets to solid Earth deformation and sea-level change (Daly, 1925; Haskell, 1935; Farrell and Clark, 1976; Peltier et al., 2015). Due to the combined effects of viscoelastic uplift and subsidence, gravitational changes, and feedbacks on Earth's rotation, GIA introduces strong geographic variability into sea-level change (e.g. Mitrovica et al., 2011). Attempts to constrain modern sea-level trends driven by climate change are further complicated by the fact that, due to the mantle's high viscosity, Earth is still responding to the retreat of major ice sheets over Fennoscandia and North America since the Last Glacial Maximum (LGM) around 20 000 years ago (e.g. Lambeck et al., 2014). All of these factors must be considered when interpreting geodetic observations such as Global Navigation Satellite System (GNSS) and GRACE data (Church et al., 2004; Tregoning et al., 2009; King et al., 2012; Ivins et al., 2013; Caron and Ivins, 2020) and a complete understanding of GIA – past, present and future – is critical in generating robust projections of future sea level.

Viscoelastic surface-loading phenomena can generally be considered as inverse problems. The rheological properties of the solid Earth are poorly known, and observations of its deformation in response to various loads must therefore be used to infer these properties. This situation is compounded in the case of GIA, since present-day observables such as crustal uplift rates or geoid changes are jointly controlled by rheological properties and the history of ice-ocean loading – both of which are poorly constrained. Progress in this field has relied on fitting combinations of ice history and Earth structure to diverse datasets including relative sea level (RSL) indicators (e.g. tide gauge records, fossil corals, paleo strand-lines), geomorphic reconstructions of past ice extent (e.g. terminal morraines, cosmogenic exposure ages), GNSS surface motions, satellite gravity data, and indirect constraints from seismology and other palaeoclimate archives (Lambeck et al., 2014; Peltier et al., 2015; Gowan et al., 2021).

Most GIA models exploring these datasets assume a radially symmetric (1D) Earth and represent its viscoelastic response with a Maxwell rheology (Wu and Peltier, 1982). However, we know from mantle convection studies that mantle rheology deviates significantly from this axisymmetric assumption, being strongly dependent on temperature and strain-rate among other factors (e.g. Moresi and Solomatov, 1998; Tackley, 2000), and resulting in lateral viscosity variations that span several orders of magnitude. Sea-level records and geodetic observations confirm that such 3D viscosity variations exert a first-order influence on GIA deformation patterns (Wu and van der Wal, 2003; Austermann et al., 2013; Nield et al., 2014; Hay et al., 2017; Barletta et al., 2018; Powell et al., 2020; Austermann et al., 2021). This recognition has driven the development of 3D GIA models using finite volume (Latychev et al., 2005), finite element (Wu, 2004; Zwinger et al., 2020; Weerdesteijn et al., 2023), and spectral approaches (Martinec, 2000; Tanaka et al., 2011). At the same time, both laboratory deformation experiments and seismic anisotropy observations indicate that non-Maxwell rheologies, including power-law creep, are likely to be important, at least in the upper mantle (Karato and Wu, 1993). In these cases, the strain rate has a non-linear dependence on driving stress, which has been explored in some more advanced GIA models (Wu, 1995; Gasperini et al., 2004; van der Wal et al., 2010). Furthermore, both laboratory experiments and observations of post-seismic deformation suggest that Earth's rheology may also depend on the frequency of loading (i.e., exhibit transient behaviour; Karato and Wu, 1993; Pollitz, 2003; Pollitz et al., 2006; Faul and Jackson, 2015; Yamauchi and Takei, 2016; Lau, 2024). In regions such as West Antarctica and coastal Greenland, recent observations of rapid GNSS-measured uplift rates, compared with longer-term post-LGM RSL curves, have sparked renewed debate, with both transient rheology and 3D viscosity variations proposed as viable explanations (e.g. Nield et al., 2014, 2018; Barletta et al., 2018; Adhikari et al., 2021; Lau et al., 2021; Pan et al., 2024; Weerdesteijn and Conrad, 2024).

While the importance of accounting for 3D Earth structure and potentially complex rheologies in GIA is now beyond doubt, 3D GIA models suffer from being computationally expensive. Conventional “brute-force” inversion strategies, whereby model sensitivities are assessed by running many forward simulations, are therefore impractical. Consequently, most 3D GIA models are still run in a forward sense and rely on LGM ice histories independently derived from 1D inversions and 3D Earth structures inferred from other datasets, such as seismic tomography. A transformative alternative to brute-force inversion is gradient-based optimisation, which is enabled by adjoint methods. This technique allows efficient computation of model sensitivities – that is, the gradient of a misfit function with respect to a large number of model parameters (or “controls”). Unlike in a finite difference approach, where cost scales with the number of control parameters (which is often tied to mesh resolution), the adjoint method computes gradients at a cost comparable to a single linearised forward solve, independent of the number of control parameters (Errico, 1997). This aspect makes it uniquely scalable and powerful for PDE-constrained optimisation and it is now widely adopted in fields such as numerical weather prediction (Kalnay, 2003), oceanography (Forget et al., 2015), seismic imaging (Tromp et al., 2005; Fichtner et al., 2006), and mantle convection (Ghelichkhan et al., 2021).

In geodynamics, the past decade has seen substantial theoretical progress in applying adjoint methods to viscoelastic loading problems. In the case of GIA, their use was first demonstrated by Al-Attar and Tromp (2014), with extension to sea-level modelling on a 1D Earth undertaken by Crawford et al. (2018) and on a 3D Earth by Lloyd et al. (2024). These developments have revealed strong diagnostic power: adjoint sensitivity kernels clearly expose the limitations of using a 1D Earth and the first-order impact of lateral viscosity variations on GIA observables. As such, they are increasingly recognised as a game-changing tool for tackling the spatial complexity and high dimensionality of the GIA inverse problem. Furthermore, recent advances in second-order adjoint techniques may offer a tantalising route to formal uncertainty quantification through efficient Hessian-based inference (e.g. Yu et al., 2025). Nevertheless, there remain obstacles to the broader use of adjoint techniques within the GIA modelling community. A central challenge lies in the complexity of deriving, implementing, and validating adjoint models for large-scale, non-linear, time-dependent problems. Compounding this is the legacy of traditional GIA workflows, which have typically relied on 1D Earth structures and relatively simple rheologies.

The shift toward 3D solid Earth structure, transient and/or non-linear rheologies, and formal data assimilation places new demands on model infrastructure, underscoring the need for automated, scalable, and efficient adjoint-capable frameworks. The automatic adjoint framework provided by the Geoscientific Adjoint Optimisation Platform (G-ADOPT; Davies et al., 2022; Ghelichkhan et al., 2024; Gibson et al., 2024) offers a compelling and timely solution. G-ADOPT combines three state-of-the-art software libraries: (i) Firedrake, a flexible and automated system for solving partial differential equations using the finite element method (Ham et al., 2023); (ii) Pyadjoint, which enables automatic derivation of discrete adjoints within Firedrake (Farrell et al., 2013; Mitusch et al., 2019); and (iii) the Rapid Optimisation Library (ROL), a high-performance optimisation engine within Trilinos (The ROL Project Team, 2022). Together, these components deliver a scalable, modular and performant modelling framework that can robustly handle complex forward models and compute their adjoints with theoretical optimal efficiency.

In this paper, we extend, benchmark, and demonstrate the applicability of G-ADOPT for solving viscoelastic surface-loading problems. The framework accommodates arbitrary rheologies, including power-law and transient behaviour, and can operate in both Cartesian and spherical geometries. We verify its accuracy against analytical solutions (Sect. 3.1), benchmark it against published 3D Cartesian models under time-dependent loading scenarios (Sect. 3.2), and extend it to incorporate compressible elasticity, power-law and Burgers rheologies (Sect. 3.3), and spherical domains with lateral viscosity variations (Sect. 3.4). Finally, we use a series of twin experiments (Sect. 3.5) to showcase the framework's adjoint-based inversion capabilities, recovering unknown ice-loading histories and mantle structure in a 2D annulus. Although self-gravitation is not yet included, these developments represent a significant step towards a next-generation, data-assimilating framework for GIA and other surface-loading problems within G-ADOPT.

2 Method

In this section we present the non-dimensional equations governing the viscoelastic loading problem, consisting of the conservation laws of mass and momentum driven by an evolving surface load. The non-dimensional mass conservation equation is given by

$$\begin{array}{}\text{(1)}& {\partial }_t\mathit{\rho }+\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{v}\right)=\mathrm{0},\end{array}$$

where v and ρ denote the non-dimensional velocity and density, respectively. In non-dimensionalising the velocity, we have used the scale $\frac{L}{\overline{\mathit{\alpha }}}$, where L is a characteristic length scale and $\overline{\mathit{\alpha }}$ is a characteristic Maxwell time that describes the transition between dominantly elastic and viscous behaviour

$$\begin{array}{}\text{(2)}& \overline{\mathit{\alpha }}={\displaystyle \frac{\overline{\mathit{\eta }}}{\overline{\mathit{\mu }}}},\end{array}$$

where $\overline{\mathit{\eta }}$ and $\overline{\mathit{\mu }}$ are characteristic dynamic viscosity and shear modulus values respectively (see Table 1).

The density field is decomposed into a time-independent background component, ρ₀, and a small time-dependent perturbation, ρ₁, such that $\mathit{\rho }={\mathit{\rho }}_{\mathrm{0}}+{\mathit{\rho }}_{\mathrm{1}}$. Substituting this decomposition into Eq. (1) yields the linearised mass conservation equation

$$\begin{array}{}\text{(3)}& {\partial }_t{\mathit{\rho }}_{\mathrm{1}}+\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{0}}\mathit{v}\right)=\mathrm{0},\end{array}$$

where we have neglected the ∇⋅(ρ₁ v) term, by assuming that perturbations in density are small compared with background values. Integrating this equation in time and using the relation ${\partial }_t\mathit{u}=\mathit{v}$, where u is the non-dimensional displacement vector, we obtain the following expression for the density perturbation

$$\begin{array}{}\text{(4)}& {\mathit{\rho }}_{\mathrm{1}}\left(t\right)=-\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{0}}\mathit{u}\left(t\right)\right),\end{array}$$

assuming a vanishing initial displacement, $\mathit{u}(t=\mathrm{0})=\mathbf{\text{0}}$.

The full GIA equations include rotational and gravitational effects arising from changes in the surface load as water is redistributed between ice and ocean, which is governed by the sea-level equation (Kendall et al., 2005). Neglecting those terms, alongside inertial terms, yields the following non-dimensional form of the momentum equation

$$\begin{array}{}\text{(5)}& \mathbf{\text{0}}=\mathrm{\nabla }\cdot \mathit{\sigma }-B_{\mathit{\mu }}({\mathit{\rho }}_{\mathrm{0}}+{\mathit{\rho }}_{\mathrm{1}})g{\widehat{\mathit{e}}}_k,\end{array}$$

where σ is the non-dimensional stress tensor, g is the non-dimensional gravity, and $B_{\mathit{\mu }}=\frac{\overline{\mathit{\rho }}\overline{g}L}{\overline{\mathit{\mu }}}$ is a non-dimensional number describing the ratio of buoyancy to elastic shear strength. Note that ${\widehat{\mathit{e}}}_k$ is aligned with either the z-axis or radial direction in Cartesian or spherical coordinates, respectively. We refer again to Table 1 for the characteristic non-dimensional scales.

Table 1 Summary of notation and model parameters.

To further simplify the momentum equation, we linearise perturbations in the stress field by expressing the total stress as a sum of a time-invariant hydrostatic background component and a perturbation: $\mathit{\sigma }={\mathit{\sigma }}_{\mathrm{0}}+{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{E}}$, where the superscript E indicates an Eulerian perturbation. Since viscoelastic constitutive relations that link stress to strain (and their time derivatives) are fundamentally defined in a Lagrangian framework, we also need to transition from an Eulerian to a Lagrangian description of the stress field. To capture the material response in a Lagrangian framework, the perturbation must account for the change in stress due to displacement of material through the background stress gradient. To first order, the incremental Lagrangian stress is ${\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{L}}={\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{E}}+\mathit{u}\cdot \mathrm{\nabla }{\mathit{\sigma }}_{\mathrm{0}}$ (see Eq. 3.16 of Dahlen and Tromp, 1998). Substituting into Eq. (5), we obtain

$$\begin{array}{}\text{(6)}& \mathbf{\text{0}}=\mathrm{\nabla }\cdot ({\mathit{\sigma }}_{\mathrm{0}}+{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{L}}-\mathit{u}\cdot \mathrm{\nabla }{\mathit{\sigma }}_{\mathrm{0}})-B_{\mathit{\mu }}({\mathit{\rho }}_{\mathrm{0}}+{\mathit{\rho }}_{\mathrm{1}})g{\widehat{\mathit{e}}}_k.\end{array}$$

We assume that time-invariant components are in equilibrium, such that the background stress field satisfies a hydrostatic balance with the background density

$$\begin{array}{}\text{(7)}& \mathrm{\nabla }\cdot {\mathit{\sigma }}_{\mathrm{0}}=B_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{0}}g{\widehat{\mathit{e}}}_k.\end{array}$$

Moreover, assuming that the background density ρ₀ varies only in the radial direction, the advection of hydrostatic pre-stress simplifies to

$$\begin{array}{}\text{(8)}& \mathit{u}\cdot \mathrm{\nabla }{\mathit{\sigma }}_{\mathrm{0}}=B_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{0}}\mathit{u}\cdot g{\widehat{\mathit{e}}}_k\mathit{I},\end{array}$$

where I is the identity tensor. Substituting these expressions leads to the final linearised form of the momentum equation

$$\begin{array}{}\text{(9)}& \mathbf{\text{0}}=\mathrm{\nabla }\cdot {\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{L}}-B_{\mathit{\mu }}\mathrm{\nabla }\left({\mathit{\rho }}_{\mathrm{0}}\mathit{u}\cdot g{\widehat{\mathit{e}}}_k\right)-B_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{1}}g{\widehat{\mathit{e}}}_k.\end{array}$$

2.1 Constitutive Equation

We follow the internal variable formulation adopted by Al-Attar and Tromp (2014) and Crawford et al. (2017, 2018), in which viscoelastic constitutive equations are expressed in integral form and reformulated using so-called internal variables. Conceptually, this approach consists of a set of elements with different shear relaxation timescales, arranged in parallel.

For a linear, compressible viscoelastic material, the non-dimensional constitutive equation takes the form

$$\begin{array}{}\text{(10)}& {\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{L}}=\mathit{\kappa }\mathrm{\nabla }\cdot \mathit{u}\left(t\right)\mathit{I}+\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}\mathit{d}\left(t\right)-\mathrm{2}\sum _i{\mathit{\mu }}_i{\mathit{m}}_i\left(t\right),\end{array}$$

where κ is the non-dimensional bulk modulus and μ₀ is the non-dimensional effective shear modulus given by

$$\begin{array}{}\text{(11)}& {\mathit{\mu }}_{\mathrm{0}}=\sum _i{\mathit{\mu }}_i\end{array}$$

where μ_i are the non-dimensional shear moduli associated with each internal variable, m_i. The deviatoric strain tensor is given as

$$\begin{array}{}\text{(12)}& \mathit{d}=\mathit{e}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\mathrm{Tr}\left(\mathit{e}\right)\mathit{I},\end{array}$$

where Tr(⋅) is the trace operator and the strain tensor, e, is

$$\begin{array}{}\text{(13)}& \mathit{e}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{\nabla }\mathit{u}+{\left(\mathrm{\nabla }\mathit{u}\right)}^T\right).\end{array}$$

We note that all non-dimensional moduli are obtained by division of terms by the characteristic shear modulus $\overline{\mathit{\mu }}$ (Table 1). Each internal variable, m_i, is defined by

$$\begin{array}{}\text{(14)}& {\mathit{m}}_i={\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_i}}\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{\mathrm{e}}^{-\frac{(t-t^{\prime })}{{\mathit{\alpha }}_i}}\mathit{d}\left(t^{\prime }\right)\mathrm{d}{t}^{\prime },\end{array}$$

where α_i is the non-dimensional Maxwell time for each element.

Equivalently, each internal variable evolves according to

$$\begin{array}{}\text{(15)}& {\partial }_t{\mathit{m}}_i+{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_i}}\left({\mathit{m}}_i-\mathit{d}\right)=\mathbf{0},{\mathit{m}}_i\left(t_{\mathrm{0}}\right)=\mathbf{0}.\end{array}$$

This formulation provides a compact, flexible and convenient means to incorporate transient rheology into viscoelastic deformation models: using a single internal variable is equivalent to a simple Maxwell material; two correspond to a Burgers model with two characteristic relaxation frequencies; and using a series of internal variables permits approximation of a continuous range of relaxation timescales for more complicated transient rheologies (e.g., Extended Burgers). This formulation is also readily extensible to non-linear constitutive equations (e.g., Crawford et al., 2017). For example, we can represent a composite rheology that combines the effects of linear diffusion creep at low stresses with non-Newtonian power-law rheology (e.g., dislocation creep) at high stresses (van der Wal et al., 2010; Kang et al., 2022) by modifying the viscosity to

$$\begin{array}{}\text{(16)}& \mathit{\eta }={\displaystyle \frac{{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{1}+{\left(\frac{\left|\right|\mathit{\tau }\left|\right|}{{\mathit{\tau }}_T}\right)}^{n-\mathrm{1}}}},\end{array}$$

where η₀ is the same spatially varying Newtonian (i.e., stress-independent) viscosity as before, $\left|\right|\mathit{\tau }\left|\right|$ is the second invariant of the deviatoric stress tensor, τ_T is a transition stress beyond which power-law deformation starts to dominate, and n is the associated power-law exponent.

2.2 Boundary conditions

To complete the problem specification, we must also define the boundary conditions. At the top of the domain (i.e., Earth's surface), normal stress is balanced by the applied surface load such that

$$\begin{array}{}\text{(17)}& \widehat{\mathit{n}}\cdot {\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{L}}=-B_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{load}}g{h}_{\mathrm{load}}{\widehat{\mathit{e}}}_k,\end{array}$$

where ρ_load and h_load are the non-dimensional density and vertical thickness of the surface load, and $\widehat{\mathit{n}}$ is the outward unit normal vector to the boundary. At the base of the domain (i.e., core–mantle boundary), we impose a no-normal-displacement condition

$$\begin{array}{}\text{(18)}& \mathit{u}\cdot \widehat{\mathit{n}}=\mathrm{0},\end{array}$$

which is consistent with the approximation of a rigid core following Weerdesteijn et al. (2023). Note that this choice departs from models that treat the core as an inviscid fluid and instead impose a free-surface condition (e.g., Zhong et al., 2003; Latychev et al., 2005). After self-gravity has been introduced into G-ADOPT, switching to a free-surface should be a trivial extension.

Finally, we assume continuity of both displacement and traction across all internal boundaries, such that

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }[\mathit{u}{]}_{-}^{+}=\mathbf{0},\text{(20)}& {\displaystyle }& {\displaystyle }[\widehat{\mathit{n}}\cdot {\mathit{\sigma }}_{L\mathrm{1}}{]}_{-}^{+}=\mathbf{0},\end{array}$$

where $[\cdot {]}_{-}^{+}$ denotes the jump in the associated property across an interface.

In summary, the governing equations for viscoelastic loading consist of the conservation of momentum (Eq. 9) driven by the evolving surface load (Eq. 17) along with the conservation of mass (Eq. 4) and a viscoelastic constitutive equation (Eqs. 10 and 15). Note that this linearised formulation is strictly relevant for small displacements with respect to the depth of the mantle and that the background stress field is assumed to be in hydrostatic equilibrium, implying no lateral variations in density (Al-Attar and Tromp, 2014). For simplicity, we have also neglected rotational and self-gravitational effects in this study. Nevertheless, they are particularly important when coupling the evolving ice and ocean loads through the sea-level equation (Kendall et al., 2005) and will be incorporated in future work.

2.3 Weak form and spatial discretisation

To derive the finite-element discretisation of these governing equations, we first translate them into their weak form. By selecting appropriate function spaces that contain both solution fields and test functions, the weak form can be obtained by multiplying the equations by their test functions and integrating over the domain, Ω. For conservation of momentum, we use the (vector) test function, ϕ, to give

$$\begin{array}{}\text{(21)}& \mathrm{0}=\underset{\mathrm{\Omega }}{\int }\mathit{\varphi }\cdot \left(\mathrm{\nabla }\cdot {\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{L}}-B_{\mathit{\mu }}\mathrm{\nabla }\left({\mathit{\rho }}_{\mathrm{0}}\mathit{u}\cdot g{\widehat{\mathit{e}}}_k\right)-B_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{1}}g{\widehat{\mathit{e}}}_k\right)\mathrm{d}x,\end{array}$$

which we then simplify by multiplying both sides by −1 and splitting into three components

$$\begin{array}{}\text{(22)}& \mathrm{0}=S+H+D,\end{array}$$

where

$$\begin{array}{}\text{(23)}& S=\underset{\mathrm{\Omega }}{\int }-\mathit{\varphi }\cdot \mathrm{\nabla }\cdot {\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{L}}\mathrm{d}x,\end{array}$$

$$\begin{array}{}\text{(24)}& H=\underset{\mathrm{\Omega }}{\int }\mathit{\varphi }\cdot B_{\mathit{\mu }}\mathrm{\nabla }\left({\mathit{\rho }}_{\mathrm{0}}\mathit{u}\cdot g{\widehat{\mathit{e}}}_k\right)\mathrm{d}x,\end{array}$$

and

$$\begin{array}{}\text{(25)}& D=\underset{\mathrm{\Omega }}{\int }\mathit{\varphi }\cdot B_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{1}}g{\widehat{\mathit{e}}}_k\mathrm{d}x.\end{array}$$

At this stage, it becomes necessary to define the finite element spaces that will be used for spatial discretisation. Generally, for each component of displacement, we use Q2 finite elements on hexahedral meshes (i.e., the piecewise continuous tri-quadratic tensor product of quadratic continuous polynomials in each direction). In Sect. 3.2, however, we instead use prismatic meshes that are unstructured in the horizontal but layered in the vertical (facilitated by Firedrake's inbuilt `extruded' mesh functionality, this space uses the piecewise continuous bi-quadratic tensor product of a quadratic polynomial in the horizontal and a quadratic polynomial in the vertical; Bercea et al., 2016; McRae et al., 2016). For the deviatoric strain tensor (and hence the internal variable), since it is proportional to the gradient of displacement, we choose the discontinuous DG1 space (i.e., linear variations within each finite element cell and discontinuous jumps between cells) for each component. For purely radial variations in density, viscosity, and shear modulus, we choose the DG0 space (i.e., constant within a finite element cell but discontinuous between cells), while for laterally varying viscosity fields, we again select DG1 finite element functions.

Once these spaces have been chosen, we can integrate the divergence of the incremental Lagrangian stress term by parts within each element, K_n, to introduce weak boundary conditions and to move derivatives from the trial function to the test function according to

$$\begin{array}{}\text{(26)}& S=\sum _n\left(\underset{K_n}{\int }\mathrm{\nabla }\mathit{\varphi }:{\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{L}}\mathrm{d}x-\underset{\partial K_n}{\int }\mathit{\varphi }\cdot \left(\widehat{\mathit{n}}\cdot {\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{L}}\right)\mathrm{d}s\right),\end{array}$$

where $\widehat{\mathit{n}}$ is the outward pointing normal vector to an element edge, represented by ∂K_n. Given that test functions are continuous across cell edges and we assume that there is no jump in stress across material discontinuities, Eq. (26) simplifies to

$$\begin{array}{}\text{(27)}& \begin{array}{rl}S& =\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }\mathit{\varphi }:\left(\mathit{\kappa }\mathrm{\nabla }\cdot \mathit{u}\mathit{I}+\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}\mathit{d}\left(\mathit{u}\right)-\mathrm{2}\sum _i{\mathit{\mu }}_i{\mathit{m}}_i\left(\mathit{u}\right)\right)\mathrm{d}x\\ & +\underset{\partial {\mathrm{\Omega }}_{\mathrm{top}}}{\int }\mathit{\varphi }\cdot \left(B_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{load}}g{h}_{\mathrm{load}}\widehat{\mathit{n}}\right)\mathrm{d}s,\end{array}\end{array}$$

where we have substituted Eq. (10) for the constitutive relation and Eq. (17) for the boundary condition at the top of the domain, ∂Ω_top. For bottom (and side) boundaries where we specify no normal displacement (i.e., Eq. 18), we utilise two different strategies depending on domain geometry. In Cartesian domains, we apply the boundary condition strongly by modifying the discrete test and trial spaces, such that the surface integral vanishes. For spherical domains, we implement that condition weakly using the Symmetric Interior Penalty Galerkin method (Epshteyn and Rivière, 2007; Hillewaert, 2013). An example showing implementation of this approach for mantle convection problems in Firedrake can be seen in Davies et al. (2022).

Following Al-Attar and Tromp (2014), we combine the hydrostatic pre-stress and buoyancy terms

$$\begin{array}{}\text{(28)}& \begin{array}{rl}H+D& =B_{\mathit{\mu }}\sum _n\left(\underset{K_n}{\int }\mathit{\varphi }\cdot \left(\mathrm{\nabla }\left({\mathit{\rho }}_{\mathrm{0}}\mathit{u}\cdot g{\widehat{\mathit{e}}}_k\right)\right.\right.\\ & \left.\left.-\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{0}}\mathit{u}\right)g{\widehat{\mathit{e}}}_k\right)\mathrm{d}x\right),\end{array}\end{array}$$

which can be written in an explicitly symmetric form

$$\begin{array}{}\text{(29)}& \begin{array}{rl}H+D& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{B}_{\mathit{\mu }}\underset{\mathrm{\Omega }}{\int }{\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{\nabla }(\mathit{u}\cdot g{\widehat{\mathit{e}}}_k)\cdot \mathit{\varphi }+\mathit{u}\cdot \mathrm{\nabla }(\mathit{\varphi }\cdot g{\widehat{\mathit{e}}}_k)\right)\mathrm{d}x\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{B}_{\mathit{\mu }}\underset{\mathrm{\Omega }}{\int }{\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{\nabla }\cdot \mathit{u}g{\widehat{\mathit{e}}}_k\cdot \mathit{\varphi }+g{\widehat{\mathit{e}}}_k\cdot \mathit{u}\mathrm{\nabla }\cdot \mathit{\varphi }\right)\mathrm{d}x,\end{array}\end{array}$$

by assuming that the background stress is in a state of hydrostatic equilibrium, such that surfaces of constant background density and gravitational potential coincide (e.g., Eqs. B22–B29 in Appendix B of Al-Attar and Tromp, 2014).

2.4 Time discretisation

We now discretise in time using the implicit Backward Euler (BE) scheme. This choice allows us to take timesteps larger than the characteristic Maxwell time without compromising numerical stability. Such flexibility is particularly advantageous when the timescale of glacial loading is substantially slower than the Maxwell time – as is often the case in low-viscosity regions – thereby avoiding having to take prohibitively small timesteps in realistic simulations of glacial cycles (cf. Bailey, 2006; Lloyd et al., 2024). Applying the BE scheme, the evolution of each internal variable in Eq. (15) becomes

$$\begin{array}{}\text{(30)}& {\mathit{m}}_i^{n+\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\displaystyle \frac{\mathrm{\Delta }t}{{\mathit{\alpha }}_i}}}}\left({\mathit{m}}_i^n+{\displaystyle \frac{\mathrm{\Delta }t}{{\mathit{\alpha }}_i}}\mathit{d}\left({\mathit{u}}^{n+\mathrm{1}}\right)\right),\end{array}$$

where the superscript n refers to the previous timestep, n+1 is the next timestep, and Δt is the timestep duration. We then substitute this result into Eq. (27), yielding

$$\begin{array}{}\text{(31)}& \begin{array}{rl}S& =\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }\mathit{\varphi }:\left(\mathit{\kappa }\mathrm{\nabla }\cdot {\mathit{u}}^{n+\mathrm{1}}\mathit{I}+\mathrm{2}\sum _i{\displaystyle \frac{{\mathit{\eta }}_i}{{\mathit{\alpha }}_i+\mathrm{\Delta }t}}\left(\mathit{d}\left({\mathit{u}}^{n+\mathrm{1}}\right)-{\mathit{m}}_i^n\right)\right)\mathrm{d}x\\ & +\underset{\partial {\mathrm{\Omega }}_{\mathrm{top}}}{\int }\mathit{\varphi }\cdot \left(B_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{load}}g{h}_{\mathrm{load}}\widehat{\mathit{n}}\right)\mathrm{d}s,\end{array}\end{array}$$

where η_i is the non-dimensional viscosity of the internal variable m_i. Recombining Eqs. (29) and (31), the final system of equations is

$$\begin{array}{}\text{(32)}& \begin{array}{rl}& \underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }\mathit{\varphi }:\left(\mathit{\kappa }\mathrm{\nabla }\cdot {\mathit{u}}^{n+\mathrm{1}}\mathit{I}+\mathrm{2}\sum _i{\displaystyle \frac{{\mathit{\eta }}_i}{{\mathit{\alpha }}_i+\mathrm{\Delta }t}}\mathit{d}\left({\mathit{u}}^{n+\mathrm{1}}\right)\right)\mathrm{d}x\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{B}_{\mathit{\mu }}\underset{\mathrm{\Omega }}{\int }{\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{\nabla }({\mathit{u}}^{n+\mathrm{1}}\cdot g{\widehat{\mathit{e}}}_k)\cdot \mathit{\varphi }+{\mathit{u}}^{n+\mathrm{1}}\cdot \mathrm{\nabla }(\mathit{\varphi }\cdot g{\widehat{\mathit{e}}}_k)\right)\mathrm{d}x\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{B}_{\mathit{\mu }}\underset{\mathrm{\Omega }}{\int }{\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{\nabla }\cdot {\mathit{u}}^{n+\mathrm{1}}g{\widehat{\mathit{e}}}_k\cdot \mathit{\varphi }+g{\widehat{\mathit{e}}}_k\cdot {\mathit{u}}^{n+\mathrm{1}}\mathrm{\nabla }\cdot \mathit{\varphi }\right)\mathrm{d}x\\ & =-\underset{\partial {\mathrm{\Omega }}_{\mathrm{top}}}{\int }\mathit{\varphi }\cdot \left(B_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{load}}g{h}_{\mathrm{load}}\widehat{\mathit{n}}\right)\mathrm{d}s\\ & +\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }\mathit{\varphi }:\left(\mathrm{2}\sum _i{\displaystyle \frac{{\mathit{\eta }}_i}{{\mathit{\alpha }}_i+\mathrm{\Delta }t}}{\mathit{m}}_i^n\right)\mathrm{d}x,\end{array}\end{array}$$

where implicit terms involving the unknown displacement at the next time step, uⁿ⁺¹, have been collected on the left-hand side and explicit terms known from the current time step are on the right-hand side. Finding the new values of uⁿ⁺¹ requires solving a linear system at each timestep, which is solved in Firedrake through PETSc's comprehensive linear algebra library (Balay et al., 1997, 2025a, b). Since the system is symmetric, we employ the Conjugate Gradient method as the Krylov method. We also choose PETSc's algebraic multigrid package, GAMG, as a preconditioner since it has proven to be a robust and effective method for other large-scale problems in geodynamics (e.g. Davies et al., 2022).

In testing, we have found that the ratio of bulk-to-deviatoric strain terms (which shows up in the first integral term of Eq. 32) plays a crucial role in determining solver performance. In particular, when the ratio of bulk-to-shear modulus is high (as is the case when approximating incompressible materials), the resulting system can become poorly conditioned and leads to an increased number of iterations. This issue is exacerbated when timesteps greatly exceed the Maxwell time, as the relative contribution of the deviatoric strain term is reduced in that scenario. Since the ratio of bulk-to-shear modulus is typically modest (i.e., ∼2) for mantle rocks, we do not expect this behaviour to be problematic when using our software to run realistic GIA simulations. Nevertheless, when approximating incompressible materials by raising this ratio, there is an increase in computational cost. One way around this issue is to reformulate the system, such that the internal variable evolution equation is solved simultaneously with the momentum equation, yielding a coupled system whereby two equations are solved for two unknowns, uⁿ⁺¹ and mⁿ⁺¹. Under this formulation, the momentum equation becomes

$$\begin{array}{}\text{(33)}& \begin{array}{rl}& \underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }\mathit{\varphi }:\left(\mathit{\kappa }\mathrm{\nabla }\cdot {\mathit{u}}^{n+\mathrm{1}}\mathit{I}+\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}\mathit{d}\left({\mathit{u}}^{n+\mathrm{1}}\right)\right)\mathrm{d}x\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{B}_{\mathit{\mu }}\underset{\mathrm{\Omega }}{\int }{\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{\nabla }({\mathit{u}}^{n+\mathrm{1}}\cdot g{\widehat{\mathit{e}}}_k)\cdot \mathit{\varphi }+{\mathit{u}}^{n+\mathrm{1}}\cdot \mathrm{\nabla }(\mathit{\varphi }\cdot g{\widehat{\mathit{e}}}_k)\right)\mathrm{d}x\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{B}_{\mathit{\mu }}\underset{\mathrm{\Omega }}{\int }{\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{\nabla }\cdot {\mathit{u}}^{n+\mathrm{1}}g{\widehat{\mathit{e}}}_k\cdot \mathit{\varphi }+g{\widehat{\mathit{e}}}_k\cdot {\mathit{u}}^{n+\mathrm{1}}\mathrm{\nabla }\cdot \mathit{\varphi }\right)\mathrm{d}x\\ & =-\underset{\partial {\mathrm{\Omega }}_{\mathrm{top}}}{\int }\mathit{\varphi }\cdot \left(B_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{load}}g{h}_{\mathrm{load}}\widehat{\mathit{n}}\right)\mathrm{d}s\\ & +\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }\mathit{\varphi }:\left(\mathrm{2}\sum _i{\mathit{\mu }}_i{\mathit{m}}^{n+\mathrm{1}}\right)\mathrm{d}x,\end{array}\end{array}$$

and Eq. (30) is converted into its weak form by integrating over the domain and testing with ψ, a (tensor) DG1 finite element test function (i.e., numerical fields vary linearly within cells but can have discontinuities across cell boundaries), yielding

$$\begin{array}{}\text{(34)}& \underset{\mathrm{\Omega }}{\int }\mathit{\psi }\cdot {\mathit{m}}_i^{n+\mathrm{1}}\mathrm{d}x=\underset{\mathrm{\Omega }}{\int }\mathit{\psi }\cdot \left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\displaystyle \frac{\mathrm{\Delta }t}{{\mathit{\alpha }}_i}}}}\left({\mathit{m}}_i^n+{\displaystyle \frac{\mathrm{\Delta }t}{{\mathit{\alpha }}_i}}\mathit{d}\left({\mathit{u}}^{n+\mathrm{1}}\right)\right)\right)\mathrm{d}x.\end{array}$$

This formulation results in a system that is no longer symmetric due to the presence of coupling terms and we therefore use GMRES as the Krylov method. A key challenge is to find an efficient preconditioner for the coupled system. Within each iteration of GMRES, we employ a block Gauss-Seidel approach (facilitated by PETSc's FieldSplit functionality) whereby we independently solve Eq. (33) for uⁿ⁺¹ and Eq. (34) for mⁿ⁺¹, using the previous value for the other unknown. The matrix for both of these individual inner solves is symmetric and, for the momentum Eq. (33), we can use the same Conjugate Gradient and GAMG strategy as before, while the internal variable Eq. (34) is computationally simpler and can be solved using the Conjugate Gradient method preconditioned with successive over-relaxation.

Another benefit of the coupled formulation is that it provides a convenient means to model non-linear composite rheologies (e.g., the power-law formulation introduced in Eq. 16). In these cases, the viscosity is a function of stress (i.e., depends on both the internal variable and displacement). Since the coupled formulation explicitly computes both of these variables at every solution step, it is therefore comparatively trivial to employ PETSc's SNES functionality to implement a non-linear solve strategy, for which we apply Newton's method with a line-search algorithm.

In addition to easy implementation of non-linear rheologies, the coupled system ensures that convergence of the momentum solve is no longer dependent on timestep size and therefore offers greater robustness in cases involving either long timesteps or large bulk-to-shear modulus ratios. Nevertheless, it comes at the computational expense of additional iterations of the outer GMRES algorithm to maintain coupling. We therefore generally prefer to use the substitution method in simulations that use compressible, linear rheologies.

We also note that an alternative, more common approach to solve the momentum equation in the incompressible limit is to reintroduce pressure as a Lagrange multiplier in order to enforce the conservation of volume constraint $\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0}$ (e.g. Zhong et al., 2003). The resulting saddle-point system can be effectively solved using a Schur complement factorisation strategy (e.g. Davies et al., 2022), which we have implemented during code development but do not discuss further here since it is not applicable to real-Earth simulations with realistic levels of compressibility.

To summarise this section, we have outlined the viscoelastic loading equations and discretisation as currently implemented in G-ADOPT. In particular, by adopting the internal variable formulation of viscoelasticity (Al-Attar and Tromp, 2014; Crawford et al., 2017, 2018), we are able to accommodate a broad class of complex rheologies with relative ease in comparison to many of the traditional, Maxwell-based approaches. Additionally, our implementation of an implicit timestepping scheme enables stable integration even when timesteps are significantly longer duration than the Maxwell time. This flexibility should prove particularly advantageous in longer glacial loading scenarios, where low-viscosity regions that have Maxwell times on the order of a year would otherwise impose prohibitively small timestep constraints.

3 Numerical Experiments and Results

We next undertake a series of benchmarking experiments to assess the accuracy, robustness, and flexibility of our code.

3.1 Analytical comparisons

In this section, we assess the accuracy of G-ADOPT in solving simple viscoelastic loading problems that have analytical solutions for both compressible and incompressible cases. Such benchmarks provide a stringent check as they also allow us to test whether, as the temporal and spatial resolution is refined, the error decreases at the expected theoretical rate.

For the first scenario, we consider the vertical surface displacement, u_r, of a 2D Maxwell viscoelastic half-space that is subjected to emplacement of an instantaneous, sinusoidal surface load. By rearranging Eq. (2b) from Cathles (2024), we can simulate a loading, rather than an unloading scenario. The analytical solution is

$$\begin{array}{}\text{(35)}& u_r=F_{\mathrm{0}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+f_e\mathit{\alpha }/\mathit{\tau }}}+{\displaystyle \frac{\mathrm{1}-{\mathrm{e}}^{-t/(\mathit{\tau }+f_e\mathit{\alpha })}}{\mathrm{1}+f_e\mathit{\alpha }/\mathit{\tau }}}\right)\cos \left(kx\right),\end{array}$$

where F₀ is the load thickness (where its density is equal to the mantle density), $\mathit{\alpha }=\frac{\mathit{\eta }}{\mathit{\mu }}$ is the Maxwell time, η is the viscosity, μ is the shear modulus, and τ is the viscous relaxation timescale given by

$$\begin{array}{}\text{(36)}& \mathit{\tau }={\displaystyle \frac{\mathrm{2}k\mathit{\eta }}{\mathit{\rho }g}},\end{array}$$

where k is the wavenumber, ρ is the density, and g is the gravitational acceleration. The compressibility parameter, f_e, is defined as

$$\begin{array}{}\text{(37)}& f_e={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{e}}+\mathrm{2}\mathit{\mu }}{{\mathit{\lambda }}_{\mathrm{e}}+\mathit{\mu }}},\end{array}$$

where λ_e is the (first) elastic Lamé parameter, given by

$$\begin{array}{}\text{(38)}& {\mathit{\lambda }}_{\mathrm{e}}=\mathit{\kappa }-{\displaystyle \frac{\mathrm{2}\mathit{\mu }}{\mathrm{3}}},\end{array}$$

where κ is the bulk modulus. When the material is incompressible, then f_e→1.

For consistency with the analytical solution, both compressible buoyancy and advection of hydrostatic pre-stress must be neglected. To ensure correct asymptotic behaviour in the viscous limit, we must include a surface integral term at Earth's surface, ${\int }_{\partial {\mathrm{\Omega }}_{\mathrm{top}}}\mathit{\varphi }\cdot \widehat{\mathit{n}}{B}_{\mathit{\mu }}{\mathit{\rho }}_{\mathrm{0}}g{u}_k\mathrm{d}s$, where u_k is the component of non-dimensional displacement in the ${\widehat{\mathit{e}}}_k$ direction. This term is equivalent to the surface term that would be obtained after integrating Eq. (24) by parts (e.g. Zhong et al., 2003), and is similar to the free-surface feedback encountered in mantle convection (Kramer et al., 2012).

We adopt a mantle viscosity of $\mathit{\eta }=\mathrm{1}\times {\mathrm{10}}^{\mathrm{21}}$ Pa s and shear modulus of $\mathit{\mu }=\mathrm{1}\times {\mathrm{10}}^{\mathrm{11}}$ Pa, yielding a Maxwell time of α≈320 years. For the compressible case, the bulk modulus is $\mathit{\kappa }=\mathrm{2}\times {\mathrm{10}}^{\mathrm{11}}$ Pa, yielding f_e=1.43. Density and gravitational acceleration are set to 4500 kg m⁻³ and 10 m s⁻², respectively. The applied load has a height of 1 km and wavelength of 375 km (equivalent to $\frac{\mathrm{1}}{\mathrm{8}}$ of the domain depth), leading to a viscous relaxation time of τ≈24 000 years.

Figure 1 Analytical comparisons for 2D viscoelastic loading problem. (a) Peak displacement vs. time for compressible (red) and incompressible (blue) cases. Highlighted boxes indicate the time-windows shown in panels (b) and (c) and the legend applies to the whole figure. (b) Elastic-limit comparison of G-ADOPT predictions (dashed) with analytical solutions (solid) at three timesteps (Δt = 0.1α, 0.05α, and 0.025α). (c) Viscous-limit comparison for Δt=16α, 8α and 4α. (d) L₂-norm of the surface displacement error in the elastic limit as a function of timestep size (Δt=0.1α to 0.0015625α), demonstrating first-order convergence. (e) Time-integrated global L₂-norm error in the viscous limit as timestep size is refined (Δt=16α to 0.5α). Diagonal lines in (d) and (e) indicate expected first-order convergence behaviour.

Figure 1a shows the temporal evolution of peak surface displacement in the analytical solutions over a duration of 160 Maxwell times (i.e., approximately 2.2 viscous relaxation timescales). Both compressible and incompressible cases exhibit similar behaviour: an immediate elastic displacement is followed by a gradual approach towards isostatic equilibrium. Including compressibility leads to a slightly larger initial elastic deformation (∼19 m, as opposed to ∼13 m; Fig. 1b), but the viscous asymptote is the same for both scenarios (Fig. 1c).

To explore numerical accuracy as a function of timestep size, we fix the mesh resolution and perform a suite of simulations over a range of timestep durations. We use a 2D domain of depth 1 and width 0.5 in non-dimensional terms that generally contains 320×320 cells, except in the incompressible, long-duration case where this resolution was doubled to 640×640 cells to ensure that the time discretisation error dominates over any spatial discretisation errors. Since the Backward Euler timestep method is unconditionally stable and first-order accurate, these tests should ideally achieve a first-order convergence rate.

At short timescales where elastic deformation is expected to dominate, we undertake a single timestep with durations of Δt=0.1α down to 0.0015625α. To approximate incompressibility within our compressible viscoelastic formulation, we set the bulk modulus to a high value of $\mathit{\kappa }=\mathrm{1}\times {\mathrm{10}}^{\mathrm{15}}$ Pa (equivalent to f_e=1.0001). For both compressible and incompressible cases, as Δt reduces, we can see that our viscoelastic response gets systematically closer to the instantaneous elastic response from the analytical solutions at time t=0 (Fig. 1b) and that the surface L₂-errors of this difference reduce with first-order accuracy (Fig. 1d).

At long timescales where viscous deformation dominates, we run numerical simulations with a timestep duration that systematically varies from 16α down to 0.5α. Since the relative contribution of elastic compressibility to total displacement is much lower over long timescales, we use a less extreme bulk modulus of $\mathit{\kappa }=\mathrm{1}\times {\mathrm{10}}^{\mathrm{13}}$ Pa for the incompressible case (i.e., f_e=1.01). Once again, we observe that, as timestep size reduces, numerical simulations approach the two analytical solutions (Fig. 1c) and that their time-integrated global L₂-error reduces with first-order convergence (Fig. 1e).

Despite the physical simplicity necessitated by these analytical solutions, the results demonstrate that G-ADOPT accurately reproduces the analytical response for both compressible and incompressible Maxwell halfspaces and that its implicit time-stepping scheme achieves the expected first-order accuracy in time. In the following benchmarks, we next evaluate G-ADOPT against other finite element and semi-analytical codes using more realistic loading scenarios in a 3D domain.

3.2 3D Cartesian benchmarks with a Maxwell rheology

Weerdesteijn et al. (2023) have recently presented a valuable suite of test simulations created using three published GIA models: two finite element codes – Abaqus (Wu, 2004) and Aspect (Weerdesteijn et al., 2023) – and one semi-analytical code that assumes an axisymmetric spherical Earth, Taboo (Spada et al., 2004). A schematic of their model setup is shown in Fig. 2, for which the domain is a 3D Cartesian box measuring 1500×1500 km laterally and 2891 km in depth (and therefore can only be approximately represented in Taboo). Two loading scenarios were considered: (i) a short-duration event that has similar spatiotemporal scales to modern ice changes, whereby a 100 m-thick ice sheet linearly grows for 100 years before being held constant for another 100 years; and (ii) a long-duration scenario that is more similar to a glacial cycle, whereby 1 km of ice linearly accumulates over 90 kyr and then fully melts by 100 kyr, followed by a further 10 kyr of evolution in the absence of further loading changes. In both cases, the ice sheet is represented by a disc with a radius of 100 km.

Figure 2 Schematic of the 3D Cartesian benchmark from Weerdesteijn et al. (2023). (a) Model domain showing rheological layers, an optional low viscosity (LV) zone, and the extent of ice at the surface. Note schematic is not drawn to scale. (b) Temporal evolution of thickness of the ice disk in the short-loading scenario. (c) Same for the long-loading scenario.

The default Earth structure follows the 1D layered model of Spada et al. (2011), composed of a lithosphere, upper mantle, transition zone, and lower mantle (Table 2). To assess sensitivity to lateral viscosity variations, a second 3D-Earth structure also incorporates a cylindrical anomaly with a radius of 100 km and a reduced viscosity of 1×10¹⁹ Pa s directly beneath the ice sheet, extending from 70–170 km depth. Since Taboo assumes spherical symmetry, it is therefore limited to modelling 1D cases only.

Table 2 Rheological parameters for 1D depth profiles following Weerdesteijn et al. (2023) and Spada et al. (2011). To approximate infinite viscosities in the lithosphere, our simulation assumes $\mathit{\eta }=\mathrm{1}\times {\mathrm{10}}^{\mathrm{40}}$ Pa s.

The configuration in G-ADOPT is set up to closely imitate Aspect's implementation of the benchmark, which includes: neglecting self-gravity; using a constant gravitational acceleration of 9.815 m s⁻²; setting viscosity in the elastic lithosphere to 1×10⁴⁰ Pa s (thereby ensuring negligible viscous deformation of this layer over the course of the simulation); and exploiting axisymmetry to model only one quarter of the domain with no-normal-flow boundary conditions imposed on its side walls and base. Ice loading is implemented using a smooth hyperbolic tangent profile with a 1 km “rollover” length-scale and the density of ice is assumed to be 931 kg m⁻³. We use Gmsh to construct a horizontally unstructured mesh, refined to 5 km under the ice and coarsened to 200 km in the far field (Geuzaine and Remacle, 2009). Vertically, each rheological layer is discretised with a fixed number of DG0 elements (e.g., 10 per layer in the default setup), providing constant values within elements but permitting sharp boundaries between layers. Timesteps are 10 and 1000 years for the short- and long-loading scenarios, respectively. For reference, Abaqus uses a similar meshing approach to G-ADOPT, with horizontal resolution of 5–200 km and 8 vertical cells per rheological layer; Aspect uses ∼6 km isotropic resolution in the top 100 km of the domain, transitioning to 50 km elsewhere; and Taboo uses a maximum spectral degree of 4096, equivalent to approximately 5 km horizontal resolution, and is semi-analytical in the vertical direction (Weerdesteijn et al., 2023). Since the benchmarks all assume incompressibility, we again approximate incompressibility by setting the bulk modulus to 1×10¹⁴ Pa (i.e., 1000 times larger than the reference shear modulus). Note that, to ensure robustness in these incompressible simulations, we solve the coupled system (i.e., Eqs. 33 and 34) for reasons outlined in Sect. 2.

Figure 3 Comparison of maximum vertical surface displacements beneath the ice load for the 3D Cartesian box benchmark scenarios of Weerdesteijn et al. (2023). (a) Short-duration loading scenario on a 1D Earth model. Dotted red = Abaqus; dash-dotted purple = Aspect; dashed orange = Taboo; and solid blue = G-ADOPT. Inset shows a zoom of the region outlined by the black-dashed box. (b) Same for the long-duration loading scenario on a 1D Earth. (c) Short-duration loading scenario on a 3D Earth containing a low-viscosity patch in the upper mantle beneath the ice sheet. (d) Same for the long-duration loading scenario on a 3D Earth; grey dashed line = G-ADOPT result on a 1D Earth (panel b) plotted to highlight impact of the low-viscosity patch. Note that results for Abaqus, Aspect and Taboo are digitised from Weerdesteijn et al. (2023) and that some of the four scenarios illustrated were not tested in those codes. Also note the different y-axis extents for the short- versus long-duration scenarios.

Figure 3 compares maximum vertical displacements beneath the ice load across all four codes in each of the four benchmark combinations of Earth model and loading scenario. The upper row shows results for 1D Earth structure, while the lower row includes the low-viscosity cylindrical anomaly. Agreement is generally excellent in all cases where cross-code comparisons are possible (Fig. 3a–c), affirming consistency with G-ADOPT's implementation. In the final panel (Fig. 3d), we present results from G-ADOPT for the long-loading scenario with 3D viscosity – a configuration not illustrated in Weerdesteijn et al. (2023). In comparison to the same loading scenario on a 1D Earth, the presence of a low-viscosity patch accelerates the surface response consistent with expectations, but only modestly: the maximum vertical displacement differs by just ∼2 % relative to the 1D case at 90 kyr, compared to a ∼60 % difference in the short-duration loading scenario. This result is consistent with the fact that, even for the 1D Earth model, the ∼100 kyr loading timescale is substantially longer than the ∼320 year Maxwell time of the upper mantle (N.B., the equivalent Maxwell time for the low-viscosity patch is only ∼3 years).

Figure 4 Sensitivity analysis for G-ADOPT using the scenario explored in Fig. 3b consisting of long-duration ice loading on a 1D Earth model. (a) Impact of horizontal resolution. Solid blue = default spacing (5 km beneath ice sheet, 200 km in far field); crosses = half of default resolution; circles = one quarter of default resolution. (b) Impact of vertical resolution. Solid blue = default case with 40 vertical cells across the model domain (i.e., 10 per layer); crosses = 20 cells; circles = eight cells; triangles = four cells. (c) Impact of timestep size. Solid blue = default of 1 kyr; crosses = 2 kyr; circles = 5 kyr; triangles = 10 kyr. (d) Impact of transitioning from incompressible to a realistic Earth-like level of compressibility by modifying the ratio of bulk-to-shear modulus. Solid blue = 1000; crosses = 100; circles = 10; triangles = 1.94.

To further assess the robustness of G-ADOPT, we conduct a resolution and parameter sensitivity analysis for the scenario involving long-duration loading of a 1D Earth model. Horizontal mesh sensitivity tests demonstrate that even as few as five elements beneath the ice-covered region (i.e., 20 km horizontal resolution) are sufficient to accurately capture the pattern of deformation (Fig. 4a). Changing the vertical resolution has a bigger effect (Fig. 4b). Halving the default values to five elements per layer (i.e., 20 vertical elements in total) results in practically indistinguishable results, but reducing to two cells per layer produces peak differences of ∼1 %, which increases to ∼6 % when only using one cell per layer. Regarding the impact of timestep size, these tests neatly demonstrate the effectiveness of our implicit time-stepping scheme. Even when the timestep is increased to 10 kyr, which is more than 30 times greater than the Maxwell time of the upper mantle in these models, the solution remains stable and accurate throughout the loading phase and produces only minor errors during the unloading phase (Fig. 4c). For this case, the computational cost per timestep is doubled because the number of outer iterations per timestep increases, but this still leads to an effective speed up of a factor of five in comparison to the 1 kyr simulation. By way of comparison, Aspect used significantly shorter timesteps (e.g., 50 years) to achieve comparable results in this benchmark (Weerdesteijn et al., 2023).

We can also use this test scenario to demonstrate the impact of realistic levels of compressibility by changing the ratio of bulk-to-shear modulus (Fig. 4d). A value of 1000 closely approximates incompressibility and turns out to be indistinguishable from a value of 100. By the time it is lowered to 1.94, there is a ∼20 % increase in maximum vertical displacement, reflecting the larger elastic deformation caused by loading compressible materials. This result further underscores the importance of including compressibility in GIA simulations (Wolf, 1985; Tanaka et al., 2011; A et al., 2013). However, it also requires taking greater care during model setup to ensure that a density inversion (with respect to depth) does not occur at any point during the simulation. Such a scenario leads to convective instabilities that can grow, depending on the timescales, and have been widely documented in GIA literature (Vermeersen and Mitrovica, 2000; Huang et al., 2023).

Table 3 Weak scaling analysis for the long-duration loading case on a 1D Earth using a structured mesh. Default iteration counts and timings use a bulk-to-shear modulus ratio of 100 (equivalent values for a ratio of 1.94 are given in brackets). DOFs = Degrees of freedom.

To test the scalability of G-ADOPT, we have undertaken a weak scaling analysis on Gadi (Australia's highest performance CPU-based supercomputer) using the long-duration loading scenario on a 1D Earth and a structured mesh. We use a number of cores ranging from 104 up to 6656 (i.e., from 1 to 64 Intel Xeon Sapphire Rapids nodes) and aim to keep the number of degrees of freedom per core constant at ∼30 000. The results show that, with a bulk-to-shear modulus ratio of 100, increasing the problem size by a factor of 64 causes the solve time per timestep to increase by ∼30 % (Table 3), which is consistent with previous scaling results for mantle convection modelling within Firedrake (Davies et al., 2022). As discussed in Sect. 2.4, for the case with a smaller bulk-to-shear modulus ratio of 1.94, conditioning of the system improves and the number of iteration counts is significantly reduced. However, the overall scaling performance is slightly worse with the solve time per timestep effectively doubling as the problem size increases by a factor of 64.

In summary, for a 3D Cartesian box with potential presence of lateral viscosity variations, G-ADOPT accurately reproduces benchmark results from published codes. It also exhibits desirable performance characteristics, including efficient timestepping and stability across a range of resolutions. These results validate the core solver and lay the groundwork for investigating more complex loading scenarios, including those involving compressibility and more complex rheologies.

3.3 Extending to complex rheologies including power-law and transience

In this section, we demonstrate the flexibility of the internal variable formulation to model both power-law and transient rheologies by extending the short- and long-duration loading 1D-Earth benchmarks of Sect. 3.2. Compressibility is implemented by adopting a bulk-to-shear modulus ratio of 1.94 (corresponding to a Poisson’s ratio of 0.28), as that value is standard in Abaqus implementations (Huang et al., 2023). As discussed before, the impact of compressibility alone is to increase the magnitude of the peak displacement through time in comparison to the incompressible Maxwell case, consistent with expectations (Fig. 5).

For composite power-law rheology, we use the same reference viscosity as before and adopt a transition stress of τ_T=0.2 MPa and exponent of n=3 (consistent with case 1 of Kang et al., 2022). All other aspects of the setup including grid resolution, timestep size, boundary conditions, and loading configuration are identical to those in Sect. 3.2. In the short-duration loading scenario, the maximum vertical displacement increases with respect to the compressible Maxwell model, becoming more obvious after the load stops growing at 100 years and reaching ∼13 % more by 200 years (Fig. 5a). In the long-duration loading scenario, the same situation is observed during loading, although the peak vertical displacement is only ∼1.7 % larger at the time of maximum loading at 90 kyr. During unloading, however, the composite power-law rheology recovers substantially faster than its Maxwell equivalent, reaching 20 % of its peak displacement by 103 kyr in comparison to 109 kyr for the Maxwell case. We also verified that setting the transition stress to a high value of τ_T=10 MPa (i.e., ensuring that the power-law component is not activated) gave equivalent results to the Maxwell case.

A Burgers rheology can be conceptually considered as a standard Maxwell body (composed of a spring and dashpot connected in series) connected in series to a single Kelvin-Voigt element (a spring and dashpot connected in parallel). It requires using two internal variables to represent the two characteristic relaxation timescales of the system (Crawford et al., 2017, 2018). To maintain the same elastic response as in the compressible Maxwell case, we halve the shear modulus of each internal variable element, thereby preserving the total effective modulus. Similarly, we halve the viscosity of the long-term viscous element, so that the overall relaxation timescale remains consistent with the Maxwell cases presented earlier. We then explore two different cases in which the viscosity of the element controlling short-term relaxation is set to either one-half or one-tenth that of the long-term element.

Figure 5 Extension of the 3D Cartesian benchmark to compressible non-Maxwell rheologies. (a) Short-duration loading scenario on a 1D Earth model. Solid blue = incompressible Maxwell rheology previously shown in Fig. 3a; solid red = equivalent for compressible Maxwell rheology (bulk-to-shear modulus ratio of 1.94) previously shown in Fig. 4d; dotted black = composite power-law rheology with transition stress τ_T=0.2 MPa and exponent n=3; dashed purple = Burgers rheology with short-term viscosity a factor of 2 lower than the long-term viscosity; dash-dotted orange = same for a factor of 10 viscosity reduction. (b) Same for the long-duration loading scenario.

In both cases, adding in a transient element with lower viscosity amplifies the response in comparison to the compressible Maxwell case (as was found for the power-law rheology), particularly where rapid changes in the load interact strongly with the shorter relaxation timescale (Fig. 5). This behaviour is most obvious in the short-duration loading scenario, but it is also visible in the unloading phase of the long-duration scenario. We also verified that setting the same viscosity for both internal variables (i.e., reducing the Burgers model to a single timescale) reproduces the original Maxwell response, confirming the internal consistency of our implementation.

It is noteworthy that, for short-duration loading, displacement curves for the compressible Maxwell case, and the compressible power-law and Burgers simulations (Fig. 5a), exhibit comparable features to those produced for the incompressible Maxwell model with a low-viscosity patch (Fig. 3c). Although the underlying deformation mechanisms differ, this similarity illustrates the difficulty of disentangling the effects of local viscosity variations, compressibility, and non-Maxwell rheologies (e.g., Pan et al., 2024). Nevertheless, it also highlights the flexibility of the G-ADOPT framework and reinforces its potential future value to these debates.

In summary, we have demonstrated G-ADOPT's ability to model compressible power-law and transient viscoelastic rheologies using the internal variable formulation and the important impact of these phenomena on GIA predictions, particularly during rapid loading changes. In the next section, we demonstrate extension of these simulations to domains with spherical geometry.

3.4 Extension to spherical geometries

A significant strength of G-ADOPT is that it is relatively trivial to change the mesh to suite a range of different domain geometries. We demonstrate this aspect by extending our previous simulations from the 3D Cartesian box to a 3D spherical Earth. We adopt a model domain from the benchmark of Spada et al. (2011), which modelled an incompressible Maxwell rheology with purely radial viscosity variations. While we cannot repeat this benchmark as we have not yet implemented self-gravity in G-ADOPT, we have extended it to a more complex rheology by using a compressible, Burgers (i.e., transient) rheology that also contains lateral viscosity variations.

To generate the mesh, we first use Firedrake's inbuilt functionality to create a cubed-sphere surface mesh with 24 576 horizontal elements (equivalent to 6 refinements). We then “extrude” the mesh in the vertical direction to create 40 vertical layers (10 per rheological layer), as in Sect. 3.2. A couple of additional changes with respect to the Cartesian setup are that we provide a rotational null space to the Firedrake solver options to improve iterative solver performance and provide a weak implementation of the no-normal-flow boundary condition on the core-mantle boundary. We retain the same radial profiles of density, shear modulus, and long-term (i.e., steady-state) viscosity as in Spada et al. (2011), but these viscosities are further modulated by a ($\mathrm{\ell }=\mathrm{4},m=\mathrm{1}$) spherical harmonic function that introduces lateral viscosity variations spanning two orders of magnitude (Fig. 6). For the short-term viscosity in the Burgers model, we assign a value that is a factor of 10 times smaller, while compressibility is set using a bulk-to-shear modulus ratio of 1.94, corresponding to a Poisson's ratio of approximately 0.28. A 1 km-thick ice disc with a radius of 10° latitude is placed on the North Pole and applied instantaneously at the start of the simulation via a Heaviside step function. Ice density is 931 kg m⁻³ and gravity is constant: g=9.815 m s⁻². A timestep of 50 years is used to resolve rapid deformation following the instantaneous ice loading event.

Figure 6 Viscoelastic deformation of a 3D spherical domain with an elastically compressible, Burgers rheology with lateral viscosity variations following instantaneous application of a disk load at the North Pole. Snapshots are at 1, 2, 5 and 10 kyr after load is applied. Red/blue = long-term, steady state viscosity (lithospheric viscosity of ∼10⁴⁰ Pa s saturated in dark blue); greyscale = magnitude of surface displacement (surface warping exaggerated by a factor of 1500 for improved visibility); black arrows = velocity vectors.

Figure 6 shows slices through the spherical domain at times of 1, 2, 5 and 10 kyr after the load is applied. Early on in the simulation, deformation velocities are large and concentrated in the upper mantle. Despite the axisymmetric nature of the load, there is a bias in their direction, with greater movement towards the neighbouring region of low viscosity. There is also a higher amplitude, shorter wavelength peripheral bulge developing on this side of the ice sheet in comparison to the direction of the highest viscosity. After 5 kyr, the remaining deformation is now concentrated on the higher viscosity side of the domain and has much lower velocity, consistent with expectations that the higher viscosity region has a longer Maxwell time and therefore takes longer to reach isostatic equilibrium. By 10 kyr, the majority of deformation has finished and the initially highly asymmetric surface displacement has gradually restored to a more symmetric pattern. At equilibrium, the surface displacement should match the symmetry of the applied load (since we are still using radially symmetric elastic structure).

This simulation is designed to showcase the relative ease in which G-ADOPT can be switched between different model domain geometries. It also simultaneously demonstrates its ability to model more complex, Earth-like problems that include compressibility, lateral viscosity variations, and transient rheologies. We believe that these aspects should render it a valuable future tool for modelling viscoelastic surface loading problems in future work, such as those routinely encountered in GIA. In the next section, we illustrate another key strength of G-ADOPT – namely its ability to compute automatic adjoints.

3.5 Adjoint twin experiments

As emphasised in the introduction, a major motivation for developing a new viscoelastic deformation code in G-ADOPT is the availability of an automatically-derived discrete adjoint that can be used to invert any observational constraints for unknown Earth structure and/or loading histories. Here, we use a synthetic 2-D annulus domain loaded with ice at its surface and demonstrate this adjoint capability by performing gradient-based optimisations for ice thickness and mantle viscosity. While this example is intentionally kept simple, our goal is to establish a clear and controlled foundation for tackling more complex and realistic inverse problems in future work.

3.5.1 Synthetic forward problem

Figure 7 shows the setup of the 2D annulus used in our adjoint twin experiments, including the configuration of viscosity fields and ice loads. In these simulations, we assume a Maxwell rheology and consider two distinct simulations that differ only in their viscosity structure. The first adopts the same purely radial viscosity profile used in our previous benchmark cases (Table 2) except that it uses a less stiff lithospheric viscosity of 1×10²⁵ Pa s. The second incorporates various lateral viscosity variations (LVV) covering two orders of magnitude that are superimposed on the same radial background. These variations are defined by Gaussian perturbations: three low-viscosity regions (i.e., 1×10²⁰ Pa s) located at equatorial parts of the core-mantle boundary and in the upper mantle near the South Pole; and two high-viscosity features (i.e., 1×10²² Pa s) in the upper mantle of the northern hemisphere. Elastic parameters and density vary only radially and follow the values in Table 2. The mesh comprises 360 horizontal cells and 80 vertical cells, resulting in 20 cells per rheological layer.

Figure 7 Setup of synthetic forward problem within a 2D annulus domain. (a) Axisymmetric viscosity variations. Blue external ring = thickness of the ice load instantaneously applied at start of simulation. (b) Same for the case with lateral viscosity variations (LVV).

The loading scenario consists of two ice sheets applied instantaneously at the start of the simulation via a Heaviside function. A larger ice sheet is centred over the South Pole and has a thickness of 2 km and half-width of 20°. A second, smaller ice sheet is offset clockwise from the north pole by 25° and has a thickness of 1 km and half-width of 10°. As in previous simulations, the density of ice is 931 kg m⁻³ and gravitational acceleration is constant, g=9.815 m s⁻². Simulations are run for 10 000 years with a 50-year timestep.

Snapshots of the surface evolution at 3 kyr are shown in Fig. 8. Looking first at the axisymmetric case, the larger South Pole ice sheet (centred at $\mathit{\theta }=-\mathrm{90}\mathit{°}$) shows a peak downwards displacement of ∼390 m, which is two and a half times larger than the displacement under the smaller northern hemisphere ice sheet (centred at $\mathit{\theta }=+\mathrm{65}\mathit{°}$; Fig. 8a). Intervening far-field locations show fairly consistent surface uplift with an amplitude of ∼40 m and ongoing rate of 4–6 mm yr⁻¹ (Fig. 8b). At this point in the simulation, the main difference for the case with LVV occurs in the low-viscosity patch beneath the large ice sheet (between −90 and −120°). Here, we see that the shape of the subsidence is markedly more asymmetric, reaching a peak downwards displacement of −420 m and accompanied by a narrower forebulge with a height of ∼40 m. The rates of subsidence and uplift are also lower in comparison to the axisymmetric case, consistent with the fact that it takes less time to reach isostatic equilibrium in the lower viscosity region. In contrast, beneath the smaller ice sheet, the locally higher upper-mantle viscosity causes a delay in the surface response relative to the axisymmetric case.

Examination of the tangential displacements tells a similar story, with the biggest difference occurring above the low-viscosity zone where the surface of the axisymmetric simulation is displaced away from the ice sheet by ∼200 m while that of the LVV simulation is displaced by ∼280 m (Fig. 8c). Similar differences occur in the rates of tangential motion (Fig. 8d). Formally examining the root-mean-squared difference in surface displacement between the two cases reveals that the radial value peaks at 0.58 m by 2.0 kyr before gradually declining, while the tangential value peaks at 0.89 m and slightly later at 2.4 kyr (Fig. 8e). Meanwhile the differences in surface velocity between the two simulations are largest at the first timestep, reaching 1.15 and 1.44 mm yr⁻¹ for radial and tangential, respectively, before rapidly falling to ∼0.1 mm yr⁻¹ by 2 kyr (Fig. 8f).

Figure 8 Comparison of surface displacement and surface velocity through time for axisymmetric and lateral viscosity variation (LVV) viscosity cases. Radial component of surface displacement (a) and surface velocity (b) at 3 kyr. (c)–(d) are corresponding tangential components. Note the large ice sheet is centred at −90° and the small ice sheet at 65°. Spatially averaged misfit between (LVV) and axisymmetric viscosity case through time for surface displacement (e) and surface velocity (f).

3.5.2 Inversion for initial load

Armed with synthetic “observations” of surface displacement through time for these two forward simulations, we can now demonstrate an adjoint-based optimisation. In our first demonstration, we retain knowledge of the true viscosity structure and only invert for the unknown initial ice thickness. It is important to note that this test and those that follow assume perfect observations (i.e., complete spatiotemporal data coverage with no uncertainties) and is therefore highly idealised with respect to real-world datasets.

The ice thickness that must be reconstructed is referred to as the (discrete) control vector, c, and is defined using a P1 finite element function (i.e., linear across elements and continuous between them) on the 1D surface of the 2D annulus domain. The misfit objective function, J(c), is written as

$$\begin{array}{}\text{(39)}& \begin{array}{rl}J\left(\mathit{c}\right)& ={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{\Delta }t}C}}\\ & \sum _n^{N_{\mathrm{\Delta }t}}\underset{\partial {\mathrm{\Omega }}_{\mathrm{top}}}{\int }\left({\displaystyle \frac{\left({\mathit{u}}^n\right(\mathit{c})-{\mathit{u}}_{\mathrm{obs}}^n{)}^{\mathrm{2}}}{{u_{\mathrm{scale}}}^{\mathrm{2}}}}+{\displaystyle \frac{\left({\mathit{v}}^n\right(\mathit{c})-{\mathit{v}}_{\mathrm{obs}}^n{)}^{\mathrm{2}}}{{v_{\mathrm{scale}}}^{\mathrm{2}}}}\right)\mathrm{d}s,\end{array}\end{array}$$

where N_Δt is the number of timesteps, C is the surface “area” (i.e., circumference) of the 2D annulus, n is the current timestep, u and v are current solution values for the (non-dimensional) displacement and velocity fields, respectively, ${\mathit{u}}_{\mathrm{obs}}^n$ and ${\mathit{v}}_{\mathrm{obs}}^n$ are the corresponding “observations” from the target forward simulation, and $u_{\mathrm{scale}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ and $v_{\mathrm{scale}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ are factors used to scale the relative sizes of the misfit coming from displacements versus velocities. J(c) is often referred to as the reduced functional, since evaluation of J also requires solution of the forward model to find u and v, which are in turn a function of c.

Once the forward problem has been specified, G-ADOPT leverages Pyadjoint to automatically derive the adjoint problem and compute the gradient of the misfit function with respect to the control vector. To verify the accuracy of these gradients, a Taylor remainder convergence test is performed (Farrell et al., 2013; Ghelichkhan et al., 2024). It is defined as

$$\begin{array}{}\text{(40)}& \left|J(\mathit{c}+h\mathit{\delta }\mathit{c})-J\left(\mathit{c}\right)-h\mathrm{\nabla }\mathit{c}\cdot \mathit{\delta }\mathit{c}\right|⟶\mathrm{0}\text{ at }O\left(h^{\mathrm{2}}\right),\end{array}$$

where c is the control vector (ice thickness), δc is a random perturbation direction, h≪1 is a small scalar, and J(c) is the reduced objective functional. The left-hand side represents the second-order Taylor remainder. As the size of the scalar is modified and if the gradient is correctly implemented, the residual should converge to zero at a rate of O(h²), which is indeed the case (Table 4).

Table 4 Taylor test results for the derivative of the misfit objective function with respect to ice thickness, for h = 0.01.

Following the approach of Ghelichkhan et al. (2024), the computed gradient is passed to ROL and we use the Lin-Moré trust-region algorithm (Lin and Moré, 1999). This algorithm permits the placing of constraints on the control variables, such as upper and lower bounds, in this case allowing enforcement of only non-negative ice thicknesses. While it would be possible to add further regularisation terms to the objective function, for example ensuring that the ice is smooth in shape, they are not considered in this experiment. The initial guess is that there is zero ice thickness across the domain. For the “observational data”, we use the results of the synthetic forward simulation with LVV from Sect. 3.5.1 and invert it with both the correct LVV viscosity structure and the incorrect axisymmetric one.

Figure 9 Adjoint-based optimisation for initial ice load. (a) Normalised ice thickness after 20 iterations. Black dashed line = true target load; blue = inversion using correct LVV mantle structure; red = inversion using incorrect axisymmetric mantle structure. (b) Objective (i.e., misfit) function value as a function of iteration number for both cases.

After only 20 iterations, the recovered ice thickness using the correct LVV viscosities is visually indistinguishable from the true target field and corresponds with the objective function (i.e., misfit in surface displacement and velocity) reducing by nearly seven orders of magnitude (Fig. 9). In contrast, misfit for the incorrect viscosity structure initially drops by two orders of magnitude before stalling. Examining the resulting ice thicknesses shows that the recovered large ice sheet is asymmetric, with ∼10 % extra ice added to compensate for absence of the low-viscosity upper mantle region near −110°. Similarly, the smaller ice sheet is ∼5 % too thin, owing to the locally reduced viscosity in the axisymmetric model. Additionally, the optimisation introduces a thin, tapering layer of spurious ice 30–50 m-thick extending either side of the small ice sheet and into the large ice sheet (i.e., between −50 and 150°). This behaviour indicates the great value of geological constraints on ice extent, such as terminal moraines and cosmogenic exposure ages, when attempting this approach with real datasets.

Regarding the computational cost of automatic adjoint calculation, we find that Pyadjoint is very efficient. Annotation of the forward model for this 2D case takes 3 % of the time required to run the forward model, while the adjoint model is only 6 % slower than the forward model with taping (Table 5). This nearly matches the expected one-to-one correspondence between the cost for the forward and adjoint model for a linear forward problem (e.g. Farrell et al., 2013). We anticipate further improvements in the adjoint-to-forward-efficiency ratio for fully 3D problems, as the cost associated with solving the linear systems increases.

Table 5 Runtimes for the forward and adjoint models.

This experiment demonstrates the benefits of automatic adjoint derivation in Pyadjoint and its effective use for recovering initial loads when the viscosity structure is known. It also highlights the biases that are introduced by incorrect rheological assumptions.

3.5.3 Inversion for viscosity structure

We now turn to the complementary inversion problem: assuming the distribution of the initial load is known, can we recover an unknown mantle viscosity structure? The same objective function (Eq. 39) is used, which quantifies misfit between simulated and target surface displacements and velocities at all timesteps. The target “observational” data is again drawn from the synthetic forward simulation using viscosity structure with lateral viscosity variations and does not contain any additional noise or uncertainties.

The viscosities in this model span five orders of magnitude. In order to improve convergence rates, we found that it is helpful to apply a logarithmic rescaling according to

$$\begin{array}{}\text{(41)}& \mathit{\eta }={\mathit{\eta }}_{\mathrm{0}}{\mathrm{10}}^{{\mathit{\eta }}_{\mathrm{control}}},\end{array}$$

where η is viscosity in the simulation, η₀ is a radial viscosity profile, and η_control is the control field of lateral variations over which the optimisation is performed. This formulation constrains the control values to a smaller, more uniform range. η₀ is initialised using the same radial viscosity profile (although it could be any other) using a DG0 finite element function. η_control is defined using a continuous P1 element and the field is initially set to zero everywhere in the domain, so that the starting guess for viscosity is equivalent to viscosities in the axisymmetric case. As before, the inversion proceeds without regularisation within the objective function. A Taylor remainder convergence test confirms the accuracy of the computed gradients (Table 6).

Table 6 Taylor test results for the derivative of the misfit objective function with respect to viscosity.

Figure 10 Optimisation of viscosity structure given a known ice load. Each row shows the state of the inversion after 0, 10, 50, and 100 iterations. Left-hand column = current inference of viscosity field; central column = misfit with respect to the true viscosity structure in control-space; right-hand column = adjoint sensitivity kernel, where positive values in red indicate that decreasing viscosity reduces the misfit and negative values in blue indicate that it should be increased.

At iteration zero where the viscosity is still axisymmetric, the largest signals in the adjoint sensitivity kernel indicate that viscosity should be reduced beneath the large southern ice sheet and increased beneath the smaller northern ice sheet (top-right panel in Fig. 10). This pattern is consistent with true features of the target model with lateral viscosity variations. By iteration 10, the inversion successfully recovers the bulk of the main low-viscosity zone in the upper mantle of the southern hemisphere, along with emerging high-viscosity features in the north. It is at this stage that the adjoint sensitivity kernel also begins to exhibit sensitivity to deep mantle structures on comparable levels to those in the shallow mantle. By iteration 50, the inversion starts resolving low-viscosity regions near the core-mantle boundary, while all of the low-viscosity features are faithfully represented by iteration 100 and only some minor blurring of the high-viscosity northern hemisphere features remains (bottom-centre panel in Fig. 10).

As expected, the magnitude of the adjoint sensitivity decreases with each iteration of the inversion and total misfit has reduced by over four orders-of-magnitude at iteration 100 (Fig. 11). This behaviour reflects the strength of adjoint-based optimisation: the algorithm first targets regions contributing most to the misfit and then gradually refines the solution in less sensitive areas. After 100 iterations, the majority of the recovered viscosity values lie within 20 %–25 % of the target $\left(\right|{\log }_{\mathrm{10}}(\mathit{\eta }/{\mathit{\eta }}_{\mathrm{target}})|<\mathrm{0.1})$. The slightly larger errors remaining in the northern hemisphere are likely due to these being high-viscosity regions with longer Maxwell times. Here, the total amount of viscous deformation over the course of the simulation is therefore much lower and our “observations” of surface deformation are less sensitive to details of the viscosity structure in this part of the domain.

Figure 11 Reduction in objective (i.e., misfit) function value as a function of iteration number during viscosity inversion with known ice load.

3.5.4 Simultaneous inversion for viscosity and load

The previous two experiments demonstrate G-ADOPT's capability to recover either the ice load or the lateral viscosity structure under the assumption that the other is known. The natural follow-up question is whether it is possible to jointly invert for both load and viscosity structure when neither is known? To investigate this question, we use a setup that mirrors that of the previous two cases but with a control vector now containing both fields. We initialise the inversion with zero ice thickness and the axisymmetric viscosity structure. Again, we assume no observational uncertainty and have not implemented any regularisation in the objective function.

Figure 12 Adjoint-based joint optimisation for both initial ice load and viscosity structure. (a) Objective (i.e., misfit) function value as a function of iteration number. (b) Normalised ice thickness. Black dashed line = true target load; pink = inversion after 2 iterations; green = 10 iterations; blue = 100 iterations.

Despite this inversion being ostensibly more challenging than the previous two experiments, the results are highly encouraging. Over the course of 100 iterations, the objective function decreases by nearly six orders of magnitude (Fig. 12a). Already by iteration 10, optimised ice thicknesses show strong similarities with the true target distribution (Fig. 12b) and key features of the viscosity structure are visible including the low-viscosity region in the southern hemisphere and emerging high-viscosity features in the north (Fig. 13c). By iteration 100, the ice load is visually indistinguishable from the target and the viscosity field is as similar to the true structure as in the viscosity-only inversion from Sect. 3.5.3.

Figure 13 Adjoint-based joint optimisation for both initial ice load and viscosity structure. (a) True target viscosity field and surface ice load. (b) Inversion results after 2 iterations. (c) Same after 10 iterations. (d) Same after 100 iterations.

A critical factor in the success of this joint inversion is finding an appropriate balance between the contribution to the gradient from the two control fields (i.e., ice and viscosity). We have used the L₂ derivative for optimisation in all experiments, as opposed to the more traditional Euclidean l₂ derivative that is defined with respect to the control vector's discrete degrees of freedom. The L₂ formulation is more appropriate because it properly accounts for variations in mesh cell size and, more importantly, for the fact that the ice is defined on a surface mesh that has different geometric dimensions to the volumetric viscosity field. Specifically, the optimisation is driven by the gradients ${\left[\frac{\partial J}{\partial h_{\text{load}}}\right]}_{L_{\mathrm{2}}(\partial {\mathrm{\Omega }}_{\text{top}})}$ and ${\left[\frac{\partial J}{\partial \mathit{\eta }}\right]}_{L_{\mathrm{2}}\left(\mathrm{\Omega }\right)}$ such that sensitivity of the objective function to perturbations δh and δη in thickness and viscosity is given by

$$\begin{array}{}\text{(42)}& \begin{array}{rl}\mathit{\delta }J& =\underset{\partial {\mathrm{\Omega }}_{\text{top}}}{\int }{\left[{\displaystyle \frac{\partial J}{\partial h_{\text{load}}}}\right]}_{L_{\mathrm{2}}(\partial {\mathrm{\Omega }}_{\text{top}})}\mathit{\delta }{h}_{\text{load}}\mathrm{d}s\\ & +\underset{\mathrm{\Omega }}{\int }{\left[{\displaystyle \frac{\partial J}{\partial \mathit{\eta }}}\right]}_{L_{\mathrm{2}}\left(\mathrm{\Omega }\right)}\mathit{\delta }\mathit{\eta }\mathrm{d}x.\end{array}\end{array}$$

We found that when instead using the Euclidean l₂ approach, which uses a simple summation of partial derivatives with respect to the discrete degrees of freedom, the inversion became dominated by sensitivity of misfit to ice thickness, leading to suboptimal results. In such cases, the inversion typically stalled, yielding a noisy ice-thickness field and failing to recover the deep mantle or northern hemisphere viscosity variations (results not shown). Importantly, the integral L₂-approach is free of these issues and is delivered automatically in G-ADOPT through combined use of Firedrake, pyadjoint and ROL.

In summary, using a suite of idealised experiments, we have demonstrated that G-ADOPT has the ability to automatically and successfully perform adjoint-based inversions for both initial ice load and/or lateral variations in mantle viscosity. This capability is powerful since it circumvents the need for lengthy derivations of the appropriate adjoint formulations for each viscoelastic loading problem.

4 Discussion, limitations and future work

A wide range of problems in geodynamics require modelling the viscoelastic response of the solid Earth to evolving surface loads. Our long-term objective is to develop a generalised framework in G-ADOPT capable of tackling this broad class of problems. Nevertheless, our more immediate motivation is to develop more advanced GIA models that incorporate lateral Earth structure, potentially complex rheologies, and use data assimilation techniques to constrain uncertain parameters such as ice histories and Earth structure. In this context, several aspects of the current implementation – its capabilities, limitations and future directions – warrant discussion.

We have demonstrated the potential of G-ADOPT to represent Earth-like viscoelastic behaviour under surface loading. The framework already supports spherical geometries, elastic compressibility, large lateral variations in viscosity, and complex rheologies, including power-law and transience via the internal variable formulation (Al-Attar and Tromp, 2014; Crawford et al., 2017). Its development is timely, as recent studies have argued that transient deformation is required to reproduce GIA observations (Simon et al., 2022; Paxman et al., 2023; Lau, 2024). In G-ADOPT, more sophisticated transient models, such as extended Burgers rheologies, can be implemented in a straightforward manner by introducing additional internal variables. This capability therefore provides a unified setting in which competing rheological scenarios can be systematically explored.

The most important physical process still absent from G-ADOPT is self-gravity, which requires solving Poisson's equation for the gravitational potential subject to boundary conditions at infinity. Implementing these conditions in discretisation-based methods is challenging. Approaches explored in previous studies include: (i) extending the computational domain with an outer mesh to approximate infinity (e.g. May and Knepley, 2011); (ii) applying a Dirichlet-to-Neumann boundary condition, whereby the analytic spectral solution outside of the domain is imposed as a weak boundary term (e.g., Chaljub and Valette, 2004); and (iii) employing semi-infinite elements with basis functions specifically chosen to reproduce the appropriate asymptotic decay (e.g., Gharti and Tromp, 2017). An alternative strategy is to couple the numerical model to an external analytic solution for gravity (Zhong et al., 2022). Our future work will focus on a hybrid of the first two approaches.

Once self-gravity is included, extending the framework to capture rotational feedbacks of load redistribution and deformation (i.e. true polar wander; Mitrovica et al., 2005) will be relatively straightforward. Together these additions will complete the physics of the sea-level equation, which links ice-sheet evolution to global sea-level change (Farrell and Clark, 1976; Kendall et al., 2005). A practical component of this equation is the continent-ocean function, which defines the ocean load. Since shoreline migration must be represented dynamically, we will adopt approaches compatible with automatic differentiation – for example, locally smoothed (e.g., sigmoidal) approximations or a level-set method that is already implemented in G-ADOPT.

Another key strength of G-ADOPT is its suitability for extension to large-scale GIA problems requiring high spatiotemporal resolution and/or long durations. Unlike explicit time-stepping approaches (e.g. Lloyd et al., 2024), implicit methods are not limited by the Maxwell relaxation time, allowing timesteps long enough to capture full glacial cycles without incurring excessive computational cost (even when viscosity contrasts span several orders of magnitude). Moreover, G-ADOPT maintains robust accuracy across a wide range of mesh resolutions and timestep sizes, underscoring its reliability in both regional and global applications using high-performance computers. These properties naturally position the framework to scale to continental- and global-scale adjoint inversions, where computational efficiency and flexibility in rheological assumptions will be essential. Scaling tests on the Gadi supercomputer show excellent parallel performance to at least ∼ 6000 cores. Preliminary testing beyond 10 000 cores also gives reasonable results, but shows a comparative increase in the time taken to assemble the preconditioner. Strategies such as replacing the algebraic with a geometric multigrid preconditioner, which is feasible using Firedrake's interface with PETSc, may prove fruitful for maintaining computational efficiency in very high spatial resolution simulations (Weerdesteijn et al., 2023).

While adjoint methods are well established in mantle convection modelling (e.g., Bunge et al., 2003; Li et al., 2017; Ghelichkhan et al., 2021), their use in GIA is more recent. Initial studies have already demonstrated their diagnostic power: adjoint sensitivity kernels clearly expose the limitations of 1D Earth assumptions and the strong influence of lateral viscosity variations on GIA observables (e.g., Al-Attar and Tromp, 2014; Crawford et al., 2018; Kim et al., 2022; Lloyd et al., 2024). G-ADOPT provides a fully automated, scalable, and extensible inversion framework that eliminates the need for manual adjoint derivations or simplifying assumptions to manage non-linearity and coupling. This functionality is achieved by combining: (i) Firedrake, which separates model formulation from implementation within a composable finite element environment (Rathgeber et al., 2016); (ii) Pyadjoint, which derives discrete adjoints directly from variational problems expressed in unified form language (Farrell et al., 2013; Mitusch et al., 2019); and (iii) ROL, which supplies high-performance optimisation algorithms with checkpointing and optional parameter constraints. Together, these components support a streamlined workflow for forward simulation, adjoint construction, and gradient-based inversion.

It is important to acknowledge that the inverse experiments presented here commit the classic “inverse crime”: the synthetic observations used for optimisation were generated by the same numerical model used for inversion and were assumed to be noise free and complete in space and time. This deliberate choice allowed us to isolate the framework's performance under idealised conditions and confirm its correctness. The next step is to apply G-ADOPT with real observational constraints and will require formal uncertainty quantification, careful treatment of data–model misfit, and potentially hybrid assimilation strategies that combine multiple observation types. As misfit functions become more complex (e.g., by incorporating discrete relative sea-level indicators or palaeo-ice extent constraints), regularisation strategies such as smoothing, damping and prior-based parameter constraints will be essential to ensure stable adjoint solutions and interpretable model updates. These cost functions can be formally weighted using the observational data uncertainties. In addition, the application of second-order adjoint techniques provides a potential route to mapping observational uncertainty into uncertainties in inferred model parameters (Yu et al., 2025).

In summary, G-ADOPT delivers a flexible, extensible, and computationally scalable framework for adjoint-enabled modelling of viscoelastic surface-loading problems, with a particular focus on GIA. With the planned inclusion of self-gravity, rotational feedbacks, and dynamic shoreline migration, it will provide a next-generation, data-assimilating system that can directly link palaeo and modern sea-level observations to the coupled evolution of ice sheets and Earth's interior.

5 Conclusions

In this study, we have extended G-ADOPT to model viscoelastic deformation of the solid Earth under evolving surface loads. We have demonstrated its ability to run simulations in 2D and 3D domains, in Cartesian and spherical geometries, and incorporating effects such as elastic compressibility, lateral viscosity variations, and non-Maxwell rheologies (including power-law and transient behaviour). Where possible, results were benchmarked against analytical solutions and published community models. We further showed that G-ADOPT scales efficiently on high-performance computing systems, enabling simulations with high spatial resolution and long time spans. A key outstanding challenge is the incorporation of self-gravity, which will be the focus of ongoing development.

Our work also lays the foundation for addressing several long-standing challenges in GIA research. As observational constraints become more diverse and models of mantle rheology more complex, there is an increasing need for inversion frameworks that can flexibly integrate these data while maintaining physical and numerical rigour. A major strength of G-ADOPT is its automatic derivation of adjoint sensitivity kernels, which enables inverse optimisation without the need for manual derivation of adjoint equations. Using synthetic twin experiments, we demonstrated simultaneous recovery of both ice-loading history and mantle structure, showing that G-ADOPT’s adjoint implementation is robust, accurate, and suitable for application to real geophysical inverse problems. At the same time, our results highlight the inherent biases that arise when inverting for ice histories with inaccurate viscosity structures.

G-ADOPT has been developed in alignment with FAIR (Findable, Accessible, Interoperable, Reusable) software principles (Wilkinson et al., 2016). We aim to lower barriers to adoption, facilitate community benchmarking, and provide a robust foundation for collaborative development. Looking forward, the composable design of G-ADOPT opens the way to incorporating new classes of observational constraints into solid-Earth deformation problems. The integration of Firedrake, Pyadjoint, and ROL enables these extensions with minimal disruption to existing workflows.

Several physical processes central to sea-level modelling remain to be incorporated including self-gravitation, rotational feedbacks, shoreline migration, and ocean loading. While these aspects are non-trivial to implement, the framework's modular design ensures that they can be added incrementally. Embedding G-ADOPT within the Firedrake ecosystem also enables coupling with other Earth system components, such as the Icepack ice-sheet model (Shapero et al., 2021), raising the possibility of coupled Earth–climate simulations capable of data assimilation over full glacial cycles.