the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
FastIsostasy v1.0 – a regional, accelerated 2D glacial isostatic adjustment (GIA) model accounting for the lateral variability of the solid Earth
Marisa Montoya
Konstantin Latychev
Alexander Robinson
Jorge AlvarezSolas
Jerry Mitrovica
The vast majority of icesheet modelling studies rely on simplified representations of the glacial isostatic adjustment (GIA), which, among other limitations, do not account for lateral variations in the lithospheric thickness and uppermantle viscosity. In studies of the last glacial cycle using 3D GIA models, this has however been shown to have major impacts on the dynamics of marinebased sectors of Antarctica, which are likely to be the greatest contributors to sealevel rise in the coming centuries. This gap in comprehensiveness is explained by the fact that 3D GIA models are computationally expensive, rarely opensource and require a complex coupling scheme. To close this gap between “best” and “tractable” GIA models, we propose FastIsostasy here, a regional GIA model capturing lateral variations in the lithospheric thickness and mantle viscosity. By means of fast Fourier transforms and a hybrid collocation scheme to solve its underlying partial differential equation, FastIsostasy can simulate 100 000 years of highresolution bedrock displacement in only minutes of singleCPU computation, including the changes in seasurface height due to mass redistribution. Despite its 2D grid, FastIsostasy parameterises the depthdependent viscosity and therefore represents the depth dimension to a certain extent. FastIsostasy is benchmarked here against analytical, as well as 1D and 3D numerical solutions, and shows good agreement with them. For a simulation of the last glacial cycle, its mean and maximal error over time and space respectively yield less than 5 % and 16 % compared to a 3D GIA model over the regional solution domain. FastIsostasy is opensource, is documented with many examples and provides a straightforward interface for coupling to an icesheet model. The model is benchmarked here based on its implementation in Julia, while a Fortran version is also provided to allow for compatibility with most existing icesheet models. The Julia version provides additional features, including a vast library of adaptive timestepping methods and GPU support.
 Article
(8272 KB)  Fulltext XML
 BibTeX
 EndNote
1.1 GIA is an important feedback on icesheet dynamics
Glacial isostatic adjustment (GIA) denotes the crustal displacement that results from changes in the ice, liquidwater and sediment columns, as well as the associated changes in Earth's gravity and rotation axis (Mitrovica et al., 2001), ultimately altering the sea level (Farrell and Clark, 1976). In the present work, we focus on the deformation and gravitational effects. For ice sheets, the former is a net negative feedback on mass balance through the lapse rate of the troposphere, and both imply additional negative feedbacks on the dynamics of marinebased regions, where enhanced melting leads not only to a groundingline retreat but also to a regional bedrock uplift and a decrease in the seasurface height (SSH), due to the reduced load applied upon the solid Earth and the lesser gravitational pull of the ice sheet on the oceans (Gomez et al., 2010, 2012, 2015, 2018; Whitehouse et al., 2019). As depicted in Fig. 1, these effects combine in a decrease in relative sea level (RSL), which is defined as the difference between the SSH and the bedrock elevation. Compared to the retreated state (Fig. 1b), this decrease ultimately leads to a readvance of the grounding line (Fig. 1c and d), therefore conditioning the marine icesheet instability along with the buttressing effect from ice shelves (Gudmundsson et al., 2012). It was shown that the representation of the deformation and gravitational feedbacks can stabilise grounding lines on retrograde slopes (Gomez et al., 2010, 2012) and that a rapid bedrock uplift can prevent the collapse of marineterminating glaciers that are transiently forced (Konrad et al., 2015, 2016). Furthermore, an uplifting bedrock might lead isolated bathymetric peaks to connect with the ice sheet, creating socalled pinning points (Adhikari et al., 2014) that further contribute to the stability of a marine ice sheet. Although the negative feedbacks are illustrated here for icesheet retreat, they conversely apply to icesheet growth.
In addition to these effects, recent work has shown that the icesheet evolution might be significantly affected by the forebulge dynamic on longer timescales, for which viscous effects become important. This process denotes the region of comparatively small bedrock uplift (subsidence) surrounding a region of pronounced bedrock subsidence (uplift), which results from a positive (negative) surface load anomaly. Albrecht et al. (2023) suggest that the forebulge formation represents a positive feedback on icesheet growth through a decrease in the RSL close to the icesheet margin, and Kreuzer et al. (2023) show that a subsiding forebulge can increase subshelf melting through the formation of oceanic gateways that ease the intrusion of warm circumpolar deep water.
1.2 Laterally variable structure of the solid Earth modulates GIA
For a given load applied to the solid Earth, the timescale of bedrock deformation is determined by the horizontal extent of the load and the mantle viscosity. The amplitude and the pattern of deformation are in turn determined by the magnitude of the load, the mantle density and the (elastic) lithospheric thickness. These properties are close to being laterally homogeneous in many regions of the solid Earth, which motivated the development of 1D GIA models, where properties are assumed to depend only on the depth coordinate (Dziewonski and Anderson, 1981). However, some regions are an exception to this and present significant lateral variability of solidEarth properties (further simply referred to as LV), even on relatively short spatial scales. Since Antarctica displays a strong dichotomy between a moderately rifting system in the west and an old craton in the east (Behrendt, 1999), it represents a prototypical example of LV. As depicted in Fig. 2, the lithospheric thickness and uppermantle viscosity are respectively as little as 50 km and 10^{18} Pa s in the west and as large as 250 km and 10^{23} Pa s in the east (Barletta et al., 2018; Heeszel et al., 2016; Ivins et al., 2022; Lloyd et al., 2015, 2020; Morelli and Danesi, 2004; Nield et al., 2014; Whitehouse et al., 2019; Wiens et al., 2022).
For simulations of the Antarctic Ice Sheet (AIS) on the timescale of glacial cycles, accounting for LV by using 3D GIA models has shown great differences compared to 1D GIA models (Albrecht et al., 2023; Gomez et al., 2018; Van Calcar et al., 2023), leading to discrepancies reaching up to 700 km for the groundingline position, more than 1 km for the ice thickness and several metres for the sealevel equivalent volume of the AIS. Although these impacts are large, they are to be expected, given that the AIS is characterised by large marinebased regions. The East Antarctic basins and the West Antarctic Ice Sheet (WAIS) respectively represent sealevel contributions from ice grounded below sea level of about 19.2 and 3.4 m at the present day (Fretwell et al., 2013), and their evolution strongly depends on the representation of the GIA feedbacks depicted in Fig. 1. While both regions are likely to present abrupt transitions to icefree conditions under warming scenarios, the WAIS is thought to have particularly low resilience, displaying a bifurcation at a mean global warming of as little as 2 °C with respect to the preindustrial era (Garbe et al., 2020). In the context of anthropogenic warming, this is very likely to result in an unprecedented rate of sealevel rise, challenging the adaptation of coastal livelihoods that represent a large number of human societies (Kulp and Strauss, 2019).
For these reasons, comprehensive projections of sealevel rise require the representation of the Antarctic LV in coupled icesheet–GIA settings. Furthermore, the mantle viscosity is uncertain and involves discrepancies of up to 2 orders of magnitude at some locations in the upper mantle (i.e. above ∼670 km depth), depending on how viscosities are inferred from sparse seismic data (Ivins et al., 2022). Such uncertainties thus need to be propagated to the solution – typically by an ensemble of simulations. This was addressed by Coulon et al. (2021) by relying on a computationally efficient 2D GIA model that allows for laterally variable relaxation times of the deformational response. However, this is not a standard practice since only a few icesheet models are coupled to a laterally variable solid Earth, typically represented by using a 3D GIA model (Albrecht et al., 2023; Gomez et al., 2018; Van Calcar et al., 2023). This is largely due to the fact that 3D GIA models are computationally too expensive for largeensemble simulations and represent a level of complexity that may be much higher than what is required to answer most of the ongoing research questions related to the evolution of ice sheets. With the exception of CitcomSVE (Zhong et al., 2022), 3D GIA models do not offer opensource implementations at this time,^{1} and some of them even require a commercial licence (e.g. Huang et al., 2023). These obstacles have prevented most icesheet modelling studies from using 3D GIA models.
1.3 FastIsostasy: reducing the misrepresentation of LV at low computational cost
The vast majority of icesheet simulations rely on greatly simplified GIA models without accounting for the parametric uncertainties of the solid Earth, thus potentially introducing biases in sealevel projections (Gomez et al., 2015). This also holds for the Ice Sheet Model Intercomparison Project (ISMIP) (Seroussi et al., 2020), used as the physical basis for the reports of the Intergovernmental Panel on Climate Change (IPCC). In summary, the icesheet modelling community faces the somewhat paradoxical situation of being increasingly aware of how important 3D GIA is without being able to represent it at a reasonable computational cost. The work of Coulon et al. (2021) partly addresses this with a computationally efficient 2D GIA model but is also characterised by an important limitation: the viscous response is parameterised by a field of relaxation times. However, the response timescale of the solid Earth depends not only on the viscosity but also on the wavelength of the load as mentioned above. For these reasons, deriving spatially coherent maps of the relaxation time and constraining them within realistic ranges is not straightforward, as pointed out by Coulon et al. (2021).
To tackle these issues, we propose FastIsostasy here, a regional 2D GIA model inspired by first principles and specially tailored for the needs of icesheet modellers. FastIsostasy (1) accounts for LV; (2) parameterises the depth dependence of the mantle viscosity; (3) captures the dependence of the response timescale on the mantle viscosity and the spatial scale of the load; (4) approximates the regional gravitation response and sealevel evolution; (5) is computationally inexpensive; (6) is extensively tested; and (7) offers a simple, opensource and extensively documented interface for a simplified coupling to an icesheet model. To illustrate its capabilities, Antarctica is used as a leitmotif of the present work since it displays (1) a high LV and depth dependence of solidEarth properties as depicted in Fig. 2; (2) a high sensitivity to GIA due to the vast marine sectors of the AIS; (3) a large impact on the future of human societies due to possible rapid sealevel rise; and (4) large uncertainties in the solidEarth parameters due, in part, to limited regional data sets. Antarctica might therefore be “the toughest test” when it comes to GIA modelling. Accordingly, the tools provided here are equally applicable to any other region covered by a past, present or future ice sheet.
1.4 FastIsostasy in the model hierarchy
The hierarchy of GIA models displays an important complexity gap between the computationally cheap models, which are largely used by the icesheet modelling community, and the computationally expensive models, developed by the GIA community. To give an impression of this, we herein present a brief overview, summarised in Tables 1 and 2, of the GIA model classes that are available to date and focus on the computation of the deformational response, with references to opensource implementations that we are aware of. The governing equations of the three first models can be found in Appendix A.
1.4.1 ELRA
The elastic lithosphere–relaxed asthenosphere (ELRA; Le Meur and Huybrechts, 1996) model conceptualises the structure of the solid Earth as two layers stacked along the depth dimension of a Cartesian coordinate system, obtained by a regional projection of the spherical Earth. The elastic lithosphere is parameterised by its constant thickness $T(x,y)=T$ and undergoes instantaneous compression under the effect of a load. It is underlain by the asthenosphere,^{2} idealised as a viscous halfspace parameterised by a constant relaxation time $\mathit{\tau}(x,y)=\mathit{\tau}$. According to this, the vertical displacement exponentially converges to the equilibrium solution, which is computed by convolving the load with a Green’s function. Parameterising the transient behaviour with a relaxation time is however simplistic, since, in reality, the response timescale of the solid Earth depends not only on the viscosity but also on the wavelength of the load, as mentioned previously. Furthermore, ELRA does not represent the depth dependence of the mantle viscosity or any LV. Due to its 2D regional domain, it ignores the gravitational and rotational feedbacks on the sea level, as well as the prestress and selfgravitation (Purcell, 1998) – the latter being partly cancelled by the lack of sphericity (Amelung and Wolf, 1994). Konrad et al. (2014) demonstrated that ELRA displays important transient differences to a 1D GIA model as well as discrepancies in the representation of the peripheral forebulge. Despite these numerous flaws, ELRA – or even simplified versions of it – remains a widespread choice among icesheet modellers (DeConto and Pollard, 2016; Lipscomb et al., 2019; Pattyn, 2017; Quiquet et al., 2018; Robinson et al., 2020; Rückamp et al., 2019), as it mimics the viscoelastic behaviour of the solid Earth with little implementation effort and at low computational cost. Its simplicity has led to a large number of implementations within opensource icesheet models (e.g. Lipscomb et al., 2019; Robinson et al., 2020), but, to our knowledge, no modular implementation is available to date. Icesheet modellers have therefore repeatedly spent time implementing ELRA, possibly with suboptimal computational performance, as it does not represent the primary focus of their work.
1.4.2 ELVA
This modelling approach was proposed by Cathles (1975), applied to icesheet modelling for the first time by Lingle and Clark (1985), and efficiently implemented by Bueler et al. (2007) through a Fourier collocation method (FCM). Although this model is sometimes named after the authors of the aforementioned work, we try to provide a unifying terminology here and therefore call it elastic lithosphere–viscous mantle (ELVA). ELVA resembles ELRA in its structure but is parameterised by the spatially homogeneous uppermantle viscosity $\mathit{\eta}(x,y)=\mathit{\eta}$. It thus avoids any conversion from viscosity to relaxation time and allows for the mechanical response to depend on the wavelength of the load (Bueler et al., 2007). Furthermore, it permits embedding more of the radial structure of the mantle viscosity by introducing a viscous channel between the elastic plate and the viscous halfspace.^{3} However, it does not address any other limitation of ELRA. It is worth mentioning that PISM (Winkelmann et al., 2011) provides an opensource implementation of ELVA, which is however embedded within a larger code base. This lack of modularity is addressed by giapy (Kachuck, 2017), a Python implementation of ELVA that might be more accessible than code traditionally written in Fortran or C$++$. ELVA was used, for example, by Kachuck et al. (2020) and Book et al. (2022) to study the stabilising potential of rapid bedrock uplift on the WAIS groundingline retreat.
1.4.3 LVELRA
The laterally variable ELRA (LVELRA) proposed by Coulon et al. (2021) is a generalisation of ELRA to include laterally variable uppermantle relaxation time τ(x,y) and lithospheric thickness T(x,y). The equilibrium displacement is obtained by solving equations derived from thinplate theory (Ventsel and Krauthammer, 2001), which requires the use of finitedifference methods (FDMs) and is computationally more expensive than ELRA, since a large system of linear equations needs to be solved. To obtain τ(x,y), Coulon et al. (2021) apply a Gaussian smoothing on a binary field, with $\mathit{\tau}(x,y)={\mathit{\tau}}_{\mathrm{1}}$ in East Antarctica and $\mathit{\tau}(x,y)={\mathit{\tau}}_{\mathrm{2}}$ in the rest of the domain. Since LVELRA does not include a lateral coupling between the transient behaviour of neighbouring cells, this smoothing ensures a certain spatial coherence when relaxing the displacement field to the equilibrium solution; i.e. neighbouring cells have similar timescales. In addition to the limitations mentioned in Sect. 1.3, this coupling prevents the representation of very localised features, depicted in Fig. 2 and inferred in many studies (e.g. Barletta et al., 2018; Heeszel et al., 2016; Nield et al., 2014). In the rest of the paper, we will refer to these limitations as the ones resulting from a relaxed rheology.
Although computing the changes in SSH resulting from changes of Earth's gravity field requires, a priori, a global domain, it can also be approximated on a regional one. This was done by Coulon et al. (2021) and combined with LVELRA, resulting in the socalled elementary GIA model. This approach represents one of the most comprehensive regional GIA models developed to date and a valuable improvement for regional modelling, as it bypasses the computational expense of more complex models. The opensource icesheet model Kori includes the effort of Coulon et al. (2021), however not in a modular way that is directly usable to other icesheet modellers.
Konrad et al. (2014)Le Meur and Huybrechts (1996)Bueler et al. (2007)Kachuck et al. (2020)Book et al. (2022)Coulon et al. (2021)1.4.4 Global 1D GIA models
Global 1D GIA models capture the radial structure of the solid Earth (down to the core–mantle boundary) but none of its lateral variability. They compute the gravitational field as well as the vertical and horizontal deformation by solving the underlying partial differential equations (PDEs) after spherical harmonic expansion of the dependent variables. The 1D GIA models typically represent the spatial heterogeneity of the SSH, the migration of shorelines and the rotational feedback (Kendall et al., 2005; Mitrovica and Milne, 2003; Spada and Melini, 2019). Most of them were crossvalidated by Spada et al. (2011) and Martinec et al. (2018), showing great agreement while presenting intermediate computational cost. However, they are incapable of rendering any LV (e.g. Klemann et al., 2008). Spada and Melini (2019) proposed an opensource implementation of 1D GIA, which remains an exception in the field.
1.4.5 Regional 3D GIA models
Based on finiteelement methods (FEMs), Nield et al. (2018) and Weerdesteijn et al. (2023) have proposed regional models to compute the solidEarth deformation in the presence of LV. Unlike ELRA, ELVA and LVELRA, these regional GIA models resolve the depth dimension (down to the core–mantle boundary) and allow for grid refinement in regions where a higher resolution is needed. In particular, the work of Weerdesteijn et al. (2023) provides an opensource implementation that is compatible with heavily parallel hardware by extending ASPECT (Advanced Solver for Problems in Earth's ConvecTion), a model originally developed to solve mantle convection problems. Despite this, ASPECT requires about an hour to compute a few hundred years of highresolution bedrock deformation on 256 CPUs. This represents a computational cost that is too high for most ongoing icesheet modelling studies, while ignoring the gravitational and rotational effects of GIA.
1.4.6 Global 3D GIA models
The 3D GIA models account for all the processes represented in 1D GIA models and are, in addition, capable of fully capturing the heterogeneity of solidEarth properties. This also results in simulations that are more complicated to set up, since the user needs to provide fields of lithospheric thickness and mantle viscosity. Unlike 1D GIA models, 3D ones have not been systematically benchmarked but can be considered to be the best technology available for cases like Antarctica. In these models, the computation of the deformational response either relies on spherical harmonics (e.g. Bagge et al., 2021), FEM (e.g. A et al., 2013; Huang et al., 2023; Martinec, 2000; Sasgen et al., 2018; van der Wal et al., 2015; Wu, 2004; Zhong et al., 2022), the finitevolume method (FVM; Latychev et al., 2005; Gomez et al., 2018) or perturbation theory (e.g. Wu and Wang, 2006). Simulations on glacial timescales typically require at least several hours and up to few days (Albrecht et al., 2023; Pan et al., 2022; Zhong et al., 2022) of computation, even with heavily parallelised code. This is particularly problematic for propagating parametric uncertainties of the solid Earth on long simulations, since the limit of computational resources is typically reached with only one or a few ensemble members. As mentioned above, Zhong et al. (2022) proposed the first opensource implementation of a 3D GIA code.
Spada and Melini (2019)Nield et al. (2018)Weerdesteijn et al. (2023)Latychev et al. (2005)Zhong et al. (2022)1.4.7 FastIsostasy
The summary above points out a gap between the elementary GIA model (Coulon et al., 2021) and regional 3D GIA models (Nield et al., 2018; Weerdesteijn et al., 2023): the former is computationally cheap but suffers the limitation of a relaxed rheology, whereas the latter accurately represent the viscoelastic response of the solid Earth but come with a computational cost that makes them impractical for most icesheet modelling studies. FastIsostasy fills this gap by relying on LVELVA, a laterally variable generalisation of ELVA, coupled to a regional sealevel model (ReSeLeM), both of which are introduced in Sect. 2, along with a discussion of the underlying limitations. In Sect. 3, we discuss the practical features of its Julia implementation, such as the adaptive time stepping used for integration and GPU support. In Sect. 4, we subsequently benchmark FastIsostasy against analytical, as well as 1D and 3D numerical, solutions. Finally, we discuss the results as well as possible future improvements.
Remarks on opensource code. Thanks to recent work, the GIA model classes listed above now present at least some opensource code, which has typically become available much later than the first equivalent piece of proprietary code. For instance, 1D and 3D GIA models have already existed for 40 and 20 years respectively, but their first opensource implementations were only published much later, in respectively Spada and Melini (2019) and Zhong et al. (2022). We specifically refer here to source code that is licensed and can be downloaded and used without request. However, all opensource code mentioned above lacks tools from modern software development, including (1) dynamically built documentations with code examples; (2) automated test suites; and (3) a transparent development, which can be eased, for instance, by the use of GitHub issues and pull requests. Although these aspects are of a technical nature, we believe that they can make GIA models more userfriendly and their development more transparent, participative and reliable. This, in turn, reduces barriers in future research, and we have therefore applied these concepts to FastIsostasy.
Since LVELVA is a generalisation of ELVA, the latter can be used in FastIsostasy by simply providing it with homogeneous parameter fields. Additionally, we implemented ELRA in an efficient way, further described in Appendix A. Thus, both of these simpler models now present a modular, optimised and documented implementation, distributed under the GNU General Public License v3.0.
2.1 Preliminary considerations
As depicted in Fig. 3, FastIsostasy assumes a rectangular domain Ω⊂ℝ^{2}, obtained from a projection of the spherical Earth onto a Cartesian plane with dimensions 2 W_{x} and 2 W_{y}, in the directions of the lateral coordinates x and y respectively. We introduce the uniform spatial discretisation step ${h}_{x}={h}_{y}=h$ such that the domain is subdivided into N_{x}×N_{y} cells, with ${N}_{x},{N}_{y}\in \mathbb{N}$. We define all variables that are not specified as scalars (cf. Table 3) to be smooth fields, such as, for instance, the vertical load ${\mathit{\sigma}}^{zz}(x,y,t):{\mathbb{R}}^{\mathrm{3}}\to \mathbb{R}$. For convenience, we will omit the space and time dependence from now on. The discretised equivalent of smooth fields is denoted by bold symbols, e.g. ${\mathit{\sigma}}^{zz}\in {\mathbb{R}}^{{N}_{x}\times {N}_{y}}$, with their entries denoted by the index notation ${\mathit{\sigma}}_{i,j}^{zz}$, with $i\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots},{N}_{x}\mathit{\}}$, $j\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots},{N}_{y}\mathit{\}}$. The vertical load field is expressed as
with g as the mean gravitational acceleration at Earth's surface and ρ^{ice}, ρ^{sw} and ρ^{sed} as the mean densities of ice, seawater and sediment respectively.^{4} The height anomalies ΔH^{ice}, ΔH^{sw} and ΔH^{sed} of the corresponding columns are defined with respect to a reference state. On this domain, the first and second spatial derivatives of an arbitrary field M can be computed with central differences:
with K as the distortion factor of the chosen projection.^{5} Furthermore, the pseudodifferential operator $\left\mathrm{\nabla}\right$ of an arbitrary matrix M is adapted from Bueler et al. (2007) to suit distorted grids:
with ⊙ as the elementwise product, ⊘ as the elementwise division, ℱ as the Fourier transform, ℱ^{−1} as its inverse and κ as the coefficient matrix derived in Bueler et al. (2007). Models that do not account for distortion underestimate the length and area of cells away from the reference latitude and therefore require a domain with restricted spatial extent, a limitation that is overcome here.
2.2 Lumping the depth dimension
As depicted in Fig. 3, the vertical structure of the solid Earth is modelled by a stack of layers along the vertical dimension z. With the layer index $l\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{\dots},L\mathrm{1},L\mathit{\}}$ going from top to bottom, the layers are the following:

l=0 – an elastic plate with a laterally constant Young modulus ${E}_{\mathrm{0}}(x,y)={E}_{\mathrm{0}}$ and Poisson ratio ${\mathit{\nu}}_{\mathrm{0}}(x,y)={\mathit{\nu}}_{\mathrm{0}}$ and laterally variable thickness T_{0}(x,y);

$l\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots},L\mathrm{2},L\mathrm{1}\mathit{\}}$ – an arbitrary number of viscous channels, each with laterally constant Young modulus ${E}_{l}(x,y)={E}_{l}$, Poisson ratio ${\mathit{\nu}}_{l}(x,y)={\mathit{\nu}}_{l}$ and laterally variable viscosity η_{l}(x,y) (as depicted in Fig. 3, the first of these layers has a laterally variable thickness T_{1}(x,y) that is complementary to T_{0}(x,y) and allows for all further layers to have a homogeneous ${T}_{l}(x,y)={T}_{l}$ for l≥2);

l=L – a viscous halfspace with a laterally constant Young modulus ${E}_{L}(x,y)={E}_{L}$, Poisson ratio ${\mathit{\nu}}_{L}(x,y)={\mathit{\nu}}_{L}$ and viscosity ${\mathit{\eta}}_{L}(x,y)={\mathit{\eta}}_{L}$.
Whereas l=0 represents the lithosphere, all further layers represent the remaining mantle. FastIsostasy lumps the latter layers into a single layer by computing a socalled effective viscosity for the whole mantle. The key to do so is provided by Cathles (1975), where a threelayer model including an elastic plate, a viscous channel and a viscous halfspace is converted into a twolayer model where the viscous channel and the viscous halfspace have been lumped into a single halfspace by introducing the following scaling factor:
with T as the channel thickness, $\stackrel{\mathrm{\u0303}}{\mathit{\eta}}$ as the channel viscosity divided by the halfspace viscosity, C=cosh (Tκ) and S=sinh (Tκ). The characteristic wavenumber $\mathit{\kappa}=\mathit{\pi}{\mathit{\lambda}}^{\mathrm{1}}$ is defined by choosing a characteristic wavelength λ for the load. Hence, solving the threelayer case can be formulated as solving the twolayer case with the halfspace viscosity scaled by R. We propose generalising this idea by performing an induction from the bottom to the top layers, i.e. with decreasing l.
Initialisation. Layer l=L is a viscous halfspace with ${\mathit{\eta}}_{l}^{\mathrm{eff}}={\mathit{\eta}}_{L}$.
Induction step. Layer l+1 can be represented as a viscous halfspace with ${\mathit{\eta}}_{l+\mathrm{1}}^{\mathrm{eff}}$ and is overlain by a viscous channel l. These can be converted in an equivalent halfspace with effective viscosity:
Thus, ${\mathit{\eta}}_{\mathrm{1}}^{\mathrm{eff}}$ is the effective viscosity of the halfspace representing the compound of layers $l\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots},L\mathrm{1},L\mathit{\}}$. In essence, this represents a nonlinear mean of the viscosity over an arbitrary number of layers, which is only computed at initialisation and can improve the parameterisation of the depth dimension compared, for instance, to Bueler et al. (2007). However, this approach presents important limitations.

Since the depth dimension is not resolved, the multimodal response of Earth to surface loading is not captured as accurately as in 1D and 3D GIA models. For instance, the larger the load, the deeper the deformation into the mantle and the more relevant the radial layering of viscosity and density, which is not captured in FastIsostasy.

The deformational response is likely to be less accurate for loads with dimensions that substantially differ from the characteristic wavelength λ, used in Eq. (4) to lump the depth dimension. We however emphasise that, for a given ice sheet, λ can be chosen such that the near field of deformation is represented well, which is typically what is required in simulations coupled to icesheet models. For all the computations presented here, we choose λ to be the mean of {W_{x},W_{y}}.

