SELEN⁴ (SELEN version 4.0): a Fortran program for solving the gravitationally and topographically self-consistent sea-level equation in glacial isostatic adjustment modeling

We present SELEN⁴ (SealEveL EquatioN solver), an open-source program written in Fortran 90 that simulates the glacial isostatic adjustment (GIA) process in response to the melting of the Late Pleistocene ice sheets. Using a pseudo-spectral approach complemented by a spatial discretization on an icosahedron-based spherical geodesic grid, SELEN⁴ solves a generalized sea-level equation (SLE) for a spherically symmetric Earth with linear viscoelastic rheology, taking the migration of the shorelines and the rotational feedback on sea level into account. The approach is gravitationally and topographically self-consistent, since it considers the gravitational interactions between the solid Earth, the cryosphere, and the oceans, and it accounts for the evolution of the Earth's topography in response to changes in sea level. The SELEN⁴ program can be employed to study a broad range of geophysical effects of GIA, including past relative sea-level variations induced by the melting of the Late Pleistocene ice sheets, the time evolution of paleogeography and of the ocean function since the Last Glacial Maximum, the history of the Earth's rotational variations, present-day geodetic signals observed by Global Navigation Satellite Systems, and gravity field variations detected by satellite gravity missions like GRACE (the Gravity Recovery and Climate Experiment). The “GIA fingerprints” constitute a standard output of SELEN⁴. Along with the source code, we provide a supplementary document with a full account of the theory, some numerical results obtained from a standard run, and a user guide. Originally, the SELEN program was conceived by Giorgio Spada (GS) in 2005 as a tool for students eager to learn about GIA, and it has been the first SLE solver made available to the community.

1 Introduction

In the last few decades, glacial isostatic adjustment (GIA) modeling has progressively gained a central role in the study of contemporary sea-level change. Sea-level variations observed at the tide gauges deployed along the world's coastlines need to be corrected for the effect of GIA to assess the impact of global warming. As discussed in the review of Spada and Galassi (2012), a precise estimate of global sea-level rise has been possible only after Peltier and Tushingham (1989) first solved the sea-level equation (SLE) using an appropriate spatial resolution, building upon the seminal papers of Farrell and Clark (1976) and Clark et al. (1978). Since then, a number of GIA models characterized by different assumptions about the Earth's rheological profile and the history of the Late Pleistocene ice sheets have been proposed, constrained by sea-level proxies dating from the Last Glacial Maximum (LGM, 21 000 years ago). For a review of the development of GIA modeling, the reader is referred to Whitehouse (2009), Spada (2017), and Whitehouse (2018). GIA models have provided increasingly accurate estimates of global mean secular sea-level rise (a summary is given in Table 1 of Spada and Galassi, 2012) but also have the potential of describing the patterns of future trends of sea level in a global change scenario (see, e.g., Bamber et al., 2009; Spada et al., 2013). Since the beginning of the “altimetry era” (1992–today) and the launch of the Gravity Recovery and Climate Experiment (GRACE; see Wahr et al., 1998) in 2002, GIA modeling has regained momentum, providing the tools for isolating the effects of global warming (i) from absolute sea-level data (Nerem et al., 2010; Cazenave and Llovel, 2010) and (ii) from the Stokes coefficients of the gravity field (see Leuliette and Miller, 2009; Cazenave et al., 2009; Chambers et al., 2010; WCRP, 2018) to infer the ocean mass variation. Despite the GIA phenomenon now being tightly integrated into the science of global change (Church et al., 2013), few efforts have been devoted so far to the development of open-source codes for the solution of the SLE, although several post-glacial rebound simulators (e.g., TABOO; see Spada et al., 2004, 2011) and Love number calculators have been made available to the community (Spada, 2008; Melini et al., 2015; Bevis et al., 2016; Kachuck and Cathles, 2019). As far as we know, the only publicly available and open-source program in which the SLE is solved in its complete form is SELEN (SealEveL EquatioN solver). The SLE solver ISSM-SESAW v1.0 (Ice Sheet System Model – Solid Earth and Sea-level Adjustment Workbench; Adhikari et al., 2016), being oriented to short-term cryosphere and climate changes, only accounts for the elastic deformation of the Earth. The open-source SLE solver giapy (Kachuck, 2017; Martinec et al., 2018), available from https://github.com/skachuck/giapy (last access: 26 November 2019), can deal with complex ice models and viscoelastic rheology but does not take rotational effects into account.

SELEN was first presented to the GIA modeling community by Spada and Stocchi (2007), who numerically implemented the SLE theory reviewed in Spada and Stocchi (2006). SELEN was fully based on the classical formulation of Farrell and Clark (1976); hence, the fixed-shoreline approximation was assumed, and no account was given for rotational effects on sea-level variations. SELEN used the Love number calculator TABOO (see Spada et al., 2011) as a subroutine and was tied to the Generic Mapping Tools (GMT; see Wessel and Smith, 1998) for the construction of the present-day ocean function. In SELEN and in all its subsequent versions, the numerical integration of the SLE over the sphere takes advantage of the icosahedron-based pixelization proposed by Tegmark (1996). Similarly, all the versions are based upon the pseudo-spectral method of Mitrovica and Peltier (1991) and Mitrovica et al. (1994) for the solution of the SLE. Originally, SELEN came without a user guide, and it was disseminated via email by the authors. After SELEN was first published in 2007, a number of improvements were made in terms of computational efficiency, portability, and versatility but left the physical ingredients of the original code unaltered. This led to a new version of the program, named SELEN 2.9 and announced by Spada et al. (2012). Since 2015, the Computational Infrastructure for Geodynamics (CIG; http://www.geodynamics.org/, last access: 26 November 2019) has been hosting SELEN 2.9 in its short-term crustal dynamics section (see http://geodynamics.org/cig/software/selen, last access: 26 November 2019), from where it can be freely downloaded along with a theory booklet and a fully detailed user guide (Spada and Melini, 2015). Since the year 2012, with the aid of Florence Colleoni and thanks to the feedback of a number of colleagues and students, Giorgio Spada (GS) and Daniele Melini (DM) have implemented new modules aimed at solving the SLE in the presence of rotational effects and taking the migration of the shorelines into account. This has progressively led to several interim versions of the program (SELEN 3.x), which have been tested intensively and validated during the years but never officially released. We note that building upon SELEN, some colleagues have independently developed other versions of the code aimed at specific tasks, such as the study of the coupling between the SLE and ice dynamics (de Boer et al., 2017).

Taking advantage of the experience developed since SELEN was first designed, we are now publishing a new version of the code named SELEN⁴ (SELEN version 4.0). With respect to previous versions, SELEN⁴ has been improved in several aspects. (i) The underlying SLE theory has been fully revised and now it accounts both for horizontal migration of shorelines and for rotational effects, resulting in a more realistic description of the GIA processes. (ii) The package has been streamlined and reorganized into two independent modules: a solver, which obtains a numerical solution of the SLE in the spectral domain, and a post-processor, which computes a full suite of observable quantities through a spherical harmonic synthesis. This new structure facilitates code portability, reusability, and customization, enabling the adaptation of SELEN⁴ to new use cases. (iii) The SELEN⁴ modules have been completely rewritten using symbol names that closely match those of the variables introduced in this paper for the ease of code readability. Particular attention has been paid to the optimization of the SLE solver, resulting in a large extent of shared-memory parallelism, which allows for an efficient scaling to high resolutions on multi-core systems. (iv) SELEN⁴ has been decoupled from the GMT software package, which is no longer strictly required to run a GIA simulation, thus facilitating code portability on high-performance systems where GMT may not be available. SELEN⁴ still takes advantage of GMT (version 4) to produce various graphical outputs through plotting scripts included in the distribution package. (v) SELEN⁴ no longer calls the post-glacial rebound solver TABOO as an internal subroutine to compute the viscoelastic loading and tidal Love numbers, which are instead supplied by the user through a data file. In this way, any set of Love numbers can be used in SELEN⁴, possibly overcoming some of the intrinsic limitations of existing Love number calculators like TABOO. (vi) Recently, a prototype version of SELEN⁴ has been successfully validated in a community benchmark of independently developed GIA codes (Martinec et al., 2018). In the benchmark tests, simplified surface loads have been employed, and the rotational effects have been ignored. In a further benchmark effort, we are planning to adopt a realistic description of the surface load and the contribution from the Earth's rotation.

The paper is organized into three main sections. In Sect. 2, we present a condensed theory background for the SLE. In Sect. 3, we describe the outputs obtained by a standard, intermediate-resolution run of SELEN⁴. In Sect. 4, we draw our conclusions.

2 Theory

Here, we obtain the SLE from first principles, leaving a number of details of the theory to the supplementary material, hereafter referred to as SSM19. The focus is on the various forms the SLE takes, which are characterized by increasing complexity; the goal is to obtain a formulation suitable for numerical discretization, which is given and analyzed in SSM19 and implemented in SELEN⁴. In an attempt to simplify the presentation, and to obtain compact expressions for all the quantities involved, we are not exactly following the traditional notation adopted in the literature since the seminal paper by Farrell and Clark (1976), hereafter referred to as FC76. The number of definitions and variables involved in the construction of the SLE is remarkable and some derivations are cumbersome; to facilitate the readers – especially those who are approaching the GIA problem for the first time – we refer to the synopsis presented in SSM19 (see Table S8). The SELEN⁴ program is written in a plain way, adopting constants and variables names that follow the same notation employed in this theory section and in SSM19, with the aim of facilitating code readability.

2.1 Surface loads

We consider the system composed by the ice and the water in the oceans at a given time t. Its mass can be expressed as

$$\begin{array}{}\text{(1)}& M\left(t\right)=\underset{e}{\int }L\mathrm{d}A,\end{array}$$

where L(γ,t) is the surface load; γ stands for (θ,λ), where θ and λ are the geocentric colatitude and longitude, respectively; the integral is over the whole Earth's surface; and dA=a²sin θdθdλ is the area element, a being the average radius of the Earth. According to Eq. (1), L represents the mass per unit area distributed over the Earth's surface:

$$\begin{array}{}\text{(2)}& L(\mathit{\gamma },t)={\displaystyle \frac{\mathrm{d}M}{\mathrm{d}A}},\end{array}$$

with $M\left(t\right)=M^{\mathrm{i}}+M^{\mathrm{w}}$, where Mⁱ and M^w are the mass of the ice and the water at time t. By obtaining explicit expressions for Mⁱ and M^w, in SSM19, we show that L can be expressed as the sum of two contributions, which account (i) for the load exerted by the grounded ice and (ii) for the load of water on the oceans' floors, respectively:

$$\begin{array}{}\text{(3)}& L(\mathit{\gamma },t)={\mathit{\rho }}^{\mathrm{i}}IC+{\mathit{\rho }}^{\mathrm{w}}BO,\end{array}$$

where ρⁱ and ρ^w are the ice and ocean water densities, I is the ice thickness, B is sea level (i.e., the offset between the sea surface and the surface of the solid Earth), $C=\mathrm{1}-O$ is the continent function (CF), where the ocean function (OF) is

$$\begin{array}{}\text{(4)}& {\displaystyle }O(\mathit{\gamma },t)\equiv \left\{\begin{array}{ll}\mathrm{1},& \text{if}T+\frac{{\mathit{\rho }}^{\mathrm{i}}}{{\mathit{\rho }}^{\mathrm{w}}}I<\mathrm{0}\\ \mathrm{0},& \text{if}T+\frac{{\mathit{\rho }}^{\mathrm{i}}}{{\mathit{\rho }}^{\mathrm{w}}}I\ge \mathrm{0}.\end{array}\right.\end{array}$$

Following Kendall et al. (2005) (see their Eq. 2), we have defined bedrock topography T as the negative of the globally defined sea level:

$$\begin{array}{}\text{(5)}& T(\mathit{\gamma },t)=-B.\end{array}$$

We remark that due to the horizontal migration of the shorelines and to the transition between floating and grounded ice, O and C are, in general, time dependent.

In the following, we are concerned with time variations of the fields involved in the SLE, which shall be denoted by calligraphic capital letters. These variations are relative to a reference equilibrium state, established for $t\le t_{\mathrm{0}}\equiv \mathrm{0}$, during which all fields have a constant value. Conventionally, time t=t₀ denotes the time at which the surface load changes (see, e.g., Kendall et al., 2005). Accordingly, using Eq. (3) in Eq. (1), the mass variation $\mathcal{M}\left(t\right)=M-M_{\mathrm{0}}$ of the system composed by ice and water is

$$\begin{array}{}\text{(6)}& \mathcal{M}\left(t\right)={\mathit{\rho }}^{\mathrm{i}}\underset{e}{\int }\left(IC-I_{\mathrm{0}}{C}_{\mathrm{0}}\right)\mathrm{d}A+{\mathit{\rho }}^{\mathrm{w}}\underset{e}{\int }\left(BO-B_{\mathrm{0}}{O}_{\mathrm{0}}\right)\mathrm{d}A,\end{array}$$

where subscript 0 denotes reference quantities (i.e., constant values attained for t≤0), and the first term on the right-hand side represents the mass variation of the grounded ice, which we denote by μ(t). Since mass must be conserved, we have

$$\begin{array}{}\text{(7)}& \mathcal{M}\left(t\right)=\underset{e}{\int }\mathcal{L}\mathrm{d}A=\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(8)}& \mathcal{L}(\mathit{\gamma },t)\equiv L-L_{\mathrm{0}}.\end{array}$$

is the surface load variation. From Eq. (7), it follows immediately that an equivalent form of the mass conservation constraint is

$$\begin{array}{}\text{(9)}& <\mathcal{L}{>}^{\mathrm{e}}\left(t\right)=\mathrm{0},\end{array}$$

where $<\mathrm{\dots }{>}^{\mathrm{e}}$ indicates the average over the whole Earth's surface. Hereinafter, we shall consider only plausible surface loads, which we define here as those loads for which the mass is conserved. A discussion on how these surface loads can be realized is given by Bevis et al. (2016), who however have adopted ad hoc assumptions on how the mass lost from the ice sheets is redistributed. The SLE, however, does not rely upon such assumptions, since it defines a gravitationally self-consistent and non-uniform mass distribution over the oceans.

In SSM19, a suitable decomposition is found for ℒ, namely

$$\begin{array}{}\text{(10)}& \mathcal{L}(\mathit{\gamma },t)={\mathcal{L}}^{\mathrm{a}}+{\mathcal{L}}^{\mathrm{b}}+{\mathcal{L}}^{\mathrm{c}},\end{array}$$

where the first term,

$$\begin{array}{}\text{(11)}& {\mathcal{L}}^{\mathrm{a}}(\mathit{\gamma },t)={\mathit{\rho }}^{\mathrm{i}}C\mathcal{I},\end{array}$$

is associated with the ice thickness variation:

$$\begin{array}{}\text{(12)}& \mathcal{I}(\mathit{\gamma },t)=I-I_{\mathrm{0}};\end{array}$$

the second,

$$\begin{array}{}\text{(13)}& {\mathcal{L}}^{\mathrm{b}}(\mathit{\gamma },t)={\mathit{\rho }}^{\mathrm{w}}O\mathcal{S},\end{array}$$

is associated with sea-level change:

$$\begin{array}{}\text{(14)}& \mathcal{S}(\mathit{\gamma },t)=B-B_{\mathrm{0}};\end{array}$$

and the third,

$$\begin{array}{}\text{(15)}& {\mathcal{L}}^{\mathrm{c}}(\mathit{\gamma },t)={\mathit{\rho }}^{\mathrm{r}}Q\mathcal{O},\end{array}$$

is associated with OF variations:

$$\begin{array}{}\text{(16)}& \mathcal{O}(\mathit{\gamma },t)=O-O_{\mathrm{0}},\end{array}$$

where ρ^r is an arbitrary reference density and Q a time-invariant auxiliary variable. We note that, by the same definition of OF, its variation 𝒪(γ,t) can only take the values $(-\mathrm{1},\mathrm{0},+\mathrm{1})$.

2.2 The sea-level equation

The sea-level equation expresses the relationship between sea-level change, the sea surface variation, and the vertical displacement of the solid Earth. In this section, various forms of this relationship are described.

Above, sea level B has been qualitatively defined as the offset between the sea surface and the surface of the solid Earth. More specifically, denoting by r^ss(γ,t) and r^se(γ,t) the radii of the sea surface and of the Earth's solid surface in a geocentric reference frame, respectively, sea level can be equivalently expressed as

$$\begin{array}{}\text{(17)}& B(\mathit{\gamma },t)=r^{\mathrm{ss}}-r^{\mathrm{se}},\end{array}$$

and in the reference state,

$$\begin{array}{}\text{(18)}& B_{\mathrm{0}}\left(\mathit{\gamma }\right)=r_{\mathrm{0}}^{\mathrm{ss}}-r_{\mathrm{0}}^{\mathrm{se}}.\end{array}$$

So, introducing the sea surface variation,

$$\begin{array}{}\text{(19)}& \mathcal{N}(\mathit{\gamma },t)=r^{\mathrm{ss}}-r_{\mathrm{0}}^{\mathrm{ss}},\end{array}$$

and the vertical displacement of the solid surface of the Earth,

$$\begin{array}{}\text{(20)}& \mathcal{U}(\mathit{\gamma },t)=r^{\mathrm{se}}-r_{\mathrm{0}}^{\mathrm{se}},\end{array}$$

using Eq. (14), we obtain the SLE in its most basic form:

$$\begin{array}{}\text{(21)}& \mathcal{S}(\mathit{\gamma },t)=\mathcal{N}-\mathcal{U},\end{array}$$

expressing the variation of the height of the water column bounded by the sea surface and the ocean floor.

The sea surface variation is tightly associated with variations in the Earth's gravity field. Indeed, FC76 have shown that

$$\begin{array}{}\text{(22)}& \mathcal{N}(\mathit{\gamma },t)=\mathcal{G}+c,\end{array}$$

where the displacement of the geoid is obtained by the Bruns formula:

$$\begin{array}{}\text{(23)}& \mathcal{G}(\mathit{\gamma },t)={\displaystyle \frac{\mathrm{\Phi }-{\mathrm{\Phi }}_{\mathrm{0}}}{g}},\end{array}$$

(Heiskanen and Moritz, 1967) where Φ is gravity potential, which includes both the effects of surface loading and those from rotational variations (Martinec and Hagedoorn, 2014), g is the reference gravity acceleration at the Earth's surface, and c is a spatially invariant term introduced by FC76 to ensure mass conservation. The c term is known within the GIA community as the FC76 constant, and its physical origin is explained in detail by Tamisiea (2011). Note that as a consequence of Eq. (22), the change in the sea surface elevation does not coincide with the variation of the geoid elevation. Hence, using Eq. (22), the SLE (21) becomes

$$\begin{array}{}\text{(24)}& \mathcal{S}(\mathit{\gamma },t)=\mathcal{R}+c,\end{array}$$

where

$$\begin{array}{}\text{(25)}& \mathcal{R}(\mathit{\gamma },t)=\mathcal{G}-\mathcal{U}\end{array}$$

shall be referred to as sea-level response function.

We now assume that the responses to surface loading and to changes in the centrifugal potential can be combined linearly. Accordingly, the SLE (Eq. 24) can be further rearranged as

$$\begin{array}{}\text{(26)}& \mathcal{S}(\mathit{\gamma },t)={\mathcal{R}}^{\mathrm{sur}}+c+{\mathcal{R}}^{\mathrm{rot}},\end{array}$$

where

$$\begin{array}{}\text{(27)}& {\mathcal{R}}^{\mathrm{sur}}(\mathit{\gamma },t)={\mathcal{G}}^{\mathrm{sur}}-{\mathcal{U}}^{\mathrm{sur}}\end{array}$$

and

$$\begin{array}{}\text{(28)}& {\mathcal{R}}^{\mathrm{rot}}(\mathit{\gamma },t)={\mathcal{G}}^{\mathrm{rot}}-{\mathcal{U}}^{\mathrm{rot}}\end{array}$$

are the surface and the rotation sea-level response functions, whereas

$$\begin{array}{}\text{(29)}& \mathcal{G}(\mathit{\gamma },t)={\mathcal{G}}^{\mathrm{sur}}+{\mathcal{G}}^{\mathrm{rot}}\end{array}$$

and

$$\begin{array}{}\text{(30)}& \mathcal{U}(\mathit{\gamma },t)={\mathcal{U}}^{\mathrm{sur}}+{\mathcal{U}}^{\mathrm{rot}}\end{array}$$

are the geoid and the vertical displacement response functions, respectively.

By the constraint of mass conservation given by Eq. (7), the FC76 constant is easily determined. Using the expression found in SSM19 (see, in particular, Eq. S52 in the Supplement), the SLE (Eq. 26) becomes

$$\begin{array}{}\text{(31)}& \begin{array}{rl}\mathcal{S}(\mathit{\gamma },t)=& {\mathcal{S}}^{\mathrm{ave}}+\left({\mathcal{R}}^{\mathrm{sur}}-<{\mathcal{R}}^{\mathrm{sur}}{>}^{\mathrm{o}}\right)\\ & +\left({\mathcal{R}}^{\mathrm{rot}}\right.\left.-<{\mathcal{R}}^{\mathrm{rot}}{>}^{\mathrm{o}}\right),\end{array}\end{array}$$

where $<\mathrm{\dots }{>}^{\mathrm{o}}$ indicates the average over the (time-dependent) ocean surface defined by O=1, and

$$\begin{array}{}\text{(32)}& {\mathcal{S}}^{\mathrm{ave}}\left(t\right)={\mathcal{S}}^{\mathrm{equ}}+{\mathcal{S}}^{\mathrm{ofu}},\end{array}$$

where 𝒮^equ and 𝒮^ofu are two spatially invariant terms. The first, referred to as equivalent sea-level change, is

$$\begin{array}{}\text{(33)}& {\mathcal{S}}^{\mathrm{equ}}\left(t\right)\equiv -{\displaystyle \frac{\mathit{\mu }}{{\mathit{\rho }}^{\mathrm{w}}{A}^{\mathrm{o}}}},\end{array}$$

where μ is the mass variation of the grounded ice and A^o is the area of the oceans. The second,

$$\begin{array}{}\text{(34)}& {\mathcal{S}}^{\mathrm{ofu}}\left(t\right)\equiv {\displaystyle \frac{\mathrm{1}}{A^{\mathrm{o}}}}\underset{e}{\int }{T}_{\mathrm{0}}\mathcal{O}\mathrm{d}A,\end{array}$$

depends explicitly upon variations of the OF, either due to the horizontal migration of the shorelines or to transitions from grounded to floating ice (or vice versa). It is important to note that 𝒮^equ(t) and 𝒮^ofu(t) are both model dependent via A^o and 𝒪. Evaluating the ocean average of both sides of Eq. (31), and observing that $<<\mathcal{R}{>}^{\mathrm{o}}{>}^{\mathrm{o}}=<\mathcal{R}{>}^{\mathrm{o}}$ and $<{\mathcal{S}}^{\mathrm{ave}}{>}^{\mathrm{o}}={\mathcal{S}}^{\mathrm{ave}}$, it is easily verified that 𝒮^ave simply represents the ocean-averaged relative sea-level change:

$$\begin{array}{}\text{(35)}& {\mathcal{S}}^{\mathrm{ave}}\left(t\right)\equiv <\mathcal{S}{>}^{\mathrm{o}}.\end{array}$$

As a consequence of that, from Eq. (31), we see that the regional imprint of GIA on relative sea-level change is totally determined by the response functions ℛ^sur and ℛ^rot. We note that according to the analysis of Mitrovica and Milne (2002), the ocean-averaged terms in Eq. (31) largely incorporate the far-field “ocean syphoning” effect, which describes the flux of water away from the equatorial regions to fill the space left by the collapse of forebulges at the periphery of previously glaciated regions.

In the classical FC76 framework, a constant OF is assumed and the effects arising from the Earth's rotation are neglected. Hence, ℛ^rot=0 and 𝒪=0, with the latter condition implying 𝒮^ofu=0 because of Eq. (34). In this context, the SLE  simply reduces to

$$\begin{array}{}\text{(36)}& {\mathcal{S}}^{\text{FC76}}(\mathit{\gamma },t)={\mathcal{S}}^{\mathrm{eus}}+\left({\mathcal{R}}^{\mathrm{sur}}-<{\mathcal{R}}^{\mathrm{sur}}{>}^{\mathrm{o}}\right),\end{array}$$

where

$$\begin{array}{}\text{(37)}& {\mathcal{S}}^{\mathrm{eus}}\left(t\right)=-{\displaystyle \frac{\mathit{\mu }}{{\mathit{\rho }}^{\mathrm{w}}{A}^{\mathrm{op}}}}\end{array}$$

is often referred to as “eustatic sea-level change”, and A^op represents the present-day area of the surface of the oceans. Hence, for a rigid and non-gravitating Earth (for which ℛ^sur=0), 𝒮^eus would represent the (spatially invariant) relative sea-level change. Note that 𝒮^eus, which only depends on the history of the past grounded ice volume (Milne and Mitrovica, 2008), should not be confused with 𝒮^ave, given by Eq. (32), since the latter is dynamically dependent upon the Earth's response through A^o(t) and 𝒪(t) (e.g., Spada, 2017). We note that in the current terminology for sea level (Gregory et al., 2019), the term “barystatic” should be preferred to “eustatic”, a term originally attributed to Suess (1906) and that has since then received many qualifications (see FC76, p. 648). Hence, it would be more correct to write

$$\begin{array}{}\text{(38)}& {\mathcal{S}}^{\mathrm{bar}}\left(t\right)=-{\displaystyle \frac{\mathit{\mu }}{{\mathit{\rho }}^{\mathrm{w}}{A}^{\mathrm{op}}}}.\end{array}$$

Following, e.g., Milne and Mitrovica (1998), the surface response function ℛ^sur in Eq. (31) is obtained by a 3-D spatiotemporal convolution:

$$\begin{array}{}\text{(39)}& {\mathcal{R}}^{\mathrm{sur}}(\mathit{\gamma },t)={\mathrm{\Gamma }}^{\mathrm{s}}\otimes \mathcal{L},\end{array}$$

where Γ^s(γ,t) is the surface sea-level Green function. The details of the expansion of ℛ^sur in series of spherical harmonics are somewhat cumbersome and can be found in SSM19. Here, we only note that using Eq. (10) we have ${\mathcal{R}}^{\mathrm{sur}}={\mathcal{R}}^{\mathrm{a}}+{\mathcal{R}}^{\mathrm{b}}+{\mathcal{R}}^{\mathrm{c}}$, where the three terms are obtained by convolving Γ^s with ℒ^a, ℒ^b, and ℒ^c, respectively. Contrary to ℛ^sur, the harmonic coefficients of the rotational response ℛ^rot(γ,t) are directly obtained by a 1-D time convolution:

$$\begin{array}{}\text{(40)}& {\mathcal{R}}_{lm}^{\mathrm{rot}}\left(t\right)={\mathrm{{\rm Y}}}_l^{\mathrm{s}}\cdot {\mathrm{\Lambda }}_{lm},\end{array}$$

where ${\mathrm{{\rm Y}}}_l^{\mathrm{s}}\left(t\right)$ is the rotation sea-level Green function and Λ_lm(t) is the coefficients of degree l and order m of the expansion of the variation of centrifugal potential Λ(γ,t) associated with changes in the Earth's angular velocity (Milne and Mitrovica, 1998). In SSM19, it is shown that Λ(γ,t) is essentially a spherical harmonic function of degree l=2 and order $m=\pm \mathrm{1}$, with the other degree l=2 terms and the degree l=0 in the expansion of Λ(γ,t) being negligible (see also Martinec and Hagedoorn, 2014).

Thus, the SLE (31) can be rearranged as

$$\begin{array}{}\text{(41)}& \mathcal{S}(\mathit{\gamma },t)={\mathcal{S}}^{\mathrm{ave}}+{\mathcal{R}}^{\prime \mathrm{a}}+{\mathcal{R}}^{\prime \mathrm{b}}+{\mathcal{R}}^{\prime \mathrm{c}}+{\mathcal{R}}^{\prime \mathrm{rot}},\end{array}$$

where the primed response functions are

$$\begin{array}{}\text{(42)}& {\displaystyle }& {\displaystyle }{\mathcal{R}}^{\prime \mathrm{abc}}(\mathit{\gamma },t)={\mathcal{R}}^{\mathrm{abc}}-<{\mathcal{R}}^{\mathrm{abc}}{>}^{\mathrm{o}},\text{(43)}& {\displaystyle }& {\displaystyle }{\mathcal{R}}^{\prime \mathrm{rot}}(\mathit{\gamma },t)={\mathcal{R}}^{\mathrm{rot}}-<{\mathcal{R}}^{\mathrm{rot}}{>}^{\mathrm{o}}.\end{array}$$

We note that ℛ^′ b depends on O𝒮 through the surface load variation ℒ^b (see Eq. 13); in SSM19, it is shown that this holds for ℛ^rot as well. Following Mitrovica and Peltier (1991), in view of the numerical solution of the SLE, it is therefore convenient to transform Eq. (41) in such a way that O𝒮 becomes the unknown in lieu of 𝒮. This is accomplished by projecting Eq. (41) on the OF (i.e., multiplying both sides of the SLE by O), which provides the final form of the SLE:

$$\begin{array}{}\text{(44)}& \mathcal{Z}(\mathit{\gamma },t)={\mathcal{Z}}^{\mathrm{ave}}+{\mathcal{K}}^{\mathrm{a}}+{\mathcal{K}}^{\mathrm{b}}\left(\mathcal{Z}\right)+{\mathcal{K}}^{\mathrm{c}}+{\mathcal{K}}^{\mathrm{rot}}\left(\mathcal{Z}\right),\end{array}$$

where

$$\begin{array}{}\text{(45)}& \mathcal{Z}(\mathit{\gamma },t)=O\mathcal{S},\end{array}$$

and

$$\begin{array}{}\text{(46)}& {\displaystyle }& {\displaystyle }{\mathcal{Z}}^{\mathrm{ave}}(\mathit{\gamma },t)=O{\mathcal{S}}^{\mathrm{ave}}=O\left({\mathcal{S}}^{\mathrm{equ}}+{\mathcal{S}}^{\mathrm{ofu}}\right),\text{(47)}& {\displaystyle }& {\displaystyle }{\mathcal{K}}^{\mathrm{abc}}(\mathit{\gamma },t)=O{\mathcal{R}}^{\prime \mathrm{abc}},\text{(48)}& {\displaystyle }& {\displaystyle }{\mathcal{K}}^{\mathrm{rot}}(\mathit{\gamma },t)=O{\mathcal{R}}^{\prime \mathrm{rot}}.\end{array}$$

The dependence of 𝒦^b and 𝒦^rot upon 𝒵 in Eq. (44) manifests the implicit nature of the SLE, which is a 3-D non-linear integral equation, similar, in some respect, to an inhomogeneous and non-linear Fredholm equation of the second kind (e.g., Jerri, 1999; Spada, 2017). For the spectral discretization of Eq. (44) and for a detailed illustration of the scheme adopted to solve the SLE, the reader is referred to Sects. S7 and S8.7 of SSM19. Here, it is only worth mentioning that the SLE is solved by a pseudo-spectral iterative approach (Mitrovica and Peltier, 1991; Mitrovica and Milne, 2003) complemented by a spatial discretization on an icosahedron-based spherical geodesic grid (Tegmark, 1996). Two nested iterations are adopted, where the external one refines the topography progressively using the solution of the SLE obtained in the internal one.

3 A test run with SELEN⁴

In the following, we illustrate some of the outputs of a standard SELEN⁴ run in which the resolution of the Tegmark grid is set to R=44 (see Sect. S8.6 in the Supplement), while the maximum harmonic degree of the spectral decomposition is $l_{max}=\mathrm{128}$. Note that condition $\mathit{\text{P}}\ge l_{max}^{\mathrm{2}}/\mathrm{3}$, where $P=\mathrm{40}R(R-\mathrm{1})+\mathrm{12}$ is the number of pixels in the grid for a given R value, must be necessarily met to preserve the properties of the spherical harmonics on the grid as discussed in Sect. S8.6 and in Tegmark (1996). In the test run, we solve the SLE by three external and three internal iterations ($n_{\mathrm{ext}}=n_{\mathrm{int}}=\mathrm{3}$; see Sect. S8.7), as customarily adopted in GIA studies as in Kendall et al. (2005). Henceforth, the notation R44/L128/I3 shall be employed to denote these fundamental SELEN⁴ settings. Of course, with increasing values of parameters (R, l_max, n_ext, n_int), more accurate results are expected, which however might come with a substantial increase in the computational burden. In the SELEN⁴ test run, we account for both loading and rotational effects.

We assume that the user has installed and executed the program following the indications given in the user guide of SELEN⁴. Most of the program outputs discussed in this section have been obtained using the same configuration file that comes with the SELEN⁴ package. However, some results shall be based on different settings, in order to appreciate the sensitivity of the outputs on key configuration parameters. We first describe the GIA model adopted in the test run, which consists of three elements, i.e., an ice melting history, a description of the present-day global relief, and a 1-D rheological model of the Earth's mantle. Then, browsing the output folders of SELEN⁴, we illustrate and discuss two distinct output sets, pertaining to the past and present effects of GIA on sea-level change and on geodetic variations, respectively. The features of the test run are summarized in Table 1. For reference, on a 12-core mid-2012 Mac Pro, the execution time of this test run of SELEN⁴ is 1 h and 15 min.

Spada et al. (2011)Spada et al. (2011)Spada et al. (2011)

Table 1 Configuration details of the SELEN⁴ test run whose results are considered in Sect. 3, with notes and links to text and figures.

3.1 GIA model

In principle, there are no restrictions on the spatial and temporal features of the ice melting history that can be employed in SELEN⁴,  provided that the model is properly discretized according to the scheme outlined in Sect. S8.7. Similarly, any linear rheological profile is a priori acceptable for the mantle and for the lithosphere, as long as the Love numbers can be cast in a normal-mode multi-exponential form (Peltier, 1974). Due to its central role in the context of contemporary GIA studies, in our test run, we have implemented an ad hoc realization of the GIA model known as ICE-6G_C (VM5a), originally introduced by Peltier et al. (2015).

3.1.1 Ice melting history

Ice thickness data for model ICE-6G_C have been downloaded from the home page of William R. Peltier in August 2016. The data span the last 26 000 years and are provided on a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ global latitude–longitude grid (the number of grid points is thus 64 800). In each of the grid cells, the time history of ice thickness is assumed to evolve in a piecewise linear manner, with variable increments of 1.0 or 0.5 kyr. Thus, to fit the SELEN⁴ default input format, we have first mapped the original thickness data on a spherical equal-area Tegmark grid described in Sect. S8.6 of SSM19. In doing that, we have chosen a resolution parameter R=44, so that the grid consists of P=75 692 pixels (or cells), each with a radius of ∼46 km. The cells' number is thus a bit larger than the number of cells in the original grid. In addition, we have transformed the original time history in a piecewise constant form with a uniform spacing of 0.5 kyr, assuming no glaciation phase prior to deglaciation. Because of the adaptations we have made, the ice model obtained is not an exact replica of ICE-6G_C but a particular realization of it. Hence, to avoid any ambiguity, in the following, it shall be referred to as I6G-T05-R44. Assuming an ice density ρ^w=931 kg m⁻³ and that the area of the oceans is fixed to the present value, I6G-T05-R44 holds 206.5 m of equivalent sea level at 26 ka and 75.1 m at present, corresponding to a total eustatic sea-level rise of 131.4 m since the inception of melting (26 ka). The ice thickness of I6G-T05-R44 for a few time frames is shown in Fig. 1.

Figure 1 Ice distribution according to model I6G-T05-R44, at six different epochs since 26 ka. The maps are obtained by direct triangulation of the pixelized ice thickness data using the GMT program pscontour. This and the following figures are drawn using GMT scripts adapted from those which are available in the output folders of SELEN⁴ after execution.

3.1.2 Present-day topography

In order to reconstruct the whole history of the Earth's topography and of relative sea level (RSL) since the inception of deglaciation, it is necessary to impose the present relief as a final condition for the sea-level equation (see Peltier, 1994). In this test run, we have utilized the “bedrock version” of the global ETOPO1 dataset (Amante and Eakins, 2009; Eakins and Sharman, 2012) as the final condition, whereas the final condition for the ice thickness is given by the last time frame of I6G-T05-R44. ETOPO1 is distributed on a longitude–latitude grid with a resolution of 1 arcmin ($\mathrm{1}/\mathrm{60}{}^{\circ }$), so that an interpolation of topography on the Tegmark grid has been performed. Note that in regions where small-scale topographic features are present, such as ocean trenches, interpolating a high-resolution relief model on a coarse Tegmark grid may yield inaccurate results. In this case, a different approach should be adopted for the realization of topography, for instance, based on a binning algorithm. Other choices of the final topography are of course possible. For example, in order to ensure the maximum accuracy in Antarctica, the original model, ICE-6G_C (VM5a), employs the Bedmap2 dataset of Fretwell et al. (2013) south of 60^∘ S latitude. In SELEN⁴, there are no restrictions on the choice of the final relief, which is left to the user. In Fig. 2, we show the realization of the modern bathymetry that we have obtained by interpolating ETOPO1 on the same Tegmark grid that we have used for the ice sheets in this test run, which shall be referred to as model ETO-R44 in the following. With SELEN⁴, other versions of this  elevation model are made available, characterized by different spatial resolutions.

Figure 2 Present-day relief according to model ETO-R44 used in the test run, obtained from model ETOPO1 by bilinear interpolation on the pixels of the Tegmark grid with resolution R=44, using the GMT program grdtrack.

3.1.3 Rheological profile

The 11-layer Maxwell rheological profile of the 1-D Earth employed in the test run is shown in Table 2. For each layer, values of density and of rigidity are obtained by volume averaging the PREM (Preliminary Reference Earth Model) of Dziewonski and Anderson (1981), while the viscosity profile is reproduced using the data available in the Supplement in Peltier et al. (2015). The 90 km thick lithosphere is elastic and the core is fluid, homogeneous, and inviscid. Note that the original VM5a profile includes elastic compressibility and reproduces the finely layered PREM structure (Peltier et al., 2015); since we are employing an incompressible rheology, in the following, we shall refer to the model in Table 2 as VM5i. As it includes N_v=9 distinct Maxwell layers in the mantle, for any given harmonic degree l, the loading Love numbers (LLNs) and tidal Love numbers (TLNs) for the VM5i model are described by a spectrum of 4N_v=36 viscoelastic normal modes (see, e.g., Spada et al., 2011). Since we are not aware of published, community-agreed sets of Love numbers for a multi-layered compressible viscoelastic model, in the test run we rest on an incompressible profile, for which agreed results have been obtained by Spada et al. (2011). We remark, however, that SELEN⁴ can work with compressible or transient rheologies as well, provided that LLNs and TLNs in normal-mode form (Wu and Peltier, 1982) are accessible to the user. Figure 3a and b show the elastic and fluid values of the LLNs in the range of harmonic degrees $\mathrm{1}\le l\le \mathrm{1024}$ for the VM5i model. The Love numbers are given in a geocentric reference frame with the origin in the Earth’s center of mass (CM). Numerical values of a few relevant LLNs and TLNs are listed in Table 3 and in its caption. They have been computed by the Love number calculator TABOO (see Spada et al., 2011) in a multi-precision environment (Smith, 1991; Spada, 2008).

Figure 3 Elastic (a) and fluid (b) LLNs as a function of harmonic degree l for the 11-layer rheological VM5i model employed in the test run (see Table 2). It is apparent that for this model asymptotic values are reached, in both cases, for l exceeding values of several hundred. Note that in panel (b), where the fluid LLN for vertical displacement is normalized by (2l+1), the relationship $h_l^{\mathrm{Lf}}\approx (\mathrm{2}l+\mathrm{1})k_l^{\mathrm{Lf}}$ is apparent.

Table 2 Density, rigidity, and viscosity profiles adopted in the rheological VM5i model, where abbreviations LT, UM, TZ, and LM stand for lithosphere, upper mantle, transition zone, and lower mantle, respectively. The radii r₋ and r₊ indicate the lower and upper radii of each layer. Some spectral properties of the VM5i model are given in Fig. 3 and Table 3.

3.2 GIA in the past

This section is devoted to the description of some outputs of the SELEN⁴ test run, concerning the effects of GIA during the whole period after the LGM. These include (i) the predictions of the history of RSL at specific sites, (ii) the time evolution of paleotopography in some regions of interest, and (iii) the excursions of the Earth's pole of rotation forced by GIA.

3.2.1 Relative sea-level curves

Figure 4 shows data (with error bars) and SELEN⁴ predictions for a small subset of the 392 sites contained in the RSL database of Tushingham and Peltier (1993), hereafter referred to as TP93. In view of its historical importance in the development of GIA studies (e.g., Tushingham and Peltier, 1991, 1992; Melini and Spada, 2019), the TP93 database is available with the SELEN⁴ package; however, there are no restrictions on the use of other datasets, or simply individual RSL records, if available to the user. Note that ages given in the TP93 database included with SELEN⁴ are expressed in radiocarbon years and have not been calibrated to calendar years. In Fig. 4, three different configurations, i.e., combinations of parameters R, l_max, n_ext, and n_int have been adopted, always assuming n_ext=n_int. Results obtained for the three configurations are shown by different colors. The black curves have been obtained using the GIA model described in Sect. 3.1, characterized by the settings R44/L128/I3 (R=44, l_max=128, $n_{\mathrm{ext}}=n_{\mathrm{int}}=\mathrm{3}$). Blue curves have been obtained by configuring SELEN⁴ with the combination of parameters R100/L512/I5, i.e., increasing the spatial resolution and the number of internal and external iterations in a significant way (a truncation degree $l_{max}=\mathrm{256}$ is often employed in GIA modeling; see, e.g., Kendall et al., 2005). With this configuration, the pixel radius is reduced to ∼20 km (see Sect. S8.6 and Table S6). Of course, the execution time of SELEN⁴ increases significantly with respect to the test run, requiring 2.5 d (∼60 h) on a 56-core Intel Xeon E7 Broadwell system. It is apparent that this high-resolution case is providing results substantially matching those of the standard run with R44/L128/I3. Minor differences can be noted in the early stages of deglaciation, which however do not exceed the typical uncertainty on the observed RSL values. These differences are likely to be caused by the significant changes that the topography undergoes in this early phase in the polar regions, which are better captured by increasing the model resolution. Finally, red curves have been obtained for a low-resolution run with R30/L64/I2, whose execution time is 15 min on a 12-core mid-2012 Mac Pro. The curves clearly indicate that computationally inexpensive runs can provide reliable results in the far field of the previously glaciated areas (e.g., at sites 639 and 525), but in the near field (e.g., site 283) they can diverge significantly from both high- and intermediate-resolution results.

Figure 4 RSL data (with error bars) at eight of the 392 sites of the TP93 RSL database. Black curves show the results of the standard test run with configuration R44/L128/I3, the blue ones are for the high-resolution run with R100/L512/I5, while those in red are for a low-resolution configuration with R30/L64/I2. The map shows the location of the RSL sites.

From a visual inspection of Fig. 4, it is apparent that at some sites the best GIA predictions fit very well the observations, like at sites 101 and 283. For others, the trend of the RSL data is captured satisfactorily (see sites 155, 209, and 328), while in some others the fit is quite poor (sites 639, 525, and 570). The identification of the possible sources of the evidenced misfits, which should be measured using rigorous statistical methods after a proper calibration of the ages, is not the purpose of this work. We only note that they do not necessarily stem from limitations of the GIA model adopted, since it is well known that at a specific site tectonic deformations can have an important role (see, e.g., Antonioli et al., 2009, for a significant example), and these are not taken into account when solving the SLE. Similarly, in our formulation of the SLE, we are neglecting the possible effects from the loads exerted by sediments (Dalca et al., 2013). Since the SELEN⁴ program is open source, the users can modify the code to account for non-glacial loads and change the configuration to determine more suitable combinations of the basic ingredients of GIA modeling, i.e., the history of deglaciation and the rheological layering of the mantle, in order to improve the fit between model predictions and any preferred RSL dataset.

3.2.2 Paleotopography

SELEN⁴ allows for a gravitationally and topographically self-consistent description of the evolution of sea level along the route highlighted by Peltier (1994) and Lambeck (2004). This implies that SELEN⁴ can iteratively reconstruct the time evolution of the coastlines and of the OF, in a fashion that is consistent with the gravitational, rotational, and deformational effects induced by deglaciation (for a full account of the theory, the reader is referred to SSM19). These features make the SLE an integral 3-D non-linear equation (Spada, 2017). The importance of the evolution of paleotopography for the development of human culture since the LGM has been pointed out in a number of works (see, e.g., Cavalli-Sforza et al., 1993; Peltier, 1994; Lambeck, 2004; Dobson, 2014, and references therein). Recently, in the context of GIA modeling, Spada and Galassi (2017) have faced the problem of the dynamic evolution of aquaterra, i.e., the land that has been inundated and exposed during the last glacial cycle (Dobson, 1999, 2014), using the same approach adopted in this work.

A standard SELEN⁴ output is shown in Fig. 5, where the Earth's relief at the LGM is reconstructed in the post-processing phase of SELEN⁴, according to the solution of the SLE. Note that the figure only shows the bedrock relief at the LGM, which is not what Peltier (1994) has called true paleotopography (PT), which also includes the contribution of ice elevation. The user can easily obtain maps of the full PT merging maps like that shown in Fig. 5 with those of Fig. 1. Light blue areas across the polar regions covered by thick ice at the LGM correspond to places where the ice was grounded below sea level in that epoch. In SELEN⁴, it is also possible to visualize the full time evolution of the OF in order to better appreciate, for example, the spatiotemporal distribution of the ice shelves, which cannot be seen in this paleotopography map. At the global scale, the major land masses that were exposed at the LGM are clearly visible, as evidenced by the low-elevation green areas; these include Beringia, Sundaland, Patagonia, Sahul, Doggerland, and the East Siberian Arctic Shelf. Using GIA models, the time evolution of these land masses has been investigated in an a number of papers, also in view of the important implications on the dispersal of modern humans. For specific case studies, the reader is referred to the works of Peltier and Drummond (2002), Lambeck (2004), Klemann et al. (2015), Spada and Galassi (2017), and references therein.

Figure 5 Sample SELEN⁴ output of the global relief at the LGM (21 ka), according to the test run with R44/L128/I3, where topography ETO-R44 of Fig. 2 has been used as a final condition for the SLE.

By increasing the spatial resolution, SELEN⁴ can also be employed to resolve the past sea-level variations on a regional scale. As a specific case study, we consider the Mediterranean Sea. The history of RSL in the Mediterranean Sea has been the subject of various investigations, stimulated by the amount of high-quality geological, geomorphological, and archeological indicators in the region (see, e.g., Lambeck and Purcell, 2005; Antonioli et al., 2009; Evelpidou et al., 2012; Vacchi et al., 2016; Roy and Peltier, 2018, and references therein). Since on the global scale of Fig. 5 the details of the paleotopography in this area are difficult to visualize, we have used the outputs of the high-resolution run with settings R100/L512/I5, already exploited in Sect. 3.2.1. The results are shown in the map of Fig. 6, where paleotopography is shown at 26 ka. The vastly exposed continental shelf of Tunisia (Mauz et al., 2015) and the northern Adriatic Sea (Lambeck and Purcell, 2005) are now clearly visible, along with other smaller-scale regions where the topography has seen significant changes during the last deglaciation (Lambeck, 2004; Purcell and Lambeck, 2007).

Figure 6 Paleotopography of the Mediterranean Sea and the Black Sea at 26 ka, obtained by a high-resolution SELEN⁴ run with configuration R100/L512/I5, in order to highlight the exposed lands in detail.

3.2.3 Polar motion

With SELEN⁴, three configurations are possible in which rotational effects on GIA are dealt with in different manners. First, these effects can be simply ignored, as it is done in the classical FC76 GIA theory. However, when rotational effects are taken into consideration, this can be done in two different ways, i.e., either following the traditional rotation theory (Milne and Mitrovica, 1998; Spada et al., 2011) or a revised rotation theory proposed by Mitrovica et al. (2005) and Mitrovica and Wahr (2011). The reader is referred to the literature for a detailed presentation of the two theories and to Sect. S5.2 for a brief account. We also refer to Martinec and Hagedoorn (2014) for a general formulation dealing with rotational effects in GIA modeling. Here, it is useful to mention that in the traditional treatment the long-term response of the Earth is evaluated assuming that the lithosphere is characterized by a finite elastic strength, while in the revised theory the equilibrium rotational shape is, more realistically, only based on the viscous properties of the planet. Furthermore, the long-term extra flattening due to mantle dynamics is properly accounted for. We remark that in both cases the fast Chandler wobble component of polar motion is filtered out from the Liouville equations, because the timescales of GIA largely exceed the Chandler wobble period (∼14 months; see, e.g., Lambeck, 1980). The implications for the GIA response of the new theory are quite significant, as illustrated in detail by Mitrovica et al. (2005) and Mitrovica and Wahr (2011) and shown in Fig. 7.

Figure 7 Cartesian components of polar motion along the axes x and y (a) and their time derivatives (b) as a function of time since the beginning of deglaciation for the test run with configuration R44/L128/I3. The x axis points along the prime meridian and the y points to 90^∘ E. Dashed and solid curves show results for the traditional and for the revised rotational theories, respectively. The steep change in the y components at ∼14 ka is forced by the inertia variations due to the occurrence of the meltwater pulse 1A (MWP-1A).

Solid curves in Fig. 7 show the evolution of the polar motion components (m_x,m_y) and their rates of change $({\dot{m}}_x,{\dot{m}}_y)$ in response to GIA, obtained by solving the Liouville equations in the test run, in which the revised rotation theory is employed. Dashed curves show results obtained with the traditional theory. As the x and y components of polar motion are conventionally measured along the Greenwich meridian and 90^∘ E, respectively, Fig. 7a indicates that since the inception of deglaciation, the displacement of the pole has been roughly in the direction of the Hudson Bay, consistent with the seminal results of Sabadini and Peltier (1981). The two theories predict similar evolutions of the pole of rotation, which according to Fig. 7a has been displaced by ∼18 km on the Earth's surface by the glacial readjustment process since 26 ka. We note that at the time of the rapid melting episode known as meltwater pulse 1A (MWP-1A, between 14.3 and 12.8 ka; see, e.g., Blanchon, 2011), a sudden variation in the m_y component of polar motion had occurred but no changes could be observed on m_x. When the rates of polar motion are considered in Fig. 7b, differences between the predictions of the two rotations theories are more apparent. In particular, the present-day (0 ka) rates of polar motion are found to be ∼1 and ∼3^∘ Ma⁻¹ for the revised and the traditional theories, respectively, which fits the predictions of Mitrovica and Wahr (2011) and confirms that the traditional theory largely overestimates the effects of GIA on polar motion. We note also that MWP-1A has caused a remarkable acceleration of polar motion, with a rate as large of ∼5^∘ Ma⁻¹ near the end of the pulse; this value exceeds the rate of polar motion observed since the year 1900 by a factor of ∼5 (see, e.g., Dickman, 1977).

3.3 GIA at present

In this section, we describe further outputs of the SELEN⁴ test run, focusing on the effects of GIA at the present time. In particular, we consider (i) the global pattern of the so-called GIA fingerprints, (ii) predictions of the rate of sea-level change at tide gauges, and (iii) the time variations of the Stokes coefficients of the Earth's gravity field induced by GIA.

3.3.1 GIA fingerprints

Figure 8 shows another standard output of SELEN⁴, i.e., the present-day rates of variation of four fundamental quantities associated with GIA, obtained for the test run. These are relative sea-level change ($\dot{\mathcal{S}}$; Fig. 8a), vertical displacement of the crust ($\dot{\mathcal{U}}$; Fig. 8b), absolute sea level ($\dot{\mathcal{N}}$; Fig. 8c), and geoid height ($\dot{\mathcal{G}}$; Fig. 8d). Generalizing the definition of sea-level fingerprint introduced by Plag and Jüettner (2001) in the context of contemporary ice mass changes, these have been referred to as GIA fingerprints by Spada and Melini (2019). The spatial variability of the GIA fingerprints reflect the effects of deformation, gravitational attraction, and rotation within the system composed by the solid Earth, the oceans, and the ice sheets (Clark et al., 1978; Mitrovica and Milne, 2002). In view of their importance on the interpretation of ground-based (King et al., 2010) or satellite geodetic observations (Peltier, 2004) and of tide gauge secular trends (e.g., Spada and Galassi, 2012; Wöppelmann and Marcos, 2016), their properties have been the subject of various investigations during last decade (see, e.g., Mitrovica et al., 2011; Tamisiea, 2011; Spada and Galassi, 2015; Spada, 2017; Husson et al., 2018; Melini and Spada, 2019; Spada and Melini, 2019).

It should be remarked that the four fingerprints shown in Fig. 8 are not independent of one another. In particular, the SLE gives $\dot{\mathcal{S}}=\dot{\mathcal{N}}-\dot{\mathcal{U}}$ according to Eq. (21). Furthermore, $\dot{\mathcal{N}}=\dot{\mathcal{G}}+\dot{c}$, where c is the FC76 constant (see Sect. S2.4 in SSM19). The two relationships above hold, regardless of the particular combination of rheology and ice model employed, and the preferred rotation theory adopted. However, the patterns of the GIA fingerprints and the numerical value of $\dot{c}$ are model dependent. Other interesting results hold for the spatial averages of the fingerprints, which reflect some physical aspects of GIA (Spada, 2017; Spada and Melini, 2019) and are useful to correct geodetic observations from the effects of deglaciation (e.g., Spada and Galassi, 2015). In Table 4, we summarize the numerical values of whole Earth surface averages (denoted by the symbol $<\mathrm{\cdots }{>}^{\mathrm{e}}$) and ocean averages ($<\mathrm{\cdots }{>}^{\mathrm{o}}$) of the GIA fingerprints in the test run. In addition, we have also executed SELEN⁴ adopting the traditional rotation theory and neglecting rotational effects, and the corresponding averages are shown in Table 4 as well. We note that by virtue of mass conservation, $<\dot{\mathcal{G}}{>}^{\mathrm{e}}=<\dot{\mathcal{U}}{>}^{\mathrm{e}}=\mathrm{0}$ (see Sects. S4.3 and S6.2), regardless of the rotation theory adopted, and as a consequence $<\dot{\mathcal{S}}{>}^{\mathrm{e}}=<\dot{\mathcal{N}}{>}^{\mathrm{e}}=\dot{c}$. We also note that the numerical value of $<\dot{\mathcal{N}}{>}^{\mathrm{o}}$, commonly employed to correct the altimetry observations of absolute sea-level change for the effects of GIA, is in fair agreement with predictions from state-of-the-art GIA models (e.g., Church et al., 2013; Spada and Galassi, 2015; Spada, 2017). Notably, $<\dot{\mathcal{N}}{>}^{\mathrm{o}}$ is not very sensitive to the choice of the rotation theory. Since model I6G-T05-R44 assumes that melting of the major ice sheets ceased ∼4000 years ago, the small value of $<\dot{\mathcal{S}}{>}^{\mathrm{o}}$ only reflects ongoing changes in the area of the oceans due to GIA. In the FC76 fixed-shoreline approximation, $<\dot{\mathcal{S}}{>}^{\mathrm{o}}$ would be identically zero by virtue of the mass conservation principle (see, e.g., Spada, 2017; Spada and Melini, 2019).

Figure 8 Present-day GIA fingerprints obtained for the test run with R44/L128/I3. Note that the color table is saturated in a narrow interval. The effects of Earth rotation, evaluated according to the revised rotation theory, can be well discerned for $\dot{\mathcal{N}}$ and $\dot{\mathcal{G}}$, with the characteristic high-amplitude lobes with a harmonic degree l=2 and order m=1 symmetry (Spada and Galassi, 2015). Spatial averages of these maps are given in Table 4. Min/max values for $\dot{\mathcal{S}}$, $\dot{\mathcal{U}}$, $\dot{\mathcal{N}}$, and $\dot{\mathcal{G}}$, after rounding the SELEN⁴ outputs to one decimal figure, are $-\mathrm{13.3}/\mathrm{4.9}$, $-\mathrm{5.2}/\mathrm{14.2}$, $-\mathrm{0.8}/\mathrm{1.2}$, and $-\mathrm{0.6}/\mathrm{1.4}$ in units of mm yr⁻¹, respectively.

Table 3 Numerical values of the LLNs for the VM5i rheological model (see Table 2) for some harmonic degrees l. Here x^e and x_e are abbreviations for x×10^e and $x\times {\mathrm{10}}^{-e}$, respectively. Note that, for this model, the elastic TLNs of degree l=2 are $(k_{\mathrm{2}}^{\mathrm{Te}},h_{\mathrm{2}}^{\mathrm{Te}},l_{\mathrm{2}}^{\mathrm{Te}})=({\mathrm{0.289}}^{\mathrm{0}},{\mathrm{0.524}}^{\mathrm{0}},{\mathrm{0.108}}^{\mathrm{0}})$, while the fluid values are $(k_{\mathrm{2}}^{\mathrm{Tf}},h_{\mathrm{2}}^{\mathrm{Tf}},l_{\mathrm{2}}^{\mathrm{Tf}})$ = $({\mathrm{0.931}}^{\mathrm{0}},{\mathrm{0.191}}^{\mathrm{1}},{\mathrm{0.514}}^{\mathrm{0}})$.

Table 4 Whole Earth and ocean averages of the GIA fingerprints according to the test run based upon the new rotation theory (column a), and the traditional theory (column b). In column (c), we also consider the case when no rotational effects are taken into account. In all the computations, we have adopted the combination R44/L128/I3. Note that $<\dot{\mathcal{N}}{>}^{\mathrm{e}}=\dot{c}$, where c is the FC76 constant. In this table, the SELEN⁴ outputs are rounded to two significant figures.

3.3.2 GIA at tide gauges

Estimating global mean sea-level rise in response to climate change requires correcting tide gauge records for the effects of GIA. Since the late 1980s, with the awareness of global warming and the availability of numerical solutions of the SLE (Peltier and Tushingham, 1989), GIA corrections have been routinely applied to the observed trends of sea level (for a review, see Spada and Galassi, 2012; Spada et al., 2015; Wöppelmann and Marcos, 2016). However, since GIA models are progressively improved to provide a better description of reality, corrections at tide gauges are not given once and for all (Kendall et al., 2006; Tamisiea, 2011; Melini and Spada, 2019). Furthermore, new constraints from past sea level or modern geodetic observations have permitted gradual refinements (either by formal inverse methods or simply by trial and error) of the two basic ingredients of GIA modeling, i.e., the Earth rheological profile and the history of deglaciation since the LGM. Uncertainties in modeling are significant (Melini and Spada, 2019), which constitutes an additional motivation to improve the approach to GIA.

In Table 5, we show SELEN⁴ predictions for $\dot{\mathcal{S}}$ at a few tide gauges for the test run and other possible configurations as well. Here, we only show results for the 23 sites that have been considered by Douglas (1997) in his redetermination of global sea-level rise, which obey specific criteria that make them suitable to represent the trend of secular sea-level rise. The sites chosen by Douglas (1997) are located in the periphery of the regions covered by thick ice sheets at the LGM, since GIA predictions at sites formerly beneath the ice sheets are expected to be more affected by uncertainties in GIA modeling. This has been recently confirmed by Melini and Spada (2019). The post-processing phase of SELEN⁴ can be configured to handle any properly formatted input dataset with coordinates of geodetic points of interest, where all the variables considered in Fig. 8 can be evaluated. This can be useful, for instance, to estimate the effects of GIA on vertical movements at specific GPS points (Serpelloni et al., 2013). Modules for the computation of horizontal displacements shall be included in future releases. Comparing column (d) with (b) and (c) (all these runs are characterized by the intermediate-resolution configuration R44/L128/I3), we note that rotational effects are important at tide gauges; however, from columns (b) and (c), we also note that differences between the revised and the traditional rotation theories generally do not exceed the 0.1 mm yr⁻¹ level at the tide gauges considered here. Comparing outputs in column (b) with the high-resolution run in column (a), we note maximum differences of 0.03 mm yr⁻¹, which further confirms the reliability of the test run.

Table 5 Present-day rates of sea-level change at the 23 Douglas (1997) tide gauges for the test run of SELEN⁴ (column b) and for other configurations. Results based upon the original implementation of model ICE-6G_C (VM5a) are reproduced in column (e). The average rate is also shown in the bottom line. The SELEN⁴ outputs have been rounded to two significant figures.

Figure 9 Present-day rate of change of the low harmonic degree (l≤6) Stokes coefficients obtained for the test run of SELEN⁴ (a) and their full spectrum extended to harmonic degree $l_{max}=\mathrm{128}$ (b). The solid line in panel (b) represents the power law that best fits the spectrum, in the least squares sense, for l≥10.

Comparing the high-resolution results in Table 5 (R100/L512/I5, column a) with those obtained using the ICE-6G_C (VM5a) model in the original implementation of W. R. Peltier (see http://www.atmosp.physics.utoronto.ca/~peltier/data.php, last access: 6 June 2019), reported in column (e), we note a fair agreement between the two model predictions. The two sets of predictions have the same sign with the only exception of Lagos. The differences between the values obtained with R100/L512/I5 and the original implementation, ICE-6G_C (VM5a), are generally close to 0.2 mm yr⁻¹ or slightly larger. The values of $\dot{\mathcal{S}}$ averaged over the tide gauges differ, but they are both <0.1 mm yr⁻¹. More significant differences in the $\dot{\mathcal{S}}$ values are however apparent when we compare model predictions for sites located in the polar regions beneath the former ice sheets, which are not considered in Table 5. At these locations, the expected $\dot{\mathcal{S}}$ values are of the order of several mm yr⁻¹, due to the large isostatic disequilibrium still associated with ice unloading. For example, at the tide gauge of Stockholm, we obtain a rate of −4.44 mm yr⁻¹ for the high-resolution run R100/L512/I5, while in the original ICE-6G_C (VM5a) implementation, the rate is −3.75 mm yr⁻¹. To better evaluate the meaning of these differences, it is worth noting that they largely exceed the typical 2σ uncertainty in the observed rates, which according to Spada and Galassi (2012) is never larger than 0.3 mm yr⁻¹ (see their Table 2). We have also ascertained that misfits of the order of 1 mm yr⁻¹ are not uncommon in other high-latitude sites of both hemispheres. Tracing the origin of the discrepancies in the two sets of GIA predictions (and therefore in the whole set of the GIA geodetic fingerprints considered in Fig. 8) is not easy at this stage and would demand a detailed model intercomparison study like those performed in the GIA community by Spada et al. (2011) and Martinec et al. (2018). We can, however, guess that the misfit between the two sets of GIA predictions stems from the different discretizations of the ice time histories, from the effects of mantle compressibility, and possibly from the different rotation theories adopted.

3.3.3 Stokes coefficients of the gravity field

In Fig. 9, we study the present-day rates of change of the variations of Stokes coefficients (${\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }c}}}_{lm},{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }s}}}_{lm}$) induced by GIA, computed in the test run with R44/L128/I3. These quantities represent the coefficients of the expansion of the gravity potential variation $\mathrm{\Phi }(\mathit{\gamma },t)-\mathrm{\Phi }(\mathit{\gamma },t_{\mathrm{0}})$ in series of spherical harmonics; hence, they contain information upon the response of the Earth to surface loading and to movements of the axis of rotation. In SELEN⁴, we use a real, fully normalized representation for the Stokes coefficients, following GRACE conventions for spherical harmonics (see, in particular, Bettadpur, 2018 and Sect. S8.9 in SSM19). However, it is important to note that in Fig. 9, the Stokes coefficients also include the direct effect of the Earth's rotation, associated with the change in centrifugal potential, on the TLNs of harmonic degree l=2 (i.e., they account for the “δ(t)” term in “$\mathit{\delta }\left(t\right)+k_{\mathrm{2}}^T\left(t\right)$”; see Eq. S173); hence, the rates we have computed are not directly comparable with the GRACE rates. In fact, since in its orbit GRACE is not physically connected with the Earth, it cannot be influenced by the direct rotational effect (the whole issue has been the subject of discussion a few years ago; see Chambers et al., 2010; Peltier et al., 2012; Chambers et al., 2012). The user of SELEN⁴, however, can supply the program with a rotation response function (𝒢^rot) that does not include the direct term in order to produce GRACE-compliant Stokes coefficients which are only indirectly affected by the Earth's rotation.

The fully normalized cosine (squares) and sine (circles) coefficients (${\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }c}}}_{lm},{\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }s}}}_{lm}$) are shown in Fig. 9a only for harmonics with degree l≤6. The dominance of the degree 2 coefficients is apparent, which reflect the symmetries of the $\dot{\mathcal{G}}$ fingerprint in Fig. 8d. We note that since $\dot{\mathcal{N}}=\dot{\mathcal{G}}+\dot{c}$, where $\dot{\mathcal{N}}$ is the absolute sea-level fingerprint in Fig. 8c and c is the FC76 constant (see Eq. 22), the Stokes coefficients for $\dot{\mathcal{G}}$ and for $\dot{\mathcal{N}}$ coincide for l≥2. For reference, the numerical values of the degree l=2 coefficients obtained in the test run are ${\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }c}}}_{\mathrm{20}}=+\mathrm{1.56}$, ${\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }c}}}_{\mathrm{21}}=-\mathrm{0.66}$, ${\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }s}}}_{\mathrm{21}}=+\mathrm{3.37}$, ${\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }c}}}_{\mathrm{22}}=-\mathrm{0.34}$, and ${\dot{\stackrel{\mathrm{‾}}{\mathit{\delta }s}}}_{\mathrm{22}}=+\mathrm{0.05}$ in units of 10⁻¹¹ yr⁻¹; the modulus of all other coefficients is <1 in these units. To better visualize the decay of the Stokes coefficients with increasing l, in the diagram of Fig. 9b, we show the harmonic spectrum defined in Sect. S8.9 (see Eq. S477). By inspection of the spectrum plot, it is now apparent that the energy contained in the degree l=2 harmonic component exceeds by at least 1 order of magnitude all those with l≥3. After a plateau that indicates a substantial power equipartition in the range of harmonic degrees $\mathrm{3}\le l\le \mathrm{7}$, the spectrum clearly shows a red character and decays very rapidly, closely following a power law $\sim l^{-\mathrm{5}}$ (solid line). This result is consistent with those obtained by Spada and Galassi (2015), although they have used a simplified GIA model with fixed shorelines and the traditional rotation theory. They confirm that, for l≥3, the power contained in the GIA-induced regional variations in absolute sea level is negligible when compared with the spectrum observed during the altimetry era.

4 Conclusions

We have presented an updated version of the SLE solver SELEN, which was originally introduced by Spada and Stocchi (2007) and principally meant as a tool for students. Along with a condensed theory background and the basic features of the new program, we have provided a step-by-step description of the outcomes of a medium-resolution test run that requires modest computing resources. However, the run accounts for an up-to-date description of the time history of melting since the LGM and a realistic rheological profile, being based upon a realization of the ICE-6G_C (VM5a) model of Peltier et al. (2015). The outputs of the test run, which cover different temporal scales, have been briefly discussed in order to appreciate some of the possible geodynamical implications. Outputs of a high-resolution test run have been also presented to illustrate the effects of spatial and harmonic resolution on some GIA predictions.

With respect to the original version of the code, two major improvements have been made in SELEN⁴. The first is represented by an increased realism in the description of the GIA process. Indeed, now the program accounts for the migration of the shorelines and for the rotational feedback on sea-level change, which enable a fully topographically and gravitationally self-consistent modelization, in the sense defined by Peltier (1994). Furthermore, SELEN⁴ can be configured assuming two different rotation theories or even excluding rotational effects. The second improvement is in terms of usability, efficiency, and versatility, and covers various aspects. First, the solution of the SLE is now performed by a single Fortran program unit, leaving to a flexible and customizable post-processor the computation of various outputs encompassing the broad spectrum of the GIA phenomenology. Second, on modern multi-core systems, SELEN⁴ can take advantage of multi-threaded parallelism to speed up the most computationally intensive portions of the code. The reader can find details about computational aspects as execution time and code parallelism in Sect. S8.7 of SSM19. Third, the user can easily customize the time evolution of the surface load and the rheological layering of the Earth providing precomputed loading and tidal Love numbers. Last, SELEN⁴ comes with a user guide and with a fully detailed theory background in a supplement, which is particularly meant to illustrate the basic concepts of GIA to young scientists or colleagues and to allow transparency and reproducibility.

In the framework of a recent benchmarking initiative (Martinec et al., 2018), a preliminary version of the new program has been successfully tested against other independently developed SLE solvers, in the particular case in which rotational effects are not taken into account and assuming simplified surface loads. The version of SELEN⁴ that is distributed with the present work reproduces the numerical results published in Martinec et al. (2018), if the code is configured with the parameters employed in the benchmark exercise. After it was progressively developed in various interim versions, SELEN⁴ has now been released to the GIA and to the global geodynamics community as an open-source tool.