Anthropocene isostatic adjustment on an anelastic mantle

Abstract

A general assumption in geodesy is that solid Earth deformation in the presence of recent hydrological and ice loading is well approximated by a purely elastic response. In cases where there is clear evidence that thermal and petrological conditions exist that favor vigorous high-temperature creep behavior, such as in the mantle beneath Iceland, Patagonia, Alaska, Japan, and Svalbard, many response models have been approximated by using a Maxwell viscoelasticity. However, non-Maxwellian transient viscoelastic rheology is required for many post-seismic relaxation studies. Here, we reconsider the solid Earth response in light of highly temperature-dependent transient viscoelastic responses currently favored in the mineral physics and seismological communities. We develop a mantle response Green’s function that accounts for the vertical isostatic motion of the mantle caused by the acceleration of ice mass loss for Greenland and Patagonia measured by spaceborne and airborne remote sensing since 1992 and 1945, respectively. The Green’s function may be used to examine how anelasticity influences the uplift associated with accelerated surface mass loss. We perform a parameter exploration of the constants that define the Extended Burgers Material model, a rheology having firm experimental and theoretical underpinnings, in order to isolate those material model parameters that have the greatest impact on anelastic-isostatic uplift over interannual and interdecadal timescales. Anelasticity corrections may be required for solid Earth vertical uplift in space gravimetric solutions for long-term hydrology and cryospheric change.

1 Introduction

A long standing issue in solid Earth geophysical geodesy is the non-unique choices available for describing rock behavior at periods longer than several hours to days (e.g., Agnew 1979; Lau and Faul 2019; Huang et al. 2021). In this paper, we take up this general issue as it applies to observations of the Earth’s changing climate that mainly come from Global Positioning System (GPS) and Gravity Recovery and Climate Experiment (GRACE) satellites.

Recent loss of ice and water storage mass on the continents is robustly recorded in both space and terrestrial geodetic data that measure the solid Earth response to the change in surface loading. Loss of near-surface water and ice correlate with increasing atmospheric and ocean temperature and the rise in 20th and 21st Century sea-level. The losses are measured as evolving interannual to interdecadal trends in time (e.g., Tapley et al. 2019; White et al. 2022). A ubiquitous observational paradigm today is the recording of accelerations in ice loss and sea-level rise (e.g., Khan et al. 2022a; Kierulf et al. 2022; Hamlington et al. 2022). Both space gravimetric observations using the GRACE and GRACE Follow-on (GFO) data releases and continuous time series of station positions from GPS, record such accelerations as gravitational and solid Earth deformational responses to the changes in surface loads.

Ice mass loss rates in both Antarctica and Greenland have been accelerating since sometime in the mid-20th Century (e.g., Rignot et al. 2019; Shepherd et al. 2020; Mouginot et al. 2019). In Greenland, this acceleration has been recorded in the evolution of GPS-determined rates of uplift (e.g., Khan et al. 2010; Jiang et al. 2010; Bevis et al. 2012, 2019). Observations of Greenland mass balance offer clues for developing a physical model for intercomparison of elastic versus viscoelastic forward model predictions of solid Earth uplift. In Patagonia, such acceleration in loss is likely to have encompassed the early era of both airborne and space-based observations (e.g., Aniya et al. 1996; Rivera et al. 2002) and near-surface water storage now appears to be in steady decline (e.g., Ivins et al. 2011; Richter et al. 2019).

The prospect for detecting the uplift style that is unique to transient viscosity in response to decadal timescale surface load loss (or gain) is challenging. This is because the time-signature of the vertical land motion (VLM) response is expected to mimic that of either elastic or Maxwell viscoelastic signatures. In this paper, we develop a simple model of a single disk over a compressible gravitating half-space with an analytic form of acceleration in time for the surface mass loss. The half-space has one of three rheological types; isotropic elasticity, Maxwell viscoelasticity or Extended Burgers Material (EBM) transient viscoelasticity. The latter also serves as a model of anelasticity (Jackson 2015). Our goal is to reveal the VLM signatures that distinguish the individual rheological types. When convolved with the homogeneous solutions of the momentum equations, the unloading scenario forms a Green’s response function for the Anthropocene (GFA). Our definition of the Anthropocene (Waters and Turner 2022) follows from the Antarctic ice core stratigraphic horizon that was formed in AD 1945, the year of the first atomic bomb blasts (e.g., Magand et al. 2004). The rate of mean sea-level rise determined after the launch of the first dedicated space ocean altimetry missions (1992) exceed the highest known geologically based rates over the last 5000 years by an order of magnitude (e.g., Onac et al. 2022), thus providing testament to the relative strength of a sustained surface mass transport during the Anthropocene.

The predictive theory requires solutions to the glacial isostatic adjustment (GIA) momentum equations (Klemann et al. 2003; Ivins et al. 2022). The mathematical-numerical technique employed involves the Hankel-Laplace transform and convolution with the load function (Wolf 1985). The EBM and Maxwell rheologies have in common a long-term viscosity, $\eta$, and an elastic-spring analog element of rigidity, $\mu$, thus allowing a common Maxwell time constant to be defined, $\tau _M \equiv \eta / \mu$. This time constant allows a systematic comparison of the different time-signatures of uplift that are predicted by the two rheologies.

In the sections which follow, we discuss experimental and observational rationale for the EBM rheology and the GFA, and briefly introduce the mathematical forward model architecture. We focus on two examples, termed geodetic model experiments GME1 and GME2. The two are analogs to the circumstances of ice loss in Greenland and southern Patagonia, respectively. In each experiment, we limit our acceleration model of the ice loss to the eras of comprehensive space- and airborne mapping by remote sensing and limit our long-term Maxwell viscosity to bounds set by both past GIA models and tectonic settings. We assume 42 and 71 year-long accelerations for GME1 and GME2, respectively. Past studies of mantle response to GIA in the two regions allows us to form a priori bounds on the long-term Maxwell viscosity appropriate to the two regions (e.g., Milne et al. 2018; Dietrich et al. 2010; Mark et al. 2022). The distinct geodynamic settings are important considerations for formulating future models and observational strategies for treating recent isostatic responses using geodesy in light of rapid warming, especially in the high-latitude Northern Hemisphere (e.g., Rantanen et al. 2022).

There are many practical limitations of the model intercomparison that is conducted here, as we assume a single disk load on a homogeneous half-space. Practical models of geodetical time series must reconcile the influences of both a finite elastic lithosphere, a finite depth asthenosphere and a surface load that is derived from detailed, spatially varying, determination of mass balance change. However, a main lesson is illuminated by our intercomparison. If we acknowledge our contemporary scientific understanding of transient rheology in olivine, then we must conclude that neither Maxwell or elastic models robustly capture the main deformation physics that drives crustal motion response when secular surface loading occurs at interdecadal timescales.

2 Background: useful 1-D analogs

Spring-dashpot viscoelastic analogs for linear stress relaxing materials, Maxwell, Kelvin, Burgers and EBM rheologies used in solid Earth geophysics are displayed in Fig. 1. Stress relaxing materials are characterized by viscous behavior for long timescales and elastic behavior at very short timescales (examples are in Fig. 1a, c and d), while strain retarding materials retain elastic strength at both short and long timescales (Fig. 1b). The Kelvin element may be coupled in series to an additional spring to form a standard linear solid (SLS).

Fig. 1

Four spring-dashpot 1-D analogs

Models with distributed SLS elements across a broad seismic frequency band form the traditional solid Earth model of anelasticity common in analysis of load tides and the Chandler wobble (e.g., Dahlen and Tromp 1998, p. 215). The traditional SLS power-law distribution model has fallen from favor in low frequency seismology and tidal analysis during the past decade (e.g., Lau and Faul 2019; Kemper et al. 2023).

The Maxwell analog forms the rheological basis for most GIA models in current use. The Burgers analog adds a Kelvin element in series, with short-term weakening influence having a viscosity, $\eta _2$, and has become useful in models of post-seismic flow following earthquakes of moment magnitude, $M_W \, > rsim$ 7.1 (Pollitz et al. 2001). Extended Burgers material analog (EBM) assumes a power-law ($\alpha$) distribution of Kelvin elements in place of the single element of the Burgers analog of Fig. 1c. The distribution is indicated by the sum on intrinsic retardation times $\tau _i$ (see the dashed box in frame of Fig. 1d). Deformation predicted by Burgers rheology can mimic those of the laboratory-motivated extended model, but in this research paper, we focus only on comparisons between the extended Burgers and Maxwell models. It is important to emphasize that while the traditional anelastic model uses an ‘$\alpha$’ power-law distribution with SLS elements, and the formalism of the extended model also uses an ‘$\alpha$’ power-law distribution, the actual values of $\alpha$ among the two models have different consequences for prediction of dispersion, energy dissipation and deformation.

3 Rationale for an EBM rheology

Over the past two decades, seismologists and geodynamicists have increasingly adopted a strategy for interpreting seismic tomographic images that rely on laboratory experimental data sets which span the full range of possible petrological conditions of the upper mantle. The experimental data allow rigorous statistical methods to be applied to the interpretation of seismic velocity structure (e.g., Abers et al. 2014; Eilon and Abers 2017; Richards et al. 2020). Lau and Holtzman (2019) and Ivins et al. (2020) have emphasized the importance of the distribution functions derived from parameterized experimental data (e.g., Jackson and Faul 2010; Qu et al. 2021), for these may allow both tidal and post-seismic deformational constraints to be reconciled. One of the advantages of the EBM rheological law formulated by Jackson and Faul (2010) is that it can be couched as a linear constitutive equation in tensorial form with isotropic viscoelastic equivalents of the elastic shear and bulk compression moduli (e.g., Ivins et al. 2022) and, hence, may be rigorously entered into the formal GIA momentum equations. The momentum equation may be solved for vertical surface displacements as a function of time.

In our previous study, we examined simple loading cases that either involved Heaviside loading or linearly increasing load after model time t = 0, followed by an unchanging load after a time $t_0$. The latter case may have application to the vertical displacement of the sea-floor during rapid melt-water pulse events of the global ocean inundation phase of the Last deglaciation. New relative sea-level (RSL) data modeled by Simon et al. (2020) demonstrate that, in many cases, the transient rheological models out-perform those that are based on Maxwellian rheology in fitting these data. Here, we develop a formal comparison of VLM that are predicted for contemporary climate change related unloading acceleration using the three different rheologies. Our goal is to intercompare the predictions of the three rheologies, as each have compelling rationale motivating their application for interpreting geodetic data. Lau (2023) has shown the importance of EBM transient rheology to the general problem of modeling the sea-level change following global inundation caused by Meltwater Pulse 1A, some 14,000 to 15,000 years ago. Here, we focus on the importance of EBM rheology to interpretations of VLM measurements from contemporary space geodesy. Chanard et al. (2018) has advocated for the importance of mantle transient rheology at seasonal timescales. Although here we focus on interannual and interdecadal timescales, EBM might also help elucidate geodetic observations at this shorter timescale. The next section introduces the notion of mass loss and rheological timescales, much in the same spirit as in recent analysis of glacial isostasy by Lau et al. (2021) using the concept of a frequency-dependent complex viscosity.

4 Timescales for load and rheology

Geodetic techniques now robustly measure changes that have interannual timescales (e.g., Fu and Freymueller 2012; Sasgen et al. 2020; Koulali et al. 2022; Lenzano et al. 2023). These do require a superimposition onto longer timescale changes. Such superposition opens the question of the mechanical behavior of the solid Earth over the longer timescales, usually assumed to be interdecadal. The transient mechanical behavior is well known to characterize post-seismic relaxation (e.g., Han et al. 2014).

In order to better understand how the transient behavior influences models of short-term isostatic uplift, it is important to turn to models having idealized ice history. In the ice history models explored in this paper, no mass loss occurs at model time $t_{-}$ = 0, while at $t_{+} \, \ge$ 0 mass loss is initiated and slowly rises to a maximum rate at $t \, = \, \tau _{dur}$. Our focus is understanding the interaction of rheology and load evolution when $\tau _{dur} \lesssim$ 100 years.

The common fundamental timescale for the Maxwell and EBM rheologies is the Maxwell time, $\tau _M$. For the Maxwell model;

$$\begin{aligned} \tau _M \, \, \equiv \, \, \eta / \mu , \end{aligned}$$

(1)

and for the EBM model;

$$\begin{aligned} \tau _M \, \, \equiv \, \, 2 J_{U1} \eta \, \, \equiv \, \, \eta / \mu _{U} \end{aligned}$$

(2)

where $J_{U1}$ is the unrelaxed shear compliance (e.g., Jackson and Faul 2010; Faul and Jackson 2015; Ivins et al. 2020), $\mu _{U}$ the unrelaxed shear modulus (or rigidity), and we take the elastic rigidity, $\mu$, to be equivalent to $\mu _{U}$. Thus, the long-time scale model behavior of Maxwell and EBM rheologies are identical and characterized by $\tau _M$, similar to the definition adopted in the simple Burgers material (SBM) model (e.g., Caron et al. 2017). EBM rheology has two additional time constants, $\tau _L$ and $\tau _H$, that bound the time frame over which the continuous distribution of Kelvin (or Voigt) elements provide additional stress relaxation. The time constants $\tau _L$ and $\tau _H$ are the lower and upper cutoff times of the distribution, respectively. This distribution effectively weakens the material over the interval $\tau _H \, - \, \tau _L$.

Following Jackson and Faul (2010), we write the compliance function as

$$\begin{aligned} J ( t^{\prime } ) / J_{U1} \, \, = \, \,1 \, \, + \, \, t^{\prime } \, \, + \, \, \Delta \int _{\tau _L^{\prime }}^{\tau _H^{\prime }} {\mathcal {F}} ( \tau ^{\prime } ) ( 1 - e^{ - t^{\prime } / \tau ^{\prime } } ) \; \textrm{d} \tau ^{\prime }.\nonumber \\ \end{aligned}$$

(3)

We have introduced $t^{\prime } \equiv t / \tau _M$ and here retain only the high-temperature background transient term, as in the development by Ivins et al. (2022). The distribution function, representing the role of Kelvin elements as in Fig. 1d, is

$$\begin{aligned} {\mathcal {F}} ( \tau ^{\prime } ) \, \, = \, \, \alpha \frac{\tau ^{\prime (\alpha - 1)}}{\tau _H^{\prime \alpha } \, - \, \tau _L^{\prime \alpha }}. \end{aligned}$$

(4)

There is special significance to the parameter $\Delta$ defined as

$$\begin{aligned} \Delta \, \, \equiv \, \, ( J_{R1} - J_{U1} ) / J_{U1}, \end{aligned}$$

(5)

as it controls the amplitude of the transient response during the time $\tau _L$ to $\tau _H$. Note in Eq. (3) and (4) we have introduced an exponent that governs the shape of the distribution, $\alpha$, and two dimensionless time parameters $\tau _L^{\prime }$ and $\tau _H^{\prime }$, also scaled to $\tau _M$. In Eq. (5), $J_{R1}$ is termed the relaxed compliance with $J_{R1} \, > \, J_{U1}$, meaning that as time passes the material becomes more compliant, or ’softer’.

The relaxation strength parameter may also be written as the ratio

$$\begin{aligned} \Delta \, \, = \, \, ( \mu _{U} - \mu _{R} ) / \mu _{R}, \end{aligned}$$

(6)

where $\mu _{U}$ and $\mu _{R}$ are the unrelaxed and relaxed shear moduli, respectively. For our purposes, $\mu _{U}$, may be referenced to a seismic frequency. SBM rheology has the pair of parameters, $\eta _1$, $\eta _2$, representing the viscous elements of the Maxwell and Kelvin elements, respectively, and $\mu _1$, $\mu _2$, the corresponding shear moduli. The relaxed SBM rigidity, $\mu _R \, \equiv \, \mu _1 \mu _2 / ( \mu _1 + \mu _2 )$, may be derived by setting to zero the 2nd degree time derivatives on both strain and stress variables in the 1-D analog constitutive equation of Reiner (1958) (see his Eq. 20.4). The SBM ratio of rigidities, ${\bar{\Delta }}$, is the complete analog to the EBM relaxation strength, $\Delta$, or

$$\begin{aligned} {\bar{\Delta }} \, \equiv \mu _{1} / \mu _2. \end{aligned}$$

(7)

Table 1 Rheological parameters

This is consistent with the definition advocated for SBM by Faul and Jackson (2015) and consistent with the amplitude factors in the creep function identified in Ivins et al. (2020). While the two strength relaxation parameters, ${\bar{\Delta }}$ and $\Delta$, are not identical, they do play analogous roles for the two viscoelastic rheologies. The comparison allows us to infer SBM parameters from of post-seismic data sets. A cautionary note is that the ’Achilles heel’ of SBM rheology is its inability to capture solid deformational phenomenon across the tidal spectrum (e.g., Lambeck and Nakiboglu 1983; Ivins et al. 2020) due to its narrow dissipation band (e.g., Cooper 2002). The rheological parameters relevant to our discussion are summarized in Table 1.

