Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: an application to Glacial Isostatic Adjustment in Antarctica and Canada

SUMMARY

We develop a 3-D finite-element model to study the viscoelastic response of a compressible Earth to surface loads. The effects of centre of mass motion, polar wander feedback, and self-consistent ocean loading are implemented. To assess the model's accuracy, we benchmark the numerical results against a semi-analytic solution for spherically symmetric structure. We force our model with the ICE-5G global ice loading history to study the effects of laterally varying viscosity structure on several glacial isostatic adjustment (GIA) observables, including relative sea-level (RSL) measurements in Canada, and present-day time-variable gravity and uplift rates in Antarctica. Canadian RSL observations have been used to determine the Earth's globally averaged viscosity profile. Antarctic GPS uplift rates have been used to constrain Antarctic GIA models. And GIA time-variable gravity and uplift signals are error sources for GRACE and altimeter estimates of present-day Antarctic ice mass loss, and must be modelled and removed from those estimates. Computing GIA results for a 3-D viscosity profile derived from a realistic seismic tomography model, and comparing with results computed for 1-D averages of that 3-D profile, we conclude that: (1) a GIA viscosity model based on Canadian relative sea-level data is more likely to represent a Canadian average than a true global average; (2) the effects of 3-D viscosity structure on GRACE estimates of present-day Antarctic mass loss are probably smaller than the difference between GIA models based on different Antarctic deglaciation histories and (3) the effects of 3-D viscosity structure on Antarctic GPS observations of present-day uplift rate can be significant, and can complicate efforts to use GPS observations to constrain 1-D GIA models.

Time variable gravity, Global change from geodesy, Dynamics of lithosphere and mantle, Antarctica

INTRODUCTION

Observations of glacial isostatic adjustment (GIA), the viscoelastic relaxation of the Earth induced by deglaciation following the last glacial maximum, have provided valuable constraints on late Pleistocene ice history and the internal viscoelastic structure of the solid Earth. Arguably, the most useful GIA observations are those provided by relative sea-level (RSL) measurements (e.g. Peltier 1998; Mitrovica & Forte 2004), particularly from near the locations of the Laurentide ice sheet in northern Canada. However, geodetic observations have also proven useful for this purpose, including those that monitor secular trends in the Earth's time-variable gravity field (e.g. from the GRACE satellite gravity mission (Tapley et al.2004)) and in surface deformation (e.g. GPS crustal motion measurements). For instance, since RSL constraints are scarce in Antarctica, GPS observations of present-day bedrock uplift rates are proving to be especially useful for assessing GIA Antarctic modelling results (see, e.g. Bevis et al.2009; Argus et al.2011; Thomas et al.2011).

The GIA signal is also a significant source of noise for other applications. For example, errors in GIA models due to errors in the assumed ice deglaciation history and mantle viscosity profile are generally assumed to be the largest source of uncertainty when using GRACE time-variable satellite gravity data to estimate present-day thinning rates of the Antarctic ice sheet (Chen et al.2006, 2008, 2009; Velicogna & Wahr 2006; Velicogna 2009). To estimate the GIA uncertainty in GRACE-derived Antarctic ice loss rate, the GIA contribution is often computed for different ice histories (e.g. ICE-5G from Peltier 2004; IJ05 from Ivins & James 2005; and W12 from Whitehouse et al.2012), and for different lower mantle and upper mantle viscosity profiles (Velicogna & Wahr 2006).

The vast majority of all these efforts employ GIA models that assume a 1-D (i.e. radially dependent) viscosity profile. And yet the Earth certainly possesses laterally inhomogeneous (3-D) structure; particularly in its viscosity, which is likely to be strongly temperature-dependent. It is thus natural to ask whether GIA results for a 3-D Earth can be adequately represented by the predictions from 1-D models. And, if so, what 1-D viscosity profile will best predict a set of GIA observations, such as RSL measurements in northern Canada? Similarly, if a 1-D viscosity profile is used to make GIA predictions for Antarctica, what errors might be introduced?

To answer these questions, comprehensive GIA models must be built to handle a realistic Earth structure. Over the last decade, various numerical schemes have been developed to address the GIA modelling problem for a 3-D spherical Earth (e.g. Latychev 2005 and Whitehouse et al.2006 for finite-volume modelling of a compressible Earth; Martinec 2000 and Tanaka et al.2011 for spectral finite-element modelling of a compressible Earth; Martinec 1999 for spectral, initial value approach for a compressible Earth; Wu & van der Wal 2003; Zhong et al.2003; Wu 2004; Paulson et al.2005 and Wang & Wu 2006 for finite-element modelling of an incompressible Earth). In this study, we adopt the methodology discussed in Zhong et al. (2003) and Paulson et al. (2005), and develop a 3-D finite-element model in which the effect of compressibility is included, the gravitational potential arising from internal density perturbations is properly treated, and an arbitrary 3-D viscosity structure can be processed. We develop this finite-element model primarily to study load-induced deformation of the Earth, such as that produced by the GIA process. However, with minor modifications, this model can also be used to study the body tide response of the Earth or any other planetary body (Zhong et al.2012).

This paper is structured into three parts. Firstly, we present the physical model for the viscoelastic response of a laterally inhomogeneous Earth. We discuss in detail the numerical methods we employ to solve the equation of motion, and to implement the effects of centre of mass motion, polar wander feedback and self-consistent ocean loading. Secondly, by computing the response of an spherically symmetric Earth to both a Heaviside loading history and the ICE-5G loading history, we produce a benchmark comparison between the finite-element result and a semi-analytic solution. We assess the performance of the numerical model and study the relationship between accuracy and spatial resolution. Using the ICE-5G benchmark results, we quantify the numerical errors associated with the finite-element model for different observation types. Finally, we apply the finite-element model to a plausible 3-D viscosity structure and discuss the possible errors that might be introduced into Antarctic mass loss and uplift rate estimates, if we approximate the 3-D viscosity profile with its 1-D average. We also consider the difference between 3-D and various 1-D model predictions of RSL observations in northern Canada, to assess whether viscosity models provided by those observations are likely to represent truly global averages, or are more apt to be regional representations.

PHYSICAL MODELS FOR VISCOELASTIC DEFORMATION

Governing equations

Our model of viscoelastic deformation assumes a compressible self-gravitational Earth. The Earth's mantle is treated as a Maxwell solid underlying an inviscid fluid core. For an spherically symmetric density distribution, the governing equations of mass and momentum conservation, along with Poisson's equation for the perturbation in gravity, can be written as (e.g. Tromp & Mitrovica 1999; Zhong et al.2003)

\begin{equation} \rho _1^E = - (\rho _0 u_i )_{,i} , \end{equation}

(1)

\begin{equation} \sigma _{ij,j} + \rho _0 \phi _{,i}^{} - \rho _1^E g_i - (\rho _0 gu_r )_{,i} = 0, \end{equation}

(2)

\begin{eqnarray} \nonumber\ \phi _{,ii}^{} = - 4\pi G\rho _1^E , \end{eqnarray}

(3)

where ρ^E₁ is the Eulerian density perturbation, u_i is the displacement, u_r is the radial component of the displacement, σ_ij is the stress tensor, ρ₀ and g are the unperturbed density and gravitational acceleration, Φ is the perturbation of the gravitational potential, G is the gravitational constant, the notation A,_i represents the derivative of the variable A with respect to x_i and the repeated index implies summation.

Boundary conditions

The boundary conditions for normal tractions at the surface and the core–mantle boundary (CMB) are given by (see eqs 52, 53 and 60 of Dahlen 1974)

\begin{equation} \sigma _{ij} n_j = - \Gamma (t,\theta ,\phi )n_i {\rm , at\,\, the \,\,outer \,\,surface,} \end{equation}

(4)

\begin{equation} \sigma _{ij} n_j = - \rho _c \phi n_i + \rho _c gn_i u_j n_j {\rm , at\,\, the\,\, CMB,} \end{equation}

(5)

where n_i represents the normal vector of the boundary surface, Γ(t, Θ, Φ) is the surface load and ρ_c is the density of the core.

Mechanical properties

The Earth's mantle is treated as a compressible Maxwell solid, and so the constitutive equation can be written as (Wu & Peltier 1982)

\begin{equation} {\mathop \sigma \limits^ \bullet} _{ij} + \frac{\mu }{\eta }\left(\sigma _{ij} - \frac{1}{3}\sigma _{kk} \delta _{ij} \right) = 2\mu {\mathop \varepsilon \limits^ \bullet} _{ij} + \lambda {\mathop \varepsilon \limits^ \bullet} _{kk} \delta _{ij} , \end{equation}

(6)

where η is the viscosity and λ and μ are the Lamé parameters. The strain | $ \mathop \varepsilon _{ij}$| is related to the deformation by

\begin{equation} \varepsilon_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i} ). \end{equation}

(7)

Our numerical model can process a layered or spherically symmetric density distribution, along with a fully 3-D structure for the mantle viscosity and Lamé parameters.

NUMERICAL ANALYSIS OF VISCOELASTIC DEFORMATION

