# Energy balance of glacial isostatic adjustment: importance of the rotational feedback

Energy balance of glacial isostatic adjustment: importance of the rotational feedback Abstract Understanding the feedback between the glacial isostatic adjustment (GIA) and the Earth’s rotation is important for an accurate prediction of sea level changes induced by climate and tectonic processes. Here we consider a simple, four-layer incompressible Earth model, recently used for a benchmark of GIA codes to estimate the accuracy of the linearized Liouville equation (LE) and to demonstrate that models with an incomplete or missing rotational feedback violate the principle of energy conservation. First, we compute GIA on a rotating Earth by solving the equation of motion coupled with LE in its full nonlinear form. By comparing the nonlinear LE solution with the traditional linearized one, we find that the radial component of the angular velocity vector is inaccurate in the latter case, with an error exceeding 10 per cent already after 1 kyr of evolution. To understand the cause of this discrepancy, we investigate the time evolution of different kinds of energy involved in the process. While the rotational, elastic and dissipative energies are straightforward to compute, the formula for the gravitational energy contains an integral that requires a careful, higher-order accurate evaluation of the gravitational potential perturbation. We circumvent this problem by transforming the integral into a different one, formulated in terms of displacement instead of potential. We find that the solution of the linearized LE equation does not conserve the energy, since, in the linearized case, the rate of change of the rotational energy is not equal to the power of the centrifugal force. We also compute the energy balance of GIA on a constantly rotating Earth, and demonstrate the importance of the rotational feedback in the equation of motion. The formalism derived in this study allows a detailed examination of the energy balance for a rotating, incompressible planetary body deformed by a surface load. As such, it may not only serve as a reliable tool for a posteriori testing of GIA numerical solutions but it can also be used in different planetary science applications. Elasticity and anelasticity, Earth rotation variations, Loading of the Earth, Sea level change, Numerical modelling 1 INTRODUCTION Glacial isostatic adjustment (GIA) has been investigated within the geophysical community for many decades, recently gaining attention due to the increasing precision of geodetic measurements obtained from the GRACE (Tapley et al. 2004) and GOCE (Floberghagen et al. 2011) satellite missions and the growing interest in sea level changes. The response of the Earth to surface loading was used primarily to constrain the viscosity profile of the mantle (Cathles 1975). However, the depth to which viscosity can be inferred by measuring postglacial rebound is limited by the lateral extent of the past glaciers, leaving the lowermost mantle poorly resolved (e.g. Mitrovica & Peltier 1991). The motion of the rotation axis of the Earth, induced by mass redistribution that accompanies the periodic accumulation and melting of ice during glacial cycles, can be used to overcome this limitation. Assuming that the ice load history is known, lower-mantle viscosity can be constrained by fitting the observed rate of the secular drift (e.g. Sabadini et al. 1982; Wu & Peltier 1984; Vermeersen et al. 1998; Lau et al. 2016). Also worth noting is the pioneering attempt of Nakada & Karato (2012) to infer the lower-mantle viscosity by comparing the Q-factor of the predicted Chandler wobble with the observed value (e.g. Benjamin et al. 2006). It was only later that the perturbation of the centrifugal force, inherent to any change in the rotation vector, was also included in modelling GIA-related suite of observables. While its effect on the sea level equation has been extensively studied (e.g. Milne & Mitrovica 1998; Mitrovica et al. 2001), its impact on the GIA-induced deformation of the Earth has only been addressed recently by Martinec & Hagedoorn (2014). Precise computation of GIA-induced polar motion (PM) thus plays an important role in numerical simulations of postglacial rebound related phenomena. Following Munk & MacDonald (1960) it has become common to use the linearized form of the Liouville equation (LE) to obtain the PM solution. The error introduced by the linearization is not straightforward to estimate and depends on the amplitudes of the wander that the rotation axis experiences. In this paper, we couple the fully nonlinear LE (NLE) with the equations governing small deformations of a radially symmetric, self-gravitating incompressible spherical shell. By comparing the NLE solution with the linearized one we quantitatively evaluate the error due to neglecting higher-order terms in the linearized LE (LLE). In particular, we estimate the error arising from decoupling the equations for PM and the length of day (LOD) variation, which is an intrinsic property of the linearized approach. The accuracy of the LLE is tested for the Earth model M3-L70-V01 loaded with a spherical ice cap, a community benchmark by Spada et al. (2011). Another way to assess the accuracy of the linearized solution, without having to solve the NLE, is to check whether the total energy of the Earth’s model is conserved. The energy balance involves the rotational, elastic and dissipative energies, which are straightforward to evaluate, and the gravitational energy, which requires a higher-order accurate evaluation of the gravitational perturbation potential. In Section 4.1, we solve this problem by deriving an equivalent integral but with the displacement in the integrand, allowing the use of a simple spectral formula for the geoid (e.g. Choblet et al. 2007). We then compute the energy balance for various test examples relevant to GIA studies. In Section 4.2, we load a non-rotating Earth with a spherical ice cap. To assess the energetic importance of free wobble, triggered by a hypothetical rapid accumulation of a surface load on a rotating Earth, we separately study this phenomenon in Section 4.3. In Section 4.4, we analyse the energy balance of a rotating Earth loaded with a spherical ice cap, that is, the benchmark case by Spada et al. (2011). We show that the LLE solution does not conserve the energy and discuss the consequences of this finding for the prediction of PM and LOD. In Section 4.5, we consider the same surface loading, but on a constantly rotating Earth. A large deviation from the energy conservation is obtained in such case, proving the importance of including rotational feedback in GIA modelling. All results discussed in Section 4 are obtained for the four-layer, incompressible Earth model M3-L70-V01 used already in the benchmark of GIA codes (Spada et al. 2011). Although this structural model is relatively simple in comparison with recent models inferred from ice age data sets (e.g. Lau et al. 2016), it allows our predictions of PM and LOD to be compared with those in Spada et al. (2011) and the energy curves presented in Section 4 to be easily reproduced. The energy balance analysis introduced in this paper does not include the effect of compressibility. As shown by Cambiotti et al. (2010), compressibility of the real Earth can have non-negligible effect on both GIA and PM. Compressible models are more deformable than the incompressible ones, which makes the readjustment of the rotational bulge quicker, effectively reducing the GIA-induced PM. This reduction, however, depends on lower mantle viscosity and is almost negligible for values higher than 1022 Pa s, which are usually needed to stabilize the rotation of the Earth on geological time scales (Ricard et al. 1993). To compute small deformations of a rotating, gravitationally pre-stressed Earth we use the Eulerian formulation of the governing equations, instead of applying the traditional Lagrangian approach (e.g. Wu & Peltier 1982). The Eulerian formulation has appeared in the literature several times in the past decade (e.g. Tobie et al. 2008; Golle et al. 2012; Soucek et al. 2016), but only for spherical shells with a constant density. We extend it for the case where the internal density is a continuous piecewise linear function and we numerically demonstrate that the solution is equivalent to the Lagrangian one. Using the Eulerian formulation we demonstrate a one-by-one correspondence between individual forces in the equation of motion and the terms describing the gravitational, rotational, elastic and dissipative energies. Such a correspondence provides a powerful tool when diagnosing a numerical solution. While the evaluation of the sum of all involved energies reveals only whether or not the solution is numerically correct, the term-by-term correspondence can be used to directly identify a potential source of energy imbalance. In the present paper, we do not use the normal mode theory (e.g. Peltier 1974) or complex contour integration (e.g. Tanaka et al. 2009; Sabadini et al. 2016), as we solve all equations in the time domain, making the coupling between the equation of motion and the LE straightforward. It is worth noting that the use of the NLE does not require the amplitudes of the modelled PM to be small and the time scales on which mass redistribution occurs can be arbitrary. Our approach thus broadens the variety of processes that can be addressed compared to the traditional techniques: the linearized approach is limited to PM with small amplitudes, and the viscous quasi-fluid approximation, introduced by Lefftz et al. (1991) and Ricard et al. (1993), is valid only for mass redistribution occurring on the time scale of a few million years (see also Cambiotti et al. 2011). The rapid formation of craters and volcanoes that can subsequently shift the rotation axis of a planetary body by tens of degrees (Runcorn 1984; Kite et al. 2009), and a number of asteroids that experience a large-amplitude wobbling (Harris 1994), are examples of the problems that could newly be tackled with the tool we present. 2 GOVERNING EQUATIONS We consider small deformations of a hydrostatically pre-stressed spherical shell with radially dependent reference density profile ρ0(r) and a homogeneous fluid core. The material of the spherical shell is assumed to behave like an incompressible Maxwell-type viscoelastic fluid. The time evolution of displacement and stress can be obtained by integrating the following set of partial differential equations (e.g. Tobie et al. 2008):   $$\nabla \cdot \boldsymbol {\tau } + \boldsymbol {f} = 0 \, ,$$ (1)  $$\nabla \cdot \boldsymbol {u}=0\, ,$$ (2)  $$\boldsymbol {\tau }^d-\mu (\nabla \boldsymbol {u}+(\nabla \boldsymbol {u})^\mathrm{T})=-\frac{\mu }{\eta }\int _0^t\boldsymbol {\tau }^d \mathrm{d}t^{\prime } \, ,$$ (3)where τ is the Cauchy stress tensor, τd is its deviatoric part, f is the body force, u is the displacement, the superscript T denotes transposition of a tensor, μ is the elastic shear modulus, η is the viscosity, and t is the time. At time t = 0 the integral on the right-hand side of eq. (3) is zero and the spherical shell behaves like an elastic solid. Throughout the paper, we strictly use the Eulerian description of the problem. All variables are expressed as functions of the instantaneous position of the particle in the deformed body. A particle occupying the position r at time t would occupy the position r − u(r, t) if all forces that cause the deformation were set to zero at all times. The body force f is given as the product of density and the negative gradient of gravity potential. The gravity potential has three components: (i) the reference gravitational potential V0, (ii) the centrifugal potential Ψ due to the rotation of the body, and (iii) the perturbation Φ of the gravitational potential V0 due to the deformation of the body. V0 is the gravitational potential of a sphere of radius a with density profile ρ0(r). Inside the sphere (r ≤ a)   $$V_0(r)=-\frac{4\pi G}{r}\int _0^{r} \rho _0(r^{\prime }){r^{\prime }}^2\, \mathrm{d}r^{\prime }- 4\pi G\int _r^{a} \rho _0(r^{\prime })r^{\prime }\, \mathrm{d}r^{\prime },$$ (4)where G is the universal gravitational constant. The centrifugal potential Ψ is given as   $$\Psi = \frac{1}{2}\left((\boldsymbol {\omega }\cdot \boldsymbol {r})^2-|\boldsymbol {\omega }|^2|\boldsymbol {r}|^2 \right),$$ (5)where ω is the angular velocity vector which has the direction of the rotation axis and the magnitude equal to the angular speed of the body. As the sphere deforms, the initial density ρ0(r) changes to density ρ(r, t) which can be expressed as a sum of density ρ0(r) and Eulerian density increment δρ(r, t). If the body is incompressible and the density ρ0 is a continuous function,   $$\delta \rho (\boldsymbol {r},t) = \rho (\boldsymbol {r},t) - \rho _0(r) = \rho _0(|\boldsymbol {r}-\boldsymbol {u}(\boldsymbol {r},t)|,t) - \rho _0(r) = -\boldsymbol {u}(\boldsymbol {r},t)\cdot \nabla \rho _0(r) + O(|\boldsymbol {u}|^2),$$ (6)where the displacement u describes the deformation of the body. The Eulerian density increment δρ is non-zero also at the outer surface and at the internal density interfaces where it is equal to the density jump at the respective interface. Using δρ we can express the perturbation potential Φ in an integral form:   $$\Phi (\boldsymbol {r},t)=-G\int _{v(t)} \frac{\delta \rho (\boldsymbol {r}^{\prime },t)}{|\boldsymbol {r}-\boldsymbol {r}^{\prime }|}\, \mathrm{d}v^{\prime },$$ (7)where v(t) is the volume occupied by the body (including the core) at time t. Neglecting the terms O(|u|2) and −δρ∇(Φ + Ψ), we obtain the following expression for the body force f:   $$\boldsymbol {f} = -\rho \nabla (V_0+\Phi +\Psi ) \cong - (\boldsymbol {u}\cdot \nabla \rho _0)\boldsymbol {g}_0 - \rho _0\nabla \left(\Phi +\Psi \right) + \rho _0\boldsymbol {g}_0,$$ (8)where g0 = −∇V0. The laterally homogeneous and static contribution ρ0g0 is counteracted by the hydrostatic pre-stress p0(r), satisfying the equation − ∇p0 + ρ0g0 = 0. Eqs (1)–(3) are solved in a spherical shell with outer radius a and inner radius b. The boundary condition at the outer surface is obtained from the force equilibrium condition, taking into account the pressure due to the deformation-induced topography (e.g. Soucek et al. 2016) and the presence of a surface load:   $$\boldsymbol {\tau }\cdot \boldsymbol {e}_r=\left(u_r[\rho _0]_a+\sigma ^\mathrm{L}\right)\boldsymbol {g}_0,$$ (9)where er is the radial unit vector, ur is the radial component of displacement u, [ρ0]a = ρ0(a) is the density jump at the surface, and σL is the surface mass density of the prescribed surface load. A similar condition can be imposed at the bottom boundary, on Earth corresponding to the core–mantle boundary:   $$-\!\boldsymbol {\tau }\cdot \boldsymbol {e}_r = u_r[\rho _0]_b\, \boldsymbol {g}_0-\rho _\mathrm{c}(V_0+\Phi +\Psi )\boldsymbol {e}_r,$$ (10)where [ρ0]b is the density jump at the bottom boundary and −ρc(V0 + Φ + Ψ) is the hydrostatic pressure acting on the boundary due to the contact with a liquid core of density ρc. To compute the temporal evolution of the angular velocity vector ω with respect to the body-fixed Tisserand frame (Munk & MacDonald 1960) we solve the LE with zero external torque,   $$-\!\boldsymbol {I}\cdot \frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t} = \frac{\mathrm{d}\boldsymbol {I}}{\mathrm{d}t}\cdot \boldsymbol {\omega }+\boldsymbol {\omega }\times (\boldsymbol {I}\cdot \boldsymbol {\omega }),$$ (11)where I is the time-dependent tensor of inertia. Eqs (1)–(10) are coupled with eq. (11) both through the displacement field u, which is needed to compute the inertia tensor I, and through the centrifugal potential Ψ, which depends on the angular velocity vector ω. 3 NUMERICAL IMPLEMENTATION The LE is a set of ordinary differential equations for three unknown components of vector ω as functions of time. It can be rewritten as   $$\frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t} = \mathcal {F}(\boldsymbol {I},\boldsymbol {\omega })\, ,$$ (12)where   $$\mathcal {F}=-\boldsymbol {I}^{-1}\cdot \left(\frac{\mathrm{d}\boldsymbol {I}}{\mathrm{d}t}\cdot \boldsymbol {\omega }+\boldsymbol {\omega }\times (\boldsymbol {I}\cdot \boldsymbol {\omega })\right).$$ (13)Eq. (12) is solved numerically using the fifth order accurate Adams-Bashforth multistep method with a time step of 10−3 yr. Each evaluation of $$\mathcal {F}$$ requires the calculation of the tensor of inertia I and its time derivative. The perturbation of the tensor of inertia is computed by MacCullagh’s formula from the gravitational potential Φ and the reference density profile ρ0 (e.g. Patočka 2013, eqs 2.12),   $$\boldsymbol {I} = I_0\, \mathcal {I} + \sqrt{\frac{5}{\pi }}\, \frac{a^3}{6\, G} \, \left[ \, \left( \begin{array}{c@{\quad}c@{\quad}c}-\Phi _{20} & 0 & 0\\ 0 & -\Phi _{20} & 0\\ 0 & 0 &2 \Phi _{20} \end{array}\right) \, - \, \sqrt{6} \, \left( \begin{array}{c@{\quad}c@{\quad}c}-\mathrm{Re}\, \Phi _{22} & \mathrm{Im}\, \Phi _{22} & \mathrm{Re}\, \Phi _{21} \\ \mathrm{Im}\, \Phi _{22} & \mathrm{Re}\, \Phi _{22} & -\mathrm{Im}\, \Phi _{21} \\ \mathrm{Re}\, \Phi _{21} & -\mathrm{Im}\, \Phi _{21} & 0 \\ \end{array}\right) \right].$$ (14)Here $$I_0:=\frac{8}{3}\pi \int ^a_{0} r^4 \rho _{0} \mathrm{d}r$$, $$\mathcal {I}$$ is the identity tensor, and Φℓm are the coefficients of the spherical harmonic expansion of potential Φ,   $$\Phi (r,\vartheta ,\varphi )=\sum _{\ell =0}^{\ell _{\rm max}}\sum _{m=-\ell }^{\ell }\Phi _{\ell m}(r)Y_{\ell m}(\vartheta ,\varphi )\, ,$$ (15)where Yℓm are the fully normalized spherical harmonics of degree ℓ and order m (e.g. Jones 1985), and ℓmax is the cut-off degree. The potential Φ depends on the deformation of the body and we describe its calculation in the next paragraph. The term $$\mathrm{d}\boldsymbol {I}/\mathrm{d}t$$ is then computed numerically using the second order accurate mid-point finite difference scheme. The spectral algorithm for the numerical solution of eqs (1)–(10) for a given centrifugal potential Ψ is outlined in Tobie et al. (2008), identities needed to construct the matrix of the discretized problem are summarized in the appendix to Golle et al. (2012). The potential Φ depends on u and vice versa. It is computed iteratively by expanding the integral (7) into a spherical harmonics series, while condensing topographies into surface mass densities (see eq. 35 in Choblet et al. 2007):   \begin{eqnarray} \Phi _{\ell m}(r) = \frac{-4\pi Gr}{2\ell +1}\Bigg ( [\rho _0]_{b}\, u_{r,\ell m}(b) \left(\frac{b}{r}\right)^{\ell +2} + [\rho _0]_{a}\, u_{r, \ell m}(a) \left(\frac{r}{a}\right)^{\ell -1}- \int _{b}^r \left(\frac{r^{\prime }}{r}\right)^{\ell +2}{u}_{r,\ell m}(r^{\prime })\frac{\mathrm{d}\rho _0}{\mathrm{d}r^{\prime }}\, \mathrm{d}r^{\prime } - \int _r^{a} \left(\frac{r}{r^{\prime }}\right)^{\ell -1}{u}_{r,\ell m}(r^{\prime })\frac{\mathrm{d}\rho _0}{\mathrm{d}r^{\prime }}\, \mathrm{d}r^{\prime } \Bigg )\, , \end{eqnarray} (16)where we substituted for δρ from eq. (6) on the second line of eq. (16). The discretization in the radial direction is performed by the finite difference method on a staggered grid with constant spacing (e.g. Gerya & Yuen 2003). Our approach differs from the algorithm described in Tobie et al. (2008) by employing a higher order Crank-Nicholson integration scheme for evaluating the time integral in eq. (3). The accuracy of the numerical method used to solve the governing eqs (1)–(10) was carefully tested against the traditional Lagrangian solution. Fig. 1 shows the relative difference between the loading and tidal degree 2 Love numbers $$k^L_e$$, $$k^T_e$$, $$k^L_f$$ and $$k^T_f$$, computed numerically using the method described above, and those computed semi-analytically by Spada et al. (2011) for a viscoelastic Earth model M3-L70-V01 (for description of the model, see table 3 therein). The difference is plotted as a function of the radial resolution considered in the numerical solution in which the density profile is approximated by a continuous piecewise linear function. Inspection of Fig. 1 shows that the relative error of the numerical method converges to zero with increasing resolution, decreasing below 10−4 when the number of equally spaced radial nodes is larger than 5000. The error does not decrease monotonically because it depends on how well the equidistant discretization matches the positions of the density interfaces. In Section 4.4, we will use 460 radial grid nodes for which the error in determining the Love numbers is less than 10−4. The number of radial grid nodes could be significantly smaller if we used a variable grid spacing. Figure 1. View largeDownload slide Relative error of the numerically computed tidal and loading degree 2 Love numbers $$k^T_e$$, $$k^T_f$$, $$k^L_e$$ and $$k^L_f$$ relative to the semi-analytically computed values by Spada et al. (2011). The subscript e denotes the instantaneous response and f is the t → ∞ ‘fluid’ limit. The shadowed region marks the domain where the relative error is smaller than 10−4. Figure 1. View largeDownload slide Relative error of the numerically computed tidal and loading degree 2 Love numbers $$k^T_e$$, $$k^T_f$$, $$k^L_e$$ and $$k^L_f$$ relative to the semi-analytically computed values by Spada et al. (2011). The subscript e denotes the instantaneous response and f is the t → ∞ ‘fluid’ limit. The shadowed region marks the domain where the relative error is smaller than 10−4. Besides the computation of the Love numbers, we have also reproduced the benchmark example in Spada et al. (2011), where model M3-L70-V01 is loaded with a spherical ice cap (Spada et al. 2011, table 4), centred at the colatitude θc = 25° and the longitude λc = 75°. Our solution consists of two steps. In the first one, the model M3-L70-V01 is rotated at a constant angular velocity ω0 = [0, 0, 7.292115] × 10−5 s−1 until it reaches hydrostatic equilibrium. Numerically we achieve this by integrating eqs (1)–(10) with the prescribed ω0 and σL ≡ 0 over a period of 200 Myr (similar procedure was used to obtain the fluid Love numbers in Fig. 1). We will denote the hydrostatic values of the polar and equatorial moments of inertia by C and A respectively ($$C\,{=}\,I_{33}^{f}$$ and $$A\,{=}\,I_{11}^f\,{=}\,I_{22}^f$$). Then, at time which we formally mark as t = 0, the body is suddenly loaded with a spherical ice cap. The loading instantaneously shifts the figure axis, and we adjust the vertical component of ω(t = 0), so that the total angular momentum is preserved in first approximation (ω = ω0(1 − c33/C), where c33 is the change of polar moment of inertia due to the loading). The further evolution of the tensor of inertia, which in turn determines the evolution of ω(t), is governed by both the viscoelastic relaxation under the load and the viscoelastic readjustment of the rotational bulge in response to the induced PM. For clarity of comparison with other studies, we introduce two quantities characterizing the PM—the PM vector m and the LOD variation ΔLOD:   $$\boldsymbol {\omega }(t)=\boldsymbol {\omega }_0 + |\boldsymbol {\omega }_0|\boldsymbol {m}(t)\, , \quad \Delta \mathrm{LOD}(t) = \frac{c_{33}(t)}{C}\frac{2\pi }{|\boldsymbol {\omega }_0|}\, ,$$ (17)where c33(t) is the time variation of I33 (polar moment of inertia) in response to the loading. In Fig. 2(a), we compare the PM components m1 and m2 computed by our method (coloured points) with those obtained using the traditional Laplace domain approach and the LLE (dashed lines, adopted from Martinec & Hagedoorn 2014). Since our solution includes the wobbling while the LLE solution by Martinec & Hagedoorn (2014) gives only the mean position of the rotation axis, the two solutions are very different in the beginning, but they match each other once the wobbling is damped. This confirms that our code gives correct results on long time scales (i.e. greater than the damping time of Chandler wobble), and conversely, that neither the wobbling nor the linearization influences the long term evolution of the PM. A good agreement is also found for ΔLOD (Fig. 2b) where the relative difference between the two solutions is smaller than 1 per cent. The study by Spada et al. (2011) provides a valuable benchmark for validation of different numerical techniques to compute the PM. Nevertheless, there are two things that are worth noting because they are not explained in the original paper in sufficient details. First, in fig. 13(c) in Spada et al. (2011), the variation ΔLOD induced by the load is plotted as a positive quantity. However, since the load is imposed close to the north pole and the total mass of the Earth is conserved, the loading reduces the polar moment of inertia and the variation ΔLOD is actually negative. Second, the hydrostatic values C and A must be corrected to reproduce the benchmark results. The ΔLOD curve (Fig. 2b) is sensitive to the value of C, which has to be 8.0394 × 1037 kg m2 to match the benchmark (i.e. the value given in table 1 of Spada et al. 2011). However, to match the benchmarked PM solution (Fig. 2a), which is sensitive to the difference C − A, this difference has to be 2.6947 ×1035 kg m2. This value corresponds to model M3-L70-V01 in the hydrostatic limit, and not to the value of C − A derived from table 1 in Spada et al. (2011). Numerically we proceeded by adding a diagonal correction to our tensor of inertia, $$\boldsymbol {I}^c(t\,{\ge }\,0)=\boldsymbol {I}(t) + \gamma \mathcal {I}$$, where $$\mathcal {I}$$ is identity tensor and Ic is the corrected tensor of inertia used in our simulations. The correction γ is given as (8.0394 − 8.17848) × 1037 kg m2, with 8.17848 × 1037 kg m2 being the hydrostatic polar moment of inertia of model M3-L70-V01. We believe this information can be helpful for the future users of the benchmark. Figure 2. View largeDownload slide (a) Polar motion induced by loading the Earth with a spherical ice cap. The coloured lines show the solution obtained by the method described in this study while the dashed black lines show the linearized solution computed in the Laplace domain and omitting the free wobble (Martinec & Hagedoorn 2014). (b) Variation of ΔLOD obtained for the same loading (blue line—our solution, dashed line—linearized solution). Figure 2. View largeDownload slide (a) Polar motion induced by loading the Earth with a spherical ice cap. The coloured lines show the solution obtained by the method described in this study while the dashed black lines show the linearized solution computed in the Laplace domain and omitting the free wobble (Martinec & Hagedoorn 2014). (b) Variation of ΔLOD obtained for the same loading (blue line—our solution, dashed line—linearized solution). 4 ENERGY BALANCE OF A ROTATING EARTH The conservation of energy for a rotating self-gravitating viscoelastic body can be expressed as   $$E_{\mathrm{rot}}+E_{\mathrm{el}}+E_{\mathrm{diss}}+E_{\mathrm{grav}}={\rm const}\, ,$$ (18)where the terms on the left-hand side are the rotational, elastic, dissipative and gravitational energies, respectively. The kinetic energy associated with the rate of deformation is not considered, since the inertia term in the equation of motion (1) is neglected. The rotational energy Erot can be computed using the standard formula   $$E_{\mathrm{rot}}=\frac{1}{2}\boldsymbol {\omega }\cdot \boldsymbol {I}\cdot \boldsymbol {\omega }.$$ (19)The elastic energy Eel stored in the body and the viscous dissipation rate $$\boldsymbol {\dot{E}}_\mathrm{diss}$$ can be expressed as (cf. Joseph 1990, p. 50)   $$E_{\mathrm{el}}=\int _{v(t)}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{4\mu }\, \mathrm{d}v \cong \int _{{\rm S}_a\!-\!{\rm S}_b}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{4\mu }\, \mathrm{d}v \, ,$$ (20)  $$\boldsymbol {\dot{E}}_{\mathrm{diss}}=\int _{v(t)}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{2\eta }\, \mathrm{d}v \cong \int _{{\rm S}_a\!-\!{\rm S}_b}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{2\eta }\, \mathrm{d}v \, ,$$ (21)where $$\boldsymbol {\tau }^d:\boldsymbol {\tau }^d \equiv \tau _{ij}^d\, \tau _{ij}^d$$ in Cartesian components and Sa − Sb denotes the spherical shell of outer radius a and inner radius b. The integration over the core is omitted since the core is assumed to be filled with inviscid liquid. Finally, the gravitational potential energy is given by   $$E_{\mathrm{grav}}=\frac{1}{2}\int _{v(t)}\rho \left(V_0 + \Phi \right)\mathrm{d}v \, .$$ (22) 4.1 Alternative formula for the gravitational energy The computation of the gravitational energy is a delicate exercise which requires a careful evaluation of topographic and volumetric contributions to integral (22). In this section we derive a formula for the gravitational energy which is easy to implement and, as we will demonstrate in Section 4.4, it is sufficiently accurate if the deformation is small. We first derive the formula for a model with a finite number of layers of constant densities, and then we extend it to a model with a continuous density profile. For the derivation we will assume that the surface load density σL is equal to zero. For a non-zero surface load density the problem can be converted into a problem with zero surface load density by using a modified radial displacement $$u_r^\mathrm{mod}(a)=u_r(a)+\sigma ^\mathrm{L}/[\rho _0]_a$$. For the sake of clarity of exposition, we consider a multilayer (N-layer) sphere. The N-layered sphere is composed of N spheres Si with radii ri, i = 1, 2, … N, where r1 = a and rN = b are respectively the radii of the outer and inner boundary of the spherical shell considered in Section 2. The density ρ0(r) of the sphere is a piecewise constant function with density jumps [ρ0]i at radii ri (taken positive if the density increases with depth). When the body is deformed, each sphere Si transforms into an aspherical object vi and the density changes by   $$\delta \rho (\boldsymbol {r},t) = \rho _0(|\boldsymbol {r}-\boldsymbol {u}(\boldsymbol {r},t)|,t) - \rho _0(r) = \left\lbrace \begin{array}{l@{\quad}l@{\quad}l}0, & |{r}{-}{r}_i|\,{\ge }\,|{u}_r(\boldsymbol {r},t)| & \forall i\in \lbrace 1,2,\ldots ,N\rbrace \\ [\rho _0]_i, & {u}_r(\boldsymbol {r},t)\,{>}\,0 \quad \wedge & |{r}\,{-}\,{r}_i|\,{<}\,|{u}_r(\boldsymbol {r},t)|\\ {-}[\rho _0]_i, & {u}_r(\boldsymbol {r},t)\,{<}\,0 \quad \wedge & |{r}\,{-}\,{r}_i|\,{<}\,|{u}_r(\boldsymbol {r},t)| \end{array} \right.$$ (23)where the three possibilities given on the right are illustrated in Fig. 3. The volume integral in eq. (22) can then be rewritten as   $$\int _{v(t)}\rho \left(V_0 + \Phi \right)\mathrm{d}v = \sum _{i=1}^{N}[\rho _0]_i\int _{v_i} \left(V_0 + \Phi \right)\mathrm{d}v.$$ (24)On the right-hand side of this equation, we obtain a sum of N integrals over N homogeneous bodies which have no internal boundaries, occupy the volumes vi and have the densities [ρ0]i. Each volume vi can be expressed in terms of the spherical volume Si and the volume hi induced by the deformation (see the yellow and green regions in Fig. 3). The volume hi is taken with a positive sign if it is above the sphere Si (yellow regions) and with a negative sign if it is inside the sphere (green regions). The total volume of hi at each interface (the sum of the yellow and green domain) is zero due to the incompressibility of the body. The gravitational energy can then be rewritten using Si and hi as follows:   $$E_{\mathrm{grav}}=\frac{1}{2}\sum _{i=1}^{N}[\rho _0]_i\int _{h_i} \left(V_0 + \Phi \right)\mathrm{d}v + \frac{1}{2}\sum _{i=1}^{N}[\rho _0]_i\int _{\mathrm{S}_i} (V_0+\Phi )\, \mathrm{d}v.$$ (25) Figure 3. View largeDownload slide Deformation of N-layer sphere. Black lines show boundaries of the reference spherical shells with constant densities, red lines show the interfaces of the deformed body. In regions where the deformation produces a positive topography, the Eulerian density increment is equal to [ρ0]i, i = 1, 2, …, N, while in regions where the topography is negative, the density increment is equal to −[ρ0]i. Figure 3. View largeDownload slide Deformation of N-layer sphere. Black lines show boundaries of the reference spherical shells with constant densities, red lines show the interfaces of the deformed body. In regions where the deformation produces a positive topography, the Eulerian density increment is equal to [ρ0]i, i = 1, 2, …, N, while in regions where the topography is negative, the density increment is equal to −[ρ0]i. The gravitational potential V0 + Φ can be split into N parts, $$V_0+\Phi =\sum _{i=1}^{N}(V_0^i+\Phi ^i)$$, where $$V_0^i$$ is the potential of the homogeneous sphere Si of density [ρ0]i and Φi is the perturbation due to the aspherical shape of the i-th interface:   $$V_0^i(r)=-G[\rho _0]_i\int _{\mathrm{S}_i} \frac{1}{|\boldsymbol {r}-\boldsymbol {r}^{\prime }|}\, \mathrm{d}v,\quad \Phi ^i(\boldsymbol {r})=-G[\rho _0]_i\int _{h_i} \frac{1}{|\boldsymbol {r}-\boldsymbol {r}^{\prime }|}\, \mathrm{d}v.$$ (26)In the following, we take advantage of the well-known fact that the gravitational energy of a body A in the gravitational field of a body B is the same as the gravitational energy of body B in the gravitational field of body A. By applying this general rule to eq. (26) we obtain the following relations:   $$[\rho _0]_i\int _{\mathrm{S}_i}\Phi ^j\, \mathrm{d}v=[\rho _0]_j\int _{h_j} V_0^i\, \mathrm{d}v,$$ (27)valid for any i, j = 1, …, N. Using eq. (27) we can link the first and the second sum on the right-hand side of eq. (25):   $$\sum _{i=1}^{N}[\rho _0]_i\int _{\mathrm{S}_i} \Phi \, \mathrm{d}v = \sum _{i=1}^{N}[\rho _0]_i \int _{\mathrm{S}_i}\sum _{j=1}^{N}\Phi ^j\, \mathrm{d}v = \sum _{j=1}^{N}[\rho _0]_j\int _{h_j} \sum _{i=1}^{N} V_0^i\, \mathrm{d}v = \sum _{i=1}^{N}[\rho _0]_i\int _{h_i} V_0\, \mathrm{d}v.$$ (28)Substituting (28) into (25) gives the alternative formula for the gravitational energy Egrav of the body:   $$E_{\mathrm{grav}}=E_{\mathrm{grav}}^0+\sum _{i=1}^N E_{\mathrm{grav}}^i,$$ (29)where Egrav0 is the gravitational energy of the undeformed N-layer sphere   $$E_{\mathrm{grav}}^0=\frac{1}{2}\sum _{i=1}^N[\rho _0]_i\int _{\mathrm{S}_i} V_0\, \mathrm{d}v=2\pi \int _0^a\rho _0(r) V_0(r) r^{\, 2}\, \mathrm{d}r\, ,$$ (30)and   $$E_{\mathrm{grav}}^i=[\rho _0]_i \int _{h_i}\left(V_0 + \frac{1}{2}\Phi \right)\mathrm{d}v\, .$$ (31)Since the deformation is small, the topography of the ith interface (defined here as the deviation from a sphere) can be approximated by the radial component of the displacement vector at the interface. Eq. (31) then takes the form   $$E_{\mathrm{grav}}^i\cong [\rho _0]_i \int _{\partial {\rm S}_i}\int _{r_i}^{r_i+u_r(r_i,\vartheta ,\varphi )}\left(V_0 + \frac{1}{2}\Phi \right)\mathrm{d}r\, \mathrm{d}s\, ,$$ (32)where ∂Si denotes the surface of the sphere of radius ri, ds is the element of ∂Si, and ϑ and φ are the spherical angular coordinates. The integral in eq. (32) can be simplified by expanding V0 in a Taylor series and setting Φ(r, ϑ, φ) ≅ Φ(ri, ϑ, φ):   $$E_{\mathrm{grav}}^i\cong [\rho _0]_i \int _{\partial {\rm S}_i}\int _{r_i}^{r_i+u_r(r_i,\vartheta ,\varphi )}\left(V_0(r_i)+g_0(r_i)(r-r_i) + \frac{1}{2}\Phi (r_i,\vartheta ,\varphi )\right)\mathrm{d}r\, \mathrm{d}s\, .$$ (33)Integration of eq. (33) over the radius gives   $$E_{\mathrm{grav}}^i\cong \frac{1}{2}[\rho _0]_i \int _{\partial {\rm S}_i}\left(g_0 u_r^2 + \Phi u_r\right)\, \mathrm{d}s\, ,$$ (34)where we used that $$\int _{\partial {\rm S}_i}u_r\mathrm{d}s=0$$ for an incompressible body. The surface integral in eq. (34) can be easily evaluated in spectral domain. Representing quantities ur and Φ in terms of spherical harmonic expansions, eq. (15), and invoking the orthonormality of the spherical harmonic basis {Yℓm}, we obtain   $$E_{\mathrm{grav}}\cong E_{\mathrm{grav}}^0+\frac{1}{2}\sum _{i=1}^N[\rho _0]_i\, r_i^2 \sum _{\ell =0}^{\ell _{\rm max}}\sum _{m=-\ell }^{\ell } \Big (g_0(r_i) |u_{r,\ell m}(r_i)|^2 + u_{r,\ell m}(r_i)\Phi ^{\ast }_{\ell m}(r_i)\Big )\, ,$$ (35)where * denotes complex conjugation. The rotational potential is described by spherical harmonics of degrees 0 and 2. Since degree 0 has no effect on the deformation of an incompressible body, we will only consider the terms with ℓ = 2 and |m| ≤ 2. The formula (35), derived for a layered density model, can be generalized to the case of a continuous density profile. Replacing the density jumps [ρ0]i in eq. (31) by (−dρ0(r)/dr) dr and the sum in eq. (29) by an integral, and following the analogous procedure as above, we get the formula   \begin{eqnarray} E_{\mathrm{grav}}&\cong & E_{\mathrm{grav}}^0 + \frac{1}{2}\bigg ( [\rho _0]_a\int _{\partial {\rm S}_a}\!\!\left(g_0u_r^2+u_r\Phi \right)\mathrm{d}s +[\rho _0]_b\int _{\partial {\rm S}_b}\!\!\left(g_0u_r^2+u_r\Phi \right)\mathrm{d}s - \int _{b}^a \frac{\mathrm{d}\rho _0}{\mathrm{d}r}\int _{\partial {\rm S}_r}\!\! \left(g_0u_r^2+u_r\Phi \right)\mathrm{d}s\, \mathrm{d}r\bigg )\, , \end{eqnarray} (36)where ∂Sr denotes the surface of the sphere with radius r. Eq. (36) allows us to compute the gravitational energy using the quantities obtained by solving the governing equations in a regular computational domain Sa − Sb, that is, in the spherical shell of outer radius a and inner radius b. The radial displacement ur needed to evaluate the integrals in eq. (36) is obtained as the solution to eqs (1)–(3), (9) and (10), with the gravitational potential Φ being computed iteratively from ρ0 and ur via eq. (16), that is, using the traditional ‘condensation’ method (e.g. Choblet et al. 2007). If the deformation is small, the Eulerian and Lagrangian approaches give the same displacement (see Fig. 7b in Section 4.4). Therefore, the formula for the gravitational energy derived here can be used regardless which of the approaches is chosen to evaluate the deformation. Note, however, that a direct application of the condensation method to the energy integral (22) leads to an incorrect result. In order to demonstrate this, we will assume, for simplicity, that the density ρ is constant throughout the body. Eq. (22) can then be expressed as the sum of two integrals over the sphere Sa and two integrals over the domain ha representing the irregular shape of the body (ha ≡ h1 in the notation used above):   $$E_{\rm grav}=\frac{\rho }{2}\left( \int _{{\rm S}_a} V_0\, {\rm d}v+ \int _{{\rm S}_a} \Phi \, {\rm d}v+ \int _{h_a} V_0\, {\rm d}v+ \int _{h_a} \Phi \, {\rm d}v \right)\, .$$ (37)The first integral corresponds to $$E_{\rm grav}^0$$ (cf. eq. 30). The second integral can be evaluated by expanding Φ into spherical harmonics. Due to the orthogonality of the spherical harmonic functions, only Φ00 has non-zero contribution to the integral. In the simplified formula (16), which is based on the condensation technique, there is no coupling among coefficients of different degrees and orders. According to that formula, Φ00 is a function of ur, 00 only, and, since ur, 00 is zero for an incompressible body, the second integral vanishes. Condensing the volumetric density ρ into surface mass density ρur(a) in the last two integrals then gives   $$E_{\rm grav}=E_{\rm grav}^0+\frac{\rho }{2}\left(V_0(a)\int _{\partial {\rm S}_a}u_r\, {\rm d}s+\int _{\partial {\rm S}_a}\Phi u_r\, {\rm d}s\right) = E_{\rm grav}^0+\frac{\rho }{2}\int _{\partial \rm {S}_a}\Phi u_r\, {\rm d}s,$$ (38)while the correct solution is (eq. (36) with [ρ0]a = ρ, [ρ0]b = 0 and dρ0/dr = 0)   $$E_{\rm grav}= E_{\rm grav}^0+ \frac{\rho }{2}\int _{\partial {\rm S}_a}\left(g_0u_r^2+\Phi u_r\right)\, {\rm d}s\, .$$ (39)In other words, the second and third integrals in eq. (37) have been incorrectly evaluated as zero, giving the wrong solution (38). In this section we have shown that both integrals are in fact non-zero, they are equally important, and together yield the contribution with $$g_0 u_r^2$$ in the integrand. The advantage of our formula for the gravitational energy is that it does not require the computation of Φ00, which is a non-zero quantity inside the sphere Sa if evaluated exactly. In appendix A, we derive other formulae that can be helpful in assessing the energy balance of a deforming body, based on evaluating the power of individual forces in the equation of motion (1). 4.2 Energy balance of GIA for a non-rotating Earth We first evaluate the energy balance for a non-rotating Earth with internal structure given by model M3-L70-V01, loaded at time t = 0 with a spherical ice cap as described in Spada et al. (2011). In this case, the body has a perfectly spherical shape prior to the loading (ρ = ρ0, Φ = Ψ = 0, u = 0, τd = 0). No elastic energy is stored and the gravitational energy reaches the absolute minimum. Since the total mass of the imposed surface load is zero, the loading is associated with a redistribution of mass at the surface of the body. This, although partially compensated by the instantaneous elastic deformation of the model due to the loading, results in shape changes and thus in an increase of the gravitational energy (Fig. 4, green curve). At the same time the elastic energy (plotted in magenta) rises to a value that is about ten times smaller than the increase in the gravitational energy. With increasing time, the system returns to equilibrium—the gravitational energy decreases and is dissipated in heat (blue curve in Fig. 4). Figure 4. View largeDownload slide Time evolution of the gravitational (green), dissipative (blue) and elastic (magenta) energy for non-rotating Earth model M3-L70-V01, loaded with a spherical ice cap (for details, see Spada et al. 2011). The variation of the gravitational energy is plotted with respect to the gravitational energy $$E^0_{\rm grav}$$ of the undeformed body. The sum of the three energies is shown by a black line. Figure 4. View largeDownload slide Time evolution of the gravitational (green), dissipative (blue) and elastic (magenta) energy for non-rotating Earth model M3-L70-V01, loaded with a spherical ice cap (for details, see Spada et al. 2011). The variation of the gravitational energy is plotted with respect to the gravitational energy $$E^0_{\rm grav}$$ of the undeformed body. The sum of the three energies is shown by a black line. Since the uppermost layer in model M3-L70-V01 is purely elastic, the final equilibrium state is non-hydrostatic and the values of the elastic and gravitational energies at time t → ∞ are non-zero. The elastic energy does not decrease monotonically, showing a local maximum at t ≅ 50 kyr. This maximum is a consequence of a complex interplay between the elastic lithosphere, in which the stored elastic energy monotonically increases in time, and the deeper mantle layers, where the elastic energy behaves non-monotonically, until it finally returns to zero in the hydrostatic equilibrium. During the whole evolution, the sum of the three considered energies is constant (black line in Fig. 4). 4.3 Energy balance of a damped Chandler wobble We consider the Earth model M3-L70-V01 rotating with constant angular velocity ω0 = [0, 0, 7.292115] × 10 − 5 s−1. The Earth is initially in equilibrium and has an oblate spheroidal shape. At time t = 0, the direction of the rotational vector is suddenly changed so that the angle between the figure axis and the rotational vector is 2.69 × 10−2 degrees, which is approximately the same as the angle by which figure axis is shifted in the benchmark example in Spada et al. (2011)—see Fig. 2(a) for comparison. This change induces a periodic motion (‘wobble’) of the axis of rotation, which is gradually damped due to dissipation. After a few thousands of years, the axis of rotation returns to its initial position, but the final rotational speed is smaller than the initial one since a part of the rotational energy is transformed to heat (note that the angular momentum is not conserved when shifting the axis at time t = 0). The time evolution of the individual types of energies associated with the above process is shown in Fig. 5. All energies are plotted with respect to the initial equilibrium state. In this state (corresponding to t < 0), the elastic energy is non-zero only in the uppermost, purely elastic layer. The abrupt change of the rotational vector at t = 0 induces an instantaneous elastic deformation of the whole body (except the liquid core) and the total elastic energy increases. Since the deformation acts to reduce the flattening, both the gravitational and rotational energies decrease at t = 0. The rotational energy (red curve in Fig. 5) further decreases with increasing time as it is dissipated through viscous deformation (blue curve). The variations in the gravitational and elastic energies (green and magenta curves, respectively) significantly contribute to the global energy budget only in the beginning of evolution. For t → ∞, both curves converge to negative values which are close to zero, suggesting that the gravitational and elastic energies in the final equilibrium state are slightly smaller than those in the initial equilibrium state. This difference is associated with the change of flattening due to the decrease in the rotational speed. Figure 5. View largeDownload slide Time evolution of different types of energy for wobbling model M3-L70-V01. The variations of energy are plotted with respect to the initial equilibrium state in which the rotational vector is parallel to the main axis of the body. The wobble is induced by changing the direction of the rotational vector by 2.69 × 10−2 degrees at time t = 0. Figure 5. View largeDownload slide Time evolution of different types of energy for wobbling model M3-L70-V01. The variations of energy are plotted with respect to the initial equilibrium state in which the rotational vector is parallel to the main axis of the body. The wobble is induced by changing the direction of the rotational vector by 2.69 × 10−2 degrees at time t = 0. 4.4 Energy balance of GIA for a rotating Earth The setting of this numerical experiment corresponds to the benchmark case adopted from Spada et al. (2011), which was already discussed in Section 3 (cf. Figs 1 and 2). The rotating Earth model M3-L70-V01, initially in equilibrium state, is suddenly loaded with an ice cap. The loading and the associated deformation induce a PM (Fig. 2a) and changes in the rotational speed (Fig. 2b), which are accompanied by large variations in the rotational and gravitational energy (Fig. 6, red and green curves, respectively). The amplitudes of these variations are of the order of 1022–1023 J and they are three orders of magnitude larger than the variations in the elastic (magenta) and dissipative (blue) energy (cf. also Figs 4 and 5). During the whole evolution, the sum of all energies is constant (full black line). Figure 6. View largeDownload slide Time evolution of different types of energy for rotating model M3-L70-V01 loaded with an ice cap at t = 0 (for details, see Spada et al. 2011). The variations of energy are plotted with respect to the initial equilibrium state in which the rotational vector is parallel to the main axis of the body. The results obtained using the NLE are represented by coloured lines while the LLE solution is plotted in black. Figure 6. View largeDownload slide Time evolution of different types of energy for rotating model M3-L70-V01 loaded with an ice cap at t = 0 (for details, see Spada et al. 2011). The variations of energy are plotted with respect to the initial equilibrium state in which the rotational vector is parallel to the main axis of the body. The results obtained using the NLE are represented by coloured lines while the LLE solution is plotted in black. So far we have discussed the results computed using the NLE. In Fig. 6 the variations of the gravitational and rotational energy obtained using the LLE are shown by dashed and dotted black lines, respectively. While the curve for the gravitational energy coincides almost exactly with that obtained using the NLE, the accuracy in determining the rotational energy decreases with increasing time, resulting in a relative error of almost 100 per cent after 5 kyr. Consequently, the total energy of the system (dashed black line) is not conserved for the LLE solution. As demonstrated in Fig. 7(a), the large error in evaluating the rotational energy is associated with determining the component m3, which is significantly larger than the value obtained using the NLE, in contrast to components m1 and m2, which are determined with high accuracy (see Fig. 2a). Fig. 7(b) shows that the predicted degree 2 shape of the Earth does not depend on whether we use the NLE or LLE, suggesting that the displacement obtained using the LLE is only slightly affected by the error in m3. This explains the good agreement obtained for the gravitational energy which is a function of the radial component of displacement (see eq. 36) and, unlike the rotational energy, does not depend directly on m3. Figure 7. View largeDownload slide (a) Time evolution of m3 for the studied benchmark case (cf. Fig. 2). Red points depict our NLE solution, dashed line is the LLE solution. Green line shows our solution when the spherical ice cap is imposed on the surface of the Earth only gradually (over the period of 5 yr). (b) Displacement components ur, 20 and ur, 21 relative to the hydrostatic equilibrium for the same case as above (continuous loading is omitted for clarity of the figure). Figure 7. View largeDownload slide (a) Time evolution of m3 for the studied benchmark case (cf. Fig. 2). Red points depict our NLE solution, dashed line is the LLE solution. Green line shows our solution when the spherical ice cap is imposed on the surface of the Earth only gradually (over the period of 5 yr). (b) Displacement components ur, 20 and ur, 21 relative to the hydrostatic equilibrium for the same case as above (continuous loading is omitted for clarity of the figure). Note that the difference between the NLE and LLE solution for m3 is not related to the free wobble, as illustrated by the green line in Fig. 7(a). The line represents a case in which the ice cap is being imposed only gradually onto the surface of the Earth, over the period of 5 yr, significantly reducing the amplitude of the induced wobble but leading to the same long term behaviour as for the instantaneous loading. In the light of the results presented in Section 3, the large error in component m3 obtained using the LLE may be viewed as surprising. As demonstrated in Fig. 2(b), both the NLE and LLE predict the same LOD variation ΔLOD given by eq. (17) (right), but m3 computed as m3 = −c33/C is still significantly affected by the neglect of the nonlinear terms in the LE (Fig. 7a). To explain this apparent paradox, we start from the formula (9.2.3) in Munk & MacDonald (1960)   $$\frac{\Delta \mathrm{LOD}}{\mathrm{LOD}}=-m_3=\frac{c_{33}}{C}\, ,$$ (40)and we describe the simplifications that were used in its derivation. The first equality in eq. (40) can be derived by linearization of the exact formula for ΔLOD:   $$\frac{\Delta \mathrm{LOD}}{\mathrm{LOD}}=\frac{2\pi \omega ^{-1}-2\pi \omega _0^{-1}}{2\pi \omega _0^{-1}}=\frac{\omega _0-\omega }{\omega }\cong \frac{\omega _0-\omega }{\omega _0} =1-\sqrt{m_1^2+m_2^2+(1+m_3)^2}\, .$$ (41)Using the binomial theorem and neglecting the higher-order terms, we obtain   $$\frac{\Delta \mathrm{LOD}}{\mathrm{LOD}}\cong -m_3+\frac{1}{2}\left(m_1^2+m_2^2+m_3^2\right)\, ,$$ (42)which, after neglecting the quadratic terms, gives the first equality in eq. (40). The second equality in eq. (40) comes from the LLE for component m3. In appendix B, we derive the following higher-order approximation:   $$-\!m_3 + \frac{1}{2}(m_1^2+m_2^2+m_3^2 ) \cong \frac{{c}_{33}}{C}\, .$$ (43)Combining eq. (43) with eq. (42) yields a more accurate version of the original equation (40):   $$\frac{\Delta \mathrm{LOD}}{\mathrm{LOD}}\cong -m_3 + \frac{1}{2}(m_1^2+m_2^2+m_3^2 ) \cong \frac{{c}_{33}}{C}\, .$$ (44)The definition of ΔLOD in eq. (17) thus includes second-order accurate evaluation of the LOD variation, while the equation m3 = −c33/C is only first-order accurate, which explains the disagreement between the LLE and NLE solutions in Fig. 7(a). The fact that the m3 component of the rotation vector is affected by a large error does not invalidate the results of previous studies using the LLE solution. As demonstrated above, the traditional equation for ▵LOD, eq. (17), used in these studies, is more accurate than the linearized equation for m3 and it gives the same LOD variation as the NLE equation (Fig. 2b), provided that excursions of the rotation axis are small. Our numerical tests (not presented in this paper) suggest that the shift of the rotation axis must be larger than 1 degree for the nonlinearity to significantly (by at least a few per cent) affect the resulting PM and LOD. The ice cap considered here is rather small and, since the characteristic time of equatorial bulge readjustment is inversely proportional to the size of the load (Ricard et al. 1993), several Myr would be needed to shift the pole by 1 degree. If we considered an unrealistic, ten times larger ice cap, the equatorial bulge readjustment would occur faster and the linearized PM and LOD solutions would differ from the nonlinear one by few per cent already after about 200 kyr. 4.5 Energy balance as a tool for testing GIA models The rotational feedback in the equation of motion has not been considered in a number of papers dealing with modelling of GIA on the Earth (for details, see Introduction in Martinec & Hagedoorn 2014). While such approximation is justifiable for the displacement component of degree 2 order 0 (ur, 20 has nearly identical evolution with and without the feedback), the displacement component of degree 2 and order 1 behaves in time completely differently for the two cases. A long term evolution of ur, 21 is governed by the balance between the viscoelastic relaxation of the model under the load and the viscoelastic readjustment of the equatorial bulge. For a homogeneous model, these two mechanisms perfectly cancel out each other for all times, as pointed out by Munk & MacDonald (1960), Section 6.4, and Wu & Peltier (1984). The effect of the rotational feedback has been recently discussed by Martinec & Hagedoorn (2014), who demonstrated that it may have a significant impact on the prediction of ur, 21 and sea level change. It should be emphasized that including the variations of the centrifugal potential into the field equations governing GIA-induced deformation of the Earth differs from including them into the gravity appearing in the sea level equation (i.e. the rotational feedback on the sea level equation), which has already been addressed in the previous studies (e.g. Milne & Mitrovica 1998; Peltier 1998) and which is not a focus of our work. Fig. 8 shows the time evolution of different types of energy when the rotation of the Earth is held constant, that is, the LE is not being solved at all. Fig. 9 compares the rotational and gravitational energies for the cases with and without the rotational feedback in the equation of motion, respectively. We can see that the gravitational energies do not differ significantly. This can be explained by the quadratic dependence of the gravitational energy on the radial component of the displacement (cf. eq. 35). While the hydrostatic value of ur, 20 reaches thousands of metres, the hydrostatic value of ur, 21 is equal to zero. Since the deformation induced by the ice cap reaches only few metres (Fig. 7b), the change of the gravitational energy of the rotating body is dominated by the change of ur, 20. In addition, as discussed above, ur, 20(t) is not affected by the rotational feedback significantly. Figure 8. View largeDownload slide Time evolution of different types of energy for a rotating model M3-L70-V01 loaded with a spherical ice cap at t = 0 (for details, see Spada et al. 2011). The rotation of the model is kept constant throughout the simulation. Figure 8. View largeDownload slide Time evolution of different types of energy for a rotating model M3-L70-V01 loaded with a spherical ice cap at t = 0 (for details, see Spada et al. 2011). The rotation of the model is kept constant throughout the simulation. Figure 9. View largeDownload slide Comparison of time evolutions of the rotational and gravitational energies from Figs 6 and 8. The energies from Fig. 6 are in colour, dashed black lines show the curves from Fig. 8. Figure 9. View largeDownload slide Comparison of time evolutions of the rotational and gravitational energies from Figs 6 and 8. The energies from Fig. 6 are in colour, dashed black lines show the curves from Fig. 8. The rotational energies, on the other hand, differ significantly for the two cases. With the rotational feedback in the equation of motion, the rotational energy evolves in time with the opposite sign than the gravitational energy in order to conserve the total energy. Without considering the rotational feedback (i.e. ω(t) = ω0), the rotational energy evolves in time similarly to the gravitational energy. This is because the increment of (ω0 · I · ω0)/2 has, similarly to the change of the gravitational energy, a positive linear dependence on the change of ur, 20 (see MacCullagh’s formula for c33 and insert for Φ20 from eq. 16). The detailed analysis of individual terms in the energy balance is given in Appendix A. Fig. 10 shows the results of such analysis for the two simulations shown in Figs 6 and 8, for which the total energy of the Earth’s model is not conserved. We can see that the energy imbalance is associated with the rotational energy, that is, the power of the centrifugal force is not equal to the rate of change of rotational energy. Figure 10. View largeDownload slide The energy integrals given in Section 4 (black lines) and their equivalents obtained by integrating the rate of mechanical work over the computational domain and time (coloured lines). Two cases are considered: (a) the rotational feedback is computed using the linearized Liouville equation (cf. Fig. 6), and (b) the rotational feedback is completely omitted (cf. Fig. 8). The numbers have the following meaning: 1—dissipated energy (dashed black and blue lines); 2—rotational energy (dotted black and red lines); 3 and 4—contributions to gravitational energy (dashed-dotted and green and yellow lines). For more details, see Appendix A and Fig. A1. Figure 10. View largeDownload slide The energy integrals given in Section 4 (black lines) and their equivalents obtained by integrating the rate of mechanical work over the computational domain and time (coloured lines). Two cases are considered: (a) the rotational feedback is computed using the linearized Liouville equation (cf. Fig. 6), and (b) the rotational feedback is completely omitted (cf. Fig. 8). The numbers have the following meaning: 1—dissipated energy (dashed black and blue lines); 2—rotational energy (dotted black and red lines); 3 and 4—contributions to gravitational energy (dashed-dotted and green and yellow lines). For more details, see Appendix A and Fig. A1. Another example of the application of the energy analysis given in Appendix A is related to the NLE solution presented in Sections 3 and 4.4. We run a numerical experiment where the original time step of 1 × 10−3 yr is made coarser as 5 × 10−2 yr. The results are shown in Fig. 11(a). Inspecting Fig. 11(a) we can see that not all the powers match the corresponding rates of energies. Again, it is the rotational energy that is not balanced. Hence, time step of 5 × 10−2 yr is too coarse in this case, which causes that the LE, but not eqs (1)–(3), is violated. In the parallel numerical experiment, the original spatial discretization of 460 radial nodes is reduced to 40 nodes only. Eqs (1)–(3) are now violated which is seen from the fact that the forces depending on the spatial derivatives of field variables do not match their corresponding energies. Figure 11. View largeDownload slide Panel (a) shows the case when the governing equations are solved using a too large time step, while panel (b) illustrates the effect of coarse radial discretization. The types of curves and the meaning of the numbers are the same as in Fig. 10. Figure 11. View largeDownload slide Panel (a) shows the case when the governing equations are solved using a too large time step, while panel (b) illustrates the effect of coarse radial discretization. The types of curves and the meaning of the numbers are the same as in Fig. 10. We hope that such a posteriori diagnostic of numerical solutions may be of interest also in other applications. For instance, in the case where the LE is employed for determining the Q-factor of the Chandler wobble, as suggested by Nakada & Karato (2012), the solution is numerically very sensitive to time-stepping. An independent test of the solution could thus prove helpful. 5 CONCLUSIONS In this paper, the temporal evolution of rotational, elastic, dissipative and gravitational energies is studied in the context of GIA. We have derived spectral formulae for computing the energies, overcoming the difficulty associated with the use of a simplified expression for the perturbation of gravitational potential. The energy evolution is then computed for several examples, demonstrating the importance of rotational and gravitational energies in the resulting balance. Viscous dissipation is proven relevant only for the damping of the Chandler wobble. By summing all energies for a rotating Earth loaded with a spherical ice cap, we show that the LLE is energetically inconsistent—the error resulting from an inaccurate evaluation of the m3 component of rotation vector. Careful evaluation of the energy balance for a rotating Earth loaded with an ice cap shows that the omission of the rotational feedback in the momentum equation leads to a significant violation of the energy conservation law, with both rotational and gravitational energy increasing in time on a time scale typical of the Earth’s glacial cycle. This finding underlines the importance of proper implementation of rotation variations in GIA modelling. Throughout the study we strictly use the Eulerian formulation for computing the deformation of the planet, and solve the governing equations directly in time domain. The equivalence with the Lagrangian formulation and Laplace domain approach (normal mode theory) is demonstrated numerically by comparing the Love numbers. The evaluation of the mechanical work of individual forces considered in the Eulerian formulation provides a versatile numerical tool, which can be used to verify solutions to the problems dealing with small deformations of planetary bodies. Acknowledgements This work has been supported by the Charles University grant SVV 260447/2017. We thank Volker Klemann and an anonymous reviewer for helpful comments that improved the manuscript. REFERENCES Benjamin D., Wahr J., Ray R.D., Egbert G.D., Desai S.D., 2006. Constraints on mantle anelasticity from geodetic observations, and implications for the J(2) anomaly, Geophys. J. Int. , 165( 1), 3– 16. Google Scholar CrossRef Search ADS   Cambiotti G., Ricard Y., Sabadini R., 2010. Ice age True Polar Wander in a compressible and non-hydrostatic Earth, Geophys. J. Int. , 183( 3), 1248– 1264. Google Scholar CrossRef Search ADS   Cambiotti G., Ricard Y., Sabadini R., 2011. New insights into mantle convection true polar wander and rotational bulge readjustment, Earth planet. Sic. Lett. , 310(3–4), 538– 543. Google Scholar CrossRef Search ADS   Cathles L., 1975. Viscosity of the Earth’s Mantle , Princeton Legacy Library, Princeton Univ. Press. Choblet G., Cadek O., Couturier F., Dumoulin C., 2007. OE DIPUS: a new tool to study the dynamics of planetary interiors, Geophys. J. Int. , 170( 1), 9– 30. Google Scholar CrossRef Search ADS   Floberghagen R., Fehringer M., Lamarre D., Muzi D., Frommknecht B., Steiger C., Pineiro J., da Costa A., 2011. Mission design, operation and exploitation of the gravity field and steady-state ocean circulation explorer mission, J. Geod. , 85( 11), 749– 758. Google Scholar CrossRef Search ADS   Gerya T., Yuen D., 2003. Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties, Phys. Earth planet. Inter. , 140( 4), 293– 318. Google Scholar CrossRef Search ADS   Golle O., Dumoulin C., Choblet G., Cadek O., 2012. Topography and geoid induced by a convecting mantle beneath an elastic lithosphere, Geophys. J. Int. , 189( 1), 55– 72. Google Scholar CrossRef Search ADS   Gurtin M.E., 1981. An Introduction to Continuum Mechanics , vol. 158 of Mathematics in Science and Engineering, Academic Press. Harris A., 1994. Tumbling asteroids, Icarus , 107( 1), 209– 211. Google Scholar CrossRef Search ADS   Jones M.N., 1985. Spherical Harmonics and Tensors for Classical Field Theory , vol. 2 of Applied and Engineering Mathematics Series, Research Studies Press. Joseph D.D., 1990. Fluid Dynamics of Viscoelastic Liquids , vol. 84 of Applied Mathematical Sciences, Springer-Verlag. Google Scholar CrossRef Search ADS   Kite E.S., Matsuyama I., Manga M., Perron J.T., Mitrovica J.X., 2009. True Polar Wander driven by late-stage volcanism and the distribution of paleopolar deposits on Mars, Phys. Earth planet. Inter. , 280(1–4), 254– 267. Google Scholar CrossRef Search ADS   Lau H.C.P., Mitrovica J.X., Austermann J., Crawford O., Al-Attar D., Latychev K., 2016. Inferences of mantle viscosity based on ice age data sets: radial structure, J. geophys. Res. , 121( 10), 6991– 7012. Google Scholar CrossRef Search ADS   Lefftz M., Legros H., Hinderer J., 1991. Non-linear equations for the rotation of a viscoelastic planet taking into account the influence of a liquid core, Celest. Mech. Dyn. Astron. , 52( 1), 13– 43. Google Scholar CrossRef Search ADS   Martinec Z., Hagedoorn J., 2014. The rotational feedback on linear-momentum balance in glacial isostatic adjustment, Geophys. J. Int. , 199( 3), 1823– 1846. Google Scholar CrossRef Search ADS   Milne G., Mitrovica J., 1998. Postglacial sea-level change on a rotating Earth, Geophys. J. Int. , 133( 1), 1– 19. Google Scholar CrossRef Search ADS   Mitrovica J., Peltier W., 1991. A complete formalism for the inversion of postglacial rebound data - resolving power analysis, Geophys. J. Int. , 104( 2), 267– 288. Google Scholar CrossRef Search ADS   Mitrovica J., Milne G., Davis J., 2001. Glacial isostatic adjustment on a rotating earth, Geophys. J. Int. , 147( 3), 562– 578. Google Scholar CrossRef Search ADS   Munk W.H., MacDonald G.J.F., 1960. Rotation of the Earth , Cambridge Univ. Press. Nakada M., Karato S.-I., 2012. Low viscosity of the bottom of the Earth’s mantle inferred from the analysis of Chandler wobble and tidal deformation, Phys. Earth planet. Inter. , 192, 68– 80. Google Scholar CrossRef Search ADS   Patočka V., 2013. Polar Wander prediction based on the solution of the Liouville equation, Master’s thesis , Charles University in Prague, Czech Republic, http://geo.mff.cuni.cz/theses/2013-Patocka-Mgr.pdf. Peltier W.R., 1974. The impulse response of a Maxwell Earth, Rev. Geophys. Space Phys. , 12, 649– 669. Google Scholar CrossRef Search ADS   Peltier W.R., 1998. Postglacial variations in the level of the sea: implications for climate dynamics and solid-earth geophysics, Rev. Geophys. , 36( 4), 603– 689. Google Scholar CrossRef Search ADS   Ricard Y., Spada G., Sabadini R., 1993. Polar wandering of a dynamic earth, Geophys. J. Int. , 113( 2), 284– 298. Google Scholar CrossRef Search ADS   Runcorn S., 1984. The primeval axis of rotation of the moon, Phil. Trans. R. Soc. , 313( 1524), 77– 83. Google Scholar CrossRef Search ADS   Sabadini R., Vermeersen B., Cambiotti G., 2016. Global Dynamics of the Earth: Applications of Viscoelastic Relaxation Theory to Solid-Earth and Planetary Geophysics , Modern Approaches in Geophysics, Springer. Google Scholar CrossRef Search ADS   Sabadini R., Yuen D., Boschi E., 1982. Polar wandering and the forced responses of a rotating, multilayered, viscoelastic planet, J. geophys. Res. , 87( NB4), 2885– 2903. Google Scholar CrossRef Search ADS   Soucek O., Hron J., Behounkova M., Cadek O., 2016. Effect of the tiger stripes on the deformation of Saturn’s moon Enceladus, Geophys. Res. Lett. , 43( 14), 7417– 7423. Google Scholar CrossRef Search ADS   Spada G. et al.  , 2011. A benchmark study for glacial isostatic adjustment codes, Geophys. J. Int. , 185( 1), 106– 132. Google Scholar CrossRef Search ADS   Tanaka Y., Klemann V., Okuno J., 2009. Application of a numerical inverse Laplace integration method to surface loading on a viscoelastic compressible earth model, Pure appl. Geophys. , 166(8–9), 1199– 1216. Google Scholar CrossRef Search ADS   Tapley B.D., Bettadpur S., Watkins M., Reigber C., 2004. The gravity recovery and climate experiment: mission overview and early results, Geophys. Res. Let. , 31, L09607, doi:10.1029/2004GL019920. Google Scholar CrossRef Search ADS   Tobie G., Cadek O., Sotin C., 2008. Solid tidal friction above a liquid water reservoir as the origin of the south pole hotspot on Enceladus, Icarus , 196( 2), 642– 652. Google Scholar CrossRef Search ADS   Vermeersen L., Sabadini R., Devoti R., Luceri V., Rutigliano P., Sciarretta C., Bianco G., 1998. Mantle viscosity inferences from joint inversions of Pleistocene deglaciation-induced changes in geopotential with a new SLR analysis and polar wander, Geophys. Res. Lett. , 25( 23), 4261– 4264. Google Scholar CrossRef Search ADS   Wu P., Peltier W., 1982. Viscous gravitational relaxation, Geophys. J. R. astr. Soc. , 70( 2), 435– 485. Google Scholar CrossRef Search ADS   Wu P., Peltier W., 1984. Pleistocene deglaciation and the Earth’s rotation—a new analysis, Geophys. J. R. astr. Soc. , 76( 3), 753– 791. Google Scholar CrossRef Search ADS   APPENDIX A: ENERGY BALANCE AND THE RATE OF MECHANICAL WORK The energies given in Section 4 can be matched with the rates of mechanical work done by individual terms in the momentum eq. (1),   $$\nabla \cdot \boldsymbol {\tau } - (\boldsymbol {u}\cdot \nabla \rho _0)\, \boldsymbol {g}_0 - \rho _0\nabla {\left(\Phi +\Psi \right)} + \rho _0\boldsymbol {g}_0 = 0 \, ,$$ (A1)where we substituted for the body force from eq. (8). Multiplication of eq. (A1) by the velocity vector, application of the tensor identity (∇ · τ) · v = ∇ · (τ · v) − τ: ∇v and subsequent integration over the regular computational domain Sa − Sb gives the balance equation of the rate of mechanical work:   $$\int _{{\rm S}_a\!-\!{\rm S}_b}\!\!\Big(\!-\boldsymbol {\tau }:\nabla \boldsymbol {v}+\nabla \cdot (\boldsymbol {\tau }\cdot \boldsymbol {v}) -(\boldsymbol {u}\cdot \nabla \rho _0)\boldsymbol {g}_0\cdot \boldsymbol {v} - (\rho _0\nabla \Phi )\cdot \boldsymbol {v}-(\rho _0\nabla \Psi )\cdot \boldsymbol {v} + \rho _0\boldsymbol {g}_0\cdot \boldsymbol {v}\Big )\, {\rm d}v = 0\, .$$ (A2)We now evaluate the integrals of the individual terms in eq. (A2) to better understand the relationship between the mechanical work and the different types of energy discussed in the main text. The resulting expressions can be used as an alternative to the energy integrals given in Section 4 and demonstrate the equivalence of the two approaches through a numerical test (Fig. A1). Figure A1. View largeDownload slide (a) Time derivatives of the energy integrals given in Section 4 (left) and their equivalents obtained by evaluating the rate of mechanical work over the computational domain (right). The numbers 1 and 2 correspond to deformation (blue) and rotational (red) powers, respectively. The gravitational power is split into individual terms derived in Section 4.1. Number 3 (green) denotes the contributions due to the deformation in a reference gravitational field of potential V0 while number 4 (yellow) corresponds to the power associated with incremental gravitational potential Φ. (b) Numerical test of the formulae in panel a performed for the benchmark case by Spada et al. (2011) (see also Sections 3 and 4.4). The coloured curves, obtained by time integration of the power formulae given in panel (a) (right), match the black curves computed directly from the corresponding energy formula. The colours and the numbers correspond to those in panel (a). Figure A1. View largeDownload slide (a) Time derivatives of the energy integrals given in Section 4 (left) and their equivalents obtained by evaluating the rate of mechanical work over the computational domain (right). The numbers 1 and 2 correspond to deformation (blue) and rotational (red) powers, respectively. The gravitational power is split into individual terms derived in Section 4.1. Number 3 (green) denotes the contributions due to the deformation in a reference gravitational field of potential V0 while number 4 (yellow) corresponds to the power associated with incremental gravitational potential Φ. (b) Numerical test of the formulae in panel a performed for the benchmark case by Spada et al. (2011) (see also Sections 3 and 4.4). The coloured curves, obtained by time integration of the power formulae given in panel (a) (right), match the black curves computed directly from the corresponding energy formula. The colours and the numbers correspond to those in panel (a). In the following, we consider that v≅∂u/∂t under the assumption of small deformations. In Appendices A4–A6, we use two modifications of the Reynolds transport theorem (e.g. Gurtin 1981, p. 91, eq. 8):   $$\frac{{\rm d}}{{\rm d} t}\int _{v(t)}\rho f\, {\rm d}v = \int _{v(t)}\rho \frac{\mathcal {D}f}{\mathcal {D}t}\, {\rm d}v\, ,$$ (A3)where v(t) is the total volume of the deformed body at time t and $$\mathcal {D}/\mathcal {D}t$$ denotes the material derivative,   $$\frac{\mathcal {D}}{\mathcal {D}t}=\frac{\partial }{\partial t}+\boldsymbol {v}\cdot \nabla \, .$$ (A4) A1 Term − τ: ∇v Since the body is incompressible (trace of ∇v is zero) and the Cauchy stress tensor is symmetric, τ: ∇v can be rewritten as   $$\boldsymbol {\tau }:\nabla \boldsymbol {v}=\boldsymbol {\tau }^d:\frac{1}{2}\!\left(\nabla \boldsymbol {v}+(\nabla \boldsymbol {v})^{\rm T}\right)\, .$$ (A5)The strain rate $$\frac{1}{2}\!\left(\nabla \boldsymbol {v}+(\nabla \boldsymbol {v})^{\rm T}\right)$$ can be expressed from the constitutive equation (3):   $$\frac{1}{2}\!\left(\nabla \boldsymbol {v}+(\nabla \boldsymbol {v})^{\rm T}\right) \cong \frac{\partial }{\partial t}\left[\frac{1}{2}\!\left(\nabla \boldsymbol {u}+(\nabla \boldsymbol {u})^{\rm T}\right)\right] = \frac{1}{2\mu }\frac{\partial \boldsymbol {\tau }^d}{\partial t}+\frac{\boldsymbol {\tau }^d}{2\eta }\, .$$ (A6)The integral of τ: ∇v then takes the form   $$\int _{{\rm S}_a\!-\!{\rm S}_b}\!\!\boldsymbol {\tau }:\nabla \boldsymbol {v}\, \mathrm{d}v \cong \int _{{\rm S}_a\!-\!{\rm S}_b}\!\!\left[\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{2\eta } + \frac{\boldsymbol {\tau }^d}{2\mu }:\frac{\partial \boldsymbol {\tau }^d}{\partial t}\right]\, \mathrm{d}v= \int _{{\rm S}_a\!-\!{\rm S}_b}\!\!\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{2\eta }\, \mathrm{d}v+ \frac{\rm d}{{\rm d}t}\int _{{\rm S}_a\!-\!{\rm S}_b}\!\! \frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{4\mu }\, \mathrm{d}v\, .$$ (A7)Assuming that the liquid core is inviscid and neglecting the dissipation and elastic energy within the small volumes corresponding to topographic anomalies, we can replace the integrals over Sa − Sb on the right-hand side of eq. (A7) by integrals over v(t), resulting in the relation:   $$\int _{{\rm S}_a\!-\!{\rm S}_b}\boldsymbol {\tau }:\nabla \boldsymbol {v}\, \mathrm{d}v \cong \int _{v(t)}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{2\eta }+ \frac{{\rm d}}{{\rm d} t}\int _{v(t)}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{4\mu }\, \mathrm{d}v=\boldsymbol {\dot{E}}_\mathrm{diss} + \boldsymbol {\dot{E}}_\mathrm{el}\, ,$$ (A8)where we used eqs (20) and (21) from Section 4. Eq. (A8) suggests that the power represented by the term τ: ∇v is dissipated as heat or stored as elastic deformation energy in the body. For numerical validation of the relationship A8, see Fig. A1(b), curve 1. A2 Term ∇ · (τ · v) The volume integral of ∇ · (τ · v) can be transformed into a surface integral using the Gauss theorem:   $$\int _{\mathrm{S}_a\!-\!\mathrm{S}_b}\nabla \cdot (\boldsymbol {\tau }\cdot \boldsymbol {v})\, \mathrm{d}v= \int _{\partial \mathrm{S}_a}\boldsymbol {v}\cdot \boldsymbol {\tau }\cdot \boldsymbol {e}_r\, \mathrm{d}s - \int _{\partial \mathrm{S}_b}\boldsymbol {v}\cdot \boldsymbol {\tau }\cdot \boldsymbol {e}_r\, \mathrm{d}s \, .$$ (A9)Evaluating τ · er from boundary conditions (9) and (10) and omitting the term σL (see the remark in the beginning of Section 4.1), we obtain   $$\int _{\mathrm{S}_a-\mathrm{S}_b}\nabla \cdot (\boldsymbol {\tau }\cdot \boldsymbol {v})\, \mathrm{d}v= \int _{\partial \mathrm{S}_a}\boldsymbol {v}\cdot \left(u_r[\rho _0]_a\boldsymbol {g}_0\right)\, \mathrm{d}s + \int _{\partial \mathrm{S}_b}\boldsymbol {v}\cdot \left(u_r[\rho _0]_b\boldsymbol {g}_0+p_{\rm c}\boldsymbol {e}_r\right)\mathrm{d}s\, ,$$ (A10)where pc is the pressure in the core, pc = −ρc(V0 + Φ + Ψ). The first surface integral on the right-hand side of eq. (A10) can be modified as   \begin{eqnarray} \int _{\partial \mathrm{S}_a}\boldsymbol {v}\cdot \left(u_r[\rho _0]_a\boldsymbol {g}_0\right)\, \mathrm{d}s &\cong & \int _{\partial \mathrm{S}_a}\frac{\partial \boldsymbol {u}}{\partial t}\cdot \left(u_r[\rho _0]_a\boldsymbol {g}_0\right)\, \mathrm{d}s \nonumber\\ &=& -[\rho _0]_a\, g_0(a) \int _{\partial \mathrm{S}_a}\frac{1}{2}\frac{\partial u^2_r}{\partial t} \, \mathrm{d}s =-\frac{1}{2}[\rho _0]_a\, g_0(a) \frac{\rm d}{{\rm d}t}\int _{\partial \mathrm{S}_a} u^2_r \, \mathrm{d}s\, . \end{eqnarray} (A11)Analogously,   $$\int _{\partial \mathrm{S}_b}\boldsymbol {v}\cdot \left(u_r[\rho _0]_b\boldsymbol {g}_0\right)\, \mathrm{d}s \cong -\frac{1}{2}[\rho _0]_b\, g_0(b) \frac{\rm d}{{\rm d}t}\int _{\partial \mathrm{S}_b} u^2_r \, \mathrm{d}s \, .$$ (A12)To summarize,   $$\int _{\mathrm{S}_a-\mathrm{S}_b}\nabla \cdot (\boldsymbol {\tau }\cdot \boldsymbol {v})\, \mathrm{d}v= -\frac{1}{2}\frac{\rm d}{{\rm d}t} \left([\rho _0]_a\, g_0(a) \int _{\partial \mathrm{S}_a} u^2_r \, \mathrm{d}s +[\rho _0]_b\, g_0(b) \int _{\partial \mathrm{S}_b} u^2_r \, \mathrm{d}s\right) +\int _{\partial \mathrm{S}_b}p_{\rm c}\boldsymbol {v}\cdot \boldsymbol {e}_r\, \mathrm{d}s\, ,$$ (A13)or, in terms of surface integrals   \begin{eqnarray} {\int _{\partial \mathrm{S}_a}\boldsymbol {v}\cdot \boldsymbol {\tau }\cdot \boldsymbol {e}_r\, \mathrm{d}s -\int _{\partial \mathrm{S}_b}\boldsymbol {v}\cdot (\boldsymbol {\tau }+p_{\rm c}\mathcal {I})\cdot \boldsymbol {e}_r\, \mathrm{d}s =- \frac{1}{2}\frac{\rm d}{{\rm d}t} \left([\rho _0]_a\, g_0(a) \int _{\partial \mathrm{S}_a} u^2_r \, \mathrm{d}s +[\rho _0]_b\, g_0(b) \int _{\partial \mathrm{S}_b} u^2_r\, \mathrm{d}s\right)}, \end{eqnarray} (A14)where $$\mathcal {I}$$ is the identity tensor. Comparison of eq. (A13) with eq. (36) suggests that the integral of ∇ · (τ · v) over the computation domain corresponds to the changes of the gravitational energy due to the deformations of boundaries in the gravitational field of potential V0. For the numerical confirmation of eq. (A14), see Fig. A1(b), curves 3a,b. It is worth noting that the integral of ∇ · (τ · v) over the whole body v(t) would be zero. Using the Gauss theorem and considering that the real surface is stress-free, we obtain   $$\int _{v(t)}\nabla \cdot (\boldsymbol {\tau }\cdot \boldsymbol {v})\, \mathrm{d}v=\int _{\partial {\rm S}}\boldsymbol {v}\cdot \boldsymbol {\tau }\cdot \boldsymbol {n}\, \mathrm{d}s=0\, ,$$ (A15)where ∂S now denotes the surface of the deformed body and n is the normal to the surface (see also Appendix A7). A3 Term − (u · ∇ρ0) g0 · v Considering that $$\boldsymbol {v}\!\cong \!\partial {\boldsymbol u}/\partial t$$, g0 = −g0er and the functions g0 and ρ0 do not depend on time, we obtain   \begin{eqnarray} {-\int _{\mathrm{S}_a-\mathrm{S}_b}(\boldsymbol {u}\cdot \nabla \rho _0)\, \boldsymbol {g}_0\cdot \boldsymbol {v}\, \mathrm{d}v \cong \int _{\mathrm{S}_a-\mathrm{S}_b}u_r\frac{\partial u_r}{\partial t}\frac{\mathrm{d}\rho _0}{\mathrm{d}r}\, g_0\, \mathrm{d}v = \frac{1}{2}\frac{\rm d}{{\rm d} t}\int _{\mathrm{S}_a-\mathrm{S}_b} u_r^2 g_0\frac{\mathrm{d}\rho _0}{\mathrm{d}r}\, \mathrm{d}v = \frac{1}{2}\frac{\rm d}{{\rm d} t}\int _{b}^a \frac{\mathrm{d}\rho _0}{\mathrm{d}r}\int _{\partial {\rm S}_r} g_0u_r^2\, \mathrm{d}s\, \mathrm{d}r}\, , \end{eqnarray} (A16)where the resulting integral can be identified with the rate of change of the gravitational energy due to volumetric density changes in the gravitational field of potential V0 (cf. the last integral in eq. (36); for the numerical comparison, see Fig. A1b, curve 3c). A4 Auxiliary formulae for gravitational and centrifugal potential Before we proceed to the remaining terms in eq. (A2) we derive two auxiliary formulae for the gravitational and centrifugal potentials valid for the whole volume v(t) of the deformed body at time t. We first consider the energy balance of a body in an inertial frame subject to gravitational body force −ρ∇V only, where V is the gravitational potential generated by the body. If the surface of the body is free,   $$\frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\left(\rho \epsilon +\frac{1}{2}\rho \boldsymbol {v}\cdot \boldsymbol {v}\right) \mathrm{d}v = -\int _{v(t)}\rho \nabla V\cdot \boldsymbol {v} \, \mathrm{d}v\, ,$$ (A17)where the expression on the left-hand side corresponds to the time change of the internal and kinetic energies, ε denotes the internal energy density, and the integral on the right-hand side is the total power of the gravitational force. Using the definition of material derivative eq. (A17) can be rearranged as   $$\frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\left(\rho \epsilon +\frac{1}{2}\rho \boldsymbol {v}\cdot \boldsymbol {v}\right)\, \mathrm{d}v + \int _{v(t)} \rho \frac{1}{2}\frac{\mathcal {D}V}{\mathcal {D}t}\mathrm{d}v = \int _{v(t)}\rho \left(\frac{\partial V}{\partial t}-\frac{1}{2}\frac{\mathcal {D}V}{\mathcal {D}t}\right)\, \mathrm{d}v\, .$$ (A18)The time derivative in the second integral on the left-hand side of eq. (A18) can be pulled in front of the integral using the modified Reynolds theorem, eq. (A3), and the integrand on the right-hand side can be reformulated by substituting for the material derivative of V:   $$\rho \left(\frac{\partial V}{\partial t}-\frac{1}{2}\frac{\mathcal {D}V}{\mathcal {D}t}\right)= \rho \left(\frac{\partial V}{\partial t}-\frac{1}{2}\frac{\partial V}{\partial t}-\frac{1}{2}\nabla V\cdot \boldsymbol {v}\right)= \frac{1}{2}\rho \left(\frac{\partial V}{\partial t}-\boldsymbol {v}\cdot \nabla V\right)\, .$$ (A19)Eq. (A18) then gives   $$\frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\left(\rho \epsilon +\frac{1}{2}\rho \boldsymbol {v}\cdot \boldsymbol {v}+\frac{1}{2}\rho V\right)\, \mathrm{d}v = \frac{1}{2}\int _{v(t)}\rho \left(\frac{\partial V}{\partial t}-\nabla V\cdot {\boldsymbol v}\right)\, \mathrm{d}v\, .$$ (A20)On the left-hand side we now have the sum of all energies in the system, including the gravitational energy. The total energy must be conserved which implies that the right-hand side of eq. (A20) is equal to zero. Consequently,   $$\int _{v(t)}\rho \nabla V\cdot \boldsymbol {v}\, \mathrm{d}v = \int _{v(t)}\rho \frac{\partial V}{\partial t}\, \mathrm{d}v\, .$$ (A21) A similar formula can be derived for the centrifugal potential Ψ. In this case, however, we must consider the energy balance in a rotating frame which is fixed with respect to the body and consider all relevant forces that act in the rotating frame. In analogy to eq. (A17), but this time not considering any gravitational forcing, we can write   $$\frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\left(\rho \epsilon +\frac{1}{2}\rho \boldsymbol {v}\cdot \boldsymbol {v}\right)\, \mathrm{d}v = \int _{v(t)}\left(-\rho \nabla \Psi -\rho \boldsymbol {\omega }\times \boldsymbol {v}-\rho \frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t}\times \boldsymbol {r}\right)\cdot \boldsymbol {v}\, \mathrm{d}v,$$ (A22)where the terms in parentheses on the right-hand side represent the centrifugal, Coriolis and Euler body force, respectively. The scalar product of the Coriolis force with the velocity vector is zero. The power of the Euler force can be rearranged using the vector algebra identity a · (b × c) = b · (c × a):   $$\int _{v(t)}\left(-\rho \frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t}\times \boldsymbol {r}\right)\cdot \boldsymbol {v}\, \mathrm{d}v=-\frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t}\cdot \int _{v(t)}\boldsymbol {r}\times (\rho \boldsymbol {v})\, \mathrm{d}v=-\frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t}\cdot \boldsymbol {h}\, ,$$ (A23)where h is the relative angular momentum. In the Tisserand frame (Munk & MacDonald 1960), h is equal to zero by definition, and the total power of the Euler force is thus zero in such frame. Following the same procedure as above (eqs A18–A20 where V is replaced by Ψ), we obtain   $$\frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\left(\rho \epsilon +\frac{1}{2}\rho \boldsymbol {v}\cdot \boldsymbol {v}+\frac{1}{2}\rho \Psi \right)\, \mathrm{d}v = \frac{1}{2}\int _{v(t)}\rho \left(\frac{\partial \Psi }{\partial t}-\nabla \Psi \cdot {\boldsymbol v}\right)\, \mathrm{d}v\, .$$ (A24)On the left-hand side we have the sum of all energies of the now considered system, when viewed from the inertial reference frame. The time derivative of this sum must again be zero, giving   $$\int _{v(t)}\rho \nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v = \int _{v(t)}\rho \frac{\partial \Psi }{\partial t}\, \mathrm{d}v\, .$$ (A25) A5 Term − (ρ0∇Φ) · v We first evaluate the integral of − (ρ0∇Φ) · v over the volume v(t) of the deformed body and then we transform it into an integral over the computational domain Sa − Sb. The potential V in auxiliary formula (A21) can be expressed as the sum of the reference potential V0, independent of time, and the increment Φ due to the deformation, V(t) = V0 + Φ(t). Considering that ∂V0/∂t = 0,   $$\int _{v(t)}\rho \nabla (V_0+\Phi )\cdot \boldsymbol {v}\, \mathrm{d}v = \int _{v(t)}\rho \frac{\partial \Phi }{\partial t}\, \mathrm{d}v\, .$$ (A26)The term ρ∇V0 · v can be moved to the right-hand side and replaced by $$\rho \mathcal {D}V_0/\mathcal {D}t$$:   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v = \int _{v(t)}\rho \left(\frac{\partial \Phi }{\partial t}-\frac{\mathcal {D}V_0}{\mathcal {D}t}\right)\, \mathrm{d}v\, .$$ (A27)Adding ∫v(t)ρ∇Φ · v dv to both sides of eq. (A27), dividing the equation by 2 and invoking the expression for the material derivative of Φ, we obtain   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v = \frac{1}{2}\int _{v(t)}\rho \frac{\mathcal {D}}{\mathcal {D}t}\left(\Phi -V_0 \right)\, \mathrm{d}v\, ,$$ (A28)where the time derivative on the right-hand can be pulled in front of the integral using the Reynolds theorem (A3):   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v = \frac{1}{2}\frac{{\rm d}}{{\rm d} t}\int _{v(t)}\rho \left(\Phi -V_0 \right)\, \mathrm{d}v\, .$$ (A29)The integral on the right-hand side of eq. (A29) is formally the same as in eq. (22) except that the sign of V0 is opposite. Thus, the integral can be evaluated in an analogous manner as described in Section 4.1, giving   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v = \frac{1}{2}\sum _{i=1}^N\frac{\mathrm{d}}{\mathrm{d}t}\int _{h_i}[\rho _0]_i\Phi \mathrm{d}v \cong \frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{1}{2}\sum _{i=1}^N[\rho _0]_i \int _{\partial {\rm S}_i}\Phi u_r\, \mathrm{d}s\right]$$ (A30)for a layered model and   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v \cong \frac{{\rm d}}{{\rm d}t}\left[\frac{1}{2}\left( [\rho _0]_a\int _{\partial {\rm S}_a}\!\!u_r\Phi \mathrm{d}s +[\rho _0]_b\int _{\partial {\rm S}_b}\!\!u_r\Phi \mathrm{d}s -\int _{b}^a \frac{\mathrm{d}\rho _0}{\mathrm{d}r}\int _{\partial {\rm S}_r}\!\! u_r\Phi \, \mathrm{d}s\, \mathrm{d}r\right)\right]$$ (A31)for a continuous density profile. The terms on the right-hand sides of eqs (A30) and (A31) correspond to the gravitational energy due to potential Φ induced by the deformation (cf. eqs 31, 34 and 35). Now we transform the integral on the left-hand side of eq. (A29) into an integral over the computational domain Sa − Sb:   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v \cong \int _{{\rm S}_a-{\rm S}_b}\rho _0\nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v+ \rho _{\rm c}\int _{{\rm S}_b}\nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v\, ,$$ (A32)where we omitted the second order term −δρ∇Φ (cf. eq. 8) and neglected the contributions of ρ∇Φ · v within the small volumes corresponding to topographic anomalies. The volume integral over Sb can be transformed into a surface integral using the Gauss theorem. Considering that ∇Φ · v = ∇ · (Φv) for an incompressible body, we get   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v \cong \int _{{\rm S}_a-{\rm S}_b}\rho _0\nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v+ \rho _{\rm c}\int _{\partial {\rm S}_b}(\Phi \boldsymbol {v})\cdot \boldsymbol {e}_r\, \mathrm{d}v\, .$$ (A33)Combining eqs (A33) and (A30) provides the sought relation between the power of ρ0∇Φ over the computational domain and the associated energy rate:   $$\int _{{\rm S}_a-{\rm S}_b}\rho _0\nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v+ \rho _{\rm c}\int _{\partial {\rm S}_b}(\Phi \boldsymbol {v})\cdot \boldsymbol {e}_r\, \mathrm{d}v \cong \frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{1}{2}\sum _{i=1}^N[\rho _0]_i \int _{\partial {\rm S}_i}\Phi u_r\, \mathrm{d}s\right].$$ (A34)An analogous expression is obtained for body with a continuous density profile by combining eqs (A33) and (A31). The validity of eq. (A34) is numerically demonstrated by curve 4 in Fig. A1(b). A6 Term − (ρ0∇Ψ) · v Using eq. (A25) and following the analogous procedure as above, we get   \begin{eqnarray} {\int _{v(t)}\rho \nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v = \int _{v(t)}\rho \frac{\partial \Psi }{\partial t}\, \mathrm{d}v = \frac{1}{2}\int _{v(t)}\rho \left(\frac{\partial \Psi }{\partial t}+\nabla \Psi \cdot \boldsymbol {v}\right)\, \mathrm{d}v = \frac{1}{2}\int _{v(t)}\rho \frac{\mathcal {D}\Psi }{\mathcal {D}t}\, \mathrm{d}v = \frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\frac{1}{2}\rho \Psi \, \mathrm{d}v = \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{2}\boldsymbol {\omega }\cdot \boldsymbol {I}\cdot \boldsymbol {\omega }\right)}\, . \end{eqnarray} (A35)The integral of ρ∇Ψ · v over v(t) can be expressed as (cf. eq. A32)   $$\int _{v(t)}\rho \nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v \cong \int _{{\rm S}_a-{\rm S}_b}\rho _0\nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v+ \rho _{\rm c}\int _{{\rm S}_b}\nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v\, .$$ (A36)Combining eqs (A35) and (A36) we obtain   $$\int _{{\rm S}_a-{\rm S}_b}\rho _0\nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v+ \rho _{\rm c}\int _{{\rm S}_b}\nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v= \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{2}\boldsymbol {\omega }\cdot \boldsymbol {I}\cdot \boldsymbol {\omega }\right)\, ,$$ (A37)see also Fig. A1(b), curve 2. Note that the derivations of eqs (A34) and (A37) rely on the auxiliary formulae from Appendix A4, and as such they are only valid if the energy of the studied system is conserved. In energetically ill-posed problems eqs (A34) and (A37) must not be satisfied, as illustrated in Section 4.5 in the main text (cf. Fig. 10). A7 Term ρ0g0 · v Since ρ0 and g0 are spherically symmetric and the body is incompressible ($$\int _{\partial {\rm S}_r}\boldsymbol {v}\cdot \boldsymbol {e}_r\, {\rm d}s\!=\!0$$), the integral of the last term in eq. (A1) over the computational domain is equal to zero:   $$\int _{{\rm S}_a-{\rm S}_b}\rho _0\boldsymbol {g}_0\cdot \boldsymbol {v}\, \mathrm{d}v=-\int _{{\rm S}_a-{\rm S}_b}\rho _0 g_0 \boldsymbol {v}\cdot \boldsymbol {e}_r\, \mathrm{d}v=-\int _b^a\rho _0g_0\int _{\partial {\rm S}_r}\boldsymbol {v}\cdot \boldsymbol {e}_r\, {\rm d}s\, {\rm d}r=0\, .$$ (A38)Note, however, that the integral of ρg0 · v over the real volume of the body v(t) is non-zero, even if the density is constant within the mantle shell. For a body with ∇ρ0 equal to zero within the shell the integral is equal to the right-hand side of eq. (A10):   $$\int _{v(t)}\rho \boldsymbol {g}_0\cdot \boldsymbol {v}\, \mathrm{d}v \cong \int _{\partial \mathrm{S}_a}\boldsymbol {v}\cdot \left(u_r[\rho _0]_a\boldsymbol {g}_0\right)\, \mathrm{d}s + \int _{\partial \mathrm{S}_b}\boldsymbol {v}\cdot \left(u_r[\rho _0]_b\boldsymbol {g}_0+p_{\rm c}\boldsymbol {e}_r\right)\mathrm{d}s\, \, .$$ (A39)This is because in this integral there are contributions from the non-zero Eulerian density increments δρ = ±[ρ0]a, b in the vicinity of the boundaries (see Fig. 3). If the computational domain is chosen regular (Sa − Sb) as in the present study, these volumetric density increments are condensed into surface densities and we encounter them in Appendix A2 in the form of boundary tractions (9) and (10). For a body with ∇ρ0 ≠ 0 inside the mantle shell the term ρg0 · v also includes a contribution from the internal Eulerian density increment − u · ∇ρ0, which was discussed in Appendix A3. A8 Remark Inspection of the left-hand side of Fig. A1(a) suggests that the formulae used to evaluate the energy balance in the main text can be replaced by the corresponding formulae for the rates of the mechanical work evaluated over the computational domain Sa − Sb. However, there are three integrals over ∂Sb on the right-hand side of Fig. A1(a) which are related to the body force −ρc∇(V0 + Φ + Ψ) in the core. The sum of these terms is equal to zero,   $$-\int _{\partial {\rm S}_b}(\rho _{\rm c}\Psi \boldsymbol {v})\cdot \boldsymbol {e}_r\, {\rm d}s -\int _{\partial {\rm S}_b}p_{\rm c}\boldsymbol {v}\cdot \boldsymbol {e}_r\, {\rm d}s -\int _{\partial {\rm S}_b}(\rho _{\rm c}\Phi \boldsymbol {v})\cdot \boldsymbol {e}_r\, {\rm d}s=0\, ,$$ (A40)since pc = −ρc(V0 + Φ + Ψ) and   $$\int _{\partial {\rm S}_b}\rho _{\rm c}V_0\boldsymbol {v}\cdot \boldsymbol {e}_r\, {\rm d}s= \rho _{\rm c}V_0\int _{\partial {\rm S}_b}\boldsymbol {v}\cdot \boldsymbol {e}_r\, {\rm d}s=0\, .$$ (A41) APPENDIX B: SECOND-ORDER CORRECTION TO LLE The third component of the NLE is   \begin{eqnarray} 0 = \frac{\mathrm{d}}{{\rm d}t}\left(\boldsymbol {I}\cdot \boldsymbol {\omega }\right)_3 + \varepsilon _{3jk}\omega _j\left(\boldsymbol {I}\cdot \boldsymbol {\omega }\right)_k\, , \end{eqnarray} (B1)where ε is the Levi-Civita symbol. Using the traditional notation,   $$\boldsymbol {\omega }=\boldsymbol {\omega }_0 + |\boldsymbol {\omega }_0|\boldsymbol {m} \quad {\rm {and}} \quad \boldsymbol {I}={\left(\begin{array}{ccc}A+c_{11} & c_{12} & c_{13} \\ c_{12} & A + c_{22} & c_{23} \\ c_{13} & c_{12} & C+c_{33} \end{array}\right)} ,$$ (B2)with ω0 equal to [0, 0, ω0] and cij, mi denoting small quantities, and inserting from eq. (B2) to eq. (B1), we obtain   \begin{eqnarray} 0 &\cong & \frac{\mathrm{d}}{{\rm d}t}[c_{33} + c_{13} m_1 + c_{23} m_2 + c_{33} m_3 + C m_3 ] + \omega _0 m_1 c_{23} - \omega _0 m_2 c_{13} \nonumber \\ &=& (\dot{c}_{33} + \dot{c}_{13} m_1 + \dot{c}_{23} m_2 + \dot{c}_{33} m_3 + c_{13} \dot{m}_1 + c_{23} \dot{m}_2 + c_{33} \dot{m}_3 + C \dot{m}_3 ) + \omega _0(m_1 c_{23} - m_2 c_{13}), \end{eqnarray} (B3)where third order terms are neglected, the dot abbreviates time derivative and $$\dot{C}$$ and $$\dot{\omega }_0$$ are zero, since both C and ω0 are constants. Recalling the LLE for m1 and m2 (Munk & MacDonald 1960)   \begin{eqnarray} \omega _0(C-A)\left(m_1 - \frac{A\, \dot{m}_2}{(C-A)\omega _0}\right) = \omega _0 c_{13} + \dot{c}_{23} \end{eqnarray} (B4)  \begin{eqnarray} \omega _0(C-A)\left(m_2 + \frac{A\, \dot{m}_1}{(C-A)\omega _0}\right) = \omega _0 c_{23} - \dot{c}_{13}, \end{eqnarray} (B5)we can insert for $$\dot{c}_{23}$$ and $$\dot{c}_{13}$$ in eq. (B3) from eqs (B4) and (B5), introducing only a third order error, because $$\dot{c}_{23}$$ and $$\dot{c}_{13}$$ are both multiplied by a small quantity in eq. (B3). This leads to   \begin{eqnarray} 0 &\cong & \dot{c}_{33} + C \dot{m}_3 - A m_1 \dot{m}_1 - A m_2 \dot{m}_2 + \dot{c}_{33} m_3 + \omega _0(m_1 c_{23} - m_2 c_{13}) \nonumber \\ &&+\,\omega _0[ c_{23}m_1 - (C-A)m_1 m_2 - c_{13}m_2 + (C-A)m_1 m_2 + c_{13} \dot{m}_1 + c_{23} \dot{m}_2 + c_{33} \dot{m}_3] \, . \end{eqnarray} (B6)The fifth term on the first line of eq. (B6) can be expressed using the LLE for m3, $$\dot{c}_{33} = - C \dot{m}_3$$, introducing only a third order error:   \begin{eqnarray} 0 &\cong & \dot{c}_{33} + C \dot{m}_3 - A m_1 \dot{m}_1 - A m_2 \dot{m}_2 - C \dot{m}_3 m_3 + \omega _0[ 2 c_{23}m_1 - 2 c_{13}m_2 + c_{13} \dot{m}_1 + c_{23} \dot{m}_2 + c_{33} \dot{m}_3] \, . \end{eqnarray} (B7)This can be rewritten as   \begin{eqnarray} 0 &\cong & \dot{c}_{33} + C {\left( m_3 - \frac{m_1^2+m_2^2+m_3^2}{2} \right)}^{\boldsymbol {\cdot }} + \omega _0[ 2 c_{23}m_1 - 2 c_{13}m_2 + c_{13} \dot{m}_1 + c_{23} \dot{m}_2 + c_{33} \dot{m}_3] + (C-A)(\dot{m_1}m_1+\dot{m_2}m_2) , \end{eqnarray} (B8)which provides a second-order correction to the vertical component of the LLE. The second equality in eq. (42) is obtained when the relative importance of the second-order terms in eq. (B8) are assessed: the terms on the second line are negligible due to the multiplication by factors ω0 and C − A, which are both much smaller than factor C. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Energy balance of glacial isostatic adjustment: importance of the rotational feedback

, Volume 212 (2) – Feb 1, 2018
21 pages

Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx469
Publisher site
See Article on Publisher Site

### Abstract

Abstract Understanding the feedback between the glacial isostatic adjustment (GIA) and the Earth’s rotation is important for an accurate prediction of sea level changes induced by climate and tectonic processes. Here we consider a simple, four-layer incompressible Earth model, recently used for a benchmark of GIA codes to estimate the accuracy of the linearized Liouville equation (LE) and to demonstrate that models with an incomplete or missing rotational feedback violate the principle of energy conservation. First, we compute GIA on a rotating Earth by solving the equation of motion coupled with LE in its full nonlinear form. By comparing the nonlinear LE solution with the traditional linearized one, we find that the radial component of the angular velocity vector is inaccurate in the latter case, with an error exceeding 10 per cent already after 1 kyr of evolution. To understand the cause of this discrepancy, we investigate the time evolution of different kinds of energy involved in the process. While the rotational, elastic and dissipative energies are straightforward to compute, the formula for the gravitational energy contains an integral that requires a careful, higher-order accurate evaluation of the gravitational potential perturbation. We circumvent this problem by transforming the integral into a different one, formulated in terms of displacement instead of potential. We find that the solution of the linearized LE equation does not conserve the energy, since, in the linearized case, the rate of change of the rotational energy is not equal to the power of the centrifugal force. We also compute the energy balance of GIA on a constantly rotating Earth, and demonstrate the importance of the rotational feedback in the equation of motion. The formalism derived in this study allows a detailed examination of the energy balance for a rotating, incompressible planetary body deformed by a surface load. As such, it may not only serve as a reliable tool for a posteriori testing of GIA numerical solutions but it can also be used in different planetary science applications. Elasticity and anelasticity, Earth rotation variations, Loading of the Earth, Sea level change, Numerical modelling 1 INTRODUCTION Glacial isostatic adjustment (GIA) has been investigated within the geophysical community for many decades, recently gaining attention due to the increasing precision of geodetic measurements obtained from the GRACE (Tapley et al. 2004) and GOCE (Floberghagen et al. 2011) satellite missions and the growing interest in sea level changes. The response of the Earth to surface loading was used primarily to constrain the viscosity profile of the mantle (Cathles 1975). However, the depth to which viscosity can be inferred by measuring postglacial rebound is limited by the lateral extent of the past glaciers, leaving the lowermost mantle poorly resolved (e.g. Mitrovica & Peltier 1991). The motion of the rotation axis of the Earth, induced by mass redistribution that accompanies the periodic accumulation and melting of ice during glacial cycles, can be used to overcome this limitation. Assuming that the ice load history is known, lower-mantle viscosity can be constrained by fitting the observed rate of the secular drift (e.g. Sabadini et al. 1982; Wu & Peltier 1984; Vermeersen et al. 1998; Lau et al. 2016). Also worth noting is the pioneering attempt of Nakada & Karato (2012) to infer the lower-mantle viscosity by comparing the Q-factor of the predicted Chandler wobble with the observed value (e.g. Benjamin et al. 2006). It was only later that the perturbation of the centrifugal force, inherent to any change in the rotation vector, was also included in modelling GIA-related suite of observables. While its effect on the sea level equation has been extensively studied (e.g. Milne & Mitrovica 1998; Mitrovica et al. 2001), its impact on the GIA-induced deformation of the Earth has only been addressed recently by Martinec & Hagedoorn (2014). Precise computation of GIA-induced polar motion (PM) thus plays an important role in numerical simulations of postglacial rebound related phenomena. Following Munk & MacDonald (1960) it has become common to use the linearized form of the Liouville equation (LE) to obtain the PM solution. The error introduced by the linearization is not straightforward to estimate and depends on the amplitudes of the wander that the rotation axis experiences. In this paper, we couple the fully nonlinear LE (NLE) with the equations governing small deformations of a radially symmetric, self-gravitating incompressible spherical shell. By comparing the NLE solution with the linearized one we quantitatively evaluate the error due to neglecting higher-order terms in the linearized LE (LLE). In particular, we estimate the error arising from decoupling the equations for PM and the length of day (LOD) variation, which is an intrinsic property of the linearized approach. The accuracy of the LLE is tested for the Earth model M3-L70-V01 loaded with a spherical ice cap, a community benchmark by Spada et al. (2011). Another way to assess the accuracy of the linearized solution, without having to solve the NLE, is to check whether the total energy of the Earth’s model is conserved. The energy balance involves the rotational, elastic and dissipative energies, which are straightforward to evaluate, and the gravitational energy, which requires a higher-order accurate evaluation of the gravitational perturbation potential. In Section 4.1, we solve this problem by deriving an equivalent integral but with the displacement in the integrand, allowing the use of a simple spectral formula for the geoid (e.g. Choblet et al. 2007). We then compute the energy balance for various test examples relevant to GIA studies. In Section 4.2, we load a non-rotating Earth with a spherical ice cap. To assess the energetic importance of free wobble, triggered by a hypothetical rapid accumulation of a surface load on a rotating Earth, we separately study this phenomenon in Section 4.3. In Section 4.4, we analyse the energy balance of a rotating Earth loaded with a spherical ice cap, that is, the benchmark case by Spada et al. (2011). We show that the LLE solution does not conserve the energy and discuss the consequences of this finding for the prediction of PM and LOD. In Section 4.5, we consider the same surface loading, but on a constantly rotating Earth. A large deviation from the energy conservation is obtained in such case, proving the importance of including rotational feedback in GIA modelling. All results discussed in Section 4 are obtained for the four-layer, incompressible Earth model M3-L70-V01 used already in the benchmark of GIA codes (Spada et al. 2011). Although this structural model is relatively simple in comparison with recent models inferred from ice age data sets (e.g. Lau et al. 2016), it allows our predictions of PM and LOD to be compared with those in Spada et al. (2011) and the energy curves presented in Section 4 to be easily reproduced. The energy balance analysis introduced in this paper does not include the effect of compressibility. As shown by Cambiotti et al. (2010), compressibility of the real Earth can have non-negligible effect on both GIA and PM. Compressible models are more deformable than the incompressible ones, which makes the readjustment of the rotational bulge quicker, effectively reducing the GIA-induced PM. This reduction, however, depends on lower mantle viscosity and is almost negligible for values higher than 1022 Pa s, which are usually needed to stabilize the rotation of the Earth on geological time scales (Ricard et al. 1993). To compute small deformations of a rotating, gravitationally pre-stressed Earth we use the Eulerian formulation of the governing equations, instead of applying the traditional Lagrangian approach (e.g. Wu & Peltier 1982). The Eulerian formulation has appeared in the literature several times in the past decade (e.g. Tobie et al. 2008; Golle et al. 2012; Soucek et al. 2016), but only for spherical shells with a constant density. We extend it for the case where the internal density is a continuous piecewise linear function and we numerically demonstrate that the solution is equivalent to the Lagrangian one. Using the Eulerian formulation we demonstrate a one-by-one correspondence between individual forces in the equation of motion and the terms describing the gravitational, rotational, elastic and dissipative energies. Such a correspondence provides a powerful tool when diagnosing a numerical solution. While the evaluation of the sum of all involved energies reveals only whether or not the solution is numerically correct, the term-by-term correspondence can be used to directly identify a potential source of energy imbalance. In the present paper, we do not use the normal mode theory (e.g. Peltier 1974) or complex contour integration (e.g. Tanaka et al. 2009; Sabadini et al. 2016), as we solve all equations in the time domain, making the coupling between the equation of motion and the LE straightforward. It is worth noting that the use of the NLE does not require the amplitudes of the modelled PM to be small and the time scales on which mass redistribution occurs can be arbitrary. Our approach thus broadens the variety of processes that can be addressed compared to the traditional techniques: the linearized approach is limited to PM with small amplitudes, and the viscous quasi-fluid approximation, introduced by Lefftz et al. (1991) and Ricard et al. (1993), is valid only for mass redistribution occurring on the time scale of a few million years (see also Cambiotti et al. 2011). The rapid formation of craters and volcanoes that can subsequently shift the rotation axis of a planetary body by tens of degrees (Runcorn 1984; Kite et al. 2009), and a number of asteroids that experience a large-amplitude wobbling (Harris 1994), are examples of the problems that could newly be tackled with the tool we present. 2 GOVERNING EQUATIONS We consider small deformations of a hydrostatically pre-stressed spherical shell with radially dependent reference density profile ρ0(r) and a homogeneous fluid core. The material of the spherical shell is assumed to behave like an incompressible Maxwell-type viscoelastic fluid. The time evolution of displacement and stress can be obtained by integrating the following set of partial differential equations (e.g. Tobie et al. 2008):   $$\nabla \cdot \boldsymbol {\tau } + \boldsymbol {f} = 0 \, ,$$ (1)  $$\nabla \cdot \boldsymbol {u}=0\, ,$$ (2)  $$\boldsymbol {\tau }^d-\mu (\nabla \boldsymbol {u}+(\nabla \boldsymbol {u})^\mathrm{T})=-\frac{\mu }{\eta }\int _0^t\boldsymbol {\tau }^d \mathrm{d}t^{\prime } \, ,$$ (3)where τ is the Cauchy stress tensor, τd is its deviatoric part, f is the body force, u is the displacement, the superscript T denotes transposition of a tensor, μ is the elastic shear modulus, η is the viscosity, and t is the time. At time t = 0 the integral on the right-hand side of eq. (3) is zero and the spherical shell behaves like an elastic solid. Throughout the paper, we strictly use the Eulerian description of the problem. All variables are expressed as functions of the instantaneous position of the particle in the deformed body. A particle occupying the position r at time t would occupy the position r − u(r, t) if all forces that cause the deformation were set to zero at all times. The body force f is given as the product of density and the negative gradient of gravity potential. The gravity potential has three components: (i) the reference gravitational potential V0, (ii) the centrifugal potential Ψ due to the rotation of the body, and (iii) the perturbation Φ of the gravitational potential V0 due to the deformation of the body. V0 is the gravitational potential of a sphere of radius a with density profile ρ0(r). Inside the sphere (r ≤ a)   $$V_0(r)=-\frac{4\pi G}{r}\int _0^{r} \rho _0(r^{\prime }){r^{\prime }}^2\, \mathrm{d}r^{\prime }- 4\pi G\int _r^{a} \rho _0(r^{\prime })r^{\prime }\, \mathrm{d}r^{\prime },$$ (4)where G is the universal gravitational constant. The centrifugal potential Ψ is given as   $$\Psi = \frac{1}{2}\left((\boldsymbol {\omega }\cdot \boldsymbol {r})^2-|\boldsymbol {\omega }|^2|\boldsymbol {r}|^2 \right),$$ (5)where ω is the angular velocity vector which has the direction of the rotation axis and the magnitude equal to the angular speed of the body. As the sphere deforms, the initial density ρ0(r) changes to density ρ(r, t) which can be expressed as a sum of density ρ0(r) and Eulerian density increment δρ(r, t). If the body is incompressible and the density ρ0 is a continuous function,   $$\delta \rho (\boldsymbol {r},t) = \rho (\boldsymbol {r},t) - \rho _0(r) = \rho _0(|\boldsymbol {r}-\boldsymbol {u}(\boldsymbol {r},t)|,t) - \rho _0(r) = -\boldsymbol {u}(\boldsymbol {r},t)\cdot \nabla \rho _0(r) + O(|\boldsymbol {u}|^2),$$ (6)where the displacement u describes the deformation of the body. The Eulerian density increment δρ is non-zero also at the outer surface and at the internal density interfaces where it is equal to the density jump at the respective interface. Using δρ we can express the perturbation potential Φ in an integral form:   $$\Phi (\boldsymbol {r},t)=-G\int _{v(t)} \frac{\delta \rho (\boldsymbol {r}^{\prime },t)}{|\boldsymbol {r}-\boldsymbol {r}^{\prime }|}\, \mathrm{d}v^{\prime },$$ (7)where v(t) is the volume occupied by the body (including the core) at time t. Neglecting the terms O(|u|2) and −δρ∇(Φ + Ψ), we obtain the following expression for the body force f:   $$\boldsymbol {f} = -\rho \nabla (V_0+\Phi +\Psi ) \cong - (\boldsymbol {u}\cdot \nabla \rho _0)\boldsymbol {g}_0 - \rho _0\nabla \left(\Phi +\Psi \right) + \rho _0\boldsymbol {g}_0,$$ (8)where g0 = −∇V0. The laterally homogeneous and static contribution ρ0g0 is counteracted by the hydrostatic pre-stress p0(r), satisfying the equation − ∇p0 + ρ0g0 = 0. Eqs (1)–(3) are solved in a spherical shell with outer radius a and inner radius b. The boundary condition at the outer surface is obtained from the force equilibrium condition, taking into account the pressure due to the deformation-induced topography (e.g. Soucek et al. 2016) and the presence of a surface load:   $$\boldsymbol {\tau }\cdot \boldsymbol {e}_r=\left(u_r[\rho _0]_a+\sigma ^\mathrm{L}\right)\boldsymbol {g}_0,$$ (9)where er is the radial unit vector, ur is the radial component of displacement u, [ρ0]a = ρ0(a) is the density jump at the surface, and σL is the surface mass density of the prescribed surface load. A similar condition can be imposed at the bottom boundary, on Earth corresponding to the core–mantle boundary:   $$-\!\boldsymbol {\tau }\cdot \boldsymbol {e}_r = u_r[\rho _0]_b\, \boldsymbol {g}_0-\rho _\mathrm{c}(V_0+\Phi +\Psi )\boldsymbol {e}_r,$$ (10)where [ρ0]b is the density jump at the bottom boundary and −ρc(V0 + Φ + Ψ) is the hydrostatic pressure acting on the boundary due to the contact with a liquid core of density ρc. To compute the temporal evolution of the angular velocity vector ω with respect to the body-fixed Tisserand frame (Munk & MacDonald 1960) we solve the LE with zero external torque,   $$-\!\boldsymbol {I}\cdot \frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t} = \frac{\mathrm{d}\boldsymbol {I}}{\mathrm{d}t}\cdot \boldsymbol {\omega }+\boldsymbol {\omega }\times (\boldsymbol {I}\cdot \boldsymbol {\omega }),$$ (11)where I is the time-dependent tensor of inertia. Eqs (1)–(10) are coupled with eq. (11) both through the displacement field u, which is needed to compute the inertia tensor I, and through the centrifugal potential Ψ, which depends on the angular velocity vector ω. 3 NUMERICAL IMPLEMENTATION The LE is a set of ordinary differential equations for three unknown components of vector ω as functions of time. It can be rewritten as   $$\frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t} = \mathcal {F}(\boldsymbol {I},\boldsymbol {\omega })\, ,$$ (12)where   $$\mathcal {F}=-\boldsymbol {I}^{-1}\cdot \left(\frac{\mathrm{d}\boldsymbol {I}}{\mathrm{d}t}\cdot \boldsymbol {\omega }+\boldsymbol {\omega }\times (\boldsymbol {I}\cdot \boldsymbol {\omega })\right).$$ (13)Eq. (12) is solved numerically using the fifth order accurate Adams-Bashforth multistep method with a time step of 10−3 yr. Each evaluation of $$\mathcal {F}$$ requires the calculation of the tensor of inertia I and its time derivative. The perturbation of the tensor of inertia is computed by MacCullagh’s formula from the gravitational potential Φ and the reference density profile ρ0 (e.g. Patočka 2013, eqs 2.12),   $$\boldsymbol {I} = I_0\, \mathcal {I} + \sqrt{\frac{5}{\pi }}\, \frac{a^3}{6\, G} \, \left[ \, \left( \begin{array}{c@{\quad}c@{\quad}c}-\Phi _{20} & 0 & 0\\ 0 & -\Phi _{20} & 0\\ 0 & 0 &2 \Phi _{20} \end{array}\right) \, - \, \sqrt{6} \, \left( \begin{array}{c@{\quad}c@{\quad}c}-\mathrm{Re}\, \Phi _{22} & \mathrm{Im}\, \Phi _{22} & \mathrm{Re}\, \Phi _{21} \\ \mathrm{Im}\, \Phi _{22} & \mathrm{Re}\, \Phi _{22} & -\mathrm{Im}\, \Phi _{21} \\ \mathrm{Re}\, \Phi _{21} & -\mathrm{Im}\, \Phi _{21} & 0 \\ \end{array}\right) \right].$$ (14)Here $$I_0:=\frac{8}{3}\pi \int ^a_{0} r^4 \rho _{0} \mathrm{d}r$$, $$\mathcal {I}$$ is the identity tensor, and Φℓm are the coefficients of the spherical harmonic expansion of potential Φ,   $$\Phi (r,\vartheta ,\varphi )=\sum _{\ell =0}^{\ell _{\rm max}}\sum _{m=-\ell }^{\ell }\Phi _{\ell m}(r)Y_{\ell m}(\vartheta ,\varphi )\, ,$$ (15)where Yℓm are the fully normalized spherical harmonics of degree ℓ and order m (e.g. Jones 1985), and ℓmax is the cut-off degree. The potential Φ depends on the deformation of the body and we describe its calculation in the next paragraph. The term $$\mathrm{d}\boldsymbol {I}/\mathrm{d}t$$ is then computed numerically using the second order accurate mid-point finite difference scheme. The spectral algorithm for the numerical solution of eqs (1)–(10) for a given centrifugal potential Ψ is outlined in Tobie et al. (2008), identities needed to construct the matrix of the discretized problem are summarized in the appendix to Golle et al. (2012). The potential Φ depends on u and vice versa. It is computed iteratively by expanding the integral (7) into a spherical harmonics series, while condensing topographies into surface mass densities (see eq. 35 in Choblet et al. 2007):   \begin{eqnarray} \Phi _{\ell m}(r) = \frac{-4\pi Gr}{2\ell +1}\Bigg ( [\rho _0]_{b}\, u_{r,\ell m}(b) \left(\frac{b}{r}\right)^{\ell +2} + [\rho _0]_{a}\, u_{r, \ell m}(a) \left(\frac{r}{a}\right)^{\ell -1}- \int _{b}^r \left(\frac{r^{\prime }}{r}\right)^{\ell +2}{u}_{r,\ell m}(r^{\prime })\frac{\mathrm{d}\rho _0}{\mathrm{d}r^{\prime }}\, \mathrm{d}r^{\prime } - \int _r^{a} \left(\frac{r}{r^{\prime }}\right)^{\ell -1}{u}_{r,\ell m}(r^{\prime })\frac{\mathrm{d}\rho _0}{\mathrm{d}r^{\prime }}\, \mathrm{d}r^{\prime } \Bigg )\, , \end{eqnarray} (16)where we substituted for δρ from eq. (6) on the second line of eq. (16). The discretization in the radial direction is performed by the finite difference method on a staggered grid with constant spacing (e.g. Gerya & Yuen 2003). Our approach differs from the algorithm described in Tobie et al. (2008) by employing a higher order Crank-Nicholson integration scheme for evaluating the time integral in eq. (3). The accuracy of the numerical method used to solve the governing eqs (1)–(10) was carefully tested against the traditional Lagrangian solution. Fig. 1 shows the relative difference between the loading and tidal degree 2 Love numbers $$k^L_e$$, $$k^T_e$$, $$k^L_f$$ and $$k^T_f$$, computed numerically using the method described above, and those computed semi-analytically by Spada et al. (2011) for a viscoelastic Earth model M3-L70-V01 (for description of the model, see table 3 therein). The difference is plotted as a function of the radial resolution considered in the numerical solution in which the density profile is approximated by a continuous piecewise linear function. Inspection of Fig. 1 shows that the relative error of the numerical method converges to zero with increasing resolution, decreasing below 10−4 when the number of equally spaced radial nodes is larger than 5000. The error does not decrease monotonically because it depends on how well the equidistant discretization matches the positions of the density interfaces. In Section 4.4, we will use 460 radial grid nodes for which the error in determining the Love numbers is less than 10−4. The number of radial grid nodes could be significantly smaller if we used a variable grid spacing. Figure 1. View largeDownload slide Relative error of the numerically computed tidal and loading degree 2 Love numbers $$k^T_e$$, $$k^T_f$$, $$k^L_e$$ and $$k^L_f$$ relative to the semi-analytically computed values by Spada et al. (2011). The subscript e denotes the instantaneous response and f is the t → ∞ ‘fluid’ limit. The shadowed region marks the domain where the relative error is smaller than 10−4. Figure 1. View largeDownload slide Relative error of the numerically computed tidal and loading degree 2 Love numbers $$k^T_e$$, $$k^T_f$$, $$k^L_e$$ and $$k^L_f$$ relative to the semi-analytically computed values by Spada et al. (2011). The subscript e denotes the instantaneous response and f is the t → ∞ ‘fluid’ limit. The shadowed region marks the domain where the relative error is smaller than 10−4. Besides the computation of the Love numbers, we have also reproduced the benchmark example in Spada et al. (2011), where model M3-L70-V01 is loaded with a spherical ice cap (Spada et al. 2011, table 4), centred at the colatitude θc = 25° and the longitude λc = 75°. Our solution consists of two steps. In the first one, the model M3-L70-V01 is rotated at a constant angular velocity ω0 = [0, 0, 7.292115] × 10−5 s−1 until it reaches hydrostatic equilibrium. Numerically we achieve this by integrating eqs (1)–(10) with the prescribed ω0 and σL ≡ 0 over a period of 200 Myr (similar procedure was used to obtain the fluid Love numbers in Fig. 1). We will denote the hydrostatic values of the polar and equatorial moments of inertia by C and A respectively ($$C\,{=}\,I_{33}^{f}$$ and $$A\,{=}\,I_{11}^f\,{=}\,I_{22}^f$$). Then, at time which we formally mark as t = 0, the body is suddenly loaded with a spherical ice cap. The loading instantaneously shifts the figure axis, and we adjust the vertical component of ω(t = 0), so that the total angular momentum is preserved in first approximation (ω = ω0(1 − c33/C), where c33 is the change of polar moment of inertia due to the loading). The further evolution of the tensor of inertia, which in turn determines the evolution of ω(t), is governed by both the viscoelastic relaxation under the load and the viscoelastic readjustment of the rotational bulge in response to the induced PM. For clarity of comparison with other studies, we introduce two quantities characterizing the PM—the PM vector m and the LOD variation ΔLOD:   $$\boldsymbol {\omega }(t)=\boldsymbol {\omega }_0 + |\boldsymbol {\omega }_0|\boldsymbol {m}(t)\, , \quad \Delta \mathrm{LOD}(t) = \frac{c_{33}(t)}{C}\frac{2\pi }{|\boldsymbol {\omega }_0|}\, ,$$ (17)where c33(t) is the time variation of I33 (polar moment of inertia) in response to the loading. In Fig. 2(a), we compare the PM components m1 and m2 computed by our method (coloured points) with those obtained using the traditional Laplace domain approach and the LLE (dashed lines, adopted from Martinec & Hagedoorn 2014). Since our solution includes the wobbling while the LLE solution by Martinec & Hagedoorn (2014) gives only the mean position of the rotation axis, the two solutions are very different in the beginning, but they match each other once the wobbling is damped. This confirms that our code gives correct results on long time scales (i.e. greater than the damping time of Chandler wobble), and conversely, that neither the wobbling nor the linearization influences the long term evolution of the PM. A good agreement is also found for ΔLOD (Fig. 2b) where the relative difference between the two solutions is smaller than 1 per cent. The study by Spada et al. (2011) provides a valuable benchmark for validation of different numerical techniques to compute the PM. Nevertheless, there are two things that are worth noting because they are not explained in the original paper in sufficient details. First, in fig. 13(c) in Spada et al. (2011), the variation ΔLOD induced by the load is plotted as a positive quantity. However, since the load is imposed close to the north pole and the total mass of the Earth is conserved, the loading reduces the polar moment of inertia and the variation ΔLOD is actually negative. Second, the hydrostatic values C and A must be corrected to reproduce the benchmark results. The ΔLOD curve (Fig. 2b) is sensitive to the value of C, which has to be 8.0394 × 1037 kg m2 to match the benchmark (i.e. the value given in table 1 of Spada et al. 2011). However, to match the benchmarked PM solution (Fig. 2a), which is sensitive to the difference C − A, this difference has to be 2.6947 ×1035 kg m2. This value corresponds to model M3-L70-V01 in the hydrostatic limit, and not to the value of C − A derived from table 1 in Spada et al. (2011). Numerically we proceeded by adding a diagonal correction to our tensor of inertia, $$\boldsymbol {I}^c(t\,{\ge }\,0)=\boldsymbol {I}(t) + \gamma \mathcal {I}$$, where $$\mathcal {I}$$ is identity tensor and Ic is the corrected tensor of inertia used in our simulations. The correction γ is given as (8.0394 − 8.17848) × 1037 kg m2, with 8.17848 × 1037 kg m2 being the hydrostatic polar moment of inertia of model M3-L70-V01. We believe this information can be helpful for the future users of the benchmark. Figure 2. View largeDownload slide (a) Polar motion induced by loading the Earth with a spherical ice cap. The coloured lines show the solution obtained by the method described in this study while the dashed black lines show the linearized solution computed in the Laplace domain and omitting the free wobble (Martinec & Hagedoorn 2014). (b) Variation of ΔLOD obtained for the same loading (blue line—our solution, dashed line—linearized solution). Figure 2. View largeDownload slide (a) Polar motion induced by loading the Earth with a spherical ice cap. The coloured lines show the solution obtained by the method described in this study while the dashed black lines show the linearized solution computed in the Laplace domain and omitting the free wobble (Martinec & Hagedoorn 2014). (b) Variation of ΔLOD obtained for the same loading (blue line—our solution, dashed line—linearized solution). 4 ENERGY BALANCE OF A ROTATING EARTH The conservation of energy for a rotating self-gravitating viscoelastic body can be expressed as   $$E_{\mathrm{rot}}+E_{\mathrm{el}}+E_{\mathrm{diss}}+E_{\mathrm{grav}}={\rm const}\, ,$$ (18)where the terms on the left-hand side are the rotational, elastic, dissipative and gravitational energies, respectively. The kinetic energy associated with the rate of deformation is not considered, since the inertia term in the equation of motion (1) is neglected. The rotational energy Erot can be computed using the standard formula   $$E_{\mathrm{rot}}=\frac{1}{2}\boldsymbol {\omega }\cdot \boldsymbol {I}\cdot \boldsymbol {\omega }.$$ (19)The elastic energy Eel stored in the body and the viscous dissipation rate $$\boldsymbol {\dot{E}}_\mathrm{diss}$$ can be expressed as (cf. Joseph 1990, p. 50)   $$E_{\mathrm{el}}=\int _{v(t)}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{4\mu }\, \mathrm{d}v \cong \int _{{\rm S}_a\!-\!{\rm S}_b}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{4\mu }\, \mathrm{d}v \, ,$$ (20)  $$\boldsymbol {\dot{E}}_{\mathrm{diss}}=\int _{v(t)}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{2\eta }\, \mathrm{d}v \cong \int _{{\rm S}_a\!-\!{\rm S}_b}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{2\eta }\, \mathrm{d}v \, ,$$ (21)where $$\boldsymbol {\tau }^d:\boldsymbol {\tau }^d \equiv \tau _{ij}^d\, \tau _{ij}^d$$ in Cartesian components and Sa − Sb denotes the spherical shell of outer radius a and inner radius b. The integration over the core is omitted since the core is assumed to be filled with inviscid liquid. Finally, the gravitational potential energy is given by   $$E_{\mathrm{grav}}=\frac{1}{2}\int _{v(t)}\rho \left(V_0 + \Phi \right)\mathrm{d}v \, .$$ (22) 4.1 Alternative formula for the gravitational energy The computation of the gravitational energy is a delicate exercise which requires a careful evaluation of topographic and volumetric contributions to integral (22). In this section we derive a formula for the gravitational energy which is easy to implement and, as we will demonstrate in Section 4.4, it is sufficiently accurate if the deformation is small. We first derive the formula for a model with a finite number of layers of constant densities, and then we extend it to a model with a continuous density profile. For the derivation we will assume that the surface load density σL is equal to zero. For a non-zero surface load density the problem can be converted into a problem with zero surface load density by using a modified radial displacement $$u_r^\mathrm{mod}(a)=u_r(a)+\sigma ^\mathrm{L}/[\rho _0]_a$$. For the sake of clarity of exposition, we consider a multilayer (N-layer) sphere. The N-layered sphere is composed of N spheres Si with radii ri, i = 1, 2, … N, where r1 = a and rN = b are respectively the radii of the outer and inner boundary of the spherical shell considered in Section 2. The density ρ0(r) of the sphere is a piecewise constant function with density jumps [ρ0]i at radii ri (taken positive if the density increases with depth). When the body is deformed, each sphere Si transforms into an aspherical object vi and the density changes by   $$\delta \rho (\boldsymbol {r},t) = \rho _0(|\boldsymbol {r}-\boldsymbol {u}(\boldsymbol {r},t)|,t) - \rho _0(r) = \left\lbrace \begin{array}{l@{\quad}l@{\quad}l}0, & |{r}{-}{r}_i|\,{\ge }\,|{u}_r(\boldsymbol {r},t)| & \forall i\in \lbrace 1,2,\ldots ,N\rbrace \\ [\rho _0]_i, & {u}_r(\boldsymbol {r},t)\,{>}\,0 \quad \wedge & |{r}\,{-}\,{r}_i|\,{<}\,|{u}_r(\boldsymbol {r},t)|\\ {-}[\rho _0]_i, & {u}_r(\boldsymbol {r},t)\,{<}\,0 \quad \wedge & |{r}\,{-}\,{r}_i|\,{<}\,|{u}_r(\boldsymbol {r},t)| \end{array} \right.$$ (23)where the three possibilities given on the right are illustrated in Fig. 3. The volume integral in eq. (22) can then be rewritten as   $$\int _{v(t)}\rho \left(V_0 + \Phi \right)\mathrm{d}v = \sum _{i=1}^{N}[\rho _0]_i\int _{v_i} \left(V_0 + \Phi \right)\mathrm{d}v.$$ (24)On the right-hand side of this equation, we obtain a sum of N integrals over N homogeneous bodies which have no internal boundaries, occupy the volumes vi and have the densities [ρ0]i. Each volume vi can be expressed in terms of the spherical volume Si and the volume hi induced by the deformation (see the yellow and green regions in Fig. 3). The volume hi is taken with a positive sign if it is above the sphere Si (yellow regions) and with a negative sign if it is inside the sphere (green regions). The total volume of hi at each interface (the sum of the yellow and green domain) is zero due to the incompressibility of the body. The gravitational energy can then be rewritten using Si and hi as follows:   $$E_{\mathrm{grav}}=\frac{1}{2}\sum _{i=1}^{N}[\rho _0]_i\int _{h_i} \left(V_0 + \Phi \right)\mathrm{d}v + \frac{1}{2}\sum _{i=1}^{N}[\rho _0]_i\int _{\mathrm{S}_i} (V_0+\Phi )\, \mathrm{d}v.$$ (25) Figure 3. View largeDownload slide Deformation of N-layer sphere. Black lines show boundaries of the reference spherical shells with constant densities, red lines show the interfaces of the deformed body. In regions where the deformation produces a positive topography, the Eulerian density increment is equal to [ρ0]i, i = 1, 2, …, N, while in regions where the topography is negative, the density increment is equal to −[ρ0]i. Figure 3. View largeDownload slide Deformation of N-layer sphere. Black lines show boundaries of the reference spherical shells with constant densities, red lines show the interfaces of the deformed body. In regions where the deformation produces a positive topography, the Eulerian density increment is equal to [ρ0]i, i = 1, 2, …, N, while in regions where the topography is negative, the density increment is equal to −[ρ0]i. The gravitational potential V0 + Φ can be split into N parts, $$V_0+\Phi =\sum _{i=1}^{N}(V_0^i+\Phi ^i)$$, where $$V_0^i$$ is the potential of the homogeneous sphere Si of density [ρ0]i and Φi is the perturbation due to the aspherical shape of the i-th interface:   $$V_0^i(r)=-G[\rho _0]_i\int _{\mathrm{S}_i} \frac{1}{|\boldsymbol {r}-\boldsymbol {r}^{\prime }|}\, \mathrm{d}v,\quad \Phi ^i(\boldsymbol {r})=-G[\rho _0]_i\int _{h_i} \frac{1}{|\boldsymbol {r}-\boldsymbol {r}^{\prime }|}\, \mathrm{d}v.$$ (26)In the following, we take advantage of the well-known fact that the gravitational energy of a body A in the gravitational field of a body B is the same as the gravitational energy of body B in the gravitational field of body A. By applying this general rule to eq. (26) we obtain the following relations:   $$[\rho _0]_i\int _{\mathrm{S}_i}\Phi ^j\, \mathrm{d}v=[\rho _0]_j\int _{h_j} V_0^i\, \mathrm{d}v,$$ (27)valid for any i, j = 1, …, N. Using eq. (27) we can link the first and the second sum on the right-hand side of eq. (25):   $$\sum _{i=1}^{N}[\rho _0]_i\int _{\mathrm{S}_i} \Phi \, \mathrm{d}v = \sum _{i=1}^{N}[\rho _0]_i \int _{\mathrm{S}_i}\sum _{j=1}^{N}\Phi ^j\, \mathrm{d}v = \sum _{j=1}^{N}[\rho _0]_j\int _{h_j} \sum _{i=1}^{N} V_0^i\, \mathrm{d}v = \sum _{i=1}^{N}[\rho _0]_i\int _{h_i} V_0\, \mathrm{d}v.$$ (28)Substituting (28) into (25) gives the alternative formula for the gravitational energy Egrav of the body:   $$E_{\mathrm{grav}}=E_{\mathrm{grav}}^0+\sum _{i=1}^N E_{\mathrm{grav}}^i,$$ (29)where Egrav0 is the gravitational energy of the undeformed N-layer sphere   $$E_{\mathrm{grav}}^0=\frac{1}{2}\sum _{i=1}^N[\rho _0]_i\int _{\mathrm{S}_i} V_0\, \mathrm{d}v=2\pi \int _0^a\rho _0(r) V_0(r) r^{\, 2}\, \mathrm{d}r\, ,$$ (30)and   $$E_{\mathrm{grav}}^i=[\rho _0]_i \int _{h_i}\left(V_0 + \frac{1}{2}\Phi \right)\mathrm{d}v\, .$$ (31)Since the deformation is small, the topography of the ith interface (defined here as the deviation from a sphere) can be approximated by the radial component of the displacement vector at the interface. Eq. (31) then takes the form   $$E_{\mathrm{grav}}^i\cong [\rho _0]_i \int _{\partial {\rm S}_i}\int _{r_i}^{r_i+u_r(r_i,\vartheta ,\varphi )}\left(V_0 + \frac{1}{2}\Phi \right)\mathrm{d}r\, \mathrm{d}s\, ,$$ (32)where ∂Si denotes the surface of the sphere of radius ri, ds is the element of ∂Si, and ϑ and φ are the spherical angular coordinates. The integral in eq. (32) can be simplified by expanding V0 in a Taylor series and setting Φ(r, ϑ, φ) ≅ Φ(ri, ϑ, φ):   $$E_{\mathrm{grav}}^i\cong [\rho _0]_i \int _{\partial {\rm S}_i}\int _{r_i}^{r_i+u_r(r_i,\vartheta ,\varphi )}\left(V_0(r_i)+g_0(r_i)(r-r_i) + \frac{1}{2}\Phi (r_i,\vartheta ,\varphi )\right)\mathrm{d}r\, \mathrm{d}s\, .$$ (33)Integration of eq. (33) over the radius gives   $$E_{\mathrm{grav}}^i\cong \frac{1}{2}[\rho _0]_i \int _{\partial {\rm S}_i}\left(g_0 u_r^2 + \Phi u_r\right)\, \mathrm{d}s\, ,$$ (34)where we used that $$\int _{\partial {\rm S}_i}u_r\mathrm{d}s=0$$ for an incompressible body. The surface integral in eq. (34) can be easily evaluated in spectral domain. Representing quantities ur and Φ in terms of spherical harmonic expansions, eq. (15), and invoking the orthonormality of the spherical harmonic basis {Yℓm}, we obtain   $$E_{\mathrm{grav}}\cong E_{\mathrm{grav}}^0+\frac{1}{2}\sum _{i=1}^N[\rho _0]_i\, r_i^2 \sum _{\ell =0}^{\ell _{\rm max}}\sum _{m=-\ell }^{\ell } \Big (g_0(r_i) |u_{r,\ell m}(r_i)|^2 + u_{r,\ell m}(r_i)\Phi ^{\ast }_{\ell m}(r_i)\Big )\, ,$$ (35)where * denotes complex conjugation. The rotational potential is described by spherical harmonics of degrees 0 and 2. Since degree 0 has no effect on the deformation of an incompressible body, we will only consider the terms with ℓ = 2 and |m| ≤ 2. The formula (35), derived for a layered density model, can be generalized to the case of a continuous density profile. Replacing the density jumps [ρ0]i in eq. (31) by (−dρ0(r)/dr) dr and the sum in eq. (29) by an integral, and following the analogous procedure as above, we get the formula   \begin{eqnarray} E_{\mathrm{grav}}&\cong & E_{\mathrm{grav}}^0 + \frac{1}{2}\bigg ( [\rho _0]_a\int _{\partial {\rm S}_a}\!\!\left(g_0u_r^2+u_r\Phi \right)\mathrm{d}s +[\rho _0]_b\int _{\partial {\rm S}_b}\!\!\left(g_0u_r^2+u_r\Phi \right)\mathrm{d}s - \int _{b}^a \frac{\mathrm{d}\rho _0}{\mathrm{d}r}\int _{\partial {\rm S}_r}\!\! \left(g_0u_r^2+u_r\Phi \right)\mathrm{d}s\, \mathrm{d}r\bigg )\, , \end{eqnarray} (36)where ∂Sr denotes the surface of the sphere with radius r. Eq. (36) allows us to compute the gravitational energy using the quantities obtained by solving the governing equations in a regular computational domain Sa − Sb, that is, in the spherical shell of outer radius a and inner radius b. The radial displacement ur needed to evaluate the integrals in eq. (36) is obtained as the solution to eqs (1)–(3), (9) and (10), with the gravitational potential Φ being computed iteratively from ρ0 and ur via eq. (16), that is, using the traditional ‘condensation’ method (e.g. Choblet et al. 2007). If the deformation is small, the Eulerian and Lagrangian approaches give the same displacement (see Fig. 7b in Section 4.4). Therefore, the formula for the gravitational energy derived here can be used regardless which of the approaches is chosen to evaluate the deformation. Note, however, that a direct application of the condensation method to the energy integral (22) leads to an incorrect result. In order to demonstrate this, we will assume, for simplicity, that the density ρ is constant throughout the body. Eq. (22) can then be expressed as the sum of two integrals over the sphere Sa and two integrals over the domain ha representing the irregular shape of the body (ha ≡ h1 in the notation used above):   $$E_{\rm grav}=\frac{\rho }{2}\left( \int _{{\rm S}_a} V_0\, {\rm d}v+ \int _{{\rm S}_a} \Phi \, {\rm d}v+ \int _{h_a} V_0\, {\rm d}v+ \int _{h_a} \Phi \, {\rm d}v \right)\, .$$ (37)The first integral corresponds to $$E_{\rm grav}^0$$ (cf. eq. 30). The second integral can be evaluated by expanding Φ into spherical harmonics. Due to the orthogonality of the spherical harmonic functions, only Φ00 has non-zero contribution to the integral. In the simplified formula (16), which is based on the condensation technique, there is no coupling among coefficients of different degrees and orders. According to that formula, Φ00 is a function of ur, 00 only, and, since ur, 00 is zero for an incompressible body, the second integral vanishes. Condensing the volumetric density ρ into surface mass density ρur(a) in the last two integrals then gives   $$E_{\rm grav}=E_{\rm grav}^0+\frac{\rho }{2}\left(V_0(a)\int _{\partial {\rm S}_a}u_r\, {\rm d}s+\int _{\partial {\rm S}_a}\Phi u_r\, {\rm d}s\right) = E_{\rm grav}^0+\frac{\rho }{2}\int _{\partial \rm {S}_a}\Phi u_r\, {\rm d}s,$$ (38)while the correct solution is (eq. (36) with [ρ0]a = ρ, [ρ0]b = 0 and dρ0/dr = 0)   $$E_{\rm grav}= E_{\rm grav}^0+ \frac{\rho }{2}\int _{\partial {\rm S}_a}\left(g_0u_r^2+\Phi u_r\right)\, {\rm d}s\, .$$ (39)In other words, the second and third integrals in eq. (37) have been incorrectly evaluated as zero, giving the wrong solution (38). In this section we have shown that both integrals are in fact non-zero, they are equally important, and together yield the contribution with $$g_0 u_r^2$$ in the integrand. The advantage of our formula for the gravitational energy is that it does not require the computation of Φ00, which is a non-zero quantity inside the sphere Sa if evaluated exactly. In appendix A, we derive other formulae that can be helpful in assessing the energy balance of a deforming body, based on evaluating the power of individual forces in the equation of motion (1). 4.2 Energy balance of GIA for a non-rotating Earth We first evaluate the energy balance for a non-rotating Earth with internal structure given by model M3-L70-V01, loaded at time t = 0 with a spherical ice cap as described in Spada et al. (2011). In this case, the body has a perfectly spherical shape prior to the loading (ρ = ρ0, Φ = Ψ = 0, u = 0, τd = 0). No elastic energy is stored and the gravitational energy reaches the absolute minimum. Since the total mass of the imposed surface load is zero, the loading is associated with a redistribution of mass at the surface of the body. This, although partially compensated by the instantaneous elastic deformation of the model due to the loading, results in shape changes and thus in an increase of the gravitational energy (Fig. 4, green curve). At the same time the elastic energy (plotted in magenta) rises to a value that is about ten times smaller than the increase in the gravitational energy. With increasing time, the system returns to equilibrium—the gravitational energy decreases and is dissipated in heat (blue curve in Fig. 4). Figure 4. View largeDownload slide Time evolution of the gravitational (green), dissipative (blue) and elastic (magenta) energy for non-rotating Earth model M3-L70-V01, loaded with a spherical ice cap (for details, see Spada et al. 2011). The variation of the gravitational energy is plotted with respect to the gravitational energy $$E^0_{\rm grav}$$ of the undeformed body. The sum of the three energies is shown by a black line. Figure 4. View largeDownload slide Time evolution of the gravitational (green), dissipative (blue) and elastic (magenta) energy for non-rotating Earth model M3-L70-V01, loaded with a spherical ice cap (for details, see Spada et al. 2011). The variation of the gravitational energy is plotted with respect to the gravitational energy $$E^0_{\rm grav}$$ of the undeformed body. The sum of the three energies is shown by a black line. Since the uppermost layer in model M3-L70-V01 is purely elastic, the final equilibrium state is non-hydrostatic and the values of the elastic and gravitational energies at time t → ∞ are non-zero. The elastic energy does not decrease monotonically, showing a local maximum at t ≅ 50 kyr. This maximum is a consequence of a complex interplay between the elastic lithosphere, in which the stored elastic energy monotonically increases in time, and the deeper mantle layers, where the elastic energy behaves non-monotonically, until it finally returns to zero in the hydrostatic equilibrium. During the whole evolution, the sum of the three considered energies is constant (black line in Fig. 4). 4.3 Energy balance of a damped Chandler wobble We consider the Earth model M3-L70-V01 rotating with constant angular velocity ω0 = [0, 0, 7.292115] × 10 − 5 s−1. The Earth is initially in equilibrium and has an oblate spheroidal shape. At time t = 0, the direction of the rotational vector is suddenly changed so that the angle between the figure axis and the rotational vector is 2.69 × 10−2 degrees, which is approximately the same as the angle by which figure axis is shifted in the benchmark example in Spada et al. (2011)—see Fig. 2(a) for comparison. This change induces a periodic motion (‘wobble’) of the axis of rotation, which is gradually damped due to dissipation. After a few thousands of years, the axis of rotation returns to its initial position, but the final rotational speed is smaller than the initial one since a part of the rotational energy is transformed to heat (note that the angular momentum is not conserved when shifting the axis at time t = 0). The time evolution of the individual types of energies associated with the above process is shown in Fig. 5. All energies are plotted with respect to the initial equilibrium state. In this state (corresponding to t < 0), the elastic energy is non-zero only in the uppermost, purely elastic layer. The abrupt change of the rotational vector at t = 0 induces an instantaneous elastic deformation of the whole body (except the liquid core) and the total elastic energy increases. Since the deformation acts to reduce the flattening, both the gravitational and rotational energies decrease at t = 0. The rotational energy (red curve in Fig. 5) further decreases with increasing time as it is dissipated through viscous deformation (blue curve). The variations in the gravitational and elastic energies (green and magenta curves, respectively) significantly contribute to the global energy budget only in the beginning of evolution. For t → ∞, both curves converge to negative values which are close to zero, suggesting that the gravitational and elastic energies in the final equilibrium state are slightly smaller than those in the initial equilibrium state. This difference is associated with the change of flattening due to the decrease in the rotational speed. Figure 5. View largeDownload slide Time evolution of different types of energy for wobbling model M3-L70-V01. The variations of energy are plotted with respect to the initial equilibrium state in which the rotational vector is parallel to the main axis of the body. The wobble is induced by changing the direction of the rotational vector by 2.69 × 10−2 degrees at time t = 0. Figure 5. View largeDownload slide Time evolution of different types of energy for wobbling model M3-L70-V01. The variations of energy are plotted with respect to the initial equilibrium state in which the rotational vector is parallel to the main axis of the body. The wobble is induced by changing the direction of the rotational vector by 2.69 × 10−2 degrees at time t = 0. 4.4 Energy balance of GIA for a rotating Earth The setting of this numerical experiment corresponds to the benchmark case adopted from Spada et al. (2011), which was already discussed in Section 3 (cf. Figs 1 and 2). The rotating Earth model M3-L70-V01, initially in equilibrium state, is suddenly loaded with an ice cap. The loading and the associated deformation induce a PM (Fig. 2a) and changes in the rotational speed (Fig. 2b), which are accompanied by large variations in the rotational and gravitational energy (Fig. 6, red and green curves, respectively). The amplitudes of these variations are of the order of 1022–1023 J and they are three orders of magnitude larger than the variations in the elastic (magenta) and dissipative (blue) energy (cf. also Figs 4 and 5). During the whole evolution, the sum of all energies is constant (full black line). Figure 6. View largeDownload slide Time evolution of different types of energy for rotating model M3-L70-V01 loaded with an ice cap at t = 0 (for details, see Spada et al. 2011). The variations of energy are plotted with respect to the initial equilibrium state in which the rotational vector is parallel to the main axis of the body. The results obtained using the NLE are represented by coloured lines while the LLE solution is plotted in black. Figure 6. View largeDownload slide Time evolution of different types of energy for rotating model M3-L70-V01 loaded with an ice cap at t = 0 (for details, see Spada et al. 2011). The variations of energy are plotted with respect to the initial equilibrium state in which the rotational vector is parallel to the main axis of the body. The results obtained using the NLE are represented by coloured lines while the LLE solution is plotted in black. So far we have discussed the results computed using the NLE. In Fig. 6 the variations of the gravitational and rotational energy obtained using the LLE are shown by dashed and dotted black lines, respectively. While the curve for the gravitational energy coincides almost exactly with that obtained using the NLE, the accuracy in determining the rotational energy decreases with increasing time, resulting in a relative error of almost 100 per cent after 5 kyr. Consequently, the total energy of the system (dashed black line) is not conserved for the LLE solution. As demonstrated in Fig. 7(a), the large error in evaluating the rotational energy is associated with determining the component m3, which is significantly larger than the value obtained using the NLE, in contrast to components m1 and m2, which are determined with high accuracy (see Fig. 2a). Fig. 7(b) shows that the predicted degree 2 shape of the Earth does not depend on whether we use the NLE or LLE, suggesting that the displacement obtained using the LLE is only slightly affected by the error in m3. This explains the good agreement obtained for the gravitational energy which is a function of the radial component of displacement (see eq. 36) and, unlike the rotational energy, does not depend directly on m3. Figure 7. View largeDownload slide (a) Time evolution of m3 for the studied benchmark case (cf. Fig. 2). Red points depict our NLE solution, dashed line is the LLE solution. Green line shows our solution when the spherical ice cap is imposed on the surface of the Earth only gradually (over the period of 5 yr). (b) Displacement components ur, 20 and ur, 21 relative to the hydrostatic equilibrium for the same case as above (continuous loading is omitted for clarity of the figure). Figure 7. View largeDownload slide (a) Time evolution of m3 for the studied benchmark case (cf. Fig. 2). Red points depict our NLE solution, dashed line is the LLE solution. Green line shows our solution when the spherical ice cap is imposed on the surface of the Earth only gradually (over the period of 5 yr). (b) Displacement components ur, 20 and ur, 21 relative to the hydrostatic equilibrium for the same case as above (continuous loading is omitted for clarity of the figure). Note that the difference between the NLE and LLE solution for m3 is not related to the free wobble, as illustrated by the green line in Fig. 7(a). The line represents a case in which the ice cap is being imposed only gradually onto the surface of the Earth, over the period of 5 yr, significantly reducing the amplitude of the induced wobble but leading to the same long term behaviour as for the instantaneous loading. In the light of the results presented in Section 3, the large error in component m3 obtained using the LLE may be viewed as surprising. As demonstrated in Fig. 2(b), both the NLE and LLE predict the same LOD variation ΔLOD given by eq. (17) (right), but m3 computed as m3 = −c33/C is still significantly affected by the neglect of the nonlinear terms in the LE (Fig. 7a). To explain this apparent paradox, we start from the formula (9.2.3) in Munk & MacDonald (1960)   $$\frac{\Delta \mathrm{LOD}}{\mathrm{LOD}}=-m_3=\frac{c_{33}}{C}\, ,$$ (40)and we describe the simplifications that were used in its derivation. The first equality in eq. (40) can be derived by linearization of the exact formula for ΔLOD:   $$\frac{\Delta \mathrm{LOD}}{\mathrm{LOD}}=\frac{2\pi \omega ^{-1}-2\pi \omega _0^{-1}}{2\pi \omega _0^{-1}}=\frac{\omega _0-\omega }{\omega }\cong \frac{\omega _0-\omega }{\omega _0} =1-\sqrt{m_1^2+m_2^2+(1+m_3)^2}\, .$$ (41)Using the binomial theorem and neglecting the higher-order terms, we obtain   $$\frac{\Delta \mathrm{LOD}}{\mathrm{LOD}}\cong -m_3+\frac{1}{2}\left(m_1^2+m_2^2+m_3^2\right)\, ,$$ (42)which, after neglecting the quadratic terms, gives the first equality in eq. (40). The second equality in eq. (40) comes from the LLE for component m3. In appendix B, we derive the following higher-order approximation:   $$-\!m_3 + \frac{1}{2}(m_1^2+m_2^2+m_3^2 ) \cong \frac{{c}_{33}}{C}\, .$$ (43)Combining eq. (43) with eq. (42) yields a more accurate version of the original equation (40):   $$\frac{\Delta \mathrm{LOD}}{\mathrm{LOD}}\cong -m_3 + \frac{1}{2}(m_1^2+m_2^2+m_3^2 ) \cong \frac{{c}_{33}}{C}\, .$$ (44)The definition of ΔLOD in eq. (17) thus includes second-order accurate evaluation of the LOD variation, while the equation m3 = −c33/C is only first-order accurate, which explains the disagreement between the LLE and NLE solutions in Fig. 7(a). The fact that the m3 component of the rotation vector is affected by a large error does not invalidate the results of previous studies using the LLE solution. As demonstrated above, the traditional equation for ▵LOD, eq. (17), used in these studies, is more accurate than the linearized equation for m3 and it gives the same LOD variation as the NLE equation (Fig. 2b), provided that excursions of the rotation axis are small. Our numerical tests (not presented in this paper) suggest that the shift of the rotation axis must be larger than 1 degree for the nonlinearity to significantly (by at least a few per cent) affect the resulting PM and LOD. The ice cap considered here is rather small and, since the characteristic time of equatorial bulge readjustment is inversely proportional to the size of the load (Ricard et al. 1993), several Myr would be needed to shift the pole by 1 degree. If we considered an unrealistic, ten times larger ice cap, the equatorial bulge readjustment would occur faster and the linearized PM and LOD solutions would differ from the nonlinear one by few per cent already after about 200 kyr. 4.5 Energy balance as a tool for testing GIA models The rotational feedback in the equation of motion has not been considered in a number of papers dealing with modelling of GIA on the Earth (for details, see Introduction in Martinec & Hagedoorn 2014). While such approximation is justifiable for the displacement component of degree 2 order 0 (ur, 20 has nearly identical evolution with and without the feedback), the displacement component of degree 2 and order 1 behaves in time completely differently for the two cases. A long term evolution of ur, 21 is governed by the balance between the viscoelastic relaxation of the model under the load and the viscoelastic readjustment of the equatorial bulge. For a homogeneous model, these two mechanisms perfectly cancel out each other for all times, as pointed out by Munk & MacDonald (1960), Section 6.4, and Wu & Peltier (1984). The effect of the rotational feedback has been recently discussed by Martinec & Hagedoorn (2014), who demonstrated that it may have a significant impact on the prediction of ur, 21 and sea level change. It should be emphasized that including the variations of the centrifugal potential into the field equations governing GIA-induced deformation of the Earth differs from including them into the gravity appearing in the sea level equation (i.e. the rotational feedback on the sea level equation), which has already been addressed in the previous studies (e.g. Milne & Mitrovica 1998; Peltier 1998) and which is not a focus of our work. Fig. 8 shows the time evolution of different types of energy when the rotation of the Earth is held constant, that is, the LE is not being solved at all. Fig. 9 compares the rotational and gravitational energies for the cases with and without the rotational feedback in the equation of motion, respectively. We can see that the gravitational energies do not differ significantly. This can be explained by the quadratic dependence of the gravitational energy on the radial component of the displacement (cf. eq. 35). While the hydrostatic value of ur, 20 reaches thousands of metres, the hydrostatic value of ur, 21 is equal to zero. Since the deformation induced by the ice cap reaches only few metres (Fig. 7b), the change of the gravitational energy of the rotating body is dominated by the change of ur, 20. In addition, as discussed above, ur, 20(t) is not affected by the rotational feedback significantly. Figure 8. View largeDownload slide Time evolution of different types of energy for a rotating model M3-L70-V01 loaded with a spherical ice cap at t = 0 (for details, see Spada et al. 2011). The rotation of the model is kept constant throughout the simulation. Figure 8. View largeDownload slide Time evolution of different types of energy for a rotating model M3-L70-V01 loaded with a spherical ice cap at t = 0 (for details, see Spada et al. 2011). The rotation of the model is kept constant throughout the simulation. Figure 9. View largeDownload slide Comparison of time evolutions of the rotational and gravitational energies from Figs 6 and 8. The energies from Fig. 6 are in colour, dashed black lines show the curves from Fig. 8. Figure 9. View largeDownload slide Comparison of time evolutions of the rotational and gravitational energies from Figs 6 and 8. The energies from Fig. 6 are in colour, dashed black lines show the curves from Fig. 8. The rotational energies, on the other hand, differ significantly for the two cases. With the rotational feedback in the equation of motion, the rotational energy evolves in time with the opposite sign than the gravitational energy in order to conserve the total energy. Without considering the rotational feedback (i.e. ω(t) = ω0), the rotational energy evolves in time similarly to the gravitational energy. This is because the increment of (ω0 · I · ω0)/2 has, similarly to the change of the gravitational energy, a positive linear dependence on the change of ur, 20 (see MacCullagh’s formula for c33 and insert for Φ20 from eq. 16). The detailed analysis of individual terms in the energy balance is given in Appendix A. Fig. 10 shows the results of such analysis for the two simulations shown in Figs 6 and 8, for which the total energy of the Earth’s model is not conserved. We can see that the energy imbalance is associated with the rotational energy, that is, the power of the centrifugal force is not equal to the rate of change of rotational energy. Figure 10. View largeDownload slide The energy integrals given in Section 4 (black lines) and their equivalents obtained by integrating the rate of mechanical work over the computational domain and time (coloured lines). Two cases are considered: (a) the rotational feedback is computed using the linearized Liouville equation (cf. Fig. 6), and (b) the rotational feedback is completely omitted (cf. Fig. 8). The numbers have the following meaning: 1—dissipated energy (dashed black and blue lines); 2—rotational energy (dotted black and red lines); 3 and 4—contributions to gravitational energy (dashed-dotted and green and yellow lines). For more details, see Appendix A and Fig. A1. Figure 10. View largeDownload slide The energy integrals given in Section 4 (black lines) and their equivalents obtained by integrating the rate of mechanical work over the computational domain and time (coloured lines). Two cases are considered: (a) the rotational feedback is computed using the linearized Liouville equation (cf. Fig. 6), and (b) the rotational feedback is completely omitted (cf. Fig. 8). The numbers have the following meaning: 1—dissipated energy (dashed black and blue lines); 2—rotational energy (dotted black and red lines); 3 and 4—contributions to gravitational energy (dashed-dotted and green and yellow lines). For more details, see Appendix A and Fig. A1. Another example of the application of the energy analysis given in Appendix A is related to the NLE solution presented in Sections 3 and 4.4. We run a numerical experiment where the original time step of 1 × 10−3 yr is made coarser as 5 × 10−2 yr. The results are shown in Fig. 11(a). Inspecting Fig. 11(a) we can see that not all the powers match the corresponding rates of energies. Again, it is the rotational energy that is not balanced. Hence, time step of 5 × 10−2 yr is too coarse in this case, which causes that the LE, but not eqs (1)–(3), is violated. In the parallel numerical experiment, the original spatial discretization of 460 radial nodes is reduced to 40 nodes only. Eqs (1)–(3) are now violated which is seen from the fact that the forces depending on the spatial derivatives of field variables do not match their corresponding energies. Figure 11. View largeDownload slide Panel (a) shows the case when the governing equations are solved using a too large time step, while panel (b) illustrates the effect of coarse radial discretization. The types of curves and the meaning of the numbers are the same as in Fig. 10. Figure 11. View largeDownload slide Panel (a) shows the case when the governing equations are solved using a too large time step, while panel (b) illustrates the effect of coarse radial discretization. The types of curves and the meaning of the numbers are the same as in Fig. 10. We hope that such a posteriori diagnostic of numerical solutions may be of interest also in other applications. For instance, in the case where the LE is employed for determining the Q-factor of the Chandler wobble, as suggested by Nakada & Karato (2012), the solution is numerically very sensitive to time-stepping. An independent test of the solution could thus prove helpful. 5 CONCLUSIONS In this paper, the temporal evolution of rotational, elastic, dissipative and gravitational energies is studied in the context of GIA. We have derived spectral formulae for computing the energies, overcoming the difficulty associated with the use of a simplified expression for the perturbation of gravitational potential. The energy evolution is then computed for several examples, demonstrating the importance of rotational and gravitational energies in the resulting balance. Viscous dissipation is proven relevant only for the damping of the Chandler wobble. By summing all energies for a rotating Earth loaded with a spherical ice cap, we show that the LLE is energetically inconsistent—the error resulting from an inaccurate evaluation of the m3 component of rotation vector. Careful evaluation of the energy balance for a rotating Earth loaded with an ice cap shows that the omission of the rotational feedback in the momentum equation leads to a significant violation of the energy conservation law, with both rotational and gravitational energy increasing in time on a time scale typical of the Earth’s glacial cycle. This finding underlines the importance of proper implementation of rotation variations in GIA modelling. Throughout the study we strictly use the Eulerian formulation for computing the deformation of the planet, and solve the governing equations directly in time domain. The equivalence with the Lagrangian formulation and Laplace domain approach (normal mode theory) is demonstrated numerically by comparing the Love numbers. The evaluation of the mechanical work of individual forces considered in the Eulerian formulation provides a versatile numerical tool, which can be used to verify solutions to the problems dealing with small deformations of planetary bodies. Acknowledgements This work has been supported by the Charles University grant SVV 260447/2017. We thank Volker Klemann and an anonymous reviewer for helpful comments that improved the manuscript. REFERENCES Benjamin D., Wahr J., Ray R.D., Egbert G.D., Desai S.D., 2006. Constraints on mantle anelasticity from geodetic observations, and implications for the J(2) anomaly, Geophys. J. Int. , 165( 1), 3– 16. Google Scholar CrossRef Search ADS   Cambiotti G., Ricard Y., Sabadini R., 2010. Ice age True Polar Wander in a compressible and non-hydrostatic Earth, Geophys. J. Int. , 183( 3), 1248– 1264. Google Scholar CrossRef Search ADS   Cambiotti G., Ricard Y., Sabadini R., 2011. New insights into mantle convection true polar wander and rotational bulge readjustment, Earth planet. Sic. Lett. , 310(3–4), 538– 543. Google Scholar CrossRef Search ADS   Cathles L., 1975. Viscosity of the Earth’s Mantle , Princeton Legacy Library, Princeton Univ. Press. Choblet G., Cadek O., Couturier F., Dumoulin C., 2007. OE DIPUS: a new tool to study the dynamics of planetary interiors, Geophys. J. Int. , 170( 1), 9– 30. Google Scholar CrossRef Search ADS   Floberghagen R., Fehringer M., Lamarre D., Muzi D., Frommknecht B., Steiger C., Pineiro J., da Costa A., 2011. Mission design, operation and exploitation of the gravity field and steady-state ocean circulation explorer mission, J. Geod. , 85( 11), 749– 758. Google Scholar CrossRef Search ADS   Gerya T., Yuen D., 2003. Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties, Phys. Earth planet. Inter. , 140( 4), 293– 318. Google Scholar CrossRef Search ADS   Golle O., Dumoulin C., Choblet G., Cadek O., 2012. Topography and geoid induced by a convecting mantle beneath an elastic lithosphere, Geophys. J. Int. , 189( 1), 55– 72. Google Scholar CrossRef Search ADS   Gurtin M.E., 1981. An Introduction to Continuum Mechanics , vol. 158 of Mathematics in Science and Engineering, Academic Press. Harris A., 1994. Tumbling asteroids, Icarus , 107( 1), 209– 211. Google Scholar CrossRef Search ADS   Jones M.N., 1985. Spherical Harmonics and Tensors for Classical Field Theory , vol. 2 of Applied and Engineering Mathematics Series, Research Studies Press. Joseph D.D., 1990. Fluid Dynamics of Viscoelastic Liquids , vol. 84 of Applied Mathematical Sciences, Springer-Verlag. Google Scholar CrossRef Search ADS   Kite E.S., Matsuyama I., Manga M., Perron J.T., Mitrovica J.X., 2009. True Polar Wander driven by late-stage volcanism and the distribution of paleopolar deposits on Mars, Phys. Earth planet. Inter. , 280(1–4), 254– 267. Google Scholar CrossRef Search ADS   Lau H.C.P., Mitrovica J.X., Austermann J., Crawford O., Al-Attar D., Latychev K., 2016. Inferences of mantle viscosity based on ice age data sets: radial structure, J. geophys. Res. , 121( 10), 6991– 7012. Google Scholar CrossRef Search ADS   Lefftz M., Legros H., Hinderer J., 1991. Non-linear equations for the rotation of a viscoelastic planet taking into account the influence of a liquid core, Celest. Mech. Dyn. Astron. , 52( 1), 13– 43. Google Scholar CrossRef Search ADS   Martinec Z., Hagedoorn J., 2014. The rotational feedback on linear-momentum balance in glacial isostatic adjustment, Geophys. J. Int. , 199( 3), 1823– 1846. Google Scholar CrossRef Search ADS   Milne G., Mitrovica J., 1998. Postglacial sea-level change on a rotating Earth, Geophys. J. Int. , 133( 1), 1– 19. Google Scholar CrossRef Search ADS   Mitrovica J., Peltier W., 1991. A complete formalism for the inversion of postglacial rebound data - resolving power analysis, Geophys. J. Int. , 104( 2), 267– 288. Google Scholar CrossRef Search ADS   Mitrovica J., Milne G., Davis J., 2001. Glacial isostatic adjustment on a rotating earth, Geophys. J. Int. , 147( 3), 562– 578. Google Scholar CrossRef Search ADS   Munk W.H., MacDonald G.J.F., 1960. Rotation of the Earth , Cambridge Univ. Press. Nakada M., Karato S.-I., 2012. Low viscosity of the bottom of the Earth’s mantle inferred from the analysis of Chandler wobble and tidal deformation, Phys. Earth planet. Inter. , 192, 68– 80. Google Scholar CrossRef Search ADS   Patočka V., 2013. Polar Wander prediction based on the solution of the Liouville equation, Master’s thesis , Charles University in Prague, Czech Republic, http://geo.mff.cuni.cz/theses/2013-Patocka-Mgr.pdf. Peltier W.R., 1974. The impulse response of a Maxwell Earth, Rev. Geophys. Space Phys. , 12, 649– 669. Google Scholar CrossRef Search ADS   Peltier W.R., 1998. Postglacial variations in the level of the sea: implications for climate dynamics and solid-earth geophysics, Rev. Geophys. , 36( 4), 603– 689. Google Scholar CrossRef Search ADS   Ricard Y., Spada G., Sabadini R., 1993. Polar wandering of a dynamic earth, Geophys. J. Int. , 113( 2), 284– 298. Google Scholar CrossRef Search ADS   Runcorn S., 1984. The primeval axis of rotation of the moon, Phil. Trans. R. Soc. , 313( 1524), 77– 83. Google Scholar CrossRef Search ADS   Sabadini R., Vermeersen B., Cambiotti G., 2016. Global Dynamics of the Earth: Applications of Viscoelastic Relaxation Theory to Solid-Earth and Planetary Geophysics , Modern Approaches in Geophysics, Springer. Google Scholar CrossRef Search ADS   Sabadini R., Yuen D., Boschi E., 1982. Polar wandering and the forced responses of a rotating, multilayered, viscoelastic planet, J. geophys. Res. , 87( NB4), 2885– 2903. Google Scholar CrossRef Search ADS   Soucek O., Hron J., Behounkova M., Cadek O., 2016. Effect of the tiger stripes on the deformation of Saturn’s moon Enceladus, Geophys. Res. Lett. , 43( 14), 7417– 7423. Google Scholar CrossRef Search ADS   Spada G. et al.  , 2011. A benchmark study for glacial isostatic adjustment codes, Geophys. J. Int. , 185( 1), 106– 132. Google Scholar CrossRef Search ADS   Tanaka Y., Klemann V., Okuno J., 2009. Application of a numerical inverse Laplace integration method to surface loading on a viscoelastic compressible earth model, Pure appl. Geophys. , 166(8–9), 1199– 1216. Google Scholar CrossRef Search ADS   Tapley B.D., Bettadpur S., Watkins M., Reigber C., 2004. The gravity recovery and climate experiment: mission overview and early results, Geophys. Res. Let. , 31, L09607, doi:10.1029/2004GL019920. Google Scholar CrossRef Search ADS   Tobie G., Cadek O., Sotin C., 2008. Solid tidal friction above a liquid water reservoir as the origin of the south pole hotspot on Enceladus, Icarus , 196( 2), 642– 652. Google Scholar CrossRef Search ADS   Vermeersen L., Sabadini R., Devoti R., Luceri V., Rutigliano P., Sciarretta C., Bianco G., 1998. Mantle viscosity inferences from joint inversions of Pleistocene deglaciation-induced changes in geopotential with a new SLR analysis and polar wander, Geophys. Res. Lett. , 25( 23), 4261– 4264. Google Scholar CrossRef Search ADS   Wu P., Peltier W., 1982. Viscous gravitational relaxation, Geophys. J. R. astr. Soc. , 70( 2), 435– 485. Google Scholar CrossRef Search ADS   Wu P., Peltier W., 1984. Pleistocene deglaciation and the Earth’s rotation—a new analysis, Geophys. J. R. astr. Soc. , 76( 3), 753– 791. Google Scholar CrossRef Search ADS   APPENDIX A: ENERGY BALANCE AND THE RATE OF MECHANICAL WORK The energies given in Section 4 can be matched with the rates of mechanical work done by individual terms in the momentum eq. (1),   $$\nabla \cdot \boldsymbol {\tau } - (\boldsymbol {u}\cdot \nabla \rho _0)\, \boldsymbol {g}_0 - \rho _0\nabla {\left(\Phi +\Psi \right)} + \rho _0\boldsymbol {g}_0 = 0 \, ,$$ (A1)where we substituted for the body force from eq. (8). Multiplication of eq. (A1) by the velocity vector, application of the tensor identity (∇ · τ) · v = ∇ · (τ · v) − τ: ∇v and subsequent integration over the regular computational domain Sa − Sb gives the balance equation of the rate of mechanical work:   $$\int _{{\rm S}_a\!-\!{\rm S}_b}\!\!\Big(\!-\boldsymbol {\tau }:\nabla \boldsymbol {v}+\nabla \cdot (\boldsymbol {\tau }\cdot \boldsymbol {v}) -(\boldsymbol {u}\cdot \nabla \rho _0)\boldsymbol {g}_0\cdot \boldsymbol {v} - (\rho _0\nabla \Phi )\cdot \boldsymbol {v}-(\rho _0\nabla \Psi )\cdot \boldsymbol {v} + \rho _0\boldsymbol {g}_0\cdot \boldsymbol {v}\Big )\, {\rm d}v = 0\, .$$ (A2)We now evaluate the integrals of the individual terms in eq. (A2) to better understand the relationship between the mechanical work and the different types of energy discussed in the main text. The resulting expressions can be used as an alternative to the energy integrals given in Section 4 and demonstrate the equivalence of the two approaches through a numerical test (Fig. A1). Figure A1. View largeDownload slide (a) Time derivatives of the energy integrals given in Section 4 (left) and their equivalents obtained by evaluating the rate of mechanical work over the computational domain (right). The numbers 1 and 2 correspond to deformation (blue) and rotational (red) powers, respectively. The gravitational power is split into individual terms derived in Section 4.1. Number 3 (green) denotes the contributions due to the deformation in a reference gravitational field of potential V0 while number 4 (yellow) corresponds to the power associated with incremental gravitational potential Φ. (b) Numerical test of the formulae in panel a performed for the benchmark case by Spada et al. (2011) (see also Sections 3 and 4.4). The coloured curves, obtained by time integration of the power formulae given in panel (a) (right), match the black curves computed directly from the corresponding energy formula. The colours and the numbers correspond to those in panel (a). Figure A1. View largeDownload slide (a) Time derivatives of the energy integrals given in Section 4 (left) and their equivalents obtained by evaluating the rate of mechanical work over the computational domain (right). The numbers 1 and 2 correspond to deformation (blue) and rotational (red) powers, respectively. The gravitational power is split into individual terms derived in Section 4.1. Number 3 (green) denotes the contributions due to the deformation in a reference gravitational field of potential V0 while number 4 (yellow) corresponds to the power associated with incremental gravitational potential Φ. (b) Numerical test of the formulae in panel a performed for the benchmark case by Spada et al. (2011) (see also Sections 3 and 4.4). The coloured curves, obtained by time integration of the power formulae given in panel (a) (right), match the black curves computed directly from the corresponding energy formula. The colours and the numbers correspond to those in panel (a). In the following, we consider that v≅∂u/∂t under the assumption of small deformations. In Appendices A4–A6, we use two modifications of the Reynolds transport theorem (e.g. Gurtin 1981, p. 91, eq. 8):   $$\frac{{\rm d}}{{\rm d} t}\int _{v(t)}\rho f\, {\rm d}v = \int _{v(t)}\rho \frac{\mathcal {D}f}{\mathcal {D}t}\, {\rm d}v\, ,$$ (A3)where v(t) is the total volume of the deformed body at time t and $$\mathcal {D}/\mathcal {D}t$$ denotes the material derivative,   $$\frac{\mathcal {D}}{\mathcal {D}t}=\frac{\partial }{\partial t}+\boldsymbol {v}\cdot \nabla \, .$$ (A4) A1 Term − τ: ∇v Since the body is incompressible (trace of ∇v is zero) and the Cauchy stress tensor is symmetric, τ: ∇v can be rewritten as   $$\boldsymbol {\tau }:\nabla \boldsymbol {v}=\boldsymbol {\tau }^d:\frac{1}{2}\!\left(\nabla \boldsymbol {v}+(\nabla \boldsymbol {v})^{\rm T}\right)\, .$$ (A5)The strain rate $$\frac{1}{2}\!\left(\nabla \boldsymbol {v}+(\nabla \boldsymbol {v})^{\rm T}\right)$$ can be expressed from the constitutive equation (3):   $$\frac{1}{2}\!\left(\nabla \boldsymbol {v}+(\nabla \boldsymbol {v})^{\rm T}\right) \cong \frac{\partial }{\partial t}\left[\frac{1}{2}\!\left(\nabla \boldsymbol {u}+(\nabla \boldsymbol {u})^{\rm T}\right)\right] = \frac{1}{2\mu }\frac{\partial \boldsymbol {\tau }^d}{\partial t}+\frac{\boldsymbol {\tau }^d}{2\eta }\, .$$ (A6)The integral of τ: ∇v then takes the form   $$\int _{{\rm S}_a\!-\!{\rm S}_b}\!\!\boldsymbol {\tau }:\nabla \boldsymbol {v}\, \mathrm{d}v \cong \int _{{\rm S}_a\!-\!{\rm S}_b}\!\!\left[\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{2\eta } + \frac{\boldsymbol {\tau }^d}{2\mu }:\frac{\partial \boldsymbol {\tau }^d}{\partial t}\right]\, \mathrm{d}v= \int _{{\rm S}_a\!-\!{\rm S}_b}\!\!\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{2\eta }\, \mathrm{d}v+ \frac{\rm d}{{\rm d}t}\int _{{\rm S}_a\!-\!{\rm S}_b}\!\! \frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{4\mu }\, \mathrm{d}v\, .$$ (A7)Assuming that the liquid core is inviscid and neglecting the dissipation and elastic energy within the small volumes corresponding to topographic anomalies, we can replace the integrals over Sa − Sb on the right-hand side of eq. (A7) by integrals over v(t), resulting in the relation:   $$\int _{{\rm S}_a\!-\!{\rm S}_b}\boldsymbol {\tau }:\nabla \boldsymbol {v}\, \mathrm{d}v \cong \int _{v(t)}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{2\eta }+ \frac{{\rm d}}{{\rm d} t}\int _{v(t)}\frac{\boldsymbol {\tau }^d:\boldsymbol {\tau }^d}{4\mu }\, \mathrm{d}v=\boldsymbol {\dot{E}}_\mathrm{diss} + \boldsymbol {\dot{E}}_\mathrm{el}\, ,$$ (A8)where we used eqs (20) and (21) from Section 4. Eq. (A8) suggests that the power represented by the term τ: ∇v is dissipated as heat or stored as elastic deformation energy in the body. For numerical validation of the relationship A8, see Fig. A1(b), curve 1. A2 Term ∇ · (τ · v) The volume integral of ∇ · (τ · v) can be transformed into a surface integral using the Gauss theorem:   $$\int _{\mathrm{S}_a\!-\!\mathrm{S}_b}\nabla \cdot (\boldsymbol {\tau }\cdot \boldsymbol {v})\, \mathrm{d}v= \int _{\partial \mathrm{S}_a}\boldsymbol {v}\cdot \boldsymbol {\tau }\cdot \boldsymbol {e}_r\, \mathrm{d}s - \int _{\partial \mathrm{S}_b}\boldsymbol {v}\cdot \boldsymbol {\tau }\cdot \boldsymbol {e}_r\, \mathrm{d}s \, .$$ (A9)Evaluating τ · er from boundary conditions (9) and (10) and omitting the term σL (see the remark in the beginning of Section 4.1), we obtain   $$\int _{\mathrm{S}_a-\mathrm{S}_b}\nabla \cdot (\boldsymbol {\tau }\cdot \boldsymbol {v})\, \mathrm{d}v= \int _{\partial \mathrm{S}_a}\boldsymbol {v}\cdot \left(u_r[\rho _0]_a\boldsymbol {g}_0\right)\, \mathrm{d}s + \int _{\partial \mathrm{S}_b}\boldsymbol {v}\cdot \left(u_r[\rho _0]_b\boldsymbol {g}_0+p_{\rm c}\boldsymbol {e}_r\right)\mathrm{d}s\, ,$$ (A10)where pc is the pressure in the core, pc = −ρc(V0 + Φ + Ψ). The first surface integral on the right-hand side of eq. (A10) can be modified as   \begin{eqnarray} \int _{\partial \mathrm{S}_a}\boldsymbol {v}\cdot \left(u_r[\rho _0]_a\boldsymbol {g}_0\right)\, \mathrm{d}s &\cong & \int _{\partial \mathrm{S}_a}\frac{\partial \boldsymbol {u}}{\partial t}\cdot \left(u_r[\rho _0]_a\boldsymbol {g}_0\right)\, \mathrm{d}s \nonumber\\ &=& -[\rho _0]_a\, g_0(a) \int _{\partial \mathrm{S}_a}\frac{1}{2}\frac{\partial u^2_r}{\partial t} \, \mathrm{d}s =-\frac{1}{2}[\rho _0]_a\, g_0(a) \frac{\rm d}{{\rm d}t}\int _{\partial \mathrm{S}_a} u^2_r \, \mathrm{d}s\, . \end{eqnarray} (A11)Analogously,   $$\int _{\partial \mathrm{S}_b}\boldsymbol {v}\cdot \left(u_r[\rho _0]_b\boldsymbol {g}_0\right)\, \mathrm{d}s \cong -\frac{1}{2}[\rho _0]_b\, g_0(b) \frac{\rm d}{{\rm d}t}\int _{\partial \mathrm{S}_b} u^2_r \, \mathrm{d}s \, .$$ (A12)To summarize,   $$\int _{\mathrm{S}_a-\mathrm{S}_b}\nabla \cdot (\boldsymbol {\tau }\cdot \boldsymbol {v})\, \mathrm{d}v= -\frac{1}{2}\frac{\rm d}{{\rm d}t} \left([\rho _0]_a\, g_0(a) \int _{\partial \mathrm{S}_a} u^2_r \, \mathrm{d}s +[\rho _0]_b\, g_0(b) \int _{\partial \mathrm{S}_b} u^2_r \, \mathrm{d}s\right) +\int _{\partial \mathrm{S}_b}p_{\rm c}\boldsymbol {v}\cdot \boldsymbol {e}_r\, \mathrm{d}s\, ,$$ (A13)or, in terms of surface integrals   \begin{eqnarray} {\int _{\partial \mathrm{S}_a}\boldsymbol {v}\cdot \boldsymbol {\tau }\cdot \boldsymbol {e}_r\, \mathrm{d}s -\int _{\partial \mathrm{S}_b}\boldsymbol {v}\cdot (\boldsymbol {\tau }+p_{\rm c}\mathcal {I})\cdot \boldsymbol {e}_r\, \mathrm{d}s =- \frac{1}{2}\frac{\rm d}{{\rm d}t} \left([\rho _0]_a\, g_0(a) \int _{\partial \mathrm{S}_a} u^2_r \, \mathrm{d}s +[\rho _0]_b\, g_0(b) \int _{\partial \mathrm{S}_b} u^2_r\, \mathrm{d}s\right)}, \end{eqnarray} (A14)where $$\mathcal {I}$$ is the identity tensor. Comparison of eq. (A13) with eq. (36) suggests that the integral of ∇ · (τ · v) over the computation domain corresponds to the changes of the gravitational energy due to the deformations of boundaries in the gravitational field of potential V0. For the numerical confirmation of eq. (A14), see Fig. A1(b), curves 3a,b. It is worth noting that the integral of ∇ · (τ · v) over the whole body v(t) would be zero. Using the Gauss theorem and considering that the real surface is stress-free, we obtain   $$\int _{v(t)}\nabla \cdot (\boldsymbol {\tau }\cdot \boldsymbol {v})\, \mathrm{d}v=\int _{\partial {\rm S}}\boldsymbol {v}\cdot \boldsymbol {\tau }\cdot \boldsymbol {n}\, \mathrm{d}s=0\, ,$$ (A15)where ∂S now denotes the surface of the deformed body and n is the normal to the surface (see also Appendix A7). A3 Term − (u · ∇ρ0) g0 · v Considering that $$\boldsymbol {v}\!\cong \!\partial {\boldsymbol u}/\partial t$$, g0 = −g0er and the functions g0 and ρ0 do not depend on time, we obtain   \begin{eqnarray} {-\int _{\mathrm{S}_a-\mathrm{S}_b}(\boldsymbol {u}\cdot \nabla \rho _0)\, \boldsymbol {g}_0\cdot \boldsymbol {v}\, \mathrm{d}v \cong \int _{\mathrm{S}_a-\mathrm{S}_b}u_r\frac{\partial u_r}{\partial t}\frac{\mathrm{d}\rho _0}{\mathrm{d}r}\, g_0\, \mathrm{d}v = \frac{1}{2}\frac{\rm d}{{\rm d} t}\int _{\mathrm{S}_a-\mathrm{S}_b} u_r^2 g_0\frac{\mathrm{d}\rho _0}{\mathrm{d}r}\, \mathrm{d}v = \frac{1}{2}\frac{\rm d}{{\rm d} t}\int _{b}^a \frac{\mathrm{d}\rho _0}{\mathrm{d}r}\int _{\partial {\rm S}_r} g_0u_r^2\, \mathrm{d}s\, \mathrm{d}r}\, , \end{eqnarray} (A16)where the resulting integral can be identified with the rate of change of the gravitational energy due to volumetric density changes in the gravitational field of potential V0 (cf. the last integral in eq. (36); for the numerical comparison, see Fig. A1b, curve 3c). A4 Auxiliary formulae for gravitational and centrifugal potential Before we proceed to the remaining terms in eq. (A2) we derive two auxiliary formulae for the gravitational and centrifugal potentials valid for the whole volume v(t) of the deformed body at time t. We first consider the energy balance of a body in an inertial frame subject to gravitational body force −ρ∇V only, where V is the gravitational potential generated by the body. If the surface of the body is free,   $$\frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\left(\rho \epsilon +\frac{1}{2}\rho \boldsymbol {v}\cdot \boldsymbol {v}\right) \mathrm{d}v = -\int _{v(t)}\rho \nabla V\cdot \boldsymbol {v} \, \mathrm{d}v\, ,$$ (A17)where the expression on the left-hand side corresponds to the time change of the internal and kinetic energies, ε denotes the internal energy density, and the integral on the right-hand side is the total power of the gravitational force. Using the definition of material derivative eq. (A17) can be rearranged as   $$\frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\left(\rho \epsilon +\frac{1}{2}\rho \boldsymbol {v}\cdot \boldsymbol {v}\right)\, \mathrm{d}v + \int _{v(t)} \rho \frac{1}{2}\frac{\mathcal {D}V}{\mathcal {D}t}\mathrm{d}v = \int _{v(t)}\rho \left(\frac{\partial V}{\partial t}-\frac{1}{2}\frac{\mathcal {D}V}{\mathcal {D}t}\right)\, \mathrm{d}v\, .$$ (A18)The time derivative in the second integral on the left-hand side of eq. (A18) can be pulled in front of the integral using the modified Reynolds theorem, eq. (A3), and the integrand on the right-hand side can be reformulated by substituting for the material derivative of V:   $$\rho \left(\frac{\partial V}{\partial t}-\frac{1}{2}\frac{\mathcal {D}V}{\mathcal {D}t}\right)= \rho \left(\frac{\partial V}{\partial t}-\frac{1}{2}\frac{\partial V}{\partial t}-\frac{1}{2}\nabla V\cdot \boldsymbol {v}\right)= \frac{1}{2}\rho \left(\frac{\partial V}{\partial t}-\boldsymbol {v}\cdot \nabla V\right)\, .$$ (A19)Eq. (A18) then gives   $$\frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\left(\rho \epsilon +\frac{1}{2}\rho \boldsymbol {v}\cdot \boldsymbol {v}+\frac{1}{2}\rho V\right)\, \mathrm{d}v = \frac{1}{2}\int _{v(t)}\rho \left(\frac{\partial V}{\partial t}-\nabla V\cdot {\boldsymbol v}\right)\, \mathrm{d}v\, .$$ (A20)On the left-hand side we now have the sum of all energies in the system, including the gravitational energy. The total energy must be conserved which implies that the right-hand side of eq. (A20) is equal to zero. Consequently,   $$\int _{v(t)}\rho \nabla V\cdot \boldsymbol {v}\, \mathrm{d}v = \int _{v(t)}\rho \frac{\partial V}{\partial t}\, \mathrm{d}v\, .$$ (A21) A similar formula can be derived for the centrifugal potential Ψ. In this case, however, we must consider the energy balance in a rotating frame which is fixed with respect to the body and consider all relevant forces that act in the rotating frame. In analogy to eq. (A17), but this time not considering any gravitational forcing, we can write   $$\frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\left(\rho \epsilon +\frac{1}{2}\rho \boldsymbol {v}\cdot \boldsymbol {v}\right)\, \mathrm{d}v = \int _{v(t)}\left(-\rho \nabla \Psi -\rho \boldsymbol {\omega }\times \boldsymbol {v}-\rho \frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t}\times \boldsymbol {r}\right)\cdot \boldsymbol {v}\, \mathrm{d}v,$$ (A22)where the terms in parentheses on the right-hand side represent the centrifugal, Coriolis and Euler body force, respectively. The scalar product of the Coriolis force with the velocity vector is zero. The power of the Euler force can be rearranged using the vector algebra identity a · (b × c) = b · (c × a):   $$\int _{v(t)}\left(-\rho \frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t}\times \boldsymbol {r}\right)\cdot \boldsymbol {v}\, \mathrm{d}v=-\frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t}\cdot \int _{v(t)}\boldsymbol {r}\times (\rho \boldsymbol {v})\, \mathrm{d}v=-\frac{\mathrm{d}\boldsymbol {\omega }}{\mathrm{d}t}\cdot \boldsymbol {h}\, ,$$ (A23)where h is the relative angular momentum. In the Tisserand frame (Munk & MacDonald 1960), h is equal to zero by definition, and the total power of the Euler force is thus zero in such frame. Following the same procedure as above (eqs A18–A20 where V is replaced by Ψ), we obtain   $$\frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\left(\rho \epsilon +\frac{1}{2}\rho \boldsymbol {v}\cdot \boldsymbol {v}+\frac{1}{2}\rho \Psi \right)\, \mathrm{d}v = \frac{1}{2}\int _{v(t)}\rho \left(\frac{\partial \Psi }{\partial t}-\nabla \Psi \cdot {\boldsymbol v}\right)\, \mathrm{d}v\, .$$ (A24)On the left-hand side we have the sum of all energies of the now considered system, when viewed from the inertial reference frame. The time derivative of this sum must again be zero, giving   $$\int _{v(t)}\rho \nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v = \int _{v(t)}\rho \frac{\partial \Psi }{\partial t}\, \mathrm{d}v\, .$$ (A25) A5 Term − (ρ0∇Φ) · v We first evaluate the integral of − (ρ0∇Φ) · v over the volume v(t) of the deformed body and then we transform it into an integral over the computational domain Sa − Sb. The potential V in auxiliary formula (A21) can be expressed as the sum of the reference potential V0, independent of time, and the increment Φ due to the deformation, V(t) = V0 + Φ(t). Considering that ∂V0/∂t = 0,   $$\int _{v(t)}\rho \nabla (V_0+\Phi )\cdot \boldsymbol {v}\, \mathrm{d}v = \int _{v(t)}\rho \frac{\partial \Phi }{\partial t}\, \mathrm{d}v\, .$$ (A26)The term ρ∇V0 · v can be moved to the right-hand side and replaced by $$\rho \mathcal {D}V_0/\mathcal {D}t$$:   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v = \int _{v(t)}\rho \left(\frac{\partial \Phi }{\partial t}-\frac{\mathcal {D}V_0}{\mathcal {D}t}\right)\, \mathrm{d}v\, .$$ (A27)Adding ∫v(t)ρ∇Φ · v dv to both sides of eq. (A27), dividing the equation by 2 and invoking the expression for the material derivative of Φ, we obtain   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v = \frac{1}{2}\int _{v(t)}\rho \frac{\mathcal {D}}{\mathcal {D}t}\left(\Phi -V_0 \right)\, \mathrm{d}v\, ,$$ (A28)where the time derivative on the right-hand can be pulled in front of the integral using the Reynolds theorem (A3):   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v = \frac{1}{2}\frac{{\rm d}}{{\rm d} t}\int _{v(t)}\rho \left(\Phi -V_0 \right)\, \mathrm{d}v\, .$$ (A29)The integral on the right-hand side of eq. (A29) is formally the same as in eq. (22) except that the sign of V0 is opposite. Thus, the integral can be evaluated in an analogous manner as described in Section 4.1, giving   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v = \frac{1}{2}\sum _{i=1}^N\frac{\mathrm{d}}{\mathrm{d}t}\int _{h_i}[\rho _0]_i\Phi \mathrm{d}v \cong \frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{1}{2}\sum _{i=1}^N[\rho _0]_i \int _{\partial {\rm S}_i}\Phi u_r\, \mathrm{d}s\right]$$ (A30)for a layered model and   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v \cong \frac{{\rm d}}{{\rm d}t}\left[\frac{1}{2}\left( [\rho _0]_a\int _{\partial {\rm S}_a}\!\!u_r\Phi \mathrm{d}s +[\rho _0]_b\int _{\partial {\rm S}_b}\!\!u_r\Phi \mathrm{d}s -\int _{b}^a \frac{\mathrm{d}\rho _0}{\mathrm{d}r}\int _{\partial {\rm S}_r}\!\! u_r\Phi \, \mathrm{d}s\, \mathrm{d}r\right)\right]$$ (A31)for a continuous density profile. The terms on the right-hand sides of eqs (A30) and (A31) correspond to the gravitational energy due to potential Φ induced by the deformation (cf. eqs 31, 34 and 35). Now we transform the integral on the left-hand side of eq. (A29) into an integral over the computational domain Sa − Sb:   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v \cong \int _{{\rm S}_a-{\rm S}_b}\rho _0\nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v+ \rho _{\rm c}\int _{{\rm S}_b}\nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v\, ,$$ (A32)where we omitted the second order term −δρ∇Φ (cf. eq. 8) and neglected the contributions of ρ∇Φ · v within the small volumes corresponding to topographic anomalies. The volume integral over Sb can be transformed into a surface integral using the Gauss theorem. Considering that ∇Φ · v = ∇ · (Φv) for an incompressible body, we get   $$\int _{v(t)}\rho \nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v \cong \int _{{\rm S}_a-{\rm S}_b}\rho _0\nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v+ \rho _{\rm c}\int _{\partial {\rm S}_b}(\Phi \boldsymbol {v})\cdot \boldsymbol {e}_r\, \mathrm{d}v\, .$$ (A33)Combining eqs (A33) and (A30) provides the sought relation between the power of ρ0∇Φ over the computational domain and the associated energy rate:   $$\int _{{\rm S}_a-{\rm S}_b}\rho _0\nabla \Phi \cdot \boldsymbol {v}\, \mathrm{d}v+ \rho _{\rm c}\int _{\partial {\rm S}_b}(\Phi \boldsymbol {v})\cdot \boldsymbol {e}_r\, \mathrm{d}v \cong \frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{1}{2}\sum _{i=1}^N[\rho _0]_i \int _{\partial {\rm S}_i}\Phi u_r\, \mathrm{d}s\right].$$ (A34)An analogous expression is obtained for body with a continuous density profile by combining eqs (A33) and (A31). The validity of eq. (A34) is numerically demonstrated by curve 4 in Fig. A1(b). A6 Term − (ρ0∇Ψ) · v Using eq. (A25) and following the analogous procedure as above, we get   \begin{eqnarray} {\int _{v(t)}\rho \nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v = \int _{v(t)}\rho \frac{\partial \Psi }{\partial t}\, \mathrm{d}v = \frac{1}{2}\int _{v(t)}\rho \left(\frac{\partial \Psi }{\partial t}+\nabla \Psi \cdot \boldsymbol {v}\right)\, \mathrm{d}v = \frac{1}{2}\int _{v(t)}\rho \frac{\mathcal {D}\Psi }{\mathcal {D}t}\, \mathrm{d}v = \frac{\mathrm{d}}{\mathrm{d}t}\int _{v(t)}\frac{1}{2}\rho \Psi \, \mathrm{d}v = \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{2}\boldsymbol {\omega }\cdot \boldsymbol {I}\cdot \boldsymbol {\omega }\right)}\, . \end{eqnarray} (A35)The integral of ρ∇Ψ · v over v(t) can be expressed as (cf. eq. A32)   $$\int _{v(t)}\rho \nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v \cong \int _{{\rm S}_a-{\rm S}_b}\rho _0\nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v+ \rho _{\rm c}\int _{{\rm S}_b}\nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v\, .$$ (A36)Combining eqs (A35) and (A36) we obtain   $$\int _{{\rm S}_a-{\rm S}_b}\rho _0\nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v+ \rho _{\rm c}\int _{{\rm S}_b}\nabla \Psi \cdot \boldsymbol {v}\, \mathrm{d}v= \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{1}{2}\boldsymbol {\omega }\cdot \boldsymbol {I}\cdot \boldsymbol {\omega }\right)\, ,$$ (A37)see also Fig. A1(b), curve 2. Note that the derivations of eqs (A34) and (A37) rely on the auxiliary formulae from Appendix A4, and as such they are only valid if the energy of the studied system is conserved. In energetically ill-posed problems eqs (A34) and (A37) must not be satisfied, as illustrated in Section 4.5 in the main text (cf. Fig. 10). A7 Term ρ0g0 · v Since ρ0 and g0 are spherically symmetric and the body is incompressible ($$\int _{\partial {\rm S}_r}\boldsymbol {v}\cdot \boldsymbol {e}_r\, {\rm d}s\!=\!0$$), the integral of the last term in eq. (A1) over the computational domain is equal to zero:   $$\int _{{\rm S}_a-{\rm S}_b}\rho _0\boldsymbol {g}_0\cdot \boldsymbol {v}\, \mathrm{d}v=-\int _{{\rm S}_a-{\rm S}_b}\rho _0 g_0 \boldsymbol {v}\cdot \boldsymbol {e}_r\, \mathrm{d}v=-\int _b^a\rho _0g_0\int _{\partial {\rm S}_r}\boldsymbol {v}\cdot \boldsymbol {e}_r\, {\rm d}s\, {\rm d}r=0\, .$$ (A38)Note, however, that the integral of ρg0 · v over the real volume of the body v(t) is non-zero, even if the density is constant within the mantle shell. For a body with ∇ρ0 equal to zero within the shell the integral is equal to the right-hand side of eq. (A10):   $$\int _{v(t)}\rho \boldsymbol {g}_0\cdot \boldsymbol {v}\, \mathrm{d}v \cong \int _{\partial \mathrm{S}_a}\boldsymbol {v}\cdot \left(u_r[\rho _0]_a\boldsymbol {g}_0\right)\, \mathrm{d}s + \int _{\partial \mathrm{S}_b}\boldsymbol {v}\cdot \left(u_r[\rho _0]_b\boldsymbol {g}_0+p_{\rm c}\boldsymbol {e}_r\right)\mathrm{d}s\, \, .$$ (A39)This is because in this integral there are contributions from the non-zero Eulerian density increments δρ = ±[ρ0]a, b in the vicinity of the boundaries (see Fig. 3). If the computational domain is chosen regular (Sa − Sb) as in the present study, these volumetric density increments are condensed into surface densities and we encounter them in Appendix A2 in the form of boundary tractions (9) and (10). For a body with ∇ρ0 ≠ 0 inside the mantle shell the term ρg0 · v also includes a contribution from the internal Eulerian density increment − u · ∇ρ0, which was discussed in Appendix A3. A8 Remark Inspection of the left-hand side of Fig. A1(a) suggests that the formulae used to evaluate the energy balance in the main text can be replaced by the corresponding formulae for the rates of the mechanical work evaluated over the computational domain Sa − Sb. However, there are three integrals over ∂Sb on the right-hand side of Fig. A1(a) which are related to the body force −ρc∇(V0 + Φ + Ψ) in the core. The sum of these terms is equal to zero,   $$-\int _{\partial {\rm S}_b}(\rho _{\rm c}\Psi \boldsymbol {v})\cdot \boldsymbol {e}_r\, {\rm d}s -\int _{\partial {\rm S}_b}p_{\rm c}\boldsymbol {v}\cdot \boldsymbol {e}_r\, {\rm d}s -\int _{\partial {\rm S}_b}(\rho _{\rm c}\Phi \boldsymbol {v})\cdot \boldsymbol {e}_r\, {\rm d}s=0\, ,$$ (A40)since pc = −ρc(V0 + Φ + Ψ) and   $$\int _{\partial {\rm S}_b}\rho _{\rm c}V_0\boldsymbol {v}\cdot \boldsymbol {e}_r\, {\rm d}s= \rho _{\rm c}V_0\int _{\partial {\rm S}_b}\boldsymbol {v}\cdot \boldsymbol {e}_r\, {\rm d}s=0\, .$$ (A41) APPENDIX B: SECOND-ORDER CORRECTION TO LLE The third component of the NLE is   \begin{eqnarray} 0 = \frac{\mathrm{d}}{{\rm d}t}\left(\boldsymbol {I}\cdot \boldsymbol {\omega }\right)_3 + \varepsilon _{3jk}\omega _j\left(\boldsymbol {I}\cdot \boldsymbol {\omega }\right)_k\, , \end{eqnarray} (B1)where ε is the Levi-Civita symbol. Using the traditional notation,   $$\boldsymbol {\omega }=\boldsymbol {\omega }_0 + |\boldsymbol {\omega }_0|\boldsymbol {m} \quad {\rm {and}} \quad \boldsymbol {I}={\left(\begin{array}{ccc}A+c_{11} & c_{12} & c_{13} \\ c_{12} & A + c_{22} & c_{23} \\ c_{13} & c_{12} & C+c_{33} \end{array}\right)} ,$$ (B2)with ω0 equal to [0, 0, ω0] and cij, mi denoting small quantities, and inserting from eq. (B2) to eq. (B1), we obtain   \begin{eqnarray} 0 &\cong & \frac{\mathrm{d}}{{\rm d}t}[c_{33} + c_{13} m_1 + c_{23} m_2 + c_{33} m_3 + C m_3 ] + \omega _0 m_1 c_{23} - \omega _0 m_2 c_{13} \nonumber \\ &=& (\dot{c}_{33} + \dot{c}_{13} m_1 + \dot{c}_{23} m_2 + \dot{c}_{33} m_3 + c_{13} \dot{m}_1 + c_{23} \dot{m}_2 + c_{33} \dot{m}_3 + C \dot{m}_3 ) + \omega _0(m_1 c_{23} - m_2 c_{13}), \end{eqnarray} (B3)where third order terms are neglected, the dot abbreviates time derivative and $$\dot{C}$$ and $$\dot{\omega }_0$$ are zero, since both C and ω0 are constants. Recalling the LLE for m1 and m2 (Munk & MacDonald 1960)   \begin{eqnarray} \omega _0(C-A)\left(m_1 - \frac{A\, \dot{m}_2}{(C-A)\omega _0}\right) = \omega _0 c_{13} + \dot{c}_{23} \end{eqnarray} (B4)  \begin{eqnarray} \omega _0(C-A)\left(m_2 + \frac{A\, \dot{m}_1}{(C-A)\omega _0}\right) = \omega _0 c_{23} - \dot{c}_{13}, \end{eqnarray} (B5)we can insert for $$\dot{c}_{23}$$ and $$\dot{c}_{13}$$ in eq. (B3) from eqs (B4) and (B5), introducing only a third order error, because $$\dot{c}_{23}$$ and $$\dot{c}_{13}$$ are both multiplied by a small quantity in eq. (B3). This leads to   \begin{eqnarray} 0 &\cong & \dot{c}_{33} + C \dot{m}_3 - A m_1 \dot{m}_1 - A m_2 \dot{m}_2 + \dot{c}_{33} m_3 + \omega _0(m_1 c_{23} - m_2 c_{13}) \nonumber \\ &&+\,\omega _0[ c_{23}m_1 - (C-A)m_1 m_2 - c_{13}m_2 + (C-A)m_1 m_2 + c_{13} \dot{m}_1 + c_{23} \dot{m}_2 + c_{33} \dot{m}_3] \, . \end{eqnarray} (B6)The fifth term on the first line of eq. (B6) can be expressed using the LLE for m3, $$\dot{c}_{33} = - C \dot{m}_3$$, introducing only a third order error:   \begin{eqnarray} 0 &\cong & \dot{c}_{33} + C \dot{m}_3 - A m_1 \dot{m}_1 - A m_2 \dot{m}_2 - C \dot{m}_3 m_3 + \omega _0[ 2 c_{23}m_1 - 2 c_{13}m_2 + c_{13} \dot{m}_1 + c_{23} \dot{m}_2 + c_{33} \dot{m}_3] \, . \end{eqnarray} (B7)This can be rewritten as   \begin{eqnarray} 0 &\cong & \dot{c}_{33} + C {\left( m_3 - \frac{m_1^2+m_2^2+m_3^2}{2} \right)}^{\boldsymbol {\cdot }} + \omega _0[ 2 c_{23}m_1 - 2 c_{13}m_2 + c_{13} \dot{m}_1 + c_{23} \dot{m}_2 + c_{33} \dot{m}_3] + (C-A)(\dot{m_1}m_1+\dot{m_2}m_2) , \end{eqnarray} (B8)which provides a second-order correction to the vertical component of the LLE. The second equality in eq. (42) is obtained when the relative importance of the second-order terms in eq. (B8) are assessed: the terms on the second line are negligible due to the multiplication by factors ω0 and C − A, which are both much smaller than factor C. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Feb 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations