*A Viscoelastic Model for the Earth*

In accordance with the previously cited works the Earth has been modeled has a linear viscoelastic spherical shell Ω ⊂ R<sup>3</sup> since the dynamics due to the glacial isostatic adjustment does not involve the core of the planet. Following the approach presented in [26] we assume that the viscoelastic stress tensor is given *σ* = *σ<sup>e</sup>* − *q*, where *σ<sup>e</sup>* is the elastic stress and tensor *q* is an internal variable used to model the effects of viscosity. Denoting with **u** the displacement field, the elastic stress tensor *σ<sup>e</sup>* is given by the sum of the deviatoric and the volumetric parts

$$
\sigma^t = 2\mu \, e(\mathbf{u}) - pI\_{,t}
$$

where *I* denotes the identity matrix. The deviatoric part is the product between the shear modulus *μ* and the deviatoric strain *e*(**u**) defined as

$$\varepsilon(\mathbf{u}) := \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^t) - \frac{1}{3}(\nabla \cdot \mathbf{u})I.$$

The volumetric part depends on the pressure *p* which is defined by the equation

$$\nabla \cdot \mathbf{u} + \frac{3}{2\mu} \frac{1 - 2\nu}{1 + \nu} p = 0,\tag{1}$$

where *ν* is the Poisson ratio.

> The internal variable *q* is defined by the evolution equation

$$\begin{cases} \dot{q} + \frac{1}{\tau} q = \frac{1}{\tau} 2\mu \,\varepsilon(\mathbf{u}),\\\\ q(0) = 0, \end{cases} \tag{2}$$

where *q*˙ denotes the derivative of the quantity *q* with respect to time and *τ* is called relaxation time and it is related to the viscosity *η* through the relation *τ* = *η* 2*μ*.

The Equation (2) can be rewritten in the integral form as:

$$q(t) = \int\_0^t \frac{1}{\pi} e^{-\frac{t-s}{\tau}} 2\mu \, e(\mathbf{u}(s)) \, ds.$$

By means of the previous expressions the viscoelastic stress tensor *σ* is evaluated through a convolution integral defined as

$$
\sigma(t) = -pI + \int\_0^t e^{-\frac{t-s}{\overline{\tau}}} 2\mu \, e(\dot{\mathbf{u}}(s)) \, ds \, \mathcal{I}
$$

This equation coupled with the Equation (1) and with the equation of conservation of linear momentum give us a system of partial differential equations that describe the motion of a viscoelastic body:

$$\begin{cases} \nabla \cdot \boldsymbol{\sigma} + \mathbf{f} = 0 & \text{in } \Omega, \\\\ \nabla \cdot \mathbf{u} + \frac{3}{2\mu} \frac{1 - 2\nu}{1 + \nu} p = 0 & \text{in } \Omega, \\\\ \boldsymbol{\sigma} = -p\boldsymbol{I} + \int\_0^t e^{-\frac{\boldsymbol{I} - \boldsymbol{s}}{\boldsymbol{\tau}}} 2\mu \, e(\dot{\mathbf{u}}(\boldsymbol{s})) \, ds & \text{in } \Omega. \end{cases} \tag{3}$$

The unknowns of this system of equations are the displacement field **u**, the pressure field *p* and the stress tensor field *σ*. The volumetric force field **f** is the gravitational force field, how this term has been modeled will be discussed later.

Since our domain is a spherical shell its boundary is the union of two connected components: the inner and the outer surfaces of the shell. The outer surface Γout is the surface of the Earth and on this portion of the boundary we assume to know the history of the load due to the presence of the ice or other type of loads, like sediments. According to these data the following boundary condition is assumed: *σ***n** = **s**load on Γout, where **n** denotes the outer normal defined on the boundary. The inner surface Γin of the shell is chosen in such a way it corresponds to the core-mantle boundary (about 2900 km of depth). On this portion of the boundary the displacement **u** is assumed to be equal to zero, since the deformation due to the glacial isostatic adjustment involves only the shallow part of the mantle, until few hundreds of kilometer, namely **u** = 0 on Γin.

The force field **f** is given by the product between the density field *ρ* and the acceleration gravity **g**. Since we are dealing with a model of the whole Earth we cannot assume a constant value for the gravity and its value must be computed in accordance with the density field using the Gauss's law for the gravity

$$\begin{cases} \mathbf{g} = -\nabla \phi, \\ \nabla \cdot \mathbf{g} + 4\pi \mathbf{G} \rho = 0, \end{cases} \tag{4}$$

where *G* is the universal gravitational constant and the boundary conditions are **g** = −*g*surface**<sup>n</sup>** on Γout and *φ* = 0 on Γin. Even though the displacement **u** is small if compared to the characteristic length of problem (1 km vs. 6371 km) the gravity acceleration acting on the point change in accordance with the displacement field, so the force field **f** is given by *ρ***g**(**x** + **u**) ≈ *ρ***g**(**x**) + *ρ*(∇**xg**)**<sup>u</sup>**. From a physical point of view the first term *ρ***g**(**x**) is a static component, it does not change in time, on the other side the second term *ρ*(∇**xg**)**u** is the buoyancy term that determine the uplift and subsidence processes. This physical interpretation can be justified from a mathematical point of view. Exploiting the linearity of the problem the stress tensor *σ* is decomposed as a sum of two components: a static stress term *σ*<sup>0</sup> (which is not important for our purposes) and a dynamic term *σ<sup>d</sup>* which is the solution of Problem Equation (3) with **f** = *ρ*(∇**xg**)**u** and the boundary conditions.