Equation (4), as derived in Cathles (1975), applies to laterally constant viscosities. In the present case, we however allow for the viscosity to be laterally variable and apply this scaling for each column.
Finally, to account for compressibility and lateral variations in the shear modulus, a scaling α of the viscosity is introduced and described in Appendix C. This yields the corrected effective viscosity η, which brings the Maxwell time of FastIsostasy close to that of a 3D GIA model:
Although converting the 3D problem into a 2D one introduces the limitations mentioned above, it also greatly reduces the computational cost. In addition, the partial differential equation (PDE) governing an elastic plate on a viscous halfspace can be transformed into an ordinary differential equation (ODE), as we describe in the next section.
2.3 LVELVA
We now assume that the aforementioned lumping of the layers has been performed and that the lithosphere and underlying mantle are represented by an elastic plate overlaying a viscous halfspace. Since the vertical extent of the plate is typically 2 orders of magnitude smaller than its horizontal one, it is considered to be thin. By assuming a Maxwell rheology, the total vertical displacement u^{el}+u of the bedrock resulting from stress σ^{zz} can be decomposed in an elastic and a viscous response respectively denoted by u^{el} and u. As in Bueler et al. (2007), the elastic response of the lithosphere is computed by a convolution of the load σ^{zz} with an appropriate Green's function Γ^{el}:
This represents the instantaneous compression of the lithosphere and accounts for the distortion resulting from the projection. In reality, this process takes place on the timescale of days, but it can be considered to be instantaneous compared to the long timescales of the viscous response and the icesheet dynamics. To construct the elastic Green's function, a “Gutenberg–Bullen A” spherical Earth model is assumed, and values from Farrell (1972, Table A3) are used, as done by Bueler et al. (2007). This treatment of the elastic response shows great agreement with a 3D GIA model, as shown in Sect. 4.
When material from the solid Earth is displaced, a hydrostatic force counteracting the load arises. We define the pressure field p as the sum of all these effects:
with ρ^{l} and ρ^{m} as mean densities of the lithosphere and the upper mantle. Since the displacement occurs in Earth's outermost layers, we assume g to be constant over these shallow depths. The 3D GIA models usually represent the elastic lithosphere as a viscous layer with very high viscosity, and the elastic displacement therefore also implies a hydrostatic force. We argue that this is closer to reality and adapt this point of view to the present context by including the elastic displacement in Eq. (8), unlike Bueler et al. (2007) and Coulon et al. (2021). The evolution of the viscous displacement is therefore coupled to the elastic one and is governed by
with M_{xx}, M_{yy} and M_{xy} as the flexural moments for a thin plate (Coulon et al., 2021; Ventsel and Krauthammer, 2001):
In these equations, $D={E}_{\mathrm{0}}\phantom{\rule{0.125em}{0ex}}{T}_{\mathrm{0}}^{\mathrm{3}}(x,y)\phantom{\rule{0.125em}{0ex}}\left(\mathrm{12}\right(\mathrm{1}{\mathit{\nu}}_{\mathrm{0}}^{\mathrm{2}}){)}^{\mathrm{1}}$ is the laterally variable lithospheric rigidity field. The PDE can be understood as an ad hoc generalisation of ELVA (Bueler et al., 2007; Cathles, 1975; Lingle and Clark, 1985) that is inspired by LVELRA (Coulon et al., 2021) and further described in Appendix A. Though we did not manage to formally derive it by generalising the work of Cathles (1975) to heterogeneous viscosities, Eq. (9) yields results that are very close to those of a 3D GIA model, as shown in Sect. 4. F on the righthand side of the PDE can be evaluated by finite differences, as defined in Eq. (2):
Subsequently, a Fourier collocation of this equation can be achieved by making use of Eq. (3) and rearranging terms:
with ε≪1 as a regularisation term to avoid division by 0. Thus, by using a hybrid FDM–FCM approach the PDE expressed in Eq. (9) is transformed into an ODE, which can be solved with explicit integration methods. In particular, we thus avoid solving a large system of linear equations, as normally done in FDM, FVM or FEM code. Note that in Bueler et al. (2007) the closed form of a Crank–Nicolson (implicit) scheme is derived, thus providing unconditional stability. Due to the complexity of the righthand side, finding such a closed form for LVELVA is more challenging and goes beyond the scope of this work. We emphasise that the smaller time steps resulting from explicit schemes might nevertheless be needed for (1) coupling purposes, as a dense output in time can be provided to the icesheet model, and (2) accurately capturing the fast dynamics that can occur in regions of low viscosity.
Far away from the ice sheet, the changes of the load are comparatively small, and we therefore require the displacement to be zero at the domain boundaries.^{6} However, FCM does not allow explicit treatment of Dirichlet boundary conditions (BCs). To enforce its approximate representation, we subtract the mean displacement of the corner vertices from the solution at each time step, which is here expressed with the common choice of notation from programming:
Note that this differs from Bueler et al. (2007), where the whole domain boundary is used for this purpose. We argue that our approach is a better representation of the required BC because corner points are (1) further away from the load and (2) equidistant from the centre of a rectangular domain.
2.4 Regional sealevel model (ReSeLeM)
In a coupled setting, a GIA model is typically expected to take the ice thickness field as an input and to return the RSL field S as an output, which can be expressed as
with z^{ss} as the SSH, z^{b} as the bedrock elevation, s as the barystatic sea level (BSL), N as the SSH perturbation due to changes in the gravitational field and c as a timedependent scalar. In a global model, c ensures mass conservation of water, which is however impossible to ensure in a regional model with open boundaries. In contrast, we use c to impose a mean zero SSH perturbation at the domain boundary, similar to Eq. (15). The sealevel terminology used here follows the definitions of Gregory et al. (2019), which we refer to for any further detail. Whereas the displacements u and u^{el} are computed as described in Sect. 2.3, the SSH perturbation N can be regionally approximated by the convolution of the mass anomaly with an appropriate Green's function Γ^{N}, as proposed in Coulon et al. (2021):
with R_{e} as Earth's radius at the Equator, M_{e} as Earth's mass and θ as the colatitude. To avoid dividing by θ=0°, we impose a minimal colatitude of the order of the resolution. By neglecting the changes in surface load from sediments and expressing ΔH^{sw} with the RSL, Eq. (1) becomes
with 𝒪 as the ocean function; 𝒞 as the continent function; and 𝒢 as the groundedice function, which can be expressed by introducing the indicator function 𝟙 and the ice thickness above flotation H^{af}:
In Eq. (21), the activation mask 𝒜 defines what we further consider the near field of displacement: it yields 1 close to the ice sheet and 0 otherwise. This approach is similar to what is done in Coulon et al. (2021) and is, by definition, somewhat arbitrary. The choice made here is illustrated in Sect. 4 for Antarctica. Most importantly, the activation mask enforces a zero change in load close to the boundary of the domain, which is necessary to fulfil the BCs expressed in Eq. (15). This is a limitation compared to a global GIA model but allows for accounting for water column changes in the near field of the ice sheet, unlike most regional models (Book et al., 2022; Bueler et al., 2007; Kachuck et al., 2020; Lingle and Clark, 1985; Le Meur and Huybrechts, 1996; Weerdesteijn et al., 2023).
In Goelzer et al. (2020), the evolution of the BSL is described as^{7}
with A_{pd} as the presentday ocean surface and V as the volume contribution of ice sheets to the ocean. The latter is decomposed into V^{af}, the contribution from ice above flotation; V^{pov}, the contribution from changes in the bedrock height; and V^{den}, the contribution from density differences between meltwater and seawater. We refer the reader to Goelzer et al. (2020) for the detailed computation of these quantities, which are defined with respect to a reference state, typically the presentday one. Assuming a fixed ocean surface to compute the evolution of the sea level can, however, lead to a bias of tens of metres over glacial cycles. To tackle this, we propose an extension of Goelzer et al. (2020) that accounts for the time dependence of the ocean surface A(t) when computing the BSL. To this end, we first introduce a time discretisation t=k Δt, with 𝒪(Δt)=10 years. We further introduce the ocean surface as a function A(s), which is nonlinear with respect to the BSL s when a realistic topography is used, as shown in Fig. 4 and described with more detail in Sect. B. The volume change $\mathrm{\Delta}{V}_{k}={V}_{k}{V}_{k\mathrm{1}}$ over a time step leads to a change $\mathrm{\Delta}{s}_{k}={\mathrm{s}}_{k}{\mathrm{s}}_{k\mathrm{1}}$ of the BSL. Since Δt is much smaller than glaciological timescales, Δs_{k} is a small number and the nonlinear relationship between the volume contribution, the ocean surface and the BSL can be approximated by a trapezoidal rule, pictured in Fig. 4 and described by the following equation:
Equation (30) is solved by using s_{k−1} as an initial guess, and the updated BSL s_{k} is typically obtained after a few iterations of the nonlinear solver. This is of course an important simplification compared to global GIA models, which typically resolve the migration of shorelines (Kendall et al., 2005; Mitrovica and Milne, 2003). Nonetheless, this is an improvement compared to fixing A(t)=A_{pd}. In particular, Fig. 4d shows that the sea level of the Last Glacial Maximum (LGM) is overestimated by about 5 m for fixed ocean boundaries compared to our trapezoidal approximation. This can lead to differences of several kilometres in the groundingline position, depending on the local bedrock slope. We emphasise that a more sophisticated approach than ours is likely to require a global domain, which we want to avoid here.
In summary, allowing for a timevariable ocean surface is the main adaptation of the sealevel treatment described in previous work (Coulon et al., 2021; Goelzer et al., 2020), and it constitutes ReSeLeM. FastIsostasy involves, as depicted in Fig. 5, a coupling between LVELVA and ReSeLeM.
2.5 Limitations
LVELVA presents limitations, since it relies on a linear PDE describing the macroscopic behaviour of the solid Earth as a Maxwell body. Therefore it does not account for transient rheologies (Caron et al., 2017; Ivins et al., 2021), nonlinear rheologies (Gasperini et al., 2004; Kang et al., 2021), composite rheologies (van der Wal et al., 2010, 2015), anisotropy (Accardo et al., 2014; Beghein et al., 2006) or microscale properties of the material (Van Calcar et al., 2023). LVELVA only computes the vertical displacement of GIA and neglects the horizontal one, which has a negligible impact on icesheet dynamics. Nonetheless, the horizontal displacement might be used to constrain GIA models through GNSS (global navigation satellite system) measurements, and its implementation is left for future versions of the model.
The governing equation of LVELVA was postulated here in an ad hoc way, as described in Sect. A, but lacks a formal derivation. Since the depth dimension is not resolved in LVELVA, Stokes flow of the mantle is not fully represented, similar to depthintegrated solvers of icesheet dynamics and shallowwater approximations used in general circulation models. In addition, the regional nature of the domain makes it inherently complicated to ensure BCs that are consistent with a global conservation of mass. Therefore, the sea level is only solved in an approximate way, and the feedback from perturbations in Earth's rotation, which will be small near the poles, is not accounted for. In particular, we emphasise that the update of the BSL described in Eq. (30) fails to give good results if most of the contribution stems from an ice sheet that is not included in the domain. For instance, the BSL computed in FastIsostasy over a glacial cycle cannot be correct if the domain only covers Antarctica, since most of the contribution stems from the Northern Hemisphere ice sheets. However, this can however be by setting up a domain for the Northern Hemisphere and coupling it to the Antarctic one or simply by taking the BSL obtained from proxies or global GIA models.
Future releases of FastIsostasy will focus on addressing some of these problems. We emphasise that, despite these limitations, the results presented in Sect. 4 show that FastIsostasy is a significant improvement compared to ELRA and ELVA, since it represents the nearfield GIA response of laterally variable Earth structures with improved accuracy and a negligible increase in computational cost. We believe that this adequately covers the needs of most icesheet modellers.
FastIsostasy has been implemented in Julia (FastIsostasy.jl) and in Fortran. Julia (Bezanson et al., 2017) is a highperformance language with a vast ecosystem, on which FastIsostasy.jl relies to offer convenient features and efficient computation.

To evaluate the righthand side of the ODE obtained in Eq. (14) and perform the convolutions used to compute the elastic and the gravitational response, FastIsostasy.jl uses forward and inverse fast Fourier transforms (FFTs), which are implemented in an optimised way in FFTW.jl (Frigo and Johnson, 2005). Evaluating the righthand side therefore scales with a computational complexity of 𝒪(N log _{2}N), for a matrix of size $N={N}_{x}\times {N}_{y}$. To achieve an even better speed increase, (1) N_{x} and N_{y} are generally chosen as powers of 2, (2) FFTs are precomputed as far as possible, and (3) the transforms are computed in place to reduce the memory allocation. FastIsostasy owes its name to fastearth, an early implementation of ELVA (referred to in the Acknowledgements) and to its reliance on FFTs to perform all the expensive computations.

To subsequently integrate the righthand side in time, FastIsostasy.jl uses OrdinaryDiffEq.jl (Rackauckas and Nie, 2017), a package that offers a wide range of optimised routines. We restrict ourselves here to explicit methods, which range from the simplest explicit Euler scheme up to schemes of order 14. For all the results presented here, we used the Runge–Kutta method proposed in Bogacki and Shampine (1996). Explicit integration schemes typically require decreasing the time step with increasing spatial resolution, which is handled by the adaptive timestepping methods of OrdinaryDiffEq.jl to prevent instabilities. This requires more evaluations of the righthand side and increases the computational complexity of the full problem to higher than 𝒪(N log _{2}N). By providing keyword arguments, the user is able to influence any option related to the integration in time, such as the scheme, the error tolerance and the minimal time step.

FastIsostasy.jl uses CUDA.jl (Besard et al., 2019) and ParallelStencil.jl to optionally run performancerelevant computations on a GPU (so far restricted to NVIDIA hardware). Due to their heavily parallelised architecture, GPUs are able to scale better than CPUs for some computations. The speed increase thus obtained will be illustrated in Test 1 of the model validation. Offering a GPUparallelised GIA code is unprecedented to our knowledge and only requires the user to set the keyword argument
use_cuda=true
. 
In FastIsostasy.jl, the nonlinearity introduced by the timedependent ocean surface is solved by using NLsolve.jl, in combination with an interpolation of A(s), which is constructed at initialisation using Interpolations.jl. Since A(s) is monotonic and initial guesses are close to the solution, the computation time associated with this step is negligible. Furthermore, whereas the adaptive time stepping is convenient for enforcing the stability of the viscous displacement, updating the diagnostics – such as the elastic displacement, the ocean surface and sea level – can be done less frequently. For instance Δt=10 years is used in the present work and can be determined by the user through a keyword argument.
As illustrated above, FastIsostasy.jl relies on numerous Julia packages. Since it is a registered package, it can however be easily installed, along with all its dependencies, by simply running add FastIsostasy
in Julia's package manager. Furthermore, it is thoroughly documented at https://janjereczek.github.io/FastIsostasy.jl/dev/ (last access: 4 July 2024), including a description of the application programming interface; a tutorial; and practical examples, which are a simplified version of the code used for the results shown in Sect. 4. Additionally, FastIsostasy.jl is designed in a modular way that facilitates its coupling to an icesheet model, and we therefore believe that the implementation burden associated with its use is very low.
Since Julia does not yet support compilation to binaries, FastIsostasy is additionally programmed in Fortran to allow for compatibility with most existing icesheet models and has already been coupled to Yelmo (Robinson et al., 2020). Since Fortran does not provide packages that are as convenient as those of the Julia ecosystem, the Fortran version (1) does not allow for computation on a GPU, (2) only provides explicit Euler schemes for integration in time and (3) does not allow for timeevolving ocean boundaries.
We now validate FastIsostasy with series of tests.

Test 1 is a comparison to an analytical solution for an idealised load on a homogeneous, flat Earth. This aims to check that the numerics are implemented well for the simplest case and that our results are comparable to Bueler et al. (2007).

Test 2 is a comparison to benchmark solutions of three different 1D GIA models, presented in Spada et al. (2011). This aims to understand the discrepancies which can arise from the lumping of depthdependent viscosity profiles and the regional approximations of the SSH perturbation.

Test 3 is a comparison to the 3D GIA model Seakon (Gomez et al., 2018; Latychev et al., 2005) for idealised cases of LV. This aims to check whether Eq. (9) and its discretisation, Eq. (14), are valid approximations of the deformational response in the presence of LV. Here we will also compare the elastic displacement of Seakon and FastIsostasy.

Test 4 is a comparison to Seakon with realistic LV and forced by the ice loading of a full glacial cycle. This aims to check whether loads and Earth structures of typical applications can be represented reasonably well.
These tests are summarised in Table 4 and aim to quantify, as independently as possible, each source of error between FastIsostasy and the baseline solutions listed above. This is measured by an absolute and a relative value respectively defined as
with the indices “FI” and “BL” respectively indicating the FastIsostasy and baseline solutions. We refer to the mean and maximal errors over space as $\overline{e}\left(t\right)$ and $\widehat{e}\left(t\right)$. In the forthcoming analysis, we will often mention a tight upper bound for these time series to quantify them in a scalar way and will emphasise the maximal error, since the spatial mean can hide important local discrepancies.
4.1 Test 1 – analytical solution for an idealised load on a homogeneous Earth
We first reproduce the test proposed in Bueler et al. (2007) by using a twolayer model with ${W}_{x}={W}_{y}=\mathrm{3000}$ km, $N={N}_{x}={N}_{y}=\mathrm{256}$ and h≃23 km. The first layer is parameterised by the lithospheric thickness $T(x,y)=\mathrm{88}$ km, and the underlying halfspace is parameterised by the mantle viscosity $\mathit{\eta}(x,y)={\mathrm{10}}^{\mathrm{21}}$ Pa s. The load is a Heaviside function in time that represents a flat cylinder of ice, with radius R=1000 km and thickness H=1 km, placed at the centre of the computation domain. For this idealised case, an analytical solution of the viscous solution is provided in Bueler et al. (2007), yielding
with J_{0} and J_{1} as the Bessel functions of the first kind and are respectively of order 0 and 1 and $\mathit{\beta}=\mathit{\beta}\left(\mathit{\kappa}\right)={\mathit{\rho}}^{\mathrm{m}}g+D{\mathit{\kappa}}^{\mathrm{4}}$. To make the solution of FastIsostasy comparable to this, we set K_{i,j} to 1 and ρ^{l} to 0, which neglects distortion and prevents the elastic displacement from contributing to the pressure term. Figure 6a shows crosssections of the domain along the x dimension, demonstrating that the numerical solution closely follows the analytical one. Complementarily, Fig. 6b shows the corresponding maximal and mean error over time. For t≥5000 years, the viscous displacement is captured with ${\overline{e}}^{\mathrm{abs}}<{\widehat{e}}^{\mathrm{abs}}<\mathrm{1}$ m. For t≤2000 years, the displacement surface is well captured in terms of shape but appears to be slightly shifted along the z dimension due to the approximate treatment of the BCs as written in Eq. (15), leading to a larger upper bound on the error ${\widehat{e}}^{\mathrm{abs}}\left(t\right)<\mathrm{5.8}$ m, which corresponds to $\widehat{e}<\mathrm{0.021}$. Whereas in Bueler et al. (2007) a correction of this effect is applied based on the knowledge of the analytical solution, we decide not to do so here, first, because such a correction only applies to this specific case and, second, because users should be informed about the potential numerical error that will arise in their experiments. When imposing Eq. (15), the unrealistically high forcing rate resulting from the Heaviside load leads to errors that are higher than what is expected from a load that is coherent in time. The discretisation error presented here should be therefore understood as an upper bound for a flat Earth.
Figure 6c shows that the maximal and mean equilibrium error respectively decrease with a slope of −0.4 and −0.33 in log _{2}−log _{10} space, showing that convergence to the analytical solution of equilibrium can be achieved relatively quickly with increasing resolution. The runtime on a CPU (Intel i710750H, 2.60 GHz) versus a GPU (NVIDIA GeForce RTX 2070) is depicted in Fig. 6d and shows that using a GPU is advantageous for N≥64, which corresponds to the typical problem size for icesheet modelling. More specifically, the CPU and GPU computation time respectively increase with a slope of 0.73 and 0.17 in log _{2}−log _{10} space, thus giving a clear advantage to GPU computation for large problems. Thanks to the hybrid FDM–FCM scheme used to evaluate Eq. (14), the scaling of computation time on both a CPU and GPU is better than is usually obtained from FDM, FVM and FEM since all of them rely on solving a large system of linear equations.
4.2 Test 2 – 1D GIA solutions of idealised loads on a layered Earth
In Spada et al. (2011), a range of 1D GIA models are benchmarked against each other and show excellent agreement on various experiments. Here, we reproduce the benchmark tests called “1/2” (geodetic quantities) and “2/2” (geodetic rates), which are similar to Test 1 but present the following differences.

The computation domain is a stereographic projection of the spherical Earth, centred at colatitude θ=0°, and the effect of distortion is therefore included, unlike for Test 1. We apply ice loads with ρ^{ice}=931 kg m^{−3}, chosen in agreement with Spada et al. (2011), and the following geometries:
 a.
an ice cap with maximal height H_{max}=1.5 km, radius θ=10° and a shape defined by a cosine function;
 b.
a cylindrical ice load of thickness H=1 km and radius θ=10°,
 a.

Earth's structure has three layers, namely (1) a lithosphere of thickness T_{0}=70 km and shear modulus ${G}_{\mathrm{0}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{11}}$ Pa, (2) an upper mantle of thickness T_{1}=600 km and viscosity η_{1}=10^{21} Pa s, and (3) a lower mantle reaching down to the core–mantle boundary with viscosity ${\mathit{\eta}}_{\mathrm{2}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{21}}$Pa s. For any further detail, we refer the reader to the M3L70V01 profile shown in Spada et al. (2011). In FastIsostasy, these layers are translated into an elastic plate, a viscous channel and a viscous halfspace.

The results of SSH perturbation provided in Spada et al. (2011) allow us to check the validity of Eq. (19), used in FastIsostasy.
Figure 7 demonstrates that the viscous displacement, its rate and the SSH computed in FastIsostasy qualitatively follow the results of Spada et al. (2011). The latter corresponds to the outputs of PMTF, VILMA and VEENT, which show such a good agreement that they are gathered into a single output. Quantitatively, the mean displacement error between FastIsostasy and Spada et al. (2011) is small, with ${\overline{e}}^{\mathrm{abs}}\left(t\right)<\mathrm{27.0}$ m for both loading cases. In addition, the equilibrium displacement is represented well by a maximal error of ${\widehat{e}}^{\mathrm{abs}}(t=\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{kyr})<\mathrm{7}$ m and ${\widehat{e}}^{\mathrm{abs}}(t=\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{kyr})<\mathrm{12}$ m and a mean error of ${\overline{e}}^{\mathrm{abs}}(t=\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{kyr})<\mathrm{5}$ m and ${\overline{e}}^{\mathrm{abs}}(t=\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{kyr})<\mathrm{7}$ m, for the cap and the disc load respectively. The maximal difference arises at t=10 kyr for the cylindrical load and yields ${\widehat{e}}^{\mathrm{abs}}<\mathrm{47}$ m, i.e. $\widehat{e}<\mathrm{0.19}$. In both cases, the difference in vertical displacement is propagated to the computation of the SSH perturbation according to Eq. (19). As shown in the last row of Fig. 7, this leads to a maximal difference between FastIsostasy and the 1D GIA models that reaches at most 6 m but less than 2 m on average, for maximal SSH perturbations around 40 m in both models.
Since the experimental setup is as similar as possible for the 1D GIA models and FastIsostasy, these differences can be largely attributed to (1) the lumping of the depth dimension as performed in Eq. 4, which leads the two approaches to solve different equations, and (2) to the regional domain used here, which only allows for an approximate treatment of the BCs as described in Eq. (15).
4.3 Test 3 – 3D GIA solution of idealised load on an idealised LV Earth
Seakon is a global 3D GIA model that includes all the processes mentioned in Sect. 1.4.6. It solves the deformational response of the solid Earth with FVM on an unstructured grid, which is typically finer at the poles, where the bedrock displacement is largest. It has been extensively used in GIA studies (e.g. Austermann et al., 2021; Mitrovica et al., 2009; Pan et al., 2021, 2022) and coupled to an icesheet model in Gomez et al. (2018). For further details about the model and the adaptive grid refinement, we refer the reader to Latychev et al. (2005) and Gomez et al. (2018). We benchmark FastIsostasy against Seakon here for idealised cases with LV similar to that estimated across Antarctica. Here again, a cylindrical ice load with H=1 km and R=1000 km is applied on a domain with ${W}_{x}={W}_{y}=\mathrm{3000}$ km and ${N}_{x}={N}_{y}=\mathrm{128}$. To isolate the error from the lumping of the depth dimension, the vertical structure of the solid Earth is kept as simple as possible, with a mantle viscosity that is constant along z.
We distinguish four cases (a–d), which are all parameterised by a Gaussianshaped anomaly that is almost 0 on the boundary and yields its largest value at the interior of the domain. For case a (case b), this anomaly represents a decrease (increase) from T=150 km down to T=50 km (up to T=250 km) of the lithospheric thickness towards the interior of the domain. For case c (case d), this anomaly represents an exponential decrease (increase) from η=10^{21} Pa s down to η=10^{20} Pa s (up to η=10^{22} Pa s) of the mantle viscosity towards the interior of the domain. The heterogeneities (a–d) are shown in Appendix D and are used to generate results that are referred to as “LVELVA”. To quantify the improvement resulting from the use of LVELVA instead of ELVA (or ELRA), we also generate results with the nominal, homogeneous parameters $T(x,y)=\mathrm{150}$ km and $\mathit{\eta}(x,y)={\mathrm{10}}^{\mathrm{21}}$ Pa s (or τ=3000 years) and index them with “ELVA” (or “ELRA”).
As can be seen in the top and middle row of Fig. 8, LVELVA closely follows Seakon on cases a–b by showing similar timescales, amplitudes and shapes of the bedrock displacement. In the bottom row of Fig. 8, the maximal and mean relative differences respectively remain at $\widehat{e}\left(t\right)<\mathrm{0.07}$ and $\overline{e}\left(t\right)<\mathrm{0.03}$ over time. In comparison, ELVA yields similar errors for case a and slightly higher ones in case b, with values of $\widehat{e}\left(t\right)<\mathrm{0.12}$ and $\overline{e}\left(t\right)<\mathrm{0.05}$. For ELRA, the maximal error over time is slightly higher, although not significantly. We recall that the lithospheric thickness is an important control on the shape of the bedrock displacement but only an indirect one on its magnitude and timescale. This can be seen in Eq. (9), where the lithospheric rigidity D is only multiplied with spatial derivatives of the displacement. Misrepresenting the LV of lithospheric thickness therefore only has a marginal effect for cases a–b, but we emphasise that its impact on the displacement magnitude can however become important when the load presents localised features, as later shown in Test 4. Furthermore, accounting for a heterogeneous lithospheric thickness can impact the bedrock slopes significantly, which are an important control on icesheet groundingline stability and therefore on the evolution of ice sheets. Finally, as shown in an additional experiment presented in Appendix D, (LV)ELVA yields equally low error values in the absence of a lithosphere. This is the extreme case of a thin lithosphere, where the absence of flexural moments effectively decouples the displacement of neighbouring cells – a behaviour that is present in both Seakon and (LV)ELVA.
The advantage of using LVELVA over ELVA and ELRA becomes significant when studying cases c–d. For these cases, ELVA yields large transient differences compared to Seakon, with $\widehat{e}\left(t\right)<\mathrm{0.37}$ and $\overline{e}\left(t\right)<\mathrm{0.11}$. Here again, ELRA shows marginally higher error values. This clearly shows that neither ELRA nor ELVA are suited to represent the typical variations of viscosity over Antarctica. In comparison, LVELVA yields errors of only $\widehat{e}<\mathrm{0.11}$ and $\overline{e}<\mathrm{0.04}$, similar to those obtained in test cases a–b. Since these values are systematically lower than what was obtained in Test 2, it appears that the error in LVELVA mainly stems from the layered Earth structure that can only be partially accounted for on a 2D grid and not from the LV generalisation presented in Eq. (9). This is further supported by an additional test, presented in Sect. D, where Seakon and LVELVA adopt a 1D Earth structure following PREM. This leads to errors of $\widehat{e}<\mathrm{0.16}$ and $\overline{e}<\mathrm{0.06}$, which are very close to what was observed in Test 2.
For cases a–b, it should be noted that ELRA, ELVA and LVELVA present the same equilibrium state because of the constant lithospheric thickness, as can be formally deduced by setting the time derivative to 0 in the equations presented in Appendix A. We stress that the higher transient error in ELVA and ELRA can therefore be easily missed when considering equilibrium states. This is the case in Le Meur and Huybrechts (1996), where the only spatial comparison across models is made for a quasiequilibrium state. We further draw attention to the fact that the transient error metrics used in the present study are stricter than plotting the mean spatial displacement of each model over time, which is done in Le Meur and Huybrechts (1996) and which potentially hides large localised differences between ELRA and the 1D GIA model used for comparison.
Throughout Test 3, LVELVA underestimates the peripheral forebulge by about [10,15] m, which is a systematic error. When considering a layered Earth, as shown in Appendix D, this value becomes as big as 40 m for early time steps but evolves to less than 20 m at equilibrium. Konrad et al. (2014) observed a similar behaviour when comparing ELRA to a 1D GIA. This comparatively large transient is the source of the aforementioned error of $\overline{e}\left(t\right)<\mathrm{0.16}$. This most likely arises because the mantle flow contributing to the amplitude of the forebulge is not resolved in LVELVA. Since the forebulge forms in the vicinity of the ice margin, this might be an important error source to keep in mind when comparing FastIsostasy to a 3D GIA model in a coupled icesheet context, especially when studying the possibility of a forebulge feedback as proposed in Albrecht et al. (2023).
When performing Test 3, we noticed that a large lithospheric thickness, a large gradient of the lithospheric thickness or low viscosity all lead to a higher computational cost. This is consistent with theoretical insights, since all of these cases lead to a larger value of the righthand side in Eq. (14), thus making the ODE stiffer and requiring smaller time steps to resolve this with sufficient accuracy in time. For quantitative values, we refer the reader to the comparative table of runtimes provided in Appendix D.
4.4 Test 4 – 3D GIA solution of the last glacial cycle on a realistic Earth
Thus far, FastIsostasy has been tested with idealised loads and parameter fields. We now consider the more realistic case of simulating the GIA response of two different Earth structures to the last glacial cycle, as reconstructed in ICE6G_D (Peltier et al., 2018), which is an updated version of ICE6G_C (Argus et al., 2014; Peltier et al., 2015) after a mismatch with the presentday uplift was pointed out in Purcell et al. (2016). The first structure is a 1D Earth that does not present any LV. The second structure is a 3D Earth with the lithospheric thickness and the mantle viscosity fields from Pan et al. (2022), which are similar to those depicted in Fig. 2. We now compare the results of five different models: ELRA, ELVA, LVELVA, Seakon 1D (SK1D) and Seakon 3D (SK3D). It should be noted that the present comparison omits LVELRA, since its implementation goes beyond the scope of the present work.
For the regional models to be comparable with each other, ELRA, ELVA and LVELVA are coupled to ReSeLeM, which uses a BSL forcing derived from SK3D instead of Eq. (30), in accordance with the comment made in Sect. 2.5. To perform the lumping of the depth dimension for LVELVA, we define the viscous halfspace as beginning at 300 km. We observed that this yields error metrics lower than 400 and 500 km, which appears coherent since, according to Eq. (4), deeper models might overestimate the contribution of the deeper layers of the mantle to the effective viscosity. We recall that the effective viscosity should primarily capture the response to loads with characteristic lengths of continental ice sheets and that deeper layers are mostly excited by loads with larger wavelengths.
The error plot shown in Fig. 9a depicts, for all time steps, the displacement of SK3D, considered to be closest to reality, against ELRA, ELVA, SK1D and LVELVA. We hereby only represent the points within the activation mask 𝒜, represented by the black contours of Fig. 9c–h and corresponding to the typical region of interest for icesheet modellers. The position around the identity shows that ELVA leads to displacements that are biased towards lower values, especially for ${u}_{\mathrm{SK}\mathrm{3}\mathrm{D}}\le \mathrm{300}$ m, where the error comes close to e^{abs}≃130 m. Although this bias is somewhat smaller for SK1D, it still reaches similar maximal values. In comparison, LVELVA is centred around the identity and presents no such bias. This can be explained by the fact that a thinner lithosphere and a less viscous mantle in West Antarctica allows for larger transient displacements around LGM. Furthermore, the spread around the identity, especially around the lower, unbiased displacement values, is an additional metric to take into account. For ${u}_{\mathrm{SK}\mathrm{3}\mathrm{D}}\ge \mathrm{300}$ m, SK1D, LVELVA and ELVA respectively present the smallest, the intermediate and the largest spread. As expected, ELRA shows a larger bias than LVELVA but, surprisingly, a smaller one than SK1D and ELVA. However, the spread of ELRA is larger than for any other model.
Figure 9b depicts the mean and maximal relative difference for ELRA, ELVA, LVELVA and SK1D with respect to SK3D. The mean and maximal value respectively relate to the spread around the identity and the bias observed in Fig. 9a. In Fig. 9b, the error metrics are computed for the full domain and therefore include the far field, where the rotational feedback dominates the displacement. Since none of the regional models are capable of capturing this, SK1D shows the smallest mean error over time, with ${\overline{e}}_{\mathrm{SK}\mathrm{1}\mathrm{D}}\left(t\right)<\mathrm{0.01}$. In comparison, ELRA, ELVA and LVELVA show larger mean error metrics that are however similar among each other, with $\overline{e}\left(t\right)<\mathrm{0.04}$. When computing the mean error only over the active mask, this difference between SK1D and the regional models vanishes, with $\overline{e}\le \mathrm{0.01}$ for all models. We believe that this latter definition of the mean error is closer than what is relevant to most icesheet modellers, but highlighting the larger error in regional models in the far field ensures that the users of FastIsostasy are aware of this limitation. The low values of the mean error that are observed regardless of the model can be explained by the fact that most of the regions, especially the bulk of East Antarctica, can be represented well by intermediate mantle viscosity and lithospheric thickness.
Interestingly, the peak values of ${\widehat{e}}_{\mathrm{SK}\mathrm{1}\mathrm{D}}$ and ${\widehat{e}}_{\mathrm{ELVA}}$ are very close to each other, yielding about 0.22, corresponding to ${\widehat{e}}^{\mathrm{abs}}\simeq \mathrm{130}$ m. In both cases, these values are observed over $t\in [\mathrm{22},\mathrm{12}]$ kyr, which corresponds to the 10 kyr of rapid deglaciation following LGM. In comparison, ELRA presents a peak maximal error with a similar timing and a marginally smaller amplitude ${\widehat{e}}_{\mathrm{ELRA}}<\mathrm{0.19}$. However, it presents large errors over the last 6 kyr of simulation. This points out that high errors in ELRA, ELVA and SK1D are to be expected when rapid changes in ice thickness occur – a situation that could be triggered by sustained anthropogenic climate warming. For ELVA, the peak difference to SK3D is observed at $t=\mathrm{18}$ kyr. The corresponding displacement fields and their difference are plotted in Fig. 9c–e, corroborating that ELVA does not allow for enough displacement in West Antarctica after the LGM. The location of the peak error points out that the use of ELVA, coupled to an icesheet model, may lead to significant errors since the WAIS is based on a retrograde bedrock below sea level and is therefore particularly sensitive to the GIA response, as pointed out in the Introduction.
Compared to ELRA, ELVA and SK1D, LVELVA reduces the maximal error down to about ${\widehat{e}}_{\mathrm{LV}\mathrm{ELVA}}\le \mathrm{0.14}$, which corresponds to ${\widehat{e}}^{\mathrm{abs}}\le \mathrm{80}$ m. The displacement fields for the time $t=\mathrm{14}$ kyr are plotted in Fig. 9f–h and show that the nearfield displacement is reasonably well matched, even for the time step showing the largest error metrics. In particular, the displacement is slightly underestimated in most of the active mask. Since this appears to be a systematic offset, it could easily be corrected by tuning the density and/or the viscosity chosen for LVELVA.^{8} We however decide not to do so to highlight the differences to SK3D without additional tuning.
Throughout Test 4, SK3D has been used as the baseline model. In Sect. D, we show a similar analysis that compares ELRA and ELVA by using SK1D as a baseline. There it appears (1) that ELVA presents smaller errors than ELRA and (2) that LVELVA presents smaller errors compared to SK3D than ELVA compared to SK1D.
Runtime analysis. Despite the errors compared to Seakon and potentially to any 3D GIA model, we believe that FastIsostasy can be a particularly appealing tool, since the 122 kyr run takes only about 9 min to compute for a horizontal resolution of h=20 km, resulting in 350×350 grid points. For ELVA, the absence of LV leads to a reduced computation time of only about 4 min and outperforms ELRA, which requires about 6 min. These computations were performed on an NVIDIA GeForce RTX 2070, a lowcost GPU with moderate performance by the standards of 2024. Although the time stepping is adaptive, no values beyond dt=10 years are used. This potentially provides the icesheet model with an input that is very dense in time, for instance as opposed to Gomez et al. (2018).
In comparison, the Seakon simulation takes about 4.5 d on 150 CPUs with a time step of $dt\in [\mathrm{125},\mathrm{1000}]$ years. Assuming an ideal parallelisation scaling of 100 %, this corresponds to about 1×10^{6} min of CPU runtime, which roughly corresponds to 5 orders of magnitude more than what FastIsostasy requires. The models used in Albrecht et al. (2023) and Zhong et al. (2022) present smaller computation times, reducing this down to 4 orders of magnitude.
Furthermore, the power consumption of the GPU used in the present study is 185 W,^{9} compared to a typical value of more than 100 W for a single, modern CPU. As the energy consumption is expressed as the product of power and computation time, FastIsostasy appears to be less energyconsuming than Seakon by several orders of magnitude. Finally, we draw attention to the fact that the acquisition price of the GPU used here is a few hundred euros, which is far less than that of a large CPU cluster.
Of course, the runtime of Seakon and FastIsostasy are not directly comparable: Seakon solves the global GIA problem, which requires a grid with many more degrees of freedom that reaches down to the core–mantle boundary. The output of Seakon is much richer, since it includes, among others, the rotational feedback, the position of migrating shorelines, the horizontal displacement of the bedrock and the relative sea level at any point on Earth. Nonetheless, these quantities tend to be less relevant for standalone icesheet modelling, and FastIsostasy therefore offers an opportunity to regionally mimic the behaviour of a 3D GIA model at very low computational, energy and financial cost.
Throughout all the tests, FastIsostasy displays a maximal and mean error over space of $\widehat{e}\left(t\right)<\mathrm{0.2}$ and $\overline{e}\left(t\right)<\mathrm{0.05}$, both being typically much smaller for most time steps. In particular, Test 1 has shown that the discretisation error of FastIsostasy is very small and that, with an increasing problem size, its computational expense scales better than what is typically obtained with FDM, FVM and FEM solvers. Test 2 has shown that FastIsostasy represents the deformational response and the SSH perturbation with a relatively low error compared to 1D GIA models, despite solving the problem on a 2D regional grid. Test 3 and Test 4 have shown that LVELVA produces greatly reduced errors with respect to SK3D, compared to ELRA and ELVA and even to SK1D for the nearfield GIA response. This means that the model uncertainty between FastIsostasy and Seakon is smaller than the upper bound on parametric uncertainty, given by the difference between a 1D and a 3D Earth structure (Albrecht et al., 2023).
In conclusion, FastIsostasy can greatly reduce the transient error in bedrock displacement compared to ELRA and ELVA and even to 1D GIA models for regions of significant LV. This was achieved by introducing LVELVA, which generalises the work of Bueler et al. (2007) and Coulon et al. (2021) and was coupled to ReSeLeM, which regionally approximates the transient changes in ocean load. Whereas the differences between FastIsostasy and global GIA models are within the range of parametric uncertainties, the computation time is typically reduced by 3 to 5 orders of magnitude. For most icesheet models, FastIsostasy can thus represent a leap in GIA comprehensiveness at very low computational cost, even for highresolution runs on the timescale of glacial cycles. This has straightforward applications, since the GIA response is particularly relevant for the many marinebased regions of the AIS, where significant LV of the solid Earth can be present. This is the case for the Amundsen sector, which could become the largest source of sealevel rise in the coming century and therefore requires projections that account, as accurately as possible, for the stabilising effect of GIA.
Since fields of the lithospheric thickness and the uppermantle viscosity can be easily found in the literature, FastIsostasy reduces the difficulty of creating meaningful ensembles compared to relaxed rheologies (Coulon et al., 2021). The very short runtime of FastIsostasy offers an efficient method of propagating the uncertainties of the solidEarth parameters to past and future climatic scenarios.
We believe that GIA modellers, as well as the few icesheet models that are coupled to a 3D GIA model, can benefit from FastIsostasy, since it can be used as a fastprototyping tool. In particular, a scheme to tune the parameters of FastIsostasy can turn it into an emulator of a 3D GIA model with better interpretability than, for instance, machine learning techniques. Nonetheless, it should be emphasised that some scientific questions can only be answered with a global 3D GIA model and that FastIsostasy is a complementary tool that does not aim to replace it. Finally, we believe that the relatively abbreviated code of FastIsostasy and its few equations compared to 1D or 3D GIA models are particularly well suited for educational purposes.
Following Le Meur and Huybrechts (1996), the governing equations of ELRA yield
with u^{eq} as the equilibrium displacement and all further quantities already introduced in this paper. Equation (A1) can be solved by convolving the load with a Green's function derived from a Kelvin function of order 0, as described in Le Meur and Huybrechts (1996). In FastIsostasy, the convolution is performed via FFTs, which is much faster than a computation in a time domain. Coulon et al. (2021) generalised the equations of ELRA to a laterally variable lithospheric thickness and relaxation time (LVELRA):
with ${M}_{xx}^{\mathrm{eq}}$, ${M}_{xy}^{\mathrm{eq}}$ and ${M}_{yy}^{\mathrm{eq}}$ as the flexural moments at equilibrium, computed by introducing u^{eq} into Eqs. (10)–(12). For ELRA as well as LVELRA, a relaxed rheology is used, which presents the limitations explained in Sect. 1.3. ELVA, as proposed by Bueler et al. (2007), Cathles (1975), and Lingle and Clark (1985), addresses this limitation because its governing equation is directly parameterised by viscosity:
However, this assumes a constant lithospheric thickness $T(x,y)=T$ and sublithospheric mantle viscosity $\mathit{\eta}(x,y)=\mathit{\eta}$ throughout the domain, which therefore prevents the representation of LV. We tried to generalise the derivation presented in Cathles (1975) to LV; however this was unsuccessful. Instead, we combined Eqs. (A3) and (A5), thus obtaining
with the flexural moments accounting for the laterally variable lithospheric thickness. By introducing a pressure term, as done in Eq. (8), and accounting for the distortion K, we obtain the governing equation of LVELVA, as written in Eq. (9). For constant parameters, LVELVA simplifies to the equation of ELVA and can be therefore understood as an ad hoc generalisation that reduces the error made compared to a 3D GIA model. Using ELVA in FastIsostasy is simply achieved by running LVELVA with constant parameters. This seamless approach offers a code that is easier to maintain and an interface that is simpler to use, at the expense of a minor increase in computational expense compared to the use of a specific solver that takes advantage of the simplifications made in ELVA.
The function $\stackrel{\mathrm{\u0303}}{A}\left(\mathrm{s}\right):\mathbb{R}\to \mathbb{R}$ is here computed by summing the surfaces of cells situated below the BSL s, based on the 1 arcmin global topography of ETOPO1 (Amante and Eakins, 2009). Note that this slightly overestimates the ocean surface, since all regions below sea level are counted as part of the ocean, including, for instance, parts of the Netherlands. To tackle this, we introduce a bias correction scaling γ, which avoids any offset for the presentday value A_{pd} and depends on the uncorrected value $\stackrel{\mathrm{\u0303}}{A}\left(s\right)$:
To reduce the runtime, we precompute A(s) as a piecewise linear interpolator for $s\in [\mathrm{150},\mathrm{70}]$ years with a discretisation of Δs=0.1 m. The resulting function is depicted in Fig. 4c and shows that, for the range of realistic sealevel contributions over glacial cycles, the trapezoidal approximation leads to variations of the ocean surface between −7 % and +4 % around the presentday value.
Two important characteristics of the mantle have to be accounted for such that the Maxwell time $\mathit{\tau}=\mathit{\eta}\phantom{\rule{0.125em}{0ex}}{E}^{\mathrm{1}}$ of FastIsostasy is comparable to that of a 3D GIA model. This is done by introducing two correction factors. First, one of the underlying assumptions made by Cathles (1975) is that the mantle is incompressible; i.e. ν^{i}=0.5 is assumed as opposed to ν_{0}=0.28 for the elastic lithosphere. In reality, the mantle is however a compressible medium with ν^{c}≃0.28. We now look for η^{i}, the viscosity that has to be used in the incompressible case in order to match the Maxwell time of the compressible case. By introducing the shear modulus $G=E\phantom{\rule{0.125em}{0ex}}\left(\mathrm{2}\phantom{\rule{0.125em}{0ex}}\right(\mathrm{1}+\mathit{\nu}){)}^{\mathrm{1}}$, we obtain
In essence, this means that compressible media have a larger Maxwell time and that we need to slightly increase the viscosity values to render this, since Eq. (A5) is used to postulate Eq. (9) and assumes an incompressible viscous flow. This is supported by Fig. C1, which shows that a 1D GIA model displays longer decay times for a compressible mantle, compared to an incompressible one.
Second, both the shear modulus and the viscosity depend on the temperature of the medium. For instance, a positive temperature anomaly in the mantle leads to a negative anomaly of both viscosity and shear modulus. This means that the decrease in the Maxwell time due to the decrease in viscosity is somewhat compensated by the decrease in the shear modulus. We have chosen to compute this scaling by calibrating FastIsostasy to results of a 3D GIA model:
with η_{0}=10^{21} Pa s the calibration constant used throughout this work. We thus obtain a relation between the viscosity η^{c}, inferred from seismic measurements, and the corrected effective viscosity η, ultimately used in FastIsostasy:
If the depth dimension is lumped according to Eq. (B1), then the viscosity field ${\mathit{\eta}}_{\mathrm{1}}^{\mathrm{eff}}$, representing the compound of layers from l=1 to l=L, is used for η^{c}.
The anomalies of lithospheric thickness and uppermantle viscosity used in Test 3 are represented in Fig. D1 and result from a scaled Gaussian distribution 𝒩(0,σ) with $\mathit{\sigma}=(W/\mathrm{4}{)}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathit{I}$ and $\mathit{I}\in {\mathbb{R}}^{\mathrm{2}\times \mathrm{2}}$ as the identity matrix.
In addition to Test 3, we perform two simulations with laterally constant Earth structures in Seakon. The first one corresponds to PREM (Dziewonski and Anderson, 1981), and the second corresponds to a singlelayer mantle without an elastic lithosphere. The results are respectively depicted in Fig. D2a and b. Unsurprisingly, the first case shows a very similar error pattern to what is obtained in Test 2 and highlights that the main error source for FastIsostasy comes from the lumped depth dimension rather than from the generalisation of ELVA to LVELVA. The second case shows that, in the absence of a lithosphere, the match between Seakon and FastIsostasy yields $\overline{e}\left(t\right)<\mathrm{0.04}$ and $\widehat{e}\left(t\right)<\mathrm{0.08}$. In particular, this example shows that, in both models, the absence of a lithosphere effectively decouples neighbouring cells due to the absence of flexural moments.
Whereas Fig. 9 uses SK3D as a baseline for the error metrics, we propose using SK1D as a baseline to compare ELRA and ELVA in Fig. D3. In Fig. D3, it appears that ELRA is more biased towards large displacements than ELVA. Furthermore, the mean error is similar for both models, but the maximal error is higher overall for ELRA, with $\widehat{e}\left(t\right)<\mathrm{0.21}$. In comparison, the maximum error in ELVA yields $\widehat{e}\left(t\right)<\mathrm{0.16}$. The middle row of Fig. D3 highlights that the higher error in ELRA stems from an overestimated displacement in Marie Byrd Land. In comparison, the bottom row shows that ELVA displays a much more homogeneously distributed error.
FastIsostasy is available under the GNU General Public License v3.0 at https://github.com/JanJereczek/FastIsostasy.jl (last access: 4 July 2024) (Julia version) and https://github.com/palmaice/FastIsostasy (last access: 4 July 2024) (Fortran version). The data used in the present work can be found at https://github.com/JanJereczek/isostasy_data (last access: 4 July 2024). The archived versions of the code and data used for this paper can be found at https://doi.org/10.5281/zenodo.10419117 (SwierczekJereczek, 2024a) and https://doi.org/10.5281/zenodo.11175418 (SwierczekJereczek, 2024b).
Animations of the results obtained by FastIsostasy on Test 4 can be found at https://github.com/JanJereczek/FastIsostasy.jl (SwierczekJereczek, 2024c).
JSJ conceptually developed FastIsostasy as well as its Julia version. The Fortran version was codeveloped by MM, AR and JSJ. The Seakon simulations were performed by KL, who, along with JM, also provided insights into 3D GIA modelling for the purposes of this paper. This paper was prepared by JSJ with contributions from all coauthors.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We would like to thank Ed Bueler and Constantine Khroulev for the helpful email exchange. Their openly accessible implementations of ELVA, available at https://github.com/pism and https://github.com/bueler/fastearth, greatly eased the initial phase of FastIsostasy's development and inspired its name. We would also like to thank Douglas Wiens, Ana Negredo and Javier Fullea for providing valuable comments and/or data. Particular thanks goes to the reviewers of this paper, who provided many helpful comments that improved this paper. Finally, Jan SwierczekJereczek wants to thank the whole PalMA team for providing their daily support and for creating a fantastic working environment.
Jan SwierczekJereczek received funding from CriticalEarth (grant no. 956170), an EU H2020 research infrastructure of the European Commission. Alexander Robinson received funding from the European Union (ERC, FORCLIMA, grant no. 101044247). Marisa Montoya received funding from the Spanish Ministry of Science, Innovation and Universities (project MARINE, grant no. PID2020117768RBI00).
This paper was edited by Danilo Mello and reviewed by Caroline van Calcar and three anonymous referees.
A, G., Wahr, J., and Zhong, S.: Computations of the viscoelastic response of a 3D compressible Earth to surface loading: an application to Glacial Isostatic Adjustment in Antarctica and Canada, Geophys. J. Int., 192, 557–572, https://doi.org/10.1093/gji/ggs030, 2013. a
Accardo, N. J., Wiens, D. A., Hernandez, S., Aster, R. C., Nyblade, A., Huerta, A., Anandakrishnan, S., Wilson, T., Heeszel, D. S., and Dalziel, I. W. D.: Upper mantle seismic anisotropy beneath the West Antarctic Rift System and surrounding region from shear wave splitting analysis, Geophys. J. Int., 198, 414–429, https://doi.org/10.1093/gji/ggu117, 2014. a
Adhikari, S., Ivins, E. R., Larour, E., Seroussi, H., Morlighem, M., and Nowicki, S.: Future Antarctic bed topography and its implications for ice sheet dynamics, Solid Earth, 5, 569–584, https://doi.org/10.5194/se55692014, 2014. a
Adhikari, S., Ivins, E. R., Larour, E., Caron, L., and Seroussi, H.: A kinematic formalism for tracking ice–ocean mass exchange on the Earth's surface and estimating sealevel change, The Cryosphere, 14, 2819–2833, https://doi.org/10.5194/tc1428192020, 2020. a
Albrecht, T., Bagge, M., and Klemann, V.: Feedback mechanisms controlling Antarctic glacial cycle dynamics simulated with a coupled ice sheet–solid Earth model, EGUsphere [preprint], https://doi.org/10.5194/egusphere20232990, 2023. a, b, c, d, e, f, g
Amante, C. and Eakins, B. E.: ETOPO1 1 arcminute global relief model: procedure, data sources and analysis, Tech. rep., NOAA, https://doi.org/10.7289/V5C8276M, 2009. a
Amelung, F. and Wolf, D.: Viscoelastic perturbations of the earth: significance of the incremental gravitational force in models of glacial isostasy, Geophys. J. Int., 117, 864–879, https://doi.org/10.1111/j.1365246X.1994.tb02476.x, 1994. a
Argus, D. F., Peltier, W. R., Drummond, R., and Moore, A. W.: The Antarctica component of postglacial rebound model ICE6G_C (VM5a) based on GPS positioning, exposure age dating of ice thicknesses, and relative sea level histories, Geophys. J. Int., 198, 537–563, https://doi.org/10.1093/gji/ggu140, 2014. a
Austermann, J., Hoggard, M. J., Latychev, K., Richards, F. D., and Mitrovica, J. X.: The effect of lateral variations in Earth structure on Last Interglacial sea level, Geophys. J. Int., 227, 1938–1960, https://doi.org/10.1093/gji/ggab289, 2021. a
Bagge, M., Klemann, V., Steinberger, B., Latinović, M., and Thomas, M.: Glacial‐Isostatic Adjustment Models Using Geodynamically Constrained 3D Earth Structures, Geochem. Geophy. Geosy., 22, e2021GC009853, https://doi.org/10.1029/2021GC009853, 2021.Please provide article number or page range. a
Barletta, V. R., Bevis, M., Smith, B. E., Wilson, T., Brown, A., Bordoni, A., Willis, M., Khan, S. A., RoviraNavarro, M., Dalziel, I., Smalley, R., Kendrick, E., Konfal, S., Caccamise, D. J., Aster, R. C., Nyblade, A., and Wiens, D. A.: Observed rapid bedrock uplift in Amundsen Sea Embayment promotes icesheet stability, Science, 360, 1335–1339, https://doi.org/10.1126/science.aao1447, 2018. a, b
Beghein, C., Trampert, J., and van Heijst, H. J.: Radial anisotropy in seismic reference models of the mantle: ANISOTROPY IN MANTLE MODELS, J. Geophys. Res.Sol. Ea., 111, B02303, https://doi.org/10.1029/2005JB003728, 2006. a
Behrendt, J. C.: Crustal and lithospheric structure of the West Antarctic Rift System from geophysical investigations – a review, Global Planet. Change, 23, 25–44, https://doi.org/10.1016/S09218181(99)000491, 1999. a
Besard, T., Foket, C., and De Sutter, B.: Effective Extensible Programming: Unleashing Julia on GPUs, IEEE T. Parall. Distr., 30, 827–841, https://doi.org/10.1109/TPDS.2018.2872064, 2019. a
Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B.: Julia: A Fresh Approach to Numerical Computing, SIAM Rev., 59, 65–98, https://doi.org/10.1137/141000671, 2017. a
Bogacki, P. and Shampine, L.: An efficient RungeKutta (4,5) pair, Comput. Math. Appl., 32, 15–28, https://doi.org/10.1016/08981221(96)001411, 1996. a
Book, C., Hoffman, M. J., Kachuck, S. B., Hillebrand, T. R., Price, S. F., Perego, M., and Bassis, J. N.: Stabilizing effect of bedrock uplift on retreat of Thwaites Glacier, Antarctica, at centennial timescales, Earth Planet. Sc. Lett., 597, 117798, https://doi.org/10.1016/j.epsl.2022.117798, 2022. a, b, c
Bueler, E., Lingle, C. S., and Brown, J.: Fast computation of a viscoelastic deformable Earth model for icesheet simulations, Ann. Glaciol., 46, 97–105, https://doi.org/10.3189/172756407782871567, 2007. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s
Caron, L., Métivier, L., GreffLefftz, M., Fleitout, L., and Rouby, H.: Inverting Glacial Isostatic Adjustment signal using Bayesian framework and two linearly relaxing rheologies, Geophys. J. Int., 209, 1126–1147, https://doi.org/10.1093/gji/ggx083, 2017. a
Cathles, L. M.: Viscosity of the Earth's Mantle, Princeton University Press, ISBN 0691081409, 1975. a, b, c, d, e, f, g, h
Cogley, J. G.: Area of the Ocean, Mar. Geod., 35, 379–388, https://doi.org/10.1080/01490419.2012.709476, 2012. a
Coulon, V., Bulthuis, K., Whitehouse, P. L., Sun, S., Haubner, K., Zipf, L., and Pattyn, F.: Contrasting Response of West and East Antarctic Ice Sheets to Glacial Isostatic Adjustment, J. Geophys. Res.Earth, 126, e2020JF006003, https://doi.org/10.1029/2020JF006003, 2021. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s
De Boer, B., Stocchi, P., Whitehouse, P. L., and Van De Wal, R. S.: Current state and future perspectives on coupled icesheet – sealevel modelling, Quaternary Sci. Rev., 169, 13–28, https://doi.org/10.1016/j.quascirev.2017.05.013, 2017. a
DeConto, R. M. and Pollard, D.: Contribution of Antarctica to past and future sealevel rise, Nature, 531, 591–597, https://doi.org/10.1038/nature17145, 2016. a
Dziewonski, A. M. and Anderson, D. L.: Preliminary Earth Model, Phys. Earth Planet. Int., 25, 297–356, 1981. a, b, c
Farrell, W. E.: Deformation of the Earth by surface loads, Rev. Geophys., 10, 761, https://doi.org/10.1029/RG010i003p00761, 1972. a
Farrell, W. E. and Clark, J. A.: On Postglacial Sea Level, Geophys. J. Roy. Astron. Soc., 46, 647–667, https://doi.org/10.1111/j.1365246X.1976.tb01252.x, 1976. a
Fretwell, P., Pritchard, H. D., Vaughan, D. G., Bamber, J. L., Barrand, N. E., Bell, R., Bianchi, C., Bingham, R. G., Blankenship, D. D., Casassa, G., Catania, G., Callens, D., Conway, H., Cook, A. J., Corr, H. F. J., Damaske, D., Damm, V., Ferraccioli, F., Forsberg, R., Fujita, S., Gim, Y., Gogineni, P., Griggs, J. A., Hindmarsh, R. C. A., Holmlund, P., Holt, J. W., Jacobel, R. W., Jenkins, A., Jokat, W., Jordan, T., King, E. C., Kohler, J., Krabill, W., RigerKusk, M., Langley, K. A., Leitchenkov, G., Leuschen, C., Luyendyk, B. P., Matsuoka, K., Mouginot, J., Nitsche, F. O., Nogi, Y., Nost, O. A., Popov, S. V., Rignot, E., Rippin, D. M., Rivera, A., Roberts, J., Ross, N., Siegert, M. J., Smith, A. M., Steinhage, D., Studinger, M., Sun, B., Tinto, B. K., Welch, B. C., Wilson, D., Young, D. A., Xiangbin, C., and Zirizzotti, A.: Bedmap2: improved ice bed, surface and thickness datasets for Antarctica, The Cryosphere, 7, 375–393, https://doi.org/10.5194/tc73752013, 2013. a
Frigo, M. and Johnson, S.: The Design and Implementation of FFTW3, P. IEEE, 93, 216–231, https://doi.org/10.1109/JPROC.2004.840301, 2005. a
Garbe, J., Albrecht, T., Levermann, A., Donges, J. F., and Winkelmann, R.: The hysteresis of the Antarctic Ice Sheet, Nature, 585, 538–544, https://doi.org/10.1038/s4158602027275, 2020. a
Gasperini, P., Dal Forno, G., and Boschi, E.: Linear or nonlinear rheology in the Earth's mantle: the prevalence of powerlaw creep in the postglacial isostatic readjustment of Laurentia, Geophys. J. Int., 157, 1297–1302, https://doi.org/10.1111/j.1365246X.2004.02319.x, 2004. a
Goelzer, H., Coulon, V., Pattyn, F., de Boer, B., and van de Wal, R.: Brief communication: On calculating the sealevel contribution in marine icesheet models, The Cryosphere, 14, 833–840, https://doi.org/10.5194/tc148332020, 2020. a, b, c, d
Gomez, N., Mitrovica, J. X., Huybers, P., and Clark, P. U.: Sea level as a stabilizing factor for marineicesheet grounding lines, Nat. Geosci., 3, 850–853, https://doi.org/10.1038/ngeo1012, 2010. a, b, c
Gomez, N., Pollard, D., Mitrovica, J. X., Huybers, P., and Clark, P. U.: Evolution of a coupled marine ice sheetsea level model, J. Geophys. Res.Earth, 117, F01013, https://doi.org/10.1029/2011JF002128, 2012. a, b
Gomez, N., Pollard, D., and Holland, D.: Sealevel feedback lowers projections of future Antarctic IceSheet mass loss, Nat. Commun., 6, 8798, https://doi.org/10.1038/ncomms9798, 2015. a, b
Gomez, N., Latychev, K., and Pollard, D.: A Coupled Ice Sheet–Sea Level Model Incorporating 3D Earth Structure: Variations in Antarctica during the Last Deglacial Retreat, J. Climate, 31, 4041–4054, https://doi.org/10.1175/JCLID170352.1, 2018. a, b, c, d, e, f, g, h
Gregory, J. M., Griffies, S. M., Hughes, C. W., Lowe, J. A., Church, J. A., Fukimori, I., Gomez, N., Kopp, R. E., Landerer, F., Cozannet, G. L., Ponte, R. M., Stammer, D., Tamisiea, M. E., and Van De Wal, R. S. W.: Concepts and Terminology for Sea Level: Mean, Variability and Change, Both Local and Global, Surv. Geophys., 40, 1251–1289, https://doi.org/10.1007/s1071201909525z, 2019. a
Gudmundsson, G. H., Krug, J., Durand, G., Favier, L., and Gagliardini, O.: The stability of grounding lines on retrograde slopes, The Cryosphere, 6, 1497–1505, https://doi.org/10.5194/tc614972012, 2012. a
Heeszel, D. S., Wiens, D. A., Anandakrishnan, S., Aster, R. C., Dalziel, I. W. D., Huerta, A. D., Nyblade, A. A., Wilson, T. J., and Winberry, J. P.: Upper mantle structure of central and West Antarctica from array analysis of Rayleigh wave phase velocities, J. Geophys. Res.Sol. Ea., 121, 1758–1775, https://doi.org/10.1002/2015JB012616, 2016. a, b
Huang, P., Steffen, R., Steffen, H., Klemann, V., Wu, P., van der Wal, W., Martinec, Z., and Tanaka, Y.: A commercial finite element approach to modelling Glacial Isostatic Adjustment on spherical selfgravitating compressible earth models, Geophys. J. Int., 235, 2231–2256, https://doi.org/10.1093/gji/ggad354, 2023. a, b
Ivins, E. R., Caron, L., Adhikari, S., and Larour, E.: Notes on a compressible extended Burgers model of rheology, Geophys. J. Int., 228, 1975–1991, https://doi.org/10.1093/gji/ggab452, 2021. a
Ivins, E. R., van der Wal, W., Wiens, D. A., Lloyd, A. J., and Caron, L.: Antarctic upper mantle rheology, Geological Society, London, Memoirs, 56, M56–2020–19, https://doi.org/10.1144/M56202019, 2022. a, b, c
Kachuck, S. B.: giapy: Glacial Isostatic Adjustment in PYthon (1.0.0), https://github.com/skachuck/giapy/ (last access: 4 July 2024), 2017. a
Kachuck, S. B., Martin, D. F., Bassis, J. N., and Price, S. F.: Rapid Viscoelastic Deformation Slows Marine Ice Sheet Instability at Pine Island Glacier, Geophys. Res. Lett., 47, e2019GL086446, https://doi.org/10.1029/2019GL086446, 2020. a, b, c
Kang, K., Zhong, S., Geruo, A., and Mao, W.: The effects of nonNewtonian rheology in the upper mantle on relative sea level change and geodetic observables induced by glacial isostatic adjustment process, Geophys. J. Int., 228, 1887–1906, https://doi.org/10.1093/gji/ggab428, 2021. a
Kendall, R. A., Mitrovica, J. X., and Milne, G. A.: On postglacial sea level – II. Numerical formulation and comparative results on spherically symmetric models, Geophys. J. Int., 161, 679–706, https://doi.org/10.1111/j.1365246X.2005.02553.x, 2005. a, b
Klemann, V., Martinec, Z., and Ivins, E. R.: Glacial isostasy and plate motion, J. Geodynam., 46, 95–103, https://doi.org/10.1016/j.jog.2008.04.005, 2008. a
Konrad, H., Thoma, M., Sasgen, I., Klemann, V., Grosfeld, K., Barbi, D., and Martinec, Z.: The Deformational Response of a Viscoelastic Solid Earth Model Coupled to a Thermomechanical Ice Sheet Model, Surv. Geophys., 35, 1441–1458, https://doi.org/10.1007/s1071201392578, 2014. a, b, c
Konrad, H., Sasgen, I., Pollard, D., and Klemann, V.: Potential of the solidEarth response for limiting longterm West Antarctic Ice Sheet retreat in a warming climate, Earth Planet. Sc. Lett., 432, 254–264, https://doi.org/10.1016/j.epsl.2015.10.008, 2015. a
Konrad, H., Sasgen, I., Klemann, V., Thoma, M., Grosfeld, K., and Martinec, Z.: Sensitivity of groundingline dynamics to viscoelastic deformation of the solidearth in an idealized scenario, Polarforschung, Bremerhaven, Alfred Wegener Institute for Polar and Marine Research & German Society of Polar Research, 85, 89–99, https://doi.org/10.2312/POLFOR.2016.005, 2016. a
Kreuzer, M., Albrecht, T., Nicola, L., Reese, R., and Winkelmann, R.: Oceanic gateways in Antarctica – Impact of relative sealevel change on subshelf melt, EGUsphere [preprint], https://doi.org/10.5194/egusphere20232737, 2023. a
Kulp, S. A. and Strauss, B. H.: New elevation data triple estimates of global vulnerability to sealevel rise and coastal flooding, Nat. Commun., 10, 4844, https://doi.org/10.1038/s4146701912808z, 2019. a
Latychev, K., Mitrovica, J. X., Tromp, J., Tamisiea, M. E., Komatitsch, D., and Christara, C. C.: Glacial isostatic adjustment on 3D Earth models: a finitevolume formulation, Geophys. J. Int., 161, 421–444, https://doi.org/10.1111/j.1365246X.2005.02536.x, 2005. a, b, c, d, e
Le Meur, E. and Huybrechts, P.: A comparison of different ways of dealing with isostasy: examples from modelling the Antarctic ice sheet during the last glacial cycle, Ann. Glaciol., 23, 309–317, https://doi.org/10.3189/S0260305500013586, 1996. a, b, c, d, e, f, g
Lingle, C. S. and Clark, J. A.: A numerical model of interactions between a marine ice sheet and the solid earth: Application to a West Antarctic ice stream, J. Geophys. Res., 90, 1100, https://doi.org/10.1029/JC090iC01p01100, 1985. a, b, c, d
Lipscomb, W. H., Price, S. F., Hoffman, M. J., Leguy, G. R., Bennett, A. R., Bradley, S. L., Evans, K. J., Fyke, J. G., Kennedy, J. H., Perego, M., Ranken, D. M., Sacks, W. J., Salinger, A. G., Vargo, L. J., and Worley, P. H.: Description and evaluation of the Community Ice Sheet Model (CISM) v2.1, Geosci. Model Dev., 12, 387–424, https://doi.org/10.5194/gmd123872019, 2019. a, b
Lloyd, A. J., Wiens, D. A., Nyblade, A. A., Anandakrishnan, S., Aster, R. C., Huerta, A. D., Wilson, T. J., Dalziel, I. W. D., Shore, P. J., and Zhao, D.: A seismic transect across West Antarctica: Evidence for mantle thermal anomalies beneath the Bentley Subglacial Trench and the Marie Byrd Land Dome, J. Geophys. Res.Sol. Ea., 120, 8439–8460, https://doi.org/10.1002/2015JB012455, 2015. a
Lloyd, A. J., Wiens, D. A., Zhu, H., Tromp, J., Nyblade, A. A., Aster, R. C., Hansen, S. E., Dalziel, I. W. D., Wilson, T., Ivins, E. R., and O’Donnell, J. P.: Seismic Structure of the Antarctic Upper Mantle Imaged with Adjoint Tomography, J. Geophys. Res.Sol. Ea., 125, 2019JB017823, https://doi.org/10.1029/2019JB017823, 2020. a
Martinec, Z.: Spectralfinite element approach to threedimensional viscoelastic relaxation in a spherical earth, Geophys. J. Int., 142, 117–141, https://doi.org/10.1046/j.1365246x.2000.00138.x, 2000. a
Martinec, Z., Klemann, V., van der Wal, W., Riva, R. E. M., Spada, G., Sun, Y., Melini, D., Kachuck, S. B., Barletta, V., Simon, K., A, G., and James, T. S.: A benchmark study of numerical implementations of the sea level equation in GIA modelling, Geophys. J. Int., 215, 389–414, https://doi.org/10.1093/gji/ggy280, 2018. a
Mitrovica, J. X. and Milne, G. A.: On postglacial sea level: I. General theory, Geophys. J. Int., 154, 253–267, https://doi.org/10.1046/j.1365246X.2003.01942.x, 2003. a, b
Mitrovica, J. X., Milne, G. A., and Davis, J. L.: Glacial isostatic adjustment on a rotating earth, Geophys. J. Int., 147, 562–578, https://doi.org/10.1046/j.1365246x.2001.01550.x, 2001. a
Mitrovica, J. X., Gomez, N., and Clark, P. U.: The SeaLevel Fingerprint of West Antarctic Collapse, Science, 323, 753–753, https://doi.org/10.1126/science.1166510, 2009. a
Morelli, A. and Danesi, S.: Seismological imaging of the Antarctic continental lithosphere: a review, Global Planet. Change, 42, 155–165, https://doi.org/10.1016/j.gloplacha.2003.12.005, 2004. a
Morlighem, M., Rignot, E., Binder, T., Blankenship, D., Drews, R., Eagles, G., Eisen, O., Ferraccioli, F., Forsberg, R., Fretwell, P., Goel, V., Greenbaum, J. S., Gudmundsson, H., Guo, J., Helm, V., Hofstede, C., Howat, I., Humbert, A., Jokat, W., Karlsson, N. B., Lee, W. S., Matsuoka, K., Millan, R., Mouginot, J., Paden, J., Pattyn, F., Roberts, J., Rosier, S., Ruppel, A., Seroussi, H., Smith, E. C., Steinhage, D., Sun, B., Broeke, M. R. V. D., Ommen, T. D. V., Wessem, M. V., and Young, D. A.: Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet, Nat. Geosci., 13, 132–137, https://doi.org/10.1038/s4156101905108, 2020. a
Nield, G. A., Barletta, V. R., Bordoni, A., King, M. A., Whitehouse, P. L., Clarke, P. J., Domack, E., Scambos, T. A., and Berthier, E.: Rapid bedrock uplift in the Antarctic Peninsula explained by viscoelastic response to recent ice unloading, Earth Planet. Sc. Lett., 397, 32–41, https://doi.org/10.1016/j.epsl.2014.04.019, 2014. a, b
Nield, G. A., Whitehouse, P. L., van der Wal, W., Blank, B., O'Donnell, J. P., and Stuart, G. W.: The impact of lateral variations in lithospheric thickness on glacial isostatic adjustment in West Antarctica, Geophys. J. Int., 214, 811–824, https://doi.org/10.1093/gji/ggy158, 2018. a, b, c
Pan, L., Powell, E. M., Latychev, K., Mitrovica, J. X., Creveling, Gomez, N., Hoggard, M. J., and Clark, P. J.: Rapid postglacial rebound amplifies global sea level rise following West Antarctic Ice Sheet collapse, Sci. Adv., 7, eabf7787, https://doi.org/10.1126/sciadv.abf7787, 2021. a
Pan, L., Milne, G. A., Latychev, K., Goldberg, S. L., Austermann, J., Hoggard, M. J., and Mitrovica, J. X.: The influence of lateral Earth structure on inferences of global ice volume during the Last Glacial Maximum, Quaternary Sci. Rev., 290, 107644, https://doi.org/10.1016/j.quascirev.2022.107644, 2022. a, b, c, d, e
Pattyn, F.: Sealevel response to melting of Antarctic ice shelves on multicentennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0), The Cryosphere, 11, 1851–1878, https://doi.org/10.5194/tc1118512017, 2017. a
Peltier, W. R., Argus, D. F., and Drummond, R.: Space geodesy constrains ice age terminal deglaciation: The global ICE‐6G_C (VM5a) model, J. Geophys. Res.Sol. Ea., 120, 450–487, https://doi.org/10.1002/2014JB011176, 2015. a
Peltier, W. R., Argus, D. F., and Drummond, R.: Comment on “An Assessment of the ICE‐6G_C (VM5a) Glacial Isostatic Adjustment Model” by Purcell et al., J. Geophys. Res.Sol. Ea., 123, 2019–2028, https://doi.org/10.1002/2016JB013844, 2018. a, b, c
Purcell, A.: The significance of prestress advection and internal buoyancy in the flatEarth formulation, in: Dynamics of the Ice Age Earth – A Modern Perspective, Trans Tech Publications LTD, 24, 105–122, ISBN 9780878498109, 1998. a
Purcell, A., Tregoning, P., and Dehecq, A.: An assessment of the ICE6G_C(VM5a) glacial isostatic adjustment model, J. Geophys. Res.Sol. Ea., 121, 3939–3950, https://doi.org/10.1002/2015JB012742, 2016. a
Quiquet, A., Dumas, C., Ritz, C., Peyaud, V., and Roche, D. M.: The GRISLI ice sheet model (version 2.0): calibration and validation for multimillennial changes of the Antarctic ice sheet, Geosci. Model Dev., 11, 5003–5025, https://doi.org/10.5194/gmd1150032018, 2018. a
Rackauckas, C. and Nie, Q.: DifferentialEquations.jl – A Performant and FeatureRich Ecosystem for Solving Differential Equations in Julia, J. Open Res. Softw., 5, 15, https://doi.org/10.5334/jors.151, 2017. a
Robinson, A., AlvarezSolas, J., Montoya, M., Goelzer, H., Greve, R., and Ritz, C.: Description and validation of the icesheet model Yelmo (version 1.0), Geosci. Model Dev., 13, 2805–2823, https://doi.org/10.5194/gmd1328052020, 2020. a, b, c
Rückamp, M., Greve, R., and Humbert, A.: Comparative simulations of the evolution of the Greenland ice sheet under simplified Paris Agreement scenarios with the models SICOPOLIS and ISSM, Polar Sci., 21, 14–25, https://doi.org/10.1016/j.polar.2018.12.003, 2019. a
Sasgen, I., MartínEspañol, A., Horvath, A., Klemann, V., Petrie, E. J., Wouters, B., Horwath, M., Pail, R., Bamber, J. L., Clarke, P. J., Konrad, H., Wilson, T., and Drinkwater, M. R.: Altimetry, gravimetry, GPS and viscoelastic modeling data for the joint inversion for glacial isostatic adjustment in Antarctica (ESA STSE Project REGINA), Earth Syst. Sci. Data, 10, 493–523, https://doi.org/10.5194/essd104932018, 2018. a
Seroussi, H., Nowicki, S., Payne, A. J., Goelzer, H., Lipscomb, W. H., AbeOuchi, A., Agosta, C., Albrecht, T., AsayDavis, X., Barthel, A., Calov, R., Cullather, R., Dumas, C., GaltonFenzi, B. K., Gladstone, R., Golledge, N. R., Gregory, J. M., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huybrechts, P., Jourdain, N. C., Kleiner, T., Larour, E., Leguy, G. R., Lowry, D. P., Little, C. M., Morlighem, M., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Reese, R., Schlegel, N.J., Shepherd, A., Simon, E., Smith, R. S., Straneo, F., Sun, S., Trusel, L. D., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., Zhao, C., Zhang, T., and Zwinger, T.: ISMIP6 Antarctica: a multimodel ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070, https://doi.org/10.5194/tc1430332020, 2020. a
Spada, G. and Melini, D.: SELEN^{4} (SELEN version 4.0): a Fortran program for solving the gravitationally and topographically selfconsistent sealevel equation in glacial isostatic adjustment modeling, Geosci. Model Dev., 12, 5055–5075, https://doi.org/10.5194/gmd1250552019, 2019. a, b, c, d
Spada, G., Barletta, V. R., Klemann, V., Riva, R. E., Martinec, Z., Gasperini, P., Lund, B., Wolf, D., Vermeersen, L. L. A., and King, M. A.: A benchmark study for glacial isostatic adjustment codes: A GIA benchmark study, Geophys. J. Int., 185, 106–132, https://doi.org/10.1111/j.1365246X.2011.04952.x, 2011. a, b, c, d, e, f, g, h, i
SwierczekJereczek, J.: JanJereczek/FastIsostasy.jl: v1.0.0, Zenodo [code], https://doi.org/10.5281/zenodo.11277347, 2024a. a
SwierczekJereczek, J.: JanJereczek/isostasy_data: Version 0.2 (v0.2.0). Zenodo [data set], https://doi.org/10.5281/zenodo.11175418, 2024b. a
SwierczekJereczek, J.: Simulating the GIA response of Antarctica during the last glacial cycle with FastIsosatsy, GitHub [video], https://github.com/JanJereczek/FastIsostasy.j, last access: 4 July 2024c. a
van Calcar, C. J., van de Wal, R. S. W., Blank, B., de Boer, B., and van der Wal, W.: Simulation of a fully coupled 3D glacial isostatic adjustment – ice sheet model for the Antarctic ice sheet over a glacial cycle, Geosci. Model Dev., 16, 5473–5492, https://doi.org/10.5194/gmd1654732023, 2023. a, b, c
van der Wal, W., Wu, P., Wang, H., and Sideris, M. G.: Sea levels and uplift rate from composite rheology in glacial isostatic adjustment modeling, J. Geodynam., 50, 38–48, https://doi.org/10.1016/j.jog.2010.01.006, 2010. a
van der Wal, W., Whitehouse, P. L., and Schrama, E. J. O.: Effect of GIA models with 3D composite mantle viscosity on GRACE mass balance estimates for Antarctica, Earth Planet. Sc. Lett., 414, 134–143, https://doi.org/10.1016/j.epsl.2015.01.001, 2015. a, b
Ventsel, E. and Krauthammer, T.: Thin Plates and Shells: Theory, Analysis, and Applications, CRC Press, Boca Raton, 1st Edn., ISBN 9780429221316, 2001. a, b
Weerdesteijn, M. F. M., Naliboff, J. B., Conrad, C. P., Reusen, J. M., Steffen, R., Heister, T., and Zhang, J.: Modeling Viscoelastic Solid Earth Deformation Due To Ice Age and Contemporary Glacial Mass Changes in ASPECT, Geochem. Geophy. Geosy., 24, e2022GC010813, https://doi.org/10.1029/2022GC010813, 2023. a, b, c, d, e
Whitehouse, P. L., Gomez, N., King, M. A., and Wiens, D. A.: Solid Earth change and the evolution of the Antarctic Ice Sheet, Nat.e Commun., 10, 503, https://doi.org/10.1038/s4146701808068y, 2019. a, b, c, d
Wiens, D. A., Shen, W., and Lloyd, A. J.: The seismic structure of the Antarctic upper mantle, Geological Society, London, Memoirs, 56, M56–2020–18, https://doi.org/10.1144/M56202018, 2022. a
Winkelmann, R., Martin, M. A., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISMPIK) – Part 1: Model description, The Cryosphere, 5, 715–726, https://doi.org/10.5194/tc57152011, 2011. a
Wu, P.: Using commercial finite element packages for the study of earth deformations, sea levels and the state of stress, Geophys. J. Int., 158, 401–408, https://doi.org/10.1111/j.1365246X.2004.02338.x, 2004. a
Wu, P. and Wang, H.: Effects of mode coupling and location of rotational axis on glacial induced rotational deformation in a laterally heterogeneous viscoelastic earth, Geophys. J. Int., 167, 853–859, https://doi.org/10.1111/j.1365246X.2006.03103.x, 2006. a
Zhong, S., Kang, K., A, G., and Qin, C.: CitcomSVE: A Three‐Dimensional Finite Element Software Package for Modeling Planetary Mantle’s Viscoelastic Deformation in Response to Surface and Tidal Loads, Geochem. Geophy. Geosy., 23, e2022GC010359, https://doi.org/10.1029/2022GC010359, 2022. a, b, c, d, e, f, g
Some of them can, however, be obtained upon request, such as Seakon (Latychev et al., 2005).
In cratonic regions with a thick lithosphere (e.g. East Antarctica) there might be no asthenosphere. This does not prevent the use of ELRA, which is suited, more generally, to approximate the sublithospheric mantle.
A viscous channel has a finite thickness, unlike a viscous halfspace, which has a infinite thickness.
FastIsostasy's interface already accepts external forcing from sediments, but they will be ignored for the present work.
The distortion K does not appear in σ^{zz} since it cancels out when computing the volumetoarea ratio.
This is not completely correct since the ocean load changes at the domain margin, which is however impossible to represent accurately in a regional model.
We refer the reader to Adhikari et al. (2020) for an alternative treatment.
This could simply be done by hand or, for instance, with an unscented Kalman inversion, as shown in the code documentation.
This tabulated value is found at https://www.nvidia.com/eses/geforce/graphicscards/compare/?section=compare20 (last access: 4 July 2024).
 Abstract
 Introduction
 Model description
 Implementation, performance and further remarks
 Model validation and benchmarks
 Conclusions
 Appendix A: From ELRA to LVELVA
 Appendix B: Ocean surface function
 Appendix C: Scaling the effective viscosity
 Appendix D: Complementary information on Test 3 and Test 4
 Code and data availability
 Video supplement
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Model description
 Implementation, performance and further remarks
 Model validation and benchmarks
 Conclusions
 Appendix A: From ELRA to LVELVA
 Appendix B: Ocean surface function
 Appendix C: Scaling the effective viscosity
 Appendix D: Complementary information on Test 3 and Test 4
 Code and data availability
 Video supplement
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References