To solve the viscoelastic deformation problem for a compressible mantle, we follow a strategy similar to that discussed in Zhong et al. (2003) for incompressible finite-element modelling. Differences from Zhong et al. (2003) arise because u_{i, i} is in general non-zero for a compressible medium, which means that the constitutive equation does not now include pressure terms. And the Eulerian density perturbation is now non-zero, which complicates our gravitational potential solution and the momentum equation. For the remainder of this section, we address these issues, build up a matrix equation and discuss how to solve the resulting equations numerically.

Constitutive equation

We use a formulation that employs incremental displacements and strains, and discretize the constitutive equation (Martinec 1999; Zhong et al.2003). Let uⁿ_i and u^{n + 1}_i be the displacements at time tⁿ and t^{n + 1} = tⁿ + Δt. An incremental displacement vⁿ_i and an incremental strain Δɛⁿ_ij can be defined as

\begin{equation} v_i^n = u_i^{n + 1} - u_i^n , \end{equation}

(8)

\begin{eqnarray} \nonumber\\ \Delta \varepsilon _{ij}^n = \frac{1}{2}(v_{i,j}^n + v_{j,i}^n ). \end{eqnarray}

(9)

Integrating eq. (6) from tⁿ to t^{n + 1} with the second-order trapezoid rule, we get

\begin{equation} \sigma _{ij}^{n + 1} = \tau _{ij}^{{\rm pre}} + 2\tilde \mu \Delta \varepsilon _{ij}^n + \tilde \lambda \Delta \varepsilon _{kk}^n \delta _{ij} , \end{equation}

(10)

where

\begin{eqnarray} \tau _{ij}^{{\rm pre}} = \displaystyle\frac{{1 - \frac{\mu }{\eta }\frac{{\Delta t}}{2}}}{{1 + \frac{\mu }{\eta }\frac{{\Delta t}}{2}}}\sigma _{ij}^n + \frac{{\frac{\mu }{{3\eta }}\Delta t}}{{1 + \frac{\mu }{\eta }\frac{{\Delta t}}{2}}}\sigma _{kk}^n \delta _{ij} , \end{eqnarray}

(11)

\begin{equation} \tilde \mu = \frac{\mu }{{1 + \frac{\mu }{\eta }\frac{{\Delta t}}{2}}}, \end{equation}

(12)

\begin{equation} \tilde \lambda = \frac{{\lambda + (\lambda + \frac{2}{3}\mu )\frac{\mu }{\eta }\frac{{\Delta t}}{2}}}{{1 + \frac{\mu }{\eta }\frac{{\Delta t}}{2}}}, \end{equation}

(13)

and where τ^pre_ij is a parameter related to the total stress, σⁿ_ij, at the previous time step. | $ \tilde \mu$| and | $ \tilde \lambda$| depend on the Lamé parameters, the viscosities and the time interval Δt. With Eqs (8)–(13), we can express the stress tensor in terms of the incremental displacements. One of the principal differences with Zhong et al.'s (2003) incompressible study is that here, we do not have pressure terms in eq. (10), and the only unknowns are the incremental displacements.

Gravitational potential

To find the gravitational potential, we assume a layered density structure and write the density as

\begin{equation} \rho _0 (r^{'} ) = \sum\limits_{i = 1}^{N} {\Delta \rho _{i} H(r_{i}^{'} - r^{'} )} , \end{equation}

(14)

where H is the Heaviside function (H(x) = 0 for x < 0, and H(x) = 1 for x ≥  0), and a density jump of Δρ_i appears at r′ = r′_i, for i = 1,…, N from the CMB to the surface of the Earth. Here, i denotes the ith layer and N is the number of density layers. Using eq. (14), we can rewrite eq. (1) as

\begin{equation} \rho _{1}^{E} = \sum\limits_{i = 1}^{N} {\Delta \rho _{i} u_{r} \delta (r_{i}^{'} - r^{'} ) - \rho _{0} u_{i,i} } . \end{equation}

(15)

In eq. (15), the first term on the right-hand side is the surface density perturbation at the boundary of each density layer, and the second term is a volumetric density perturbation due to non-zero u_{i, i} in the compressible mantle.

Expanded into the spherical harmonic domain, the general solution to eq. (3) can be written as

\begin{equation} \phi = \sum\limits_{l,m} {Y_l^m \phi _{lm} } , \end{equation}

(16)

where

\begin{equation} {\phi _{lm}} (\vec r) = \frac{4\pi G}{2l + 1}\int {\frac{r_ {<} ^{l} }{r_ {>} ^{l + 1} }Y_{l}^{m *}} {(\theta ^{\prime},\varphi ^{\prime})\rho _{1}^{E} (r^{\prime})d^{3} \vec r^{\prime}} , \end{equation}

(17)

and where r_< and r_> are the smaller and larger of r and r′, respectively. Y^m_l is given by

\begin{equation} Y_l^m (\theta ,\varphi ) = p_{lm} (\theta )(\cos (m\phi ) + i\sin (m\phi )), \end{equation}

(18)

where p_lm(Θ) is a normalized associated Legendre polynomial as described in Zhong et al. (2008). Later, we may reference to the cosine and sine part of Y^m_l, which are computed through the linear combination of Y^m_l and Y^{− m}_l.

With ρ^E₁ given in eq. (15), and taking into account the load induced surface mass density Γ/g, we can express eq. (17) in terms of surface and volumetric contributions Φ_s and Φ_v (see Wu 2004 for a similar representation of eq. 19a),

\begin{eqnarray} &&\!\!\!\!\phi _s = \frac{{4\pi G}}{{2l + 1}} \nonumber\\ &&\times\, \left[ {\sum\limits_{i = 1}^j {\frac{{r_i^{l + 2} }}{{r^{l + 1} }}} \Delta \rho _i u_{rlm} (r_i ) + \sum\limits_{i = j + 1}^N {\frac{{r_{}^l }}{{r_i^{l - 1} }}} \Delta \rho _i u_{rlm} (r_i ) + \frac{{r^l }}{{r_N^{l - 1} }}\frac{{\Gamma _{lm} }}{g}} \right], \nonumber\\ \end{eqnarray}

(19a)

\begin{equation} \varphi _\nu = \frac{{4\pi G}}{{2l + 1}}\int {\frac{{r_ < ^1 }}{{r_ > ^{l + 1} }}Y_l^{m*} (\theta ',\phi ')\rho _0 (r')\nabla \cdot \vec ud^3 \mathord{\buildrel{$\rightharpoonup$} \over r} '} , \end{equation}

(19b)

where r_i denotes the radius of the ith boundary and r_N is the surface radius, u_rlm and Γ_lm are the spherical harmonic coefficients of the radial displacements and the surface loads, respectively, and without losing generality, we assume that the field point r satisfies r_j < r ≤ r_{j + 1}. It is notable that the surface contribution is caused by the radial displacements at density discontinuities, and the volumetric contribution is due to the compressibility. Using eq. (19a) and eq. (19b), we can compute the gravitational potential as long as we know the displacements.

Momentum equation

We follow the usual strategy in finite-element methods (Hughes 2000; Zhong et al.2003) and reformulate the momentum equation into its weak form. Eq. (2) can be written as

\begin{equation} (\sigma _{ij} + \rho _0 \phi \delta _{ij} - \rho _0 gu_r \delta _{ij} )_{,j} = \rho _{0,i} \phi + \rho _1^E g_i . \end{equation}

(20)

Multiplying eq. (20) by a weighting function w_i, integrating over the entire mantle, and substituting eqs (4), (5), (8), (9), (10) and (15) and the derivative of eq. (14) into the integral, leads to the weak form of the momentum equation

\begin{eqnarray} \begin{array}{lll} {\displaystyle\int {w_{i,j}} [\tilde \mu (v_{i,j} + v_{j,i} ) + \tilde \lambda v_{k,k} \delta _{ij} ]} dV\! -\!\! \displaystyle\int (w_{i,i} \rho _{0} gv_{r} + w_{r} \rho _{0} gv_{i,i} )dV \nonumber\\ \quad+ \displaystyle\sum\limits_l {\int {w_r \Delta \rho _l gv_r dS_l } } = { - \int {w_{i,j} \tau _{ij}^{\rm pre} dV}} \nonumber\\ \quad+ \displaystyle\int {(w_{i,i} \rho _0 gU_r + w_r \rho _0 gU_{i,i} )dV} \nonumber\\ \quad + \displaystyle\sum\limits_l \int w_r (\Delta \rho _l \phi - \Delta \rho _l gU_r )dS_l\nonumber\\ \quad - \displaystyle\int\limits_{{\rm Earth^{\prime}s}\,\,{\rm surface}}\!\!\!\!\! {w_r \Gamma _{{\rm ice + ocean}} dS - \int {w_{i,i} \rho _0 \phi dV,} }\quad\quad\quad\quad\quad\quad\quad(21)\\\end{array}\ , \end{eqnarray}

(21)

where v_i is the incremental displacement that we solve for, and U_i is the total displacement from the previous time step. In the summations, l denotes the lth density boundary.

Eq. (21) can be converted into a matrix equation, using standard finite-element methods (Hughes 2000). The second and third terms on the left-hand side of eq. (21) lead to terms in the stiffness matrix that depend on gravity, which is not usually included in non-geophysical finite-element formulations. The right-hand side of the equation contributes the forcing terms for the matrix equation. Most of the forcing terms enter either through the boundary conditions or through the stress and displacement resulting from the previous time step.

There are terms on the right-hand side of eq. (21) that depend on the gravitational potential Φ. When we solve eq. (21), we decompose Φ into Φ = Φ⁰ + ΔΦ(v_i), where Φ⁰ is the initial potential (total potential from the previous time step plus the potential induced by the load itself) and ΔΦ(v_i) is the incremental potential that depends on the unknown incremental displacements. We find ΔΦ(v_i) by using Eqs (19a) and (19b), replacing u_i in those equations by v_i.

Matrix equation and its solution

Matrix equation

To construct a matrix equation from eq. (21), we use brick elements, with displacements defined on the eight nodes at the corners, and with the elements arranged into the same finite-element grid as in Zhong et al. (2003). After including compressibility, significant changes appear in the stiffness matrix and the forcing terms in the momentum equation, compared with those in Zhong et al. (2003) for incompressible media. Following the procedure described in Appendix A, we can convert eq. (21) into

\begin{equation} KV = F_0 + F(\Delta \phi (V)). \end{equation}

(22)

where K is the total stiffness matrix, V is the incremental displacement vector containing | $ \vec v$| at all the nodes and F₀ is the force vector that depends on the surface load Γ, the pre-stresses τ^pre_ij, the displacements from the previous time step U_i and the initial gravitational potential Φ⁰ (i.e. the total gravitational potential at the previous time step plus the gravitational potential induced by the incremental load itself). F is the force vector that depends on the incremental gravitational potential ΔΦ(V), which, in turn, depends on the incremental displacements V.

Solution to the matrix equation

Eq. (22) needs to be solved iteratively, because the incremental gravitational potential ΔΦ(V) depends on the unknown displacements. For a given time step, an initial guess V₀ (chosen as the solution from the previous time step) is assigned to V, and the incremental gravitational potential is computed as ΔΦ(V₀). So the force vector can be obtained as F₀ + F(ΔΦ(V₀)). Eq. (22) is solved with the guess force vector, and the result V₁ is assigned to V. This process is repeated until we obtain a convergent solution for V. For later use, we call this process the self-gravity iteration (see also Zhong et al.2003; Wu 2004).

SOLUTION METHODS FOR DEGREE 1 DEFORMATION, POLAR WANDER FEEDBACK AND OCEAN LOADING

Our compressible model includes the effects of centre of mass motion, polar wander feedback, and self-consistent ocean loading, similar to those in Paulson et al. (2005).

Degree 1 deformation

The CM frame is the coordinate system with origin at the centre of mass of the Earth-plus-load system. Different from the study for incompressible media in Paulson et al. (2005) where the CM is only affected by the load and displacements at the surface and CMB, the CM in the current model formulation for compressible media is also affected by displacements and compressibility within the mantle. We compute deformation in the CM frame, by implementing the following two steps.

1. Centre of mass change induced by the load itself. Each time a new incremental load is applied to the Earth's surface, and the centre of mass position changes. We transform our old frame (the CM frame for the previous time step) to the new CM frame, by fixing the finite-element grids and changing the values of the physical quantities defined on those grids. The physical quantities that concern us are the incremental surface loads at the current time step, the total displacements and the stress tensors obtained from the previous time step. Suppose we find the centre of mass change induced by the incremental load is | $ \vec r_{{\rm cm}}$|⁠. We then shift the total displacement field by | $ - \vec r_{{\rm cm}}$|⁠. Correspondingly, we add an additional surface mass of − Δρr_cmcos ϑ at each boundary, where Δρ is the density jump at that boundary and ϑ is the angular distance between | $ \vec r_{{\rm cm}}$| and the position vectors at the boundary. Now the loads are identical to those that would be experienced in the new CM frame, and the degree 1 term of the gravitational potential induced by the new loads vanishes at the Earth's surface. The stress tensors are invariant under the degree 1 translation, and so we do not modify them.

2. Centre of mass change induced by the deformation. Using the total displacements, the loads, and the stress tensors, all written in the new CM frame, we solve for the incremental displacements for the current time step. In general, the incremental displacements will shift the centre of mass position, so we have a new CM frame again. We perform the degree 1 translation as described in step (1), except that we shift the incremental displacements rather than the total displacements. Adding the incremental displacements to the total displacements, we get the degree 1 deformation in the CM frame for the current time step.

Polar wander feedback

Deformation of the Earth and changes in surface loading can perturb the Earth's moment of inertia, and consequently the orientation of the Earth's rotation axis. Changes in Earth rotation, in turn, deform the Earth via the resulting perturbation of the centrifugal potential. To model this process, we write the Earth's angular velocity vector as Ω(m₁, m₂, 1 + m₃), where Ω is the unperturbed magnitude of the angular velocity and m_i are the dimensionless Cartesian components of the perturbation. We combine the equatorial components of the perturbation into the complex form, m_± = m₁ ± im₂. Similarly, for the perturbation to the Earth's inertia tensor, we use I_± = I₁₃ ± iI₂₃. For periods much longer than the Chandler Wobble, we have (Mitrovica et al.2005; Paulson et al.2005):

\begin{equation} m_ \pm = \frac{I_ \pm }{(C - A)_{\rm hyd} (1 + \delta )}, \end{equation}

(23)

where C and A are the unperturbed principal polar and equatorial moments of inertia, and (C − A)_hyd denotes the purely hydrostatic oblateness, the value of which is determined by the Earth's response to rotation, in the fluid limit. δ is the parameter that describes the non-hydrostatic oblateness and is set here to be 0.8 per cent, as described by Mitrovica et al. (2005). The perturbation in the rotation axis will induce a perturbation in the centrifugal potential

\begin{equation} \Delta \phi _c = \sqrt {\frac{{2\pi }}{{15}}} \Omega ^2 r_s^2 \left( {m_ - Y_2^1 - m_ + Y_2^{ - 1} } \right), \end{equation}

(24)

where r_s is the Earth's mean surface radius. ΔΦ_c will deform the Earth and change the Earth's rotation axis again. This feedback process is implemented through the following two steps:

1. When a new incremental load is applied to the Earth's surface at each time step, we compute the change in the Earth's moments of inertia, which are proportional to the Y¹₂ components of the load-induced gravitational potential. Using eq. (23), we find m_±. The perturbation in the centrifugal potential induced by the load itself can be obtained using these m_± in eq. (24). We denote this potential as ΔΦ⁰_c, and we add it to the initial gravitational potential Φ⁰. We are now ready to start the self-gravity iteration.

2. F₀ in eq. (22) is built up using Φ⁰ + ΔΦ⁰_c, Γ, U_i and τ^pre_ij. An initial guess V₀ is assigned to V. To take polar wander feedback into account, we compute the change in the moments of inertia induced by the deformationV₀. The m_± are then computed using eq. (23). With eq. (24), we get the deformation-induced perturbation of the centrifugal potential ΔΦ_c(V₀), and the force vector can be computed as F₀ + F(ΔΦ(V₀) + ΔΦ_c(V₀)). Eq. (22) is solved with the guess force vector and the result V₁ is assigned to V. This process is repeated until we obtain converged results for both V and m_±.

Ocean loading

Changes in ocean loading have two sources: increased volume from melted ice, and the response of the fluid ocean to changes in the topography and in the geoid. Using the same method discussed in Paulson et al. (2005), we include ocean loading via the sea-level equation

\begin{equation} L_0 (\theta ,\varphi ,t) = [N(\theta ,\varphi ,t) - U(\theta ,\varphi ,t) + c(t)]O(\theta ,\varphi ,t), \end{equation}

(25)

where L₀ is the change in height of the ocean load, N and U are the changes in the geoid and in the surface topography, c is a spatial constant needed to conserve mass and O is the ocean function (1 over the ocean, and 0 elsewhere). Integrating eq. (25) over the ocean surface, we have

\begin{equation} c(t) = \frac{1}{{A_0 (t)}}\left( { - \frac{{M_{{\rm ice}} (t)}}{{\rho _w }} - \int {(N - U)Od\Omega } } \right), \end{equation}

(26)

where A₀(t) is the area of the oceans at time t, M_ice is the change in ice mass, ρ_w is the water density and dΩ is the differential element of solid angle. In principle, A₀(t) depends on time because of three effects: (i) as the ice melts off regions that lie below sea level (e.g. Hudson Bay), those regions become ocean and add to the ocean area; (ii) as meltwater is added to the oceans, the area will increase as the water flows up over land in regions with sloping bathymetry adjacent to shore; (iii) as land uplifts (or subsides) adjacent to the coast, the adjacent ocean coverage decreases (or increases). In this study, we include only effect (i).

Combining eqs (25) and (26), we have

\begin{eqnarray} &&\!\!\! L_0 (\theta ,\varphi ,t) = - \displaystyle\frac{{M_{{\rm ice}} (t)}}{{\rho _w A_0 (t)}}O(\theta ,\varphi ,t) \nonumber\\ &&\!\!\!\!\!\!\quad+ \left( {N(\theta ,\varphi ,t) - U(\theta ,\varphi ,t) - \frac{1}{{A_o }}\int {(N - U)Od\Omega } } \right)O(\theta ,\varphi ,t). \nonumber\\ \end{eqnarray}

(27)

The first term on the right-hand side of eq. (27) represents the static ocean load, the value of which can be determined from the knowledge of the ice loading history and the ocean function. The second term represents the dynamic ocean load. In our finite-element model, these two terms are treated as follows:

• Static ocean load. Each time an incremental ice load is applied to the surface of the Earth, and the corresponding change in the static ocean load is computed according to the first term on the right-hand side of eq. (27). That change in ocean load is added to the ice load to build the complete F₀ term.

• Dynamic ocean load. The dynamic ocean load is implemented through the self-gravity iteration (see also Wu 2004; Paulson et al.2005). Within an arbitrary iteration, the incremental gravitational potential ΔΦ is computed based on V and the dynamic ocean load from the previous iteration using eq. (19a) and (19b). U is the surface value of V and N is proportional to the surface value of ΔΦ. We use the second term on the right-hand side of eq. (27) to compute the updated dynamic ocean load for this iteration. The contribution of this load is added to the surface terms of F. Eq. (22) is solved and a new displacement field is obtained. This process is repeated until we obtain convergence. For the ICE-5G ice history and VM2 viscosity profile, it takes six to eight self-gravity iterations to reach convergence.

RESULTS

Benchmarks

Using the finite-element model, we generate numerical solutions and benchmark them against a semi-analytic solution (Appendix B) for a spherically symmetric Earth model. To build the finite-element grids, we divide the Earth's mantle into 12 caps that have approximately equal size, and each cap is further divided into p cells in each of the horizontal direction and q cells in the radial direction. So, the total element number is 12×p×p×q (see, e.g. Zhong et al.2008). The Earth model we use is available online (Peliter's website), where the viscosity profile is based on VM2 (Peltier 2004), and the elastic parameters are derived from PREM (Dziewonski & Anderson 1981). Two types of loading history are used in the benchmark runs: (1) single harmonic loads with Heaviside loading history and (2) the realistic ice loading history ICE-5G (Peltier 2004). The details are shown in the next two subsections.

Single harmonic loads with Heaviside loading history

We apply a single harmonic load with Heaviside time dependence

\begin{equation} \sigma (\theta ,\varphi ,t) = Y_{l_0 }^{m_0 } (\theta ,\varphi )H(t), \end{equation}

(28)

where | $ Y_{l_0 }^{m_0 } (\theta ,\varphi )$| is the spherical harmonic function for degree l₀ and order m₀, and H(t) is the Heaviside function. For a spherically symmetric Earth, we expect to see the following two features in the Earth's response: (i) the response should have the same | $ Y_{l_0 }^{m_0 } (\theta ,\varphi )$| angular dependence as the load and (ii) the ratio of the response to the surface load should be independent of the order m₀ chosen for the load (though it should depend on the degree, l₀). To verify that the finite-element solution displays these features and matches the semi-analytic solution, we adopt two quantitative measures of success: the amplitude error ɛ_a, which measures the deviation of the finite-element solution from the semi-analytic solution for degree l₀ and order m₀; and the dispersion error ɛ_d, which measures the combined contributions of all harmonics other than (l₀, m₀). These errors are defined as

\begin{equation} \varepsilon _a = \frac{1}{T}\int\limits_0^T {\left| {\frac{{S_n (l_0 ,m_0 ,t) - S_g (l_0 ,m_0 ,t)}}{{S_g (l_0 ,m_0 ,t)}}} \right|dt} , \end{equation}

(29)

\begin{equation} \varepsilon _d = \sqrt {\frac{1}{T}\int\limits_0^T {\frac{{\sum\limits_{l \ne l_0 ,m \ne m_0 } {\left| {S_n (l,m,t)} \right|^2 } }}{{\sum\limits_{l,m} {\left| {S_n (l,m,t)} \right|^2 } }}dt} } , \end{equation}

(30)

where S_n and S_g are the response from the finite-element model and the semi-analytic solution, respectively. The | $ Y_{l_0 }^{m_0 } (\theta ,\varphi )$| used in eq. (28) has cos (m₀Φ) and sin (m₀Φ) parts. In our case, the cosine part of eq. (28) is taken as the input load, and correspondingly, the cosine coefficient of the surface uplift is used as S_n and S_g. T is the time duration for which we compute the solutions. Table 1 shows the amplitude and dispersion errors for the vertical surface displacement for benchmark calculations with loads at different harmonics and grid size, and Fig. 1 compares the Love numbers derived from the vertical surface displacement for the numerical and semi-analytic solutions.

Figure 1.

(a) Time-dependent Love number from semi-analytic solutions (dashed lines) and from the finite-element solution (solid lines) in response to single harmonic Heaviside loading. The elapsed time is normalized byτ₀ = 443 yr (a typical value of the Maxwell relaxation times for the viscosity model). The figure only includes data points with elapsed time ≤ 50, because the nearly flat viscous tails show consistent agreement between the two solutions. (b) Convergence of the degree 2 order 1 numerical results to the semi-analytic result with increasing spatial resolution in the finite-element model.

Table 1.

Amplitude and dispersion errors* for cases with Heaviside loading history.

Case	(l,m)	Grid	Amplitude error (per cent)	Dispersion error (per cent)
A1	(1,0)	12×48×48×48	0.0639	0.0508
A2	(1,1)	12×48×48×48	0.0726	0.0430
A3	(2,0)	12×48×48×48	0.2431	0.0691
A4	(2,2)	12×48×48×48	0.2644	0.0656
A5	(3,0)	12×48×48×48	0.3104	0.0994
A6	(3,3)	12×48×48×48	0.3427	0.0847
A7	(4,0)	12×48×48×48	0.4775	0.1468
A8	(2,1)	12×48×48×48	2.8895	0.1630
A9	(2,1)	12×64×64×48	1.1255	0.0936
A10	(2,1)	12×80×80×48	0.4703	0.0781

*The error results are computed for T = 300τ₀, where τ₀ is 443 yr, for both the finite-element model and the semi-analytic method. The incremental time step used in the finite-element model is set to be Δt = 0.2τ₀.

For cases A1–A7 with loading harmonics from degrees 1–4 on a 12×48×48×48 grid, the dispersion errors over 300 Maxwell times are less than 0.2 per cent, which means that the surface displacement field from the finite-element model has the same angular dependence as the surface load. The amplitude errors for these cases are smaller than 1.0 per cent. This shows a good match between the finite-element solutions and the semi-analytic solutions, which is evident in Fig. 1a. To examine the order dependence in the response, we include m ≠ 0 cases (A2, A4 and A6) in our calculations. Their errors are consistent with the corresponding m = 0 cases (A1, A3 and A5). Similar to Zhong et al. (2003), we find that the errors increase with increasing spherical harmonic degree (i.e. with decreasing spatial scales), because resolving shorter wavelengths requires a finer mesh.

As shown in cases A8–A10, for the degree 2 order 1 term, the error increases with time (Fig. 1b), and higher resolution is required to reach an accuracy of better than 1.0 per cent. This special feature is because of polar wander feedback, which, for a spherical Earth, affects only the degree 2 order 1 term, and can be understood with the aid of the analytic expression given in eqs (B9) and (B10) in Appendix B. At long periods, k^T₂ approaches k_f, and the denominators of eqs (B9) and (B10) are close to 0. Therefore, to obtain reliable results, k^T₂ needs to be resolved to high accuracy. Using the language of the finite-element model, this means we need to resolve the Earth's response induced by the perturbed centrifugal potential to high accuracy. As shown in Table 1, by increasing the resolution, we correspondingly decrease the amplitude error for the degree 2 order 1 response.

ICE-5G loading history

The ICE-5G loading history is a realistic ice model that describes the temporal and spatial distribution of ice on the Earth's surface during the last 122 000 yr. We apply the ICE-5G load to the same Earth model used in the last section. We use the 12×80×80×48

grids for our computation (Table 2). To compare with the semi-analytic solution, the surface displacement from our finite-element model is expanded into a set of spherical harmonic coefficients. And the benchmark results are obtained for l ≤ 32. For each harmonic, the amplitude errors are computed using eq. (29) either for the cosine coefficients or for the sine coefficients, whichever has the dominant amplitude. Fig. 2 shows the magnitudes of the vertical displacement at different spherical harmonic degrees and Table 2 shows the amplitude errors.

Figure 2.

Magnitudes of the surface displacement from the semi-analytic solution (dashed lines) and from the finite-element solution (solid lines) in response to the ICE-5G loading. The resolution in the numerical code is set to be 12 × 80 × 80 × 48. (a) The cosine coefficients for m = 0 terms. (b) The sine coefficients for m = 1 terms.

Table 2.

Amplitude errors* for different spherical harmonics for the ICE-5G loading history.

(l,m)	Coefficients	Amplitude error (per cent)
(1,0)	Cosine	0.6112
(2,0)	Cosine	0.2908
(7,0)	Cosine	1.108
(9,0)	Cosine	1.866
(12,0)	Cosine	4.7758
(2,1)	Sine	2.9156
(4,1)	Sine	0.4551
(5,1)	Sine	0.3752
(8,1)	Sine	2.3207
(10,1)	Sine	2.8064

* The error results are computed using the cosine coefficients for the m = 0 terms and the sine coefficients for the m = 1 terms. The time duration is T = 122 000yr. In the finite-element model, the incremental time stepΔt is 200yr from 122 000 yr ago to 17 000 yr ago, and is 25 yr from 17 000 yr ago to the present day.

For the longest wavelength terms, the amplitude errors are less than 1.0 per cent (except for the degree 2 order 1 term). At shorter wavelengths, for instance l = 10 or 12, larger errors can be observed. To further investigate these errors, we compute the following root mean square values in the spatial domain:

\begin{equation} {\rm rms}\_\delta u(\theta ,\varphi ) = \sqrt {\frac{{\int_0^T {(u_n (\theta ,\varphi ,t) - u_g (\theta ,\varphi ,t))^2 dt} }}{T}} , \end{equation}

(31)

\begin{equation} {\rm rms}\_u(\theta ,\varphi ) = \sqrt {\frac{{\int_0^T {u_g (\theta ,\varphi ,t)^2 dt} }}{T}} , \end{equation}

(32)

where u_n and u_gdenote the vertical surface displacement from the finite-element model and the semi-analytic solution, respectively. | $ {\rm rms}\_u$| is the RMS value of the surface uplift, and | $ {\rm rms}\_\delta u$| is the RMS value of the difference between the finite-element and semi-analytic solutions for the surface uplift. As shown in Fig. 3, the errors are localized to the major loading regions. If we include spherical harmonics up to degree 32, comparing the top left panel with the bottom left panel in Fig. 3, an error as large as 3.0 per cent can be observed in Canada. On the other hand, if we truncate both solutions to degree 8, comparing the top right panel and the bottom right panel in Fig. 3, we have an error of less than 1.0 per cent even in Canada. To understand the source of the relatively large error for the high-degree terms, we compute the degree amplitude for each l:

\begin{equation} a_l = \frac{1}{T}\int_0^T {\sqrt {\frac{1}{{2l + 1}}\sum\limits_{m = 0}^l {\left( {C(l,m,t)^2 + S(l,m,t)^2 } \right)} } \;{\kern 1pt} } dt, \end{equation}

(33)

where C and S are the cosine and sine coefficients for the vertical surface displacement. As shown in Fig. 4, the degree amplitude decreases rapidly with increasing l. Because the short-wavelength components have relatively small amplitudes, their high-amplitude errors do not translate to the same percentage error in Fig. 3. On the other hand, also due to their smaller amplitude, the short-wavelength components in the finite-element approach are more susceptible to leakage from the components with larger amplitudes. Our numerical experiments suggest that this error becomes nearly two times smaller when increasing the resolution from 12×48×48×48 to 12×80×80×48.

Figure 3.

The left two panels show the root mean square values of the surface uplift and the difference in surface uplift between the finite-element and semi-analytic solutions, for degree l ≤ 32; the right two panels show the rms value of the uplift and the difference in uplift for degree l ≤ 8.

Figure 4.

Degree amplitude for each degree l.

Similar to the Heaviside loading case, the degree 2 order 1 term has a larger error (2.9 per cent) than other large-scale terms (see Table 2). For the ICE-5G case, the forcing includes harmonics of all degrees and orders; and, in fact, the (2, 1) term has a generally smaller amplitude than other large-scale terms (see Fig. 2). This suggests that the relatively large (2, 1) error could be due to leakage into the (2, 1) term caused by the response to the larger amplitude harmonics, compounded by the difficulties introduced by polar wander feedback as described above for the Heaviside case.

We investigate this through the following numerical experiments. Since the calculation of the dynamic ocean load involves the integration of the Earth's response over the ocean, it couples responses at different spherical harmonics. We remove this physical coupling by turning off the sea-level iteration and considering only the static ocean load (we consider the static load alone, only for this (2,1) test; the inclusion of the dynamic ocean load is the default setup unless otherwise specified). We set up three cases with resolution of 12×48×48×48, 12×64×64×48 and 12×80×80×48, respectively (Table 2), and run the numerical code to examine its performance. As shown in Fig. 5a, the result converges to the semi-analytic solution with increasing resolution. This suggests that the relatively large (2, 1) error, as shown in Table 2, is probably a resolution issue, and could be improved if we had sufficient memory resources to compute on a finer grid. To determine the effects of leakage caused by the other harmonics in the forcing, we extract the (2, 1) component from the ICE-5G ice load and compute the response to this single harmonic load. Fig. 5b shows that the error is reduced significantly compared to the 12×80×80×48 full ICE-5G case, which does suggest that the large (2, 1) error is mainly due to leakage from components with larger amplitudes.

Figure 5.

Convergence of the degree 2 order 1 numerical results to the semi-analytic solution with (a) increasing resolution, and (b) using only the single harmonic (2,1) load extracted from the ICE-5G loading history. The figures show the elapsed time from 60 000 yr ago to 15 000 yr ago.

So far, we have discussed benchmark results only for surface uplift. By averaging the difference between the finite-element and semi-analytic results over the entire ice loading history (eq. 31), we obtain a modelling error of 3.0 per cent when including l ≤ 32 (less than 1.0 per cent for l ≤ 8) (see Fig. 3). This also provides a reasonable estimate of the modelling error for any GIA observation that depends mainly on surface uplift. For example, we thus expect a modelling error of about 3 per cent for our finite-element predictions (based on a 12×80×80×48 mesh) of historic RSL measurements, where RSL is defined as the geoid minus the surface topography.

Similar benchmark results can be used to quantify modelling errors related to present-day variability. For example, suppose we are interested in using time-variable gravity measurements from the GRACE satellite gravity mission, or the ice elevation measurements from a satellite altimeter mission, to determine ongoing mass changes of the Antarctic ice sheet. This requires that we model and remove GIA-induced secular trends from the total mass change. To model the GIA process on an spherically symmetric Earth, we can use either the semi-analytic method or the finite-element method, and the finite-element modelling error can be defined as the difference between the results from these two methods. For the GRACE

estimates, we expand our predicted present-day GIA geoid rates into spherical harmonics. The harmonic coefficients are processed using the same spatial averaging technique used by Velicogna & Wahr (2006) to determine Antarctic mass loss from GRACE data. The semi-analytic result shows that the GIA effect for ICE-5G contributes 145.40 Gton yr⁻¹ ice mass gain over Antarctica. Using the finite-element model, we obtain 147.74 Gton yr⁻¹. Our estimate of the finite-element modelling error for the GRACE GIA correction is obtained by differencing these two numbers: 2.34 Gton yr⁻¹ (an error of 1.6 per cent relative to the semi-analytic result). For the altimeter estimates, the present-day GIA uplift rates are integrated over the Antarctic ice sheet, and multiplied by the density of ice to obtain the apparent change in ice mass. The results from the semi-analytic model and the finite-element model are 25.02 and 24.46 Gton yr⁻¹, respectively. The implied finite-element modelling error is then 0.56 Gton yr⁻¹ (2.2 per cent) for the Antarctic altimeter estimate. These estimates are summarized in Table 3.

Table 3.

The impact of GIA on GRACE and altimeter estimates, of the present-day rate of increase in Antarctic mass (results in Gton yr⁻¹) for a spherically symmetric Earth and the ICE-5G deglaciation history and VM2 viscosity profile (Peltier 2004). Results computed using both the semi-analytic and finite-element models are shown. The difference is interpreted as the error in the finite-element model.

	Semi-analytic result	Finite-element result	Modelling error
GRACE	145.40	147.74	2.34
Altimeter	25.02	24.46	0.56

Another common geodetic data type used in GIA studies is GPS monitoring of crustal uplift rates. For example, in Antarctica, where RSL observations are scarce, GPS observations provide a means to assess and improve the accuracy of GIA models. In these studies, the elastic response to the contemporaneous ice mass change is modelled and removed from the GPS uplift rates, and the trend of the residuals is interpreted as an estimate of the GIA signal (see, e.g. Thomas et al.2011). This trend is compared with the appropriate GIA model result, to assess that model. To determine the computational accuracy of our finite-element predictions of Antarctic uplift rates, we compute those uplifts using ICE-5G and a spherical Earth, for both the semi-analytic and finite-element methods. The finite-element modelling error is taken to be the difference between the results of the two models. The present-day uplift rates are computed on a spatial grid and the difference between the two models is shown in Fig. 6. The modelling error increases with the amplitude of the uplift rate. The largest error we observe for the uplift rate is 0.98 mm yr⁻¹ (8.0 per cent). And the largest error observed at an existing GPS station is 0.73 mm yr⁻¹ (7.0 per cent), the details of which are shown in the next section.

Figure 6.

The left panel shows the present-day uplift rates in Antarctica from the semi-analytic method. The right panel shows the difference in the present-day uplift rates between the finite-element model and the semi-analytic model.

Three-dimensional viscosity structure, with implications for Canadian RSL observations and Antarctic geodetic measurements

Most GIA modelling efforts employ 1-D viscosity models. Viscoelastic Green functions (Love number solutions) are computed using the 1-D models and are convolved with an ice model to predict present-day observables, such as RSL observations, GRACE mass estimates, altimeter elevation estimates, etc. What errors might this introduce, given that there is almost certain to be lateral variability in the Earth's viscosity structure? The answer to this question depends both on the true viscosity structure of the Earth and on the 1-D model used to approximate the Earth.

For the former, we choose a single plausible 3-D viscosity structure, determined from a global seismic tomography model (S20RTS shear wave model of Ritsema 1999), using assumed values of the activation energy (Paulson et al.2005). We include an elastic lithosphere that has lateral variations in thickness derived from a compilation of elastic plate thicknesses from Watts (2001) (see Zhong et al.2003). For our 1-D approximation, we initially employ a 1-D model that best represents the average of the 3-D viscosities under Canada, where the most prominent deglaciation takes place and where most of the RSL observations are recorded. The rationale is that it is likely that a 1-D viscosity profile determined by fitting to real GIA observations is most likely a representative of the average viscosity structure beneath Canada. When GIA modellers use the VM2 (Peltier 2004) 1-D viscosity model, for instance, it is likely that they are using the Canadian average, rather than the global average, of the Earth's true 3-D structure. To build our 1-D model, at each depth, we average together the logarithms of the 3-D viscosities under Canada to get the 1-D value at that depth (Paulson et al.2005).

To help verify our hypothesis that a 1-D viscosity model like VM2 is likely to be a Canadian average, we use our 3-D and 1-D models to predict RSL observations for the last 8000 yr at four prominent RSL observing sites: Churchill, Cape Henrietta Maria, Ottawa Island and Ungava Peninsula. As a comparison, we also include results computed for a globally averaged 1-D viscosity profile. The 1-D viscosity models are plotted in Fig. 7, and the RSL results are shown in Fig. 8. As anticipated the Canadian average results are in good agreement with the 3-D results. The global average does not do as good a job. Our result is consistent with previous studies that considered other major loading areas (e.g. see Kaufmann & Wu 2002 for a 3-D inversion study for Fennoscandia).

Figure 7.

1-D viscosity structures.

Figure 8.

Relative sea-level estimates at Churchill (top left), Cape Henrietta Maria (top right), Ottawa Island (bottom left) and Ungava Peninsula (bottom right).

We next look at the possible impact of using a 1-D viscosity profile when interpreting geodetic observations of Antarctica. Considering that there are significant lateral variability in the viscosity and the lithosphere thickness beneath Antarctica (Fig. 9), we suspect that the effect of 3-D structure might be important for Antarctic observations (see Kaufmann et al.2005 for a 3-D study for Antarctica). We compute present-day GIA signals for the 3-D model and the 1-D models using the finite-element method and the semi-analytic method, respectively. The GRACE and altimeter Antarctic estimates are shown in Table 5 (the case IDs are summarized in Table 4). These numbers should not be used to correct real GRACE or altimeter measurements, because none of the 1-D averages of our assumed viscosity structure are in particularly close agreement with a real GIA-based viscosity profile, such as VM2. For the GRACE estimates, the ‘Canada’ result (i.e. the calculation using 1-D viscosity averaged for the Canadian region) is 11.1 Gton yr⁻¹ (5.5 per cent) larger than the 3-D result. For the altimeter estimates, the ‘Canada’ result is 0.6 Gton yr⁻¹ (1.8 per cent) larger. Because the finite-element modelling error given in the last section is 1.6 per cent for the GRACE estimates and 2.2 per cent for the altimeter estimates, we conclude that the use of the 1-D Canadian average viscosity model rather than the 3-D profile introduces an error of 5.5 ± 1.6 percent for the GRACE estimates. For the altimeter estimates, the 1-D model reproduces the 3-D result quite well, and the error we obtain (1.8 per cent) is not large enough to be significant. As a comparison, Table 5 also shows the result from the 1-D model using a globally averaged viscosity. The ‘Global’ result is 7.1 per cent larger than the 3-D result for the GRACE estimates and 8.0 per cent larger for the altimeter estimates, implying that the 1-D Canadian average provides a better match to the 3-D result.

Figure 9.

The viscosity structure at 372 km depth beneath Antarctica (left) and the lithosphere thickness for Antarctica (right).

Table 4.

The viscosity profiles considered for the 3-D comparisons. For each 1-D case, we find the lateral average of the logarithms of the 3-D viscosities, at each depth.

Case ID	Viscosity profile
3-D	Full 3-D viscosity profile
Canada	1D viscosity profile derived from the average under Canada.
Global	1D viscosity profile derived from the global average
Antarctica	1D viscosity profile derived from the average under Antarctica

Table 5.

Our ICE-5G Antarctic GRACE and altimeter estimates using different viscosity profiles. Note: these numbers should not be used to correct real GRACE or altimeter estimates (see text).

		3-D	Canada	Global	Antarctica
GRACE	Entire Ant	201.02	212.07	215.29	221.29
	West Ant	91.34	88.10	90.97	93.40
	East Ant	108.26	122.90	123.54	133.37
Altimeter	Entire Ant	33.10	33.69	35.73	36.06
	West Ant	13.75	14.01	15.23	15.24
	East Ant	19.35	19.68	20.50	20.82

The Antarctic GPS uplift rate results are shown in Table 6 and Fig. 10. Comparing the 3-D result and the 1-D ‘Canada’ result, differences of up to 0.5–1.6 mm yr⁻¹ (a relative error ranging from 7 to 60 per cent) in uplift rate can be observed at selected stations, which are significantly larger than the corresponding finite-element modelling errors at those stations. Although the ‘Canada’ case produces good agreement with the 3-D results for the altimeter estimates (an integration of the uplift rates), it does not provide a good prediction of the 3-D result at each GPS station. Similarly, for the ‘Global’ 1-D case, the 3-D/1-D differences are as large as 1.78 mm yr⁻¹ (46 per cent). This suggests that the effects of 3-D viscosity structure can be important for localized measurements.

Figure 10.

The Antarctic uplift rate results for: the 3-D case (top left), the ‘global’ case result (bottom left), the ‘Canada’ case result (top right) and the ‘Antarctica’ case result (bottom right). The GPS stations are labelled by name.

Table 6.

The ICE-5G Antarctic GPS uplift rate predictions: columns 1–3 specify the station information; columns 4–7 show the uplift rate results from four different cases; columns 8–10 present the difference between the three 1-D case results and the 3-D case result; the 11th column shows the modelling error. Note that these numbers should not be used to correct or interpret real GPS measurements (see text).

lon	lat	Station	3-D	Canada	Global	Antarctica	Canada –3-D	Global 3-D	Antarctica 3-D	Modelling error
162.57	−78.93	FTP1	2.75	2.26	2.10	2.19	0.50	0.65	0.57	0.12
163.19	−77.03	ROB1	2.99	1.32	1.53	1.46	1.67	1.46	1.52	0.12
164.1	−74.7	TNB1	2.65	1.05	1.88	1.66	1.61	0.77	1.00	0.12
166.67	−77.85	CRAR	2.77	1.63	1.42	1.48	1.14	1.35	1.28	0.12
204.98	−78.03	MBL1	5.06	6.16	6.68	6.68	−1.10	−1.62	−1.62	0.20
215.7	−76.32	MBL2	2.47	3.04	2.49	2.74	−0.57	−0.03	−0.27	0.11
218.13	−77.34	MBL3	3.86	5.47	5.64	5.75	−1.61	−1.78	−1.89	0.15
279.44	−80.04	W05A	8.93	7.94	8.93	8.82	0.98	−0.01	0.11	0.53
291.45	−85.61	W02B	7.38	6.82	6.75	6.89	0.55	0.62	0.49	0.37
302.1	−63.32	OHIG	2.64	1.50	1.99	1.80	1.15	0.65	0.84	0.12
296.97	−72.67	BREN	6.09	4.72	6.12	5.78	1.38	−0.03	0.31	0.26
306.8	−82.86	W04A	10.46	8.95	11.05	10.67	1.51	−0.59	−0.21	0.73

One more useful 1-D viscosity average is the average of the logarithms of the 3-D viscosities under Antarctica, referred to in the tables and figures as ‘Antarctica’. Given that the ‘Canada’ case does a good job reproducing 3-D results in Canada, it is natural to ask whether the 1-D Antarctic model provides a good match to the 3-D results in Antarctica. The GRACE and altimeter results are presented in the last column of Table 5, and the GPS uplift results are shown in Table 6 and Fig. 10. Compared to the 3-D results, the ‘Antarctica’ results are 10.0 per cent larger for the GRACE estimates and 9.0 per cent larger for the altimeter estimates. The error in GPS uplift rate estimates can be as large as 1.89 mm yr⁻¹ (49 per cent). So, we do not get a better match using the Antarctic average. We suspect that the reason is related to the complex viscosity structure under Antarctica. Unlike in the mantle beneath Canada, where viscosities are relatively uniform across a large region, the mantle under Antarctica exhibits considerable lateral variations in viscosity. As a result, the 1-D viscosity model averaged under Antarctica might not be expected to produce the best match to the 3-D results, suggesting that it is important to use 3-D models in Antarctica.

However, we should note that we cannot formulate this as a general conclusion, because we only consider one possible 3-D viscosity profile in this study. We cannot rule out the possibility that for the Earth's true viscosity structure, the use of an Antarctic average might better match the 3-D results.

DISCUSSION AND CONCLUSIONS

A 3-D finite-element model is developed to study the viscoelastic response of a compressible Earth to surface loads. The mantle in the model can have a layered density distribution, along with 3-D structure for the mantle viscosity and Lamé parameters. Though here, we consider 3-D inhomogeneities only in the viscosity and assume that the Lamé parameters are spherically symmetric. The effects of centre of mass motion, polar wander feedback, and self-consistent ocean loading are implemented. The numerical method is discussed in detail. Solutions for spherically symmetric structure are benchmarked against a semi-analytic solution for both a Heaviside loading history and the ICE-5G loading history. In general, accuracies of better than 1 per cent can be obtained in the surface uplift solution for a Heaviside loading history. For ICE-5G, we use surface uplift as a proxy for RSL, which is a key observation type when constructing ice deglaciation and mantle viscosity models, and we obtain a modelling error of better than 1 per cent for long-wavelength terms (l ≤ 8), and of about 3 per cent when including degrees up to 32. Our numerical experiments suggest that the error for the short-wavelength terms is mainly induced by leakage from terms with larger amplitudes, which can be reduced if we increase the resolution by decreasing the grid size. When the effects of polar wander are included, we find that for ICE-5G, the degree 2 order 1 solution shows a smaller amplitude and a larger error than other global-scale terms. The error is mainly induced by leakage from components with larger amplitudes, and it can be reduced significantly by increasing the resolution. Throughout these discussions, we quantify the numerical errors in the finite-element model by computing the differences between the finite-element results and the semi-analytic solutions. There can, of course, also be errors in the semi-analytic results. However, given that the two solutions agree with each other to high accuracy, it is unlikely that the uncertainties in the finite-element solutions could notably exceed the error estimates we have provided here.

The GIA signal is a source of noise for GRACE gravity and altimeter elevation estimates of present-day Antarctic ice mass loss. It is important to be able to accurately model and remove GIA effects from those observations and to assess the error in those GIA model results. Using the ICE-5G benchmark results, we investigate the numerical error in our Antarctic GIA estimates. When computing the GIA contribution in Antarctic mass trend using the finite-element model given a 12×80×80×48 resolution setup, we obtain a 1.6 per cent relative error for the GRACE results, and a 2.2 per cent relative error for ice altimeter results.

Antarctic GPS observations have been used to constrain Antarctic GIA models, and it is useful to know how large the GIA modelling errors are when comparing with GPS observations. Our ICE-5G benchmark results show that if we use the finite-element model to compute the GIA-induced present-day uplift rates at existing Antarctic GPS sites, we obtain a maximum numerical error of 7 per cent. This is a significantly larger relative error than is obtained for the altimeter estimates, despite the fact that both are based on uplift rate results. The reason that the GPS errors tend to be larger than the altimeter errors is that the altimeter error is an integrated value over the entire ice sheet, while the uplift rate at a specific GPS station is more sensitive to short-wavelength features. This issue can be mitigated in the near future if we can further refine both the spatial and the temporal resolution in our calculation, given increased capacity in parallel computing.

Since there are a large number of GIA studies that employ incompressible Earth models, it is useful to briefly discuss the effect of compressibility on the GIA estimates. We use our semi-analytic method to compute two sets of GIA estimates, one using the original PREM structure and the other using an incompressible version of PREM. We find that making the Earth model incompressible reduces the GIA-induced Antarctic GRACE mass gain estimates by roughly 2 per cent, and reduces the present-day Antarctic uplift rates by about 5 per cent.

We apply the finite-element model to a plausible 3-D viscosity structure. The 3-D viscosity is determined from a global seismic tomography model along with a reasonable model of elastic lithosphere thickness. Using the model results, we study the effects of 3-D viscosity structure on various GIA observables. RSL results from northern Canada show that the 3-D case is better predicted using the 1-D Canadian average than the global average. This suggests that a GIA viscosity model based on RSL data is more likely to represent a Canadian average than a true global average.

We investigate the error that might be introduced into GRACE and altimeter estimates of total Antarctic ice mass loss by approximating the 3-D viscosity structure using its 1-D Canadian average. The 1-D model produces a small difference compared to the 3-D results for both the GRACE and the altimeter estimates (5.5 and 1.8 per cent, respectively). For uplift rates computed at the GPS stations, the 3-D results and the results for the 1-D Canadian average show differences ranging from 7 to 60 per cent, which indicates that 3-D effects might be more important for localized measurements. We also compute the GIA estimates using the Antarctic average of the 3-D viscosity profile. The results suggest that the Antarctic average does not provide a better match to the 3-D case. We suspect that this may be due to the complex viscosity structure beneath Antarctica.

In this study, we explore the effects of only one 3-D viscosity profile. So, it is difficult to make specific, quantitative estimates of the likely effects of realistic 3-D structure. However, it is still useful to compare our results for the effects of 3-D structure based on this one plausible 3-D viscosity model, with the results from a recent Antarctic mass balance study (Shepherd et al. 2012) that are based on GIA models with different ice deglaciation histories. Shepherd et al. (2012) use two newly developed GIA models, W12a (Whitehouse et al.2012) and IJ05_R2 (Ivins et al.2013), to compute GIA corrections for GRACE Antarctic mass loss estimates. They find that those two corrections differ by roughly 20 per cent. This relative difference is much larger than the relative effects of 3-D structure we obtain here (5.5 per cent for GRACE estimates). Shepherd et al. (2012) also compare the W12a and IJ05_R2 corrections with corrections based on ICE-5G, and find absolute differences that are 5–10 times larger than the differences between the W12a and IJ05_R2 corrections. So, in general, we conclude that the effects of 3-D viscosity structure on GRACE estimates of present-day Antarctic mass loss are probably smaller than the differences between GIA models based on different Antarctic deglaciation histories. On the other hand, the effects of 3-D viscosity structure on Antarctic GPS observations of present-day uplift rate can be significant, and can complicate efforts to use GPS observations to constrain 1-D GIA models.

This work is partially supported by NASA grants NNX08AF02G and NNXI0AR66G, by NASA's ‘Making Earth Science Data Records for Use in Research Environments (MEaSUREs) Program’, and by NSF grant NSF-1114168, all to the University of Colorado. We thank Mark Tamisiea for helpful discussions, and Patrick Wu and an anonymous reviewer for their careful reviews.

APPENDIX A

The stiffness matrix

The first term on the left-hand side of eq. (21) represents the ordinary stiffness matrix that does not depend on gravity. It can be written as

\begin{equation} \int {w_{i,j} [\tilde \mu (v_{i,j} + v_{j,i} ) + \tilde \lambda v_{k,k} \delta _{ij} ]} dV = \int {\vec \varepsilon (\vec w)} ^T D\vec \varepsilon (\vec v)dV, \end{equation}

(A1)

where | $ \vec \varepsilon$| and D are given by

\begin{equation} \vec \varepsilon (\vec w) = \left[ {\begin{array}{*{20}c} {\varepsilon _{\theta \theta } } \\ {\varepsilon _{\varphi \varphi } } \\ {\varepsilon _{rr} } \\ {\varepsilon _{\theta \varphi } } \\ {\varepsilon _{\theta r} } \\ {\varepsilon _{r\varphi } } \\ \end{array}} \right] = \left[ {\begin{array}{*{20}c} {\displaystyle\frac{{w_{\theta ,\theta } }}{r} + \frac{{w_r }}{r}} \\ {\displaystyle\frac{{w_{\varphi ,\varphi } }}{{r\sin \theta }} + \frac{{w_\theta \cot \theta }}{r} + \frac{{w_r }}{r}} \\ {w_{r,r} } \\ { - \displaystyle\frac{{\cot \theta w_\varphi }}{r} + \frac{{w_{\varphi ,\theta } }}{r} + \frac{{w_{\theta ,\varphi } }}{{r\sin \theta }}} \\ {w_{\theta ,r} - \displaystyle\frac{{w_\theta }}{r} + \frac{{w_{r,\theta } }}{r}} \\ {w_{\varphi ,r} - \displaystyle\frac{{w_\varphi }}{r} + \frac{{w_{r,\varphi } }}{{r\sin \theta }}} \\ \end{array}} \right], \end{equation}

(A2)

\begin{equation} D = \left[ {\begin{array}{*{20}c} {\lambda + 2\mu } &\quad \lambda &\quad \lambda &\quad {} &\quad {} &\quad {} \\ \lambda &\quad {\lambda + 2\mu } &\quad \lambda &\quad {} &\quad 0 &\quad {} \\ \lambda &\quad \lambda &\quad {\lambda + 2\mu } &\quad {} &\quad {} &\quad {} \\ {} &\quad {} &\quad {} &\quad \mu &\quad {} &\quad {} \\ {} &\quad 0 &\quad {} &\quad {} &\quad \mu &\quad {} \\ {} &\quad {} &\quad {} &\quad {} &\quad {} &\quad \mu \\ \end{array}} \right], \end{equation}

(A3)

and | $ \vec \varepsilon$| can be represented by the nodal displacements

\begin{equation} \varepsilon _i (\vec w) = \sum\limits_{a,j} {B_{ij}^a w_j^a } , \end{equation}

(A4)

where a goes through all the nodes and i = 1, ⋅⋅⋅, 6, j = 1, 2, 3. w^a_j denotes the jth component of | $ \vec w$| at the ath node, and B^a is a 6 × 3 matrix defined at the ath node. The form of B^a depends on the coordinate system, which, in this case, is spherical coordinates, and B^a is found through the calculation of the nodal positions and the shape functions (see Zhong et al.2000). Eq. (A1) can be rewritten as

\begin{equation} \int {\vec \varepsilon (\vec w)} ^T D\vec \varepsilon (\vec \nu )dV = \sum\limits_{a,b} {(\vec w^a )} ^T K_{ab}^1 \vec \nu ^b , \end{equation}

where

\begin{equation} K_{ab}^1 = \int {(B^a )^T DB^b dV} . \end{equation}

(A5)

The volume integration over the Earth's entire mantle can be decomposed into sums of elemental contributions, which can be computed directly. The second term on the left-hand side of eq. (21) represents another symmetric matrix. The radial displacement at a given point can be expressed as the linear combination of the nodal displacements

\begin{equation} v_r (\vec x) = v_3 (\vec x) = \sum\limits_{b,j} {N_b (\vec x)C_{3j}^b (\vec x)v_j^b } , \end{equation}

(A6)

where we denote the radial component as the third component, and | $ N_b (\vec x)$| is the shape function at the bth node. C^b is a 3 × 3 matrix, whose form also depends on the coordinate system (Zhong et al.2000). Using eqs (A2) and (A4), we have

\begin{equation} \nabla \cdot \vec w = \sum\limits_{a,i} {(B_{1i}^a + B_{2i}^a + B_{3i}^a )w_i^a } . \end{equation}

(A7)

Using eqs (A6) and (A7), the second term on the left-hand side of eq. (21) can be written as

\begin{equation} (w_{i,i} \rho _0 g\nu _r + w_r \rho _0 g\nu _{i,i} )dV = \sum\limits_{a,b} {(\vec w^a )} ^T K_{ab}^2 \vec \nu ^b , \end{equation}

where

\begin{eqnarray} (K_{ab}^2 )_{ij} \!&=&\! \displaystyle\int {\left[ {(B_{1i}^a + B_{2i}^a + B_{3i}^a )N_b C_{3j}^b + (B_{1j}^b + B_{2j}^b + B_{3j}^b )N_a C_{3i}^a } \right]}\nonumber\\ &&\times\rho _0 gdV . \end{eqnarray}

(A8)

It is evident that

\begin{equation} (K_{ab}^2 )_{ij} = (K_{ba}^2 )_{ji} . \end{equation}

(A9)

 The third term on the left-hand side of eq. (21) represents the restoring forces induced by displacements at every boundary that has a density discontinuity. These can be computed through surface integration.

The forcing terms

The right-hand side of eq. (21) is composed of the forcing terms. Following the same procedure that we use for the stiffness matrix, the forcing terms can be expressed as the sum of two force vectors, and eq. (21) becomes

\begin{equation} KV = F_0 + F(\Delta \phi (V)), \end{equation}

(A10)

where K is the total stiffness matrix, V is the incremental displacement vector containing | $ \vec v$| at all the nodes and F₀ is the force vector that depends on the surface load Γ, the pre-stresses τ^pre_ij, the displacements from the previous time step U_i and the initial gravitational potential Φ⁰ (i.e. the total gravitational potential at the previous time step plus the gravitational potential induced by the incremental load itself). F is the force vector that depends on the incremental gravitational potential ΔΦ(V), which, in turn, depends on the incremental displacements V.

APPENDIX B: SEMI-ANALYTIC METHOD

For an spherically symmetric, compressible Earth, we use a spectral method (Han & Wahr 1995) to find the Love numbers in the Laplace transformation domain for each spherical harmonic degree l and order m. Following Wu & Peltier 1982, we expand the Love numbers as a Laurent series of first-order poles:

\begin{equation} \tilde h_l (s) = h_l^E + \sum\limits_j {\frac{{r_j^h }}{{s + s_j }}} , \end{equation}

(B1)

\begin{equation} \tilde k_l (s) = k_l^E + \sum\limits_j {\frac{{r_j^k }}{{s + s_j }}} , \end{equation}

(B2)

where | $ \tilde h_l (s)$| is the Love number for the surface vertical displacement, | $ \tilde k_l (s)$| is the Love number for the surface potential, h^E_land k^E_l are the elastic Love numbers, s is the Laplace transform variable, − s_j denotes the poles and r^h_j and r^k_j are the residues for | $ \tilde h_l$| and | $ \tilde k_l$|⁠, respectively. Instead of determining the modes − s_j, and the residues r^h_j and r^k_j directly, we apply the collocation technique (Peltier 1974). Evenly spaced points { s_j|j = 1, 2, …, N} are chosen so that they cover the expected range of the negatives of the poles, and the Earth's response | $ \tilde h_l (s_j )$| and | $ \tilde k_l (s_j )$| are computed at those points. Using least-square fitting, we solve for r^h_j and r^k_j such that

\begin{equation} \tilde h_l (s_k ) = h_l^E + \sum\limits_{j = 1}^N {\frac{{r_j^h }}{{s_k + s_j }}} , \end{equation}

(B3)

\begin{equation} \tilde k_l (s_k ) = k_l^E + \sum\limits_{j = 1}^N {\frac{{r_j^k }}{{s_k + s_j }}} . \end{equation}

(B4)

The robustness of this application of the collocation technique has been verified in Zhong et al. (2003) for an incompressible Earth model. For a compressible Earth, because we are approximating an infinite series (since we have an infinite number of normal modes) by a finite summation, care needs to be taken when choosing s_j. We define

\begin{equation} \delta _l^h = \frac{{\tilde h_l (s \to 0) - h_l^E }}{{\sum\limits_{j = 1}^N {\frac{{r_j^h }}{{s_j }}} }}, \end{equation}

(B5)

\begin{equation} \delta _l^{k^{\prime}} = \frac{{\tilde k_l (s \to 0) - k_l^E }}{{\sum\limits_{j = 1}^N {\frac{{r_j^k }}{{s_j }}} }}. \end{equation}

(B6)

To obtain reliable results, we adjust the spacing and the range of the s_j's so that the numerical results for δ^h_l and δ^k_l are close to 1 (Mitrovica & Peltier 1992). We compute the numerators in eqs (B5) and (B6) directly through our code, by finding | $ \tilde h_l$| and | $ \tilde k_l$| at the smallests_j; we compute the denominators by summing the residues. Once we have chosen the s_j's so that the results for δ^h_l and δ^k_l differ from 1 by less than 3 per cent, we use the model results computed with those values to obtain Love number solutions for impulse forcing in the time domain:

\begin{equation} h_l (t) = h_l^E \delta (t) + \sum\limits_{j = 1}^N {r_j^h \exp ( - s_j t)} , \end{equation}

(B7)

\begin{equation} k_l (t) = k_l^E \delta (t) + \sum\limits_{j = 1}^N {r_j^k \exp ( - s_j t)} , \end{equation}

(B8)

where h_l(t) and k_l(t) are the time-dependent Love number solution and δ(t) is the Dirac delta function. Since the time dependence of the loading we apply is either a Heaviside function or a piecewise linear function (ICE-5G), it can be analytically convolved with the impulse response in the time domain. For the semi-analytic model, the treatments for degree-1 deformation, polar wander feedback, and ocean loading have been discussed in Paulson et al. (2005). To find the correction to the degree 2 order 1 Love numbers induced by polar wander feedback, we compute

\begin{equation} h_2^{ \pm 1} \left( s \right) = h_2^L \left( s \right) + h_2^T \left( s \right)\left( {\frac{{1 + k_2^L \left( s \right)}}{{k_f - k_2^T \left( s \right)}}} \right), \end{equation}

(B9)

\begin{equation} k_2^{ \pm 1} \left( s \right) = k_2^L \left( s \right) + k_2^T \left( s \right)\left( {\frac{{1 + k_2^L \left( s \right)}}{{k_f - k_2^T \left( s \right)}}} \right), \end{equation}

(B10)

where the superscripts L and Tdenote the load and tide love numbers, respectively. Here, the fluid love number is defined as

\begin{equation} k_f = (1 + \delta ) \cdot k_{{\rm hyd}}^{} , \end{equation}

(B11)

where δ is set to be 0.8 per cent (Mitrovica et al.2005), and k_hyd is the degree-2 tide Love number, k^T₂, computed in the fluid limit. That is, k_hyd represents the tidal response of an Earth that has the same density profile as the model Earth, but where the material is all assumed to be fluid (Mitrovica et al.2005 refer to k_hyd as k^T₂(s → 0; LT = 0)). At long temporal scales (i.e. as s → 0), k^T₂ approaches k_f, which makes the denominators in eqs (B9) and (B10) close to 0. In fact, if the lithosphere were modelled as viscoelastic (presumably with a very large viscosity), then k^T₂(s → 0) → k_hyd exactly, and so the dominators would converge exactly to 0 if δ had been chosen to be 0. So, the estimation of these two equations is very sensitive to the value of k_f. Numerical errors that occur at very small values of s prohibit us from computing k_hyd by taking the s → 0 limit of k^T₂(s). Instead, to obtain an accurate solution, we calculate k_hyd by computing k^T₂ for a fluid Earth using the boundary conditions for the gravitational potential and its gradient (Dahlen 1974). The same fluid love number is also used to obtain (C − A)_hyd in the finite-element model.