Fig. 2

Load function applied to IMBIE-2 reports for Greenland ice mass balance. The inset shows the formal dimensionless functional form used for the mathematical architecture of the GFA for the case of R = 10. The case of R = 10 may be interpreted as the duration of the acceleration ($\tau _{dur}$) being shorter than the Maxwell time ($\tau _M$) by a factor of 1/10. For example, in the IMBIE-2 case shown in the inset, with R = 10, $\tau _{dur}$ = 22 years, and $\mu _{U}$ = 67 GPa, corresponds to a Maxwell viscosity, $\eta \, \approx \, 4.65 \, \times \, 10^{20}$, roughly consistent with the Maxwell estimate for upper mantle viscosity beneath Greenland

4.1 Time-dependence of the load function

As space and airborne observations have been systematically monitoring change in water/ice mass since the mid-20th Century, acceleration in change has been measured on interdecadal time scales. As one example, in Fig. 2 we show the subset of trends (blue dots) that were published in the intercomparison summary for the 2nd Ice Mass Balance Intercomparison Exercise (IMBIE-2) (IMBIE-Team 2020). For this case, it may easily be seen that an acceleration function captures the essence of this evolution in mass for the whole of the ice sheet and peripheral glaciers of Greenland. It should be carefully pointed out that if we were to examine the individual drainage basins, we would not find such a simple story, and this owes to the fact that the losses additionally have a geographical evolution (e.g., Khan et al. 2010; Mouginot et al. 2019).

Table 2 Geodetic model experiments

The detailed specification of a 42-year-long acceleration for GME1 (Greenland loss) in our study is not particularly relevant to the intercomparison of rheologies, nor to contrasts with the presumed weaker mantle of the GME2 case. Yang et al. (2016), for example, has advocated a similar acceleration model for Greenland that is initiated with zero-loss rates at 1990. Neither the assumption of zero-loss rates at 1980, or at 1990, and finite loss rates thereafter, are inconsistent with the 1-$\sigma$ values during the first few years of the 1990s in the latest IMBIE report (e.g., Otosaka et al. 2023, ). The important point is that the proxy surface load model in both the experiments is constantly accelerating, driving a constantly evolving GIA response as vertical surface uplift even in the elastic case.

Our goal is to simplify the physics to the extent that we may perform a set of experiments that are decoupled from the nuances of the time evolution as observed in space and time, and that capture the interdecadal evolution in continental bound ice mass. In the case of mass loss in Greenland and Patagonia, there has been an acceleration phase followed by years of constant mass loss rate. The loss evolution in the two cases may be grossly captured by an approximate, but analytical, function, $\vartheta _{o} ( t^{\prime } )$, as follows:

$$\begin{aligned} \vartheta _{o} ( t^{\prime } )= & {} \left[ - 1 \, + \, Cos \left( \frac{\pi R t^{\prime } }{2} \right) \right] H [ R t^{\prime } \,, \, 1 - R t^{\prime } ] \, \nonumber \\{} & {} + \, \frac{1}{2} ( - 2 \, - \, \pi R \, t^{\prime } \, + \, \pi ) H [ R t^{\prime } \, - \, 1 ]. \end{aligned}$$

(8)

wherein $t^{\prime }$ is dimensionless (scaled to $\tau _M$). We note that a factor R is related to a transformed time, $t^{\prime \prime } \, \equiv \, \frac{\tau _M }{ \tau _{\textrm{dur}}} t^{\prime }$, with $\tau _{\textrm{dur}}$ representing the total duration of the acceleration phase of the surface load. The function also accounts for the linear change in load that occurs after $t^{\prime } \ge \frac{1}{R}$, with $R \equiv \, \frac{\tau _M}{ \tau _{\textrm{dur}}}$. The function $\vartheta _{o} ( t^{\prime } )$ and first derivative are continuous at $t^{\prime } = \frac{1}{R}$. The first term is associated with the Heaviside of two arguments; $t^{\prime } \cdot R$ and $1 - t^{\prime } \cdot R$, is associated with the accelerated unloading. (When either of the two arguments is less than zero, the function is zero). The inset in Fig. 2 shows the function $\vartheta _{o} ( t^{\prime } )$ for the case R = 10. In Fig. 2, the complete acceleration phase is shown. The second term associated with the Heaviside function in Eq. (8) represents a longer-term linear ramp (not shown in Fig. 2). To streamline the later discussion, we may separate the two loading phases for $t^{\prime } \le 1/R$ and $t^{\prime } \ge 1/R$, as $\vartheta ^{(-)}_{o} ( t^{\prime } )$ and $\vartheta ^{(+)}_{o} ( t^{\prime } )$, respectively.

An important facet of the load function time-dependence is that the embedded parameter R is linked to the assumed model $\tau _M$. For example, a case with $R \, > \, 1$ corresponds to a situation where the acceleration phase is shorter than the Maxwell time, and hence, we anticipate a rich exhibition of transient-like geodetic signal.

In our two experiments, GME1 and GME2, we explore predictions in the contrasting cases $R \, > \, 1$ and $R \, < \, 1$, respectively. The two situations correspond to what we generally understand about mantle viscosity beneath the two ice regimes of Greenland and southern Patagonia, one having upper mantle $\eta \, \ge \, 3 \times 10^{20}$ Pa s, and $\eta \, \le \, 1 \times 10^{19}$ Pa s, respectively (e.g., Lecavalier et al. 2014; Mark et al. 2022). We summarize the contrasting values in Table 2.

The last column in Table 2 gives a notion of the uplift rates (${\dot{w}}$) as observed in Greenland and southern Patagonia, respectively. In Greenland, there are more than 57 stations distributed all along the coastlines peripheral to the ice sheet. Residuals to GIA models with elastic correction are generally at the level of 2 – 5 mm/yr (e.g., Adhikari et al. 2021) when applying the an optimal GIA model (HUY3) with respect to a large number of regional RSL data (e.g., Lecavalier et al. 2014; Milne et al. 2018). It should be noted that Kappelsberger et al. (2021) have shown that HUY3 is the best preforming model for northeastern Greenland with respect to geodetically determined uplift trends. Generally, we wish to demonstrate that a residual of about 2 – 5 mm/yr can be explained by anelasticity.

The potential impact of EBM rheology on the interpretation of LIA and shorter unloading scenarios for Greenland is uncertain. However, a recent model analysis of the combination of Greenland upper mantle seismic data and parameterizations of both EBM and a semi-empirical broadband rheology have led Paxman et al. (2023) to conclude that the most significant reductions in apparent mantle strength, or effective viscosity, occur on time scales of $10^2$ – $10^3$ years. This inference helps to motivate the simple investigation of GIA responses that we undertake in the EBM broadband rheological experiment GME1.

In southern Patagonia, 43 GPS stations have recorded uplift trends. Some stations lie very close to the the Northern Patagonian (NPI) and Southern Patagonian Icefields (SPI) (e.g., Richter et al. 2016). The largest recorded uplift rate at one station is roughly 41 mm/yr (also see Lenzano et al. (2023)) which can be explained by a Maxwell GIA model with viscosity as low as $\mathrm {1.6 \times 10^{18}}$ Pa s and a load history that includes the LIA and 20th Century ice loss (e.g., Lange et al. 2014; Richter et al. 2016). The simple anelastic GIA model may, in fact, demonstrate that part of the uplift signal is caused by the ’softness’ induced by a spectrum of mantle relaxation processes over the timescale $\tau _H \, - \, \tau _L$.

Fig. 3

Load function and estimates published by Rivera for southern Patagonia using combined airborne and spaceborne mass change and spaceborne rate-based estimates (a). Corresponding load function represented evolution of mass loss rates for SPI + NPI + Cordillera Darwin and regional adjacent water storage (b). Here, this function is used explicitly for the GME2 case of uplift intercomparisons

Figure 3 shows the form of $\vartheta _{o} ( t^{\prime } )$ that is used to represent changes since 1945 as reported by Rivera et al. (2002) for the ice mass of Patagonia. These include the NPI and SPI. Here the mass change (a) and rate of change (b) are with respect to time in years, and hence there is no reference to a value for R in Fig. 3. The functional form also accommodates the mean rate of loss (Fig. 3b) over the GRACE observing period (e.g., Ivins et al. 2011; Richter et al. 2019).

Table 3 Static GFA parameters (see Fig. 4)

4.2 Spatial-dependence of the load function

Figure 4 shows the geometry assumed for loading of our experiments GME1 and GME2. At time $t_{-} \, \le \, 0$ the half-space is in gravitational equilibrium with the square edged axisymmetric disk load (of radius, $\alpha _D$), while at $t_{+} \, \ge \, 0$ mass begins to be lost, and a finite stress is imparted to the half-space. Our goal is then to predict the time evolution of the vertical displacement, w, at the top ($z = 0$) of the half-space as a function of time. We select a site for evaluation of $w(z = 0, r_0, t)$ that is close to the center of the disk, for it can be shown that the structure of the solution is quite uniform when $r_0 \,<< \, \alpha _D$ (e.g., Wolf 1985; Ivins et al. 2022). It will be especially instructive to predict the evolution of w(t) and first derivative ${\dot{w}} ( t )$ as $t \rightarrow \tau _{dur}$ for the assumptions of elasticity, Maxwell viscoelasticity and our prototype of an anelastic material, the EBM viscoelasticity. In terms of the time scaled to the Maxwell time, $t^\prime$, this is also the evolution as $t^\prime \rightarrow 1 / R$.

Fig. 4

Square-edged disk load on gravitating half-space. The fourth-order differential equation system uses stress boundary conditions at the surface (as shown) and regularity of the solutions at $z \rightarrow \infty$. The static parameters of the model are given in Table 3

Here, we briefly introduce the theory for the homogeneous half-space case by introducing the Laplace transform equations of motion for GIA. The latter may be written as

$$\begin{aligned} \partial _i {\tilde{\sigma }}_{ij} ( s ) \, \, + \, \, \partial _i {\tilde{u}}_{j} ( s ) \partial _j p \, \, - \, \, \rho g_i \partial _j {\tilde{u}}_{j} ( s ) \, \, \, = \, \, \, 0. \end{aligned}$$

(9)

The Laplace transform is defined as

$$\begin{aligned} {\mathcal {L}} \{ F ( t ); t, s \} \, \, \equiv \, \, \int _{0}^{\infty } e^{-st} F ( t ) \, d t \, \, \equiv \, \, {\tilde{f}} ( s ) \end{aligned}$$

(10)

where we use the over-mark tilde ($\tilde{...}$) to represent parameters that are to be interpreted as having been transformed from the domain of t to s. Here, the Cauchy stress ${\tilde{\sigma }}_{ij} ( s )$ is expressed by the compressible form of the constitutive equation for EBM rheology as given in Sect. 3 of the recent paper by Ivins et al. (2022). In Eq. (9), ${\tilde{u}}_{j}$ is the displacement vector. The two terms to the right of the divergence of the stress tensor are, respectively, the advection of the hydrostatic pre-stress, p, with $\hat{{\textbf{k}}} \cdot \nabla p \, = \, \rho g$, and the buoyancy force owing to changes in density caused by finite local dilatation, $\partial _j {\tilde{u}}_{j}$, and $\hat{{\textbf{k}}}$ defined as the unit dyad vector in the z direction. Klemann et al. (2003) studied the character of the system as written in Eq. (9) with solutions for the homogeneous half-space problem. Following the guidance provided by Klemann et al. (2003), we set the last term in Eq. (9) to zero, thus reducing the system to that studied by Wolf (1985). With a homogeneous half-space as an assumption, we also take the downward directed gravity force $\rho g_i$ to be constant.

4.3 Hankel transform

The solution method for the time-dependence may be conducted analytically using the normal mode formalism of a viscoelastic transfer function that was developed by Wolf (1985) and applied to SBM rheology in Rümpker and Wolf (1996). Ivins et al. (2022) showed how the function may be used to setup the numerical integrations to solve for the inversion of the Hankel-Laplace transformed vertical displacement variable at the surface, $\tilde{\hat{w}}_0 (k, s, z = 0)$ in the case of EBM rheology. Here, we use the symbol $\hat{ }$ to indicate the Hankel transform domain, with transform variable k. In this paper, we relegate these details of the convolution with the load to Appendix A with examination of solutions in the Hankel transform domain to Appendix B. In the latter, we examine the dependencies of short-term responses ($t^\prime \le 1/R$) on $\Delta$, $\alpha$, and the ratios $\tau _L / \tau _H$, $\tau _H / \tau _M$, identifying $\Delta$ as the most critical parameter for geodetic applications.

The method by which we can develop the load response at the half-space surface to the forcing function $\vartheta _0 ( t )$ is to employ the natural basis functions for the axial symmetry of a disk load. The basis functions are Bessel’s functions, $J_n ( k_m r )$ where r is the radial coordinate of cylindrical coordinate frame as in Fig. 4. Here, n is the order of the Bessel function of accompanying wave number $k_m$, having units of inverse length. The orthogonality of the functions allow them to be used to represent any well-behaved function f(r) in the interval 0 to $\alpha _D$. For a function f(r) which is 0 for $r > \alpha _D$, then

$$\begin{aligned} f ( r ) = \sum _{n = 0}^{\infty } c_m J_n ( k_m r ), \end{aligned}$$

(11)

with $k_m$ chosen to such that $J_n ( k_m \alpha _D ) = 0$. The orthogonality condition of Bessel functions allows each $c_m$ to be determined as

$$\begin{aligned} c_m = \frac{\int _{0}^{\alpha _D} f ( r ) J_n ( k_m r ) r dr}{\frac{{\alpha _D}^2}{2} [ J_{n+1} ( k_m \alpha _D ) ]^2}. \end{aligned}$$

(12)

If the function f(r) is finite at $r = 0$, then all functions except those with order $n = 0$ are eliminated from this expansion (e.g., Mathews and Walker 1970, p. 178). The space we consider has coordinate r extending to $\infty$, hence the wave numbers $k_m$ are continuous, and we may approach the formal boundary value problem with the formalism of the Hankel transform. This approach is described by Farrell (1972) and Wolf (1985). Ivins et al. (2022) discusses the transform as it applies to our disk load shown in Fig. 4. It is also discussed here in Appendix A (see Eqs. A1-A9) and the reader is recommended to consult Wolf (1985) and Klemann et al. (2003) who detail the forward Hankel-Laplace transform of Eq. (9). While elegant in its formalism, the full inverse Hankel-Laplace transform must be computed numerically for EBM and Maxwell materials.

The solution for the vertical displacement at the surface of the half-space is

(13)

where the integral represents the reduced inverse Hankel transform. Eq. (A9) of Appendix A, and its derivation, detail the time-dependence, wherein (in units $m^3 / N$) represents the inverse Laplace transform Green’s function that convolves both solid Earth structure parameters and ice load history. Subscripts $\nu = 1 - 3$ are for elastic (E), Maxwell viscoelastic (M) and extended Burgers Material (EBM), respectively. The inversions are discussed in Appendix A. Much of the inter-comparative physics of the rheologies are revealed in the Hankel transform domain (see Appendix B). Henceforth, the discussion assumes that the load disk density ($\rho _D$) is that of ice, $\rho _{ice}$. In Eq. (13), $\delta h$ is the change in total height over the time span $\tau _{dur}$. A summary of basic elements needed for deriving Eq. (13) are provided in Table 4.

Table 4 GFA Theory

4.4 Amplitude of the transient

Let us now return to some basic questions concerning the intercomparison of rheologies. A fundamental amplification factor is $\Delta$ in Eq. (3). There are two avenues to pursue an a priori value, or at minimum some bound on geophysically reasonable values for $\Delta$. A first avenue is found by using values that are reported in laboratory torsion experiments and we give a summary in Table 5.

Table 5 Experimental ‘optima’ $\Delta$ for EBM

In Table 6, the definition of ${\bar{\Delta }}$ that Faul and Jackson (2015) and Jackson (2015) define for the SBM rheology is employed, thus allowing values of the relaxation strength from SBM studies to play a role in assessing viable a priori estimates of $\Delta$ for EBM models. SBM rheology is assumed in numerical models of GPS-determined post-seismic relaxation (e.g., Pollitz 2003; Broerse et al. 2015, P03 & BRSGV15). In addition, RSL data-constrained GIA models (e.g., Simon et al. 2020, SRB22) also serve as a basis for estimating a bound on ${\bar{\Delta }}$. In addition, a retrospective study has used laboratory creep data with naturally occurring mantle rock samples (dunites) (e.g., Chopra 1997) to better characterize the value of ${\bar{\Delta }}$ for the SBM rheology (e.g., Masuti and Barbot 2021, MB21). The values relating to SBM studies are summarized in Table 6.

Experimental values for $\Delta$ reported by Jackson and Faul (2010) (JF10) range from 0.2 to 1.9, with $\Delta = 1.4$ reported as an optimum value that is determined from a large number of experiments. More recent experiments by Qu et al. (2021) (QFJ21) using polycrystalline olivine and pyroxene suggest values slightly less than 1 and earlier experiments by McCarthy et al. (2011) (MTH11) also report smaller $\Delta$. Considerably larger values for EBM and SBM-based modeling are reported using geodetic and RSL data, ranging from 2 to 10 (see Table 6). Experiments using SBM are not considered, as there are no viable proposals that SBM forms a basis for solid Earth anelasticity that bridge seismic, tidal and GIA timescales.

Table 6 Geodesy and reanalysis of SBM relaxation strength

Water content is known to strongly reduce olivine steady creep strength above 50–100 ppm. Upper mantle rock generally has water fugacity levels that exceed this. It is generally understood that water content enhances seismic attenuation and dispersion (e.g., Liu et al. 2023, ). Currently, experimental quantification of EBM parameters relies primarily on dry olivine rock samples. Rigorous determination of the relaxation strength, $\Delta$, as a function of water fugacity must await future experimentation.

In the section which follows we demonstrate the significant role of dimensionless parameters R and $\Delta$ for predicting the magnitude of a geodetically detectable anelastic-EBM-related uplift for the 21st century monitoring of surface load changes.

5 Acceleration responses

In the case of the analog acceleration of ice mass loss to Patagonia, Fig. 3 shows the detailed model mass loss and rate over the full duration of the acceleration in the experiment GME2. The full spaceborne record has been summarized for Greenland by Mouginot et al. (2019) and this is used as a guide for constructing a model experiment. The GME1 evolution of unloading rate increase is best represented by the reconstruction reported for the Southeast Greenland sector (Fig. 3 g of Mouginot et al. (2019)) for which we select $\tau _{dur}$ = 42 years, a total mass loss of $-1.5 \times 10^3$ Gt, and a disk radius $\alpha _D$ = 375 km. Using a constant ice density of $\rho _{ice}$ = $\mathrm {0.89 \, Gt / km^3 \, = \, 890 \, kg/m^3}$, the mean rate of height change is over the model 42 years is $\mathrm {-4.5 \, cm/yr}$ for the model selected disk radius. The mean model rate of loss ${\dot{M}} \, = \, \mathrm {- 35 \, Gt/yr}$, roughly consistent with changes estimated by Hurkmans et al. (2014) for the southeastern basin of Greenland. These details are for clarification of the model experiment performed here. The main point is that the acceleration function applied corresponds to Eq. (8) and Fig. 2. The experimental setup is unlike Greenland ice loss in that the GME1 loss occurs uniformly over a disk, rather that being concentrated along the coastline (e.g., Hurkmans et al. 2014; IMBIE-Team 2020; Mouginot et al. 2019). Again, the model seeks only to characterize the different physics involved in the prediction of elastic, Maxwell viscoelastic and EBM viscoelastic vertical crustal rebound.

5.1 Uplift in the GME1 experiment

In Fig. 5, the uplift and uplift rate prediction is shown for GME1 over the course of the accelerating model mass loss. The history of uplift and uplift rate shown indicates that the uplift rate for EBM model is roughly a factor of 2 larger than both the elastic and Maxwell predictions.

Fig. 5

Near-center uplift (w) and rate (${\dot{w}}$) for the GME1 experiment (solid and dashed lines, respectively). $\tau _{dur}$ = 42 yr. Half-space unrelaxed elastic moduli are $\kappa \, = \, 124.44 \,$ GPa, $\mu \, = \, 67.0 \,$ GPa for compression and shear, respectively. Density is assumed $\rho \, = \, 3380 \, kg/m^3$ and gravity constant at $g \, = \, 9.8 \, m/sec^2$. The density of ice assumed is $\rho _{ice} \, = 890 \, kg/m^3$. EBM parameters are as indicated and $\alpha = \,$0.5 and $\tau _L = \,$ 0.0

In the case of GME1, we assume that the half-space viscosity, $\eta$ is well understood and constrained by non-geodetic data. This is similar to the case of Greenland and we assume a value of the viscosity close to those assumed by Lecavalier et al. (2014) and Milne et al. (2018). With a value of viscosity well above $1.0 \times 10^{20}$ Pa s, the ratio of Maxwell time to load acceleration time (R) is well above 1. In the example shown in Fig. 5, $R \, =$ 3.3. The importance of this dimensionless parameter R will become more obvious when we examine the experiment GME2. Using the fixed ratio, R, as in the GME1 experimental result shown in Fig. 5, we may examine the prediction of the difference of the EBM and elastic uplift rates ($\delta {\dot{w}}$) when evaluated at the time $t \, = \, \tau _{dur}$. The results are shown in Fig. 6 for varying EBM rheology parameters $\tau _H$ and $\Delta$.

Fig. 6

$\tau _H$ and $\Delta$ sensitivity for predicting uplift rates at $t = \tau _{dur}$ in excess of elastic rate. $\alpha$ and $\tau _L$, solid Earth and ice parameters are the same as Fig. 5

The contours of constant excess uplift $\delta {\dot{w}}$, relative to the elastic uplift rates, shown in Fig. 6 reveal the dominant dependence upon the relaxation strength, $\Delta$, as discussed in Sect. 4.4. For example, a selection of $\Delta \, = \, \mathrm {0.87}$ predicts about 1.7 mm/yr greater uplift than elastic, quite decoupled from the dependency on $\tau _H$. Note that as $\Delta$ increases to value above those reported in torsion experiments at $\sim 2.0$ (e.g., Jackson and Faul 2010), the dependence on $\tau _H$ increases. In the example shown in Fig. 6, if the EBM parameters are set to $\Delta \, = \, \mathrm {4.4}$ and $\tau _H \, = \, \mathrm { 9 \, yr}$, the excess uplift is predicted more than 4 mm/yr, or three times the elastic rate shown in Fig. 5. We attribute the stronger uplift for smaller $\tau _H$ to the fact that the ensemble of micro-deformational mechanisms that are captured in the integral from $\tau _L$ to $\tau _H$ have more effective weakening, as they are concentrated into a shorter period of time. This dependency is not easily revealed for the case of a linearly varying load (see Ivins et al. (2022)), wherein the external load change has no higher order derivatives. In the accelerating load case, stress and stress rate in constantly changing, so the action of the ensemble of weakening effects are omnipresent during the acceleration.

Fig. 7

Uplift and rate over 71 years in the GME2 experiment for the minimum torsion experimental value of $\Delta$. Note $t \, = \,\tau _{dur}$ = 71 yr. Predictions using the three rheologies are shown. The half-space elastic parameters and other EBM parameters are as in Fig. 5, $\rho _{ice} \, = \, \mathrm {850 \, kg/m^3}$. Solid lines are for uplift, dashed for uplift rate

Fig. 8

Uplift and rate over 71 years in the GME2 experiment for the ’optimum’ torsion experimental value of $\Delta$. Solid lines are for uplift, dashed for uplift rate. Other parameters as in Fig. 7

5.2 Uplift in the GME2 experiment

We now turn to the case of experiment GME2, patterned after Patagonian glacier loss (e.g., Aniya et al. 1996; Rivera et al. 2002; Richter et al. 2019) and tectonic setting (e.g., Mark et al. 2022). Our best estimates for the Patagonian Icefields is that the upper mantle viscosity is weak and the accelerating loss longer than in southeastern Greenland. Hence, our setup for GME2 will assume a half-space $\eta \, \le \, 1 \times 10^{20}$ Pa s, and $\tau _{dur}$ is as long as the combined actual airborne and spaceborne observations which we set to 71 years. Hence, we will examine cases in which $R \, \le \, 1$. In this case, a corresponding upper mantle viscosity in southern Patagonia may be very low, and it possibly is characteristic of a youthful slab window environment (e.g., Lange et al. 2014; Mark et al. 2022). The EBM sensitivities in the case of our second model experiment (GME2) are more nuanced than in GME1.

In Fig. 7, we show uplift history and rate for the case of the dimensionless R and $\Delta$ having relatively small values (0.2 in both cases). Here, the experiment shows that EBM and Maxwell predictions closely track one another (magenta and red lines, respectively). This is a somewhat predictable result for two reasons: the amplitude of EBM weakness is suppressed (low $\Delta$) and the Maxwell time is exceeded during the time evolution in the experiment ($\tau _M \, < \, \tau _{dur}$), yielding a more robust Maxwellian rebound than in GME1.

Fig. 9

Uplift and rate over 71 years in the GME2 experiment using a high torsion experimental value of $\Delta$ and a high end value for the long-term mantle viscosity in a slab-wedge tectonic environment. Solid lines are for uplift, dashed for uplift rate. Other parameters as in Fig. 7

Figure 8 assumes a value of $\Delta$ increased to the optimum value reported by Jackson and Faul (2010) with the same value of R as in Fig. 7. The experiments shown in Figs. 8 and 9 demonstrate that as $\Delta$ increases there is large separation of Maxwellian and EBM ${\dot{w}}$. The higher values of $\Delta$ offer a prediction of nearly 40 mm/yr uplift rate, about 8 mm/yr faster than the Maxwell case (Fig. 8), and close to the maximum uplift observed geodetically (e.g., Richter et al. 2016). In Fig. 9, the value of R is raised to 1.5 (by increasing the assumed viscosity to $\eta \, = \, 1.0 \times 10^{20}$ Pa s), simultaneously raising the relaxation strength to $\Delta \, = \, 1.9$, the highest value found in the torsion experiments of Jackson and Faul (2010). In this case, in spite of the increase in $\Delta$ relative to the case in Fig. 8, the experiment has smaller ${\dot{w}}$, lowering by more that 10 mm/yr, at model time t = 71 yr. This is because the strengthening due to higher long-term viscosity now has a compensating (and relatively stronger) influence on the solution than does the increase in $\Delta$. The three cases of GME2 uplift shown in Figs. 7, 8 and 9 suggest that the tradeoff of long-term $\eta$ and relaxation strength $\Delta$ have greatest influence on the rates of crustal rebound throughout the history.

Fig. 10

Residual uplift rates sensitivity to $\eta$ and $\Delta$ sensitivity at $t = \tau _{dur}$ for GME2. A residual value of 0 means that the viscoelastic model explains the peak uplift of 40 mm/yr of the GME2 experiment. Other parameters as in Fig. 7

For the purpose of studying the sensitivity to combinations of $\eta$ and $\Delta$ in experiment GME2, we test the model rate ${\dot{w}} - max\{\textrm{observation}\}$ at $t \, = \, \tau _{dur}$ in Fig. 10. This investigation of the tradeoff in assumed half-space viscosity and relaxation strength in experiment GME2 reveals the more important role of long-term viscosity than EBM parameter $\Delta$. In Fig. 10, we show the residual of the uplift rate at t = 71 yr to that maximum value reported by Lange et al. (2014) and Richter et al. (2016) (i.e., predicted rate minus 40 mm/yr), thus a contour value of 0.0 meets the observational criteria.

Clearly, the case of $R \, \le \, 1$ requires consideration of the strength of relaxation and the value of the viscosity, or the Maxwell time of the mantle layer for which EBM rheology is assumed. In the specific case shown in Fig. 10, sensitivity to $\Delta$, relative to the Maxwell viscosity, is relatively mute for the $\Delta$-$\eta$ combinations that satisfy the uplift requirement of experiment GME2 at the zero contour curve.

Table 7 Comparison of VLM rate evolution with $\Delta = 1.9$

Model computations shown in Figs. 5 and 9 with $\Delta = 1.9$ offer a comparison of the model VLM rate evolution for GME1 vs GME2. In Table 7, we give the ratios of the predicted rates of VLM at 10, 20 and 30 years after the initiation of the accelerating ice unloading for EBM vs both elastic and Maxwell constitutive assumptions. At this value of $\Delta$, VLM rate evolution over the first 10–30 years of surface load acceleration predicts robust EBM enhancement by 50 - 100 percent over both elastic and Maxwell predictions.

6 Summary and conclusions

The two idealized experiments with accelerating mass loss are guided by the mass loss and regional VLM observed in southeast Greenland and southern Patagonia. In order for the geodesy community to generalize some lessons learned from the experiments to other regions of recent rapid surface mass change, the tectonic setting is an important consideration. Neotectonic regions (tectonic activity during the past 20 million years) generally have much weaker shallow mantle/crust than do regions of stable crust and lithosphere. Our experiments offer some guidance for assessing how transient EBM viscoelasticity will impact models of isostatically driven VLM. Table 8 summarizes generalizations that we infer from the two experiments conducted in this paper. In the stable continental environment, wherein estimates of upper mantle viscosity exceed $1 - 2 \times 10^{20}$ Pa s, transient rheology will invigorate any VLM that is induced by decadal scale ice/water loss. By contrast, in the Neotectonic environment, the details of the Maxwellian long-term viscosity and the strength of the transient rheology (here measured by $\Delta$) have a more nuanced interaction: at the lowest end of a priori estimates of $\eta$, Maxwell solutions to the GIA equations adequately explain the evolution of VLM, though the introduction of EBM transient rheology will increase any geodetically constrained GIA model estimate of $\eta$. As the a priori estimate of $\eta$ increases, transient viscosity should be imposed in the geodetic-GIA model in order to accurately predict the time-dependency of VLM. In all cases, accounting for EBM rheology will correct a bias of Maxwellian models. The essential bias is that models having only Maxwellian viscosity may predict too low a long-term mantle viscosity.

Table 8 Lessons for surface load geodesy in the anthropocene: summary

The most important conclusion deduced from the two geodetic experiments presented here pertains to the interpretations of solid Earth signal in geodetic observations over the duration of the Anthropocene. The time frame spans both the space and airborne observational record, including the GRACE and GRACE-FO missions (e.g., Scanlon et al. 2018; Liu et al. 2022) and the period of robust climate reconstruction that began after 1980 (e.g., Cowtan and Way 2014). The implications of EBM viscoelasticity are important in all aspects of geodesy that involve interannual to centennial timescale mantle deformation, including the surface load models of center-of-mass-center-of-figure offset rate (e.g., Lambert et al. 2023) which are highly sensitive to lower mantle rheology.

A clear message from experiment GME1 is that the concept of an ‘instantaneous elastic response’ must be viewed with some skepticism when considering geodetic observables on interannual, or longer, timescales. Both homogeneous half-space experiments GME1 and GME2, on the other hand, suggest that forward model strategies for interpreting geodetic observations of rapid uplift should consider transient rheology. Our approach has been to tie viable transient rheologies to formulations that come from the mineral physics community and which find application to both tidal and post-seismic relaxation deformation data (e.g., Lau and Faul 2019; Boulze et al. 2022). The transient EBM rheology, using a range of plausible constants for the material behavior, promotes the anelastic part of the solid Earth response as the observing period increases. Anelastic corrections may be needed for inversions for mass loading changes and incorporated into future gravity mission plans (e.g., Wiese et al. 2022) and GPS observations adjacent to ice systems that are projected to accelerate through the 21st Century (e.g., Khan et al. 2022b). We may anticipate an interesting future in geodesy as the computational software for solving Love’s classical boundary value problem with laboratory-based transient flow laws begin to flourish.

Some important caveats should be considered when integrating what we have learned from our two experiments with geodetic data. First, the replacement of a simple homogeneous half-space with one having a lithosphere and finite asthenospheric depth will broaden the longer timescale response character of the solid Earth in ways that will depend on the choices for structure and parameter values. For example, we have not explored how a radially structured EBM model might add new complexities to the spatiotemporal evolution of bulge migration and horizontal crustal velocities. Secondly, while the EBM formulation of the rheology is convenient from a modeling perspective, it is not the only one being considered in the mineral physics community. Alternatives might include microscopic-model physics that is controlled by dislocation dynamics (e.g., Wallis et al. 2021), and consequently, the corresponding model for isostatic motions would require a nonlinear transient viscoelasticity. Transient rheology may also be caused by smaller scale creep heterogeneity (Ivins and Sammis 1996). Currently, however, there is a dearth of constraints on the relative volume and distribution of weak viscoelastic components in the sub-lithospheric environment. Thirdly, there is the question of fully understanding the error budget associated with observing the elastic-gravitational load response. These include non-periodic atmospheric and oceanic loads (e.g., Adhikari et al. 2017; Pan et al. 2022; Gómez et al. 2022) and shallow crustal deformational responses to changes in upper crustal water storage that are non-gravitational in nature (e.g., Neely et al. 2021; Razeghi et al. 2022).

The main purpose of this paper is to stimulate geodesists to reconsider recent isostatic motions in response to secular changes in the climate and the associated water loads, in light of the recent advances in our understanding of transient rheology and anelasticity within the solid Earth community. One goal of future model research with transient viscoelasticity in geodesy could be to better quantify the uncertainties in quantitatively relating VLM observations to contemporary surface mass change.

Appendices

Appendix A GIA Momentum Equation: Half-space Solutions

The geometry of the IVBVP is shown in Fig. 4. In the normal mode treatment by Wolf (1985), the z-dependency is solved as an eigenvalue problem upon separation of the partial differential equation into r-dependent and z-dependent parts. The displacement field in cylindrical coordinates is ${\textbf{u}} ( r, z ) = u \hat{{\textbf{r}}} + w \hat{{\textbf{k}}}$, where $\hat{{\textbf{r}}}$ is the unit dyad vector in the radial coordinate (see Fig. 4). The radial part of the field is solved using the first order Hankel transform,

$$\begin{aligned} \hat{f}_1 ( k ) \, \, \equiv \, \, {\mathcal {H}}_1 \{ f ( r ); k \} \, \, \, = \, \, \, \int _{0}^\infty \, r f ( r ) J_1 ( k r ) dr, \end{aligned}$$

(A1)

while the z-component of the momentum balance employs the zeroth order Hankel transform (Wolf 1985, Eqs. 7 and 8). Using the transform on the governing equations naturally generates both the 0th and 1st order Hankel transforms, as, for example,

$$\begin{aligned} {\mathcal {H}}_1 \{ d f ( r ) / d r; k \} \, \, \, = \, \, \, - k {\mathcal {H}}_0 \{ f ( r ); k \}, \end{aligned}$$

(A2)

where the order is defined by the orders of the Bessel function $J_\upsilon$ with $\upsilon = 0, 1$ for zeroth and 1st order Hankel transforms, respectively. The convolution at the boundary is over the disk load geometry. This is straightforwardly represented as

$$\begin{aligned} \hat{q}_0 ( k )= & {} \int _{0}^\infty \, r H ( \alpha _D - r ) J_0 ( k r ) dr \, \, \, \nonumber \\= & {} \frac{\alpha _D}{k} J_1 ( k \alpha _D ), \end{aligned}$$

(A3)

where $\alpha _D$ is the radius of the disk and k plays the role of a wave number having units $\mathrm {m^{-1}}$ and $H( \alpha _D - r )$ is the Heaviside function. The disk places a stress on the surface,

$$\begin{aligned} q (r, t ) \, \, = \, \, \rho _{D} g h \vartheta _{o} ( t ) H( \alpha _D - r ), \end{aligned}$$

(A4)

where h and $\rho _{D}$ are the disk height and density, respectively. Although four Hankel transformed dependent variables enter the fully transformed equations, only the boundary condition for the vertical component of displacement need be considered to generate time-dependent solutions for the vertical displacement history w(r, t) after inverting the Laplace-Hankel transform (e.g., Wolf 1985; Klemann et al. 2003; Ivins et al. 2022). This boundary condition for the stress component ${\tilde{\sigma }}_{zz}$ (see Fig. 4) allows the transformed vertical displacement to take the form

$$\begin{aligned} \tilde{\hat{w}}_0 ( k, s, z = 0 )= & {} \tilde{\hat{q}}_0 ( k, s ) \frac{ 1 \, - \, s {\tilde{\nu }} ( s ) }{ s [ {\tilde{\mu }} ( s ) k + \rho g {\tilde{\nu }} ( s ) ] } \, \, \, \nonumber \\= & {} \tilde{\hat{q}}_0 ( k, s ) \hat{{\tilde{T}}}_{[\nu ]}^{(e v)} ( k, s ). \end{aligned}$$

(A5)

In Eq. (A5), $\hat{{\tilde{T}}}_{[\nu ]}^{(e v)} ( k, s )$ is the viscoelastic transfer function, with ${\tilde{\nu }}$ and ${\tilde{\mu }}$ the Laplace transformed Poisson’s ratio and shear modulus, respectively.

Here, the $\, \tilde{\hat{ }} \,$ symbol on dependent variable w and boundary parameter q indicate that the parameter is considered in the Laplace-Hankel transform space. We employ the definition $s {\tilde{\mu }}$ and $s {\tilde{\nu }}$ in place of $^{\star }{\tilde{\mu }}$ and $^{\star }{\tilde{\nu }}$ that are used by Wolf (1985) (compare Ivins et al. (2022), for example, with Eq. (23b) of Wolf (1985)). For a compressible Maxwell material, the second term on the right-hand-side of Eq. (A5), $\hat{{\tilde{T}}}_M^{(e v)} ( k, s )$ may be written as

$$\begin{aligned} \hat{{\tilde{T}}}_M^{(e v)} ( k, s ) = \frac{( 1 \, + \, s ) [ 3 ( 1 \, + \, s ) \kappa \, + \, 4 s \mu ]}{2 k \mu s [ 3 ( 1 \, + \, s ) \kappa \, + \, s \mu ] \, + \, \rho g [ ( 1 \, + \, s ) \kappa \, - \, 2 s \mu ] }.\nonumber \\ \end{aligned}$$

(A6)

For each of the rheological assumptions, we take the Laplace transform of $\vartheta _{o} ( t )$ to be ${\tilde{\vartheta }}_{o} ( s )$. For Maxwell viscoelasticity, we assume the transfer function, as derived in Wolf (1985) and Ivins et al. (2022), $\hat{{\tilde{T}}}_M^{(e v)} ( k, s )$. The analytic convolution into the time-domain is possible in the compressible case. The resulting expression yields efficient computations, but the expression is unwieldy to write here. By assuming incompressibility, the Laplace transform inversion of the convolved load function and viscoelastic transfer function is easily expressed as follows:

$$\begin{aligned}{} & {} {\mathcal {L}}^{-1} \{ {\tilde{\vartheta }}_{o} ( s ) \, \hat{{\tilde{T}}}_M^{(ev)} (k, s), \, s, \, t \}\nonumber \\{} & {} \quad = \left\{ - 2 k \mu (2 k \mu + \rho g ) [ 1 - e^{ - \rho g t / (2 k \mu + \rho g )} ] \right. \nonumber \\{} & {} \qquad + \frac{4 \rho g k \mu }{\pi R} \, Sin \left( \frac{ \pi R t}{2} \right) \, \nonumber \\{} & {} \qquad + \rho g \, \left[ 2 k \mu \, + \, \rho g \left( 1 + \frac{4}{\pi ^2 R^2} \right) \right] \, \nonumber \\{} & {} \qquad \left. \left[ - 1 \, + \, Cos \left( \frac{\pi R t}{2} \right) \right] \, \right\} \, / \nonumber \\{} & {} \qquad \rho g \, \left\{ 4 k \mu ( k \mu \, + \, \rho g ) \, + \, (\rho g)^2 \left[ 1 \, + \, \frac{4}{\pi ^2 R^2} \right] \right\} . \end{aligned}$$

(A7)

when $t \le 1/ R$. Note that the independent variable t is dimensionless (scaled to the Maxwell time), while the overall expression has units $\mathrm {m^3 / N}$. A noteworthy property of Eq. (A7) is that there is a term that is exactly in phase with the time-dependence of the load function. In order to isolate the rheological source for this term, we may take the limit as $\mu \, \rightarrow \, 0$ as this will isolate the viscous part of the solution. For $t \le 1/R$, this limit yields,

$$\begin{aligned} \lim _{\mu \rightarrow 0} {\mathcal {L}}^{-1} \{ {\tilde{\vartheta }}_{o} ( s ) \,, \, s, \, t \} \, \, = \, \, \frac{1}{\rho g} \, \left[ \, - 1 \, + \, Cos \left( \frac{\pi R t}{2} \right) \, \right] .\nonumber \\ \end{aligned}$$

(A8)

To simplify notation when we introduce the EBM rheology, we let

(A9)

where we employ the subscripts $[ \nu ]$ having either M, EBM or E for referring to the rheology-specific inversions. More generally, assuming $t \ge 0$,

(A10)

with the focus of this paper on the shorter timescale solutions . The time character of is essentially identical to the solutions for ramp forcing discussed in Ivins et al. (2022). The solution is simply the right-hand-side of Eq. (A7). We give the incompressible form of the convolution with the load for $t > 1/R$ here

(A11)

Clearly, the solution has a term that descends linearly in t and one that has slow exponential decay with advancing time t, as we note $R t > 1$. The constant term, , is

(A12)

with

(A13)

(A14)

(A15)

Clearly, these analytic expressions are complicated in the case of incompressibility. They are much more so in the case of compressibility, and we provide guidance for obtaining these at the end of this paper.

It is also of interest, quite naturally, to examine the purely elastic convolution. The elastic solution is recovered by taking the limit $s \rightarrow \infty$ of Eq. (A6) and then convolving with the Laplace transformed load function, $\tilde{\vartheta }( s )$. The two elastic solutions are as follows:

(A16)

and

(A17)

with

(A18)

The reader can readily observe that the solution has continuity at $t = 1 / R$

(A19)

similar to the continuity of the viscoelastic solutions. In the elastic case .

Appendix B Time evolution in wave number domain

Fig. 11

Inverse Laplace transformed Maxwellian response versus $t^{\prime }$ for $t^{\prime } \le 1 / R$ with $R = 10$. In frame (a), the ordinates on the left and right are for the wave number (k) cases in blue and green, respectively. Solid and dashed curves are for compressible and incompressible solutions, respectively. In frame (b), the two wave number cases (solid blue and green curves) are shown again for the compressible cases, with an additional wave number case shown by the red curve. The case of a Newtonian viscous fluid solution, for comparison, is shown in the black-dashed curve

Fig. 12

Inverse Laplace transformed response versus $t^{\prime }$ for $t^{\prime } \le 1 / R$ with $R = 2$. The same wave number cases shown in solid green ($k = 7.85 \times 10^{-7} m^{-1}$) and blue ($k = 7.85 \times 10^{-6} m^{-1}$) of Fig. 11a. The larger amplitudes correspond to the smaller wave number. Note the slight differences between Maxwell viscoelastic and elastic response evolution

Fig. 13

Inverse Laplace transformed response versus $t^{\prime }$ for $t^{\prime } \le 1 / R$ with $R = 10$ for three different ratios $\tau _H / \tau _M$ (solid lines) and corresponding Maxwell (red-dashed line) and elastic (black-dashed line). Frames (a) and (b) are for the same low and high wave numbers (k), respectively, as in Figure (11). For EBM, $\Delta$ and $\alpha$ are held fixed at 1.2 and 0.26, respectively. For this example, $\tau _L$, $\tau _M$ are also fixed to 57 min. and 37.2 years, respectively

In this Appendix, we examine the inverse Laplace transform solutions for elastic, Maxwell and EBM assumptions. Inversion provides the function , a k-dependent proxy for vertical motion as a function of time, thus allowing extensive examination the parameter-dependencies of the EBM rheology in contrast to elasticity and Maxwell counterparts in during the early phases of the accelerating surface unloading. In the figures that follow in the Appendix, we employ the $\prime$ affixed to s and t to remind the readers that these are non-dimensional variables.

The equivalent expressions for Laplace transform inversion cannot be derived analytically and we resort to numerical methods as discussed in Ivins et al. (2022). Comparisons of the k-dependent , and reveal a great deal about the relative robustness of uplift in response to the unloading function $\vartheta _{o} ( t )$ (i.e., Eq. 8). Physically, is a response to the negative downward stress that follows the removal of the surface load.

In Fig. 11, the proxy function is shown for the case $R \, = \, \textrm{10}$ from 0 to $t^{\prime }$ = 0.1, spanning the time over which the load is accelerating in load loss. In this case, only the Maxwell solutions are shown. Figure 11a presents a comparison of incompressible versus compressible solutions are shown at two wave numbers k. The two values of k correspond to long-wavelength and short-wavelength time evolutions. These two wave-lengths correspond to roughly a spherical harmonic degree $n = 5$ and $n = 300$, respectively. Also shown is the purely viscous solution which has no k-dependency. In this Appendix, all computations are shown assuming a half-space density, $\rho$ of $3380 \, \mathrm {kg / m^3}$, and elastic moduli, $\mu \, = \, 67.0$ GPa (unrelaxed) and $\kappa \, = \, 124.44$ GPa, with gravity fixed to $g \, = \, \mathrm {9.8 \, m/sec^2}$. The comparison to the viscous solution reminds us that a property of elasticity provides strength to the material, or effectively, a resistance to deformation. Figure 12 shows a comparison on elastic and Maxwell viscoelastic solutions. Figure 11, and subsequent figures in this Appendix only show the solutions for the compressible cases of the three rheologies under study.

The proxy evolution for EBM is shown in Fig. 13 for the same two values of k for compressible cases as shown in Fig. 11a and b. Here, the results for varying ratios of high cutoff period ($\tau _H$) to Maxwell time ($\tau _M$) for fixed low cutoff period ($\tau _L$) are shown with the two wave number cases separated into frames (a and b). Interestingly, the shorter ratio, $\tau _H/\tau _M$, reveals the more robust response. In the case of EBM rheology, all of the relaxation processes that occur over the time interval $\tau _L$ to $\tau _H$ will be operative over the entire time of evolution, 0 to $\tau _{dur}$, or in non-dimensional time, $t^{\prime }$, 0 to 1/R.

Fig. 14

Inverse Laplace transformed response versus $t^{\prime }$ for $t^{\prime } \le 1 / R$ with $R = 10$. As in Figure (13), also shown are the Maxwell and elastic cases (red- and black-dashed lines, respectively). Here, the anelastic relaxation strength, $\Delta$, is varied within values typical of laboratory experiments. The ratio of high cutoff time to Maxwell time, $\tau _H / \tau _M$, is held fixed to 0.5. Also shown are cases $[\nu ] = EBM$ (solid lines), M (red-dashed) and E (black-dashed), using differing values for $\Delta$, holding the ratio $\tau _H / \tau _M$ as indicated. Other EBM parameters, $\alpha$, $\tau _L$, and $\tau _M$, fixed to 0.26, 57 min., 37.2 years, respectively

Figures 14 and 15 reveal the $\Delta$-dependency of the solutions with long-wavelength k value in time for fixed $\alpha$, $\tau _L$, $\tau _L$ and $\tau _M$, again in the case R = 10. The important dependency of uplift solutions on the value of $\Delta$ is discussed at length in the space-time domain in the main text.

Fig. 15

Inverse Laplace transformed response versus $t^{\prime }$ for $t^{\prime } \le 1 / R$ with $R = 3$ and EBM parameters consistent with shear modulus dispersion observed in the GPS analysis of the $M_2$ ocean load tide by Bos et al. (2015) and Wang et al. (2020). EBM parameters are as indicated, with $\tau _M$ selected to be appropriate to continental Europe and the East China Sea. Elastic and Maxwell responses are also shown. The wave number is $k = 7.85 \times 10^{-7} m^{-1}$

Appendix C Potential Role of Power-law Exponent $\alpha$

The viscoelastic transfer function for EBM, $\hat{{\tilde{T}}}_{EBM}^{(ev)}$, is used for the convolution with the load function to obtain the z-component of the surface displacement as given in Eq. (13). In the experiments preformed, we have not explored the power-law- dependency of the solution space in any detail ($\alpha$-dependency, see Eq. 4). Experience with developing solutions for Boussinesq’s problem for EBM suggests that this dependency is weak (e.g., Ivins et al. 2020).

Fig. 16

The EBM transfer function for the half-space and dependence on $\alpha$ across 8 decades of $1/s^{\prime }$. (Units are $10^{-5}$ $\textrm{m}^3$ / Newton). Solid curves: $\alpha = 0.5$, dashed: $\alpha = 0.25$, dotted: $\alpha = 0.05$. Two Hankel wave numbers, k, are tested. The transfer function samples small $1/s^{\prime }$ at the elastic limit and and viscous limit to the extreme right on the abscissa. The high and low wave numbers sample shallow and deep portions of the half-space, respectively. All other parameters are as in the model for Fig. 5

However, closer examination of $\hat{{\tilde{T}}}_{EBM}^{(ev)}$ as a function of wave number, k, suggests some sensitivity to $\alpha$ when varied by an order of magnitude ($0.05 \le \alpha \le 0.5$). We show this dependency in Fig. 16. This dependency appears to be at about a 20% level at high wave number and should be examined in realistic models of GIA that have both upper mantle and lithospheric structure. Note that the dependency shown in Fig. 16 is independent of the load convolution. Also, note that at the lower wave number the curves are essentially identical, indicating a lack of sensitivity to $\alpha$. (We employ $s^{\prime }$ to remind the reader that it is scaled by the Maxwell time). The effect of lowering $\alpha$ below $\frac{1}{2}$ is to increase the amplitude of the transfer function in a part of the $1 / s^{\prime }$ space that corresponds to short wavelengths and short timescales. As one might suspect, smaller $\alpha$ should enhance the amplitude of short-timescale deformation response.

For the decadal timescale load modeling considered in this paper, we assume a value of $\tau _L$ that is always small, and when $\tau _L / \tau _H \, \ll \, 10^{-3}$ sensitivity to this parameter is negligible.