THM Thermo-Hydro-Mechanical

### **Appendix A. The Numerical Solver for the Isostatic Glacial Rebound Model**

The numerical solver for the glacial isostatic adjustment is implemented using the deal.II library [40,41]. This library provides the tools for an efficient implementation of a parallel solver based on the Discontinuous Galerkin method. Let T*h* be a subdivision of the geometric domain Ω consisting of non-overlapping hexahedra with characteristic mesh size *h*, we introduce the finite dimensional spaces

$$\begin{aligned} \mathbf{V}\_h^k &= \{ \mathbf{v} \in L^2(\Omega; \mathbb{R}^3) \, : \, \mathbf{v}|\_K \in \mathbb{P}^k(K; \mathbb{R}^3) \, \forall K \in \mathcal{T}\_h \}, \\ Q\_h^k &= \{ q \in L^2(\Omega) \, : \, q|\_K \in \mathbb{P}^k(K) \, \forall K \in \mathcal{T}\_h \}. \end{aligned}$$

Let E*h* = *<sup>K</sup>*∈T*h ∂K* denote the set of the faces of T*h* then E0*h* = E*h*\*∂*<sup>Ω</sup> is the set of internal faces. Let *e* ∈ E0*h* be a face shared by two elements *K*<sup>+</sup> and *K*<sup>−</sup>, define the unit normal vector **n**<sup>+</sup> and **n**<sup>−</sup> on the *e* pointing exterior to *K*<sup>+</sup> and *K*<sup>−</sup>, respectively. With *φ*<sup>±</sup> = *φ*|*K*<sup>±</sup> we set the average operators as

$$\{\mathbf{v}\} = \frac{\mathbf{v}^+ + \mathbf{v}^-}{2} \qquad \{q\} = \frac{q^+ + q^-}{2} \mathbf{v}$$

together with the jump operators:

$$[\mathbf{v}] = \mathbf{v}^+ \otimes \mathbf{n}^+ + \mathbf{v}^- \otimes \mathbf{n}^- \qquad [q] = q^+ \mathbf{n}^+ + q^- \mathbf{n}^-.$$

Equation (4) is discretized using the standard SIPG method [42]: find *φh* ∈ *Qkh*such that

$$\begin{split} \int\_{\Omega} \nabla\_{h} \Phi\_{\boldsymbol{h}} \cdot \nabla\_{h} \boldsymbol{\upmu}\_{\boldsymbol{h}} + 4 \pi \mathsf{G} \boldsymbol{\uprho} \boldsymbol{\upmu}\_{\boldsymbol{h}} \, d \mathbf{x} + \\ & \quad + \int\_{\mathcal{E}\_{\boldsymbol{h}}^{0}} \boldsymbol{\uprho} [\boldsymbol{\upmu}\_{\boldsymbol{h}}] \cdot [\boldsymbol{\upmu}\_{\boldsymbol{h}}] - \{\nabla\_{h} \boldsymbol{\upphi}\_{\boldsymbol{h}}\} \cdot [\boldsymbol{\upmu}\_{\boldsymbol{h}}] - [\boldsymbol{\upphi}\_{\boldsymbol{h}}] \{\nabla\_{h} \boldsymbol{\upmu}\_{\boldsymbol{h}}\} \, ds + \mathrm{B.C.} = 0 \qquad \forall \boldsymbol{\upmu}\_{\boldsymbol{h}} \in Q\_{\boldsymbol{h}}^{k}. \end{split}$$

where B.C. is the term related to the boundary conditions and *γ* is a proper penalization factor. The algebraic system obtained from this weak formulation is solved using the standard conjugate gradient method preconditioned with a geometric multigrid method, available in deal.II library [43].

The Equations (3) are discretized using a standard second order accurate, one-step and unconditionally stable scheme [26] and the integral is approximated using the mid-point formula. Using this approach the equations describing the motion of the viscoelastic Earth model can be rewritten in semi-discrete form using only the variables **u**, *p* and *h*:

$$\begin{cases} \nabla \cdot \left( 2\mu e^{-\frac{\Delta h}{2\tau}} \varepsilon (\mathbf{u}\_n - \mathbf{u}\_{n-1}) - p\_n I + e^{-\frac{\Delta h}{\tau}} h\_{n-1} \right) + \rho (\nabla\_\mathbf{x} \mathbf{g}) \mathbf{u}\_n = 0 & \text{in } \Omega, \\\\ \nabla \cdot \mathbf{u}\_n + \frac{3}{2\mu} \frac{1 - 2\nu}{1 + \nu} p\_n = 0 & \text{in } \Omega, \\\\ h\_n = 2\mu e^{-\frac{\Delta h}{2\tau}} \varepsilon (\mathbf{u}\_n - \mathbf{u}\_{n-1}) + e^{-\frac{\Delta h}{\tau}} h\_{n-1} & \text{in } \Omega. \end{cases}$$

$$\begin{cases} h\_n = 2\mu e^{-\frac{\mathbf{u}\cdot\mathbf{x}}{2\tau}}e(\mathbf{u}\_n - \mathbf{u}\_{n-1}) + e^{-\frac{\mathbf{u}\cdot\mathbf{x}}{\tau}}h\_{n-1} & \text{in } \Omega. \end{cases}$$

The first two equations of this system have the same structure of the linear elastic problem and they are discretized using the Discontinuous Galerkin scheme. In order to keep the notation simpler

ˆ

the shear modulus will be redefined as *μ*ˆ*n* = *μe*<sup>−</sup> Δ*tn* 2*τ* and all the known terms in the first equation are written as a unique term **F***n*:

$$\mathbf{F}\_{\rm ll} = \nabla \cdot \left( 2\hat{\mu}\_{\rm tr} e(\mathbf{u}\_{n-1}) - e^{-\frac{\Delta t\_{\rm tr}}{\bar{\tau}}} h\_{n-1} \right) \cdot \mathbf{f}$$

These equations can be rewritten in a full discrete form using the standard SIPG-method for the linear elasticity [44]: find (**<sup>u</sup>***h*,*n*, *ph*,*<sup>n</sup>*) ∈ **V***k*+<sup>1</sup> *h* × *Qkh* such that

$$\begin{cases} a(\mathbf{u}\_{h,\boldsymbol{\nu}}, \mathbf{v}\_{h}) + b(\mathbf{v}\_{h\boldsymbol{\nu}}, p\_{h,\boldsymbol{\nu}}) = \int\_{\Omega} \mathbf{F}\_{\boldsymbol{\nu}} \cdot \mathbf{v}\_{h} \, d\mathbf{x} + \text{B.C.} & \forall \mathbf{v}\_{h} \in \mathbf{V}\_{h}^{k+1}, \\ b(\mathbf{u}\_{h,\boldsymbol{\nu}'}q\_{h}) - c(p\_{h,\boldsymbol{\nu}'}q\_{h}) = \text{B.C.} & \forall q\_{h} \in Q\_{h}^{k}. \end{cases}$$

where the bilinear forms are defined by

$$\begin{split} a(\mathbf{u}\_{h\prime},\mathbf{v}\_{h}) &= \int\_{\Omega} 2\hat{\mu}\_{h} e(\mathbf{u}\_{h}) : e(\mathbf{v}\_{h}) + \rho(\nabla\_{\mathbf{x}}\mathbf{g}) \mathbf{u}\_{h} \cdot \mathbf{v}\_{h} \, d\mathbf{x} + \\ &\quad + \int\_{\mathcal{E}\_{h}^{0}} \gamma[\mathbf{u}\_{h}] : [\mathbf{v}\_{h}] - \{2\hat{\mu}\_{h} e(\mathbf{u}\_{h})\} : [\mathbf{v}\_{h}] - [\mathbf{u}\_{h}] : \{2\hat{\mu}\_{h} e(\mathbf{v}\_{h})\} \, ds\_{\prime} \\ b(\mathbf{u}\_{h\prime}, q\_{h}) &= -\int\_{\Omega} \nabla \cdot \mathbf{u}\_{h} q\_{h} \, d\mathbf{x} + \int\_{\mathcal{E}\_{h}^{0}} [\mathbf{u}\_{h}] : \{q\_{h}I\} \, ds\_{\prime} \\ c(p\_{h\prime}, q\_{h}) &= \int\_{\Omega} \frac{3}{2\mu} \frac{1 - 2\nu}{1 + \nu} p\_{h} q\_{h} \, d\mathbf{x}\_{\prime} \end{split}$$

The terms B.C. are related to the boundary conditions and *γ* is a proper penalization factor. This problem is equivalent to the algebraic block structured linear system

$$
\begin{bmatrix} 2\hat{\mu}\_{\text{nl}}A & B^T \\\\ B & -\frac{3}{2\mu} \frac{1-2\nu}{1+\nu} \mathcal{C} \end{bmatrix} \begin{bmatrix} \mathcal{U}\_n \\\\ P\_n \end{bmatrix} = \begin{bmatrix} F\_n \\ G\_n \end{bmatrix}
$$

This system is solved using the GMRES method and it is preconditioned with the block-triangular preconditioner

$$P = \begin{bmatrix} 2\hat{\mu}\_n \hat{A} & -B^T \\\\ 0 & \frac{3}{2\mu}M \end{bmatrix}$$

where *A* is the matrix associated to the Discontinuous Galerkin approximation of the vector-valued Laplace operator in **V***k*+<sup>1</sup> *h* and *M* is the mass matrix in *Qkh*. The implementation of *P*−<sup>1</sup> requires the computation of the inverse matrix for *A* ˆ and *M*. The mass matrix *M* can be inverted easily, because a proper choice of the base functions for the space *Qkh* and the quadrature rule to evaluate the integral leads to a diagonal mass matrix. The inverse of matrix *A* ˆ is replaced with a geometric multigrid preconditioner [45].

Even if the numerical solver is general with respect to the polynomial order *k*, the results presented in this work are obtained using quadratic polynomials for the gravity potential field (in order to evaluate the gradient of gravity acceleration) and the stable pair of spaces **V**2*h* × *Q*1*h* is used for the displacement and the pressure unknowns in the viscoelastic problem.

### **Appendix B. The Numerical Solver for the THM Model**

The surface erosion on top of the basin is taken into account as a prescribed evolution of the upper part of the sedimentary basin. To model erosion, we use the cut finite element method, briefly CutFEM, in which the boundary of the physical domain is represented on a background grid using a level set function see for example [46]. This approach requires that the computational domain must embed all possible eroded configurations. As a consequence, we immerse the physical domain that describes the basin <sup>Ω</sup>(*t*), into a larger computational domain Ω(*t*) ∪ Ω*out* as shown in Figure A1. The background or computational grid is also used to approximate the solution of the governing problem. We ideally divide the computational grid into three regions, <sup>Ω</sup>(*t*), Ω*out* and Γ, where the level set is lower, grater and equal to zero, respectively; Ω(*t*) is the physical domain where the poromechanical problem is solved; Ω*out* is a dummy zone that does not affect the solution, while Γ is the top surface of the basin where the top boundary conditions are enforced.

**Figure A1.** On the left we show the physical domain. On the right the physical domain (gray) is embedded in a larger computational domain (Ω ∪ Ω*out*). Γ marks the top surface of the basin.

**Figure A2.** A schematic of the complete algorithm for the solution of the thermal poro mechanical problem. Matrices *Bu*, *Bp*, *BT* correspond to the discretization of Equations (A1) *a*, *b*, *c* respectively.

Since Equations (5)–(7) are solved numerically, we briefly introduce the corresponding finite element discretization, which is based on the weak formulation of the problem.

Let T*h* := {*K*} denote a triangulation of ΩT = Ω(*t*) ∪ Ω*out* that does not necessarily conform to the surface Γ and let us introduce the discrete spaces

$$\begin{aligned} \mathbf{V}\_{\mathcal{h}} &:= \left\{ \mathbf{v}\_{\mathcal{h}} \in H^{1}(\Omega\_{\mathcal{T}}, \mathbb{R}^{3}) : \mathbf{v}\_{\mathcal{h}}|\_{T} \in \mathbb{P}^{1}(K, \mathbb{R}^{3}), \forall K \in \mathcal{T}\_{\mathcal{h}} \right\}, \\ Q\_{\mathcal{h}} &:= \left\{ q\_{\mathcal{h}} \in H^{1}(\Omega\_{\mathcal{T}}) : q\_{\mathcal{h}}|\_{T} \in \mathbb{P}^{1}(K), \forall K \in \mathcal{T}\_{\mathcal{h}} \right\}, \\ W\_{\mathcal{h}} &:= \left\{ w\_{\mathcal{h}} \in H^{1}(\Omega\_{\mathcal{T}}) : w\_{\mathcal{h}}|\_{T} \in \mathbb{P}^{1}(K), \forall K \in \mathcal{T}\_{\mathcal{h}} \right\}, \end{aligned}$$

where P<sup>1</sup> denotes the space of scalar piecewise linear polynomials on T*<sup>h</sup>*. We introduce the following bilinear forms:*a*(**<sup>u</sup>**, **v**) := 2 Ω(*t*) *με*(**u**) : *ε*(**v**)*d*<sup>Ω</sup> + Ω(*t*) *λ*(∇ · **u**)(∇ · **v**)*d*<sup>Ω</sup> and *<sup>c</sup>*(*p*, *q*, *K*) := Ω*in K*∇*p* · ∇*qd*Ω. Let us also introduce the operators *Dp*(*p*, *q*, **u**) = 1/*M*(*p*, *q*) + *α*(∇ · **u**, *q*) , *DT*(*<sup>T</sup>*, *w*, **v**) = *ρbcb* (*<sup>T</sup>*, *w*) + *ρf c f* (**v** · ∇*T*, *w*) and *<sup>Θ</sup> p*, *<sup>q</sup> <sup>K</sup>*Γ = − Γ *K*∇*p* · **n***qd*<sup>Γ</sup> + *γh*−<sup>1</sup> Γ *pqd*<sup>Γ</sup> − Γ *K*∇*q* · **<sup>n</sup>***pd*Γ, where (·, ·) is the standard inner product in the space *<sup>L</sup>*<sup>2</sup>(Ω(*t*)), *h* is the characteristic size of the quasi-uniform computational mesh and *γ* > 0 denotes a penalization parameter. For the imposition of the pressure and temperature boundary conditions on the internal unfitted interface Γ, we rely to the Nitsche's method following the approaches proposed in [47–50]. This technique allows to weakly enforce interface conditions at the discrete level by adding to the variational formulation of the problem appropriate penalization terms (*γ*). Finally, considering a backward-Euler time discretization scheme, the fully discretized problem at time *t<sup>n</sup>*, *n* = 1, 2, ..., *N* can be written as follows: find (**<sup>u</sup>***h*, *ph*, *Th*) ∈ **V***h*× *Qh*× *Wh*such that:

$$\begin{cases} \begin{aligned} &a(\mathbf{u}\_{h}^{n},\mathbf{v}\_{h}) - a(p\_{h\prime}^{n}\nabla\cdot\mathbf{v}\_{h}) = (\mathbf{f},\mathbf{v}\_{h}) \quad \forall \mathbf{v}\_{h} \in \mathbf{V}\_{h\prime} \\ &D\_{\mathcal{V}}(p\_{h\prime}^{n},q\_{h\prime}\mathbf{u}\_{h}^{n}) + \tau\Lambda(p\_{h\prime}^{n},q\_{h\prime}K) = D\_{\mathcal{V}}(p\_{h\!\!\!h\prime}^{n-1},q\_{h\prime}\mathbf{u}\_{h}^{n}) \quad \forall q\_{h\!\!\!h\prime} \in Q\_{h\prime} \\ &D\_{\mathcal{T}}(T\_{h\!\!\!h\prime}^{n}w\_{h\prime}\mathbf{v}\_{D\!\!h\prime}^{n}) + \tau\Lambda(T\_{h\prime}^{n},w\_{h\prime}b) = D\_{\mathcal{T}}(T\_{h\!\!\!h\prime}^{n-1},w\_{h\prime}\mathbf{0}) + (q\_{\prime\!\!r}^{n},w\_{h}) \quad \forall w\_{h} \in \mathcal{W}\_{h} .\end{aligned} \end{cases} \tag{A1}$$

where <sup>Λ</sup>(*p*, *q*, *k*) = *<sup>c</sup>*(*p*, *q*, *K*) + *<sup>Θ</sup> p*, *<sup>q</sup> <sup>K</sup>*Γ and *τ* is the computational time step.

The solution of Equations (A1) is not trivial. We decouple the solution of poromechanical problem (i.e., the first two Equations of (A1)) from the solution of the thermal problem (last Equation of (A1)). However the solution of the poromechanical problem is still challenging due to the tight coupling between deformation and flow. To solve such problem we adopt the fixed stress iterative scheme as proposed in [31,51]. This algorithm is a sequential procedure where the flow is solved first followed by the solution of the mechanical problem. In particular, in every time steps, the algorithm is iterated until the solution converges within an acceptable tolerance. In the first step the fixed stress algorithm, given (**u***n*,*<sup>k</sup> h* , *p<sup>n</sup>*,*<sup>k</sup> h* ) ∈ **V***rh* × *Qsh* we evaluate *p<sup>n</sup>*,*k*+<sup>1</sup> *h* ∈ *Qsh* where ·*k* is the index of the fixed stress iteration. In the second step, given the new pressure *p<sup>n</sup>*,*k*+<sup>1</sup> *h* ∈ *Qsh* we find the new **u***n*,*k*+<sup>1</sup> *h* ∈ **<sup>V</sup>***rh*. The iterative steps are performed until the following convergence criterion is fulfilled

$$\frac{p\_h^{n,k+1} - p\_h^{n,k}}{p\_h^{n,0}} < \eta\_{p^\*} \qquad \frac{\mathbf{u}\_h^{n,k+1} - \mathbf{u}\_h^{n,k}}{\mathbf{u}\_h^{n,0}} < \eta\_{u\prime} \tag{A2}$$

where *ηp* and *ηu* are the desired tolerances (here chose equal to <sup>10</sup>−7). A schematic of the complete algorithm for the solution of the thermal poro mechanical problem in a compact form is shown in Figure A2. For further details about the performance of this method we remand to [52].
