**3. The Thermo-Hydro-Mechanical Model Including Erosion**

In this section we focus the attention on the mathematical framework describing the mechanical and thermal evolution of the basin. First, we show how the basin model is built and how we basin tilting is extracted from the large scale isostasy model. Second, we use mathematical models to describe the thermal and mechanical evolution of the basin under the effect of a glaciation cycle. The combination of these models consist in the multiscale modeling of glacial loading illustrated in Figure 1.

In the next sections the THMC model is introduced; the details of the numerical method can be found in Appendix B.

### *3.1. The Geological Model of a Sedimentary Basin*

We assume that information about sedimentary units thicknesses, porosity and/or mineral composition are available at selected locations across a sedimentary system, typically from well logs, i.e., along the vertical direction (as sketched in Figure 1). In this framework we employ the stochastic inverse modelling procedure implemented in [7] to interpret vertical distributions of system properties with a one-dimensional model, which was developed in [5] starting from classical approaches to vertical compaction modeling (e.g., [2]). At each location such one-dimensional model provides an approximation of layers interface locations, whose characterization under uncertainty is investigated in detail in [27]. These interface locations can then be approximated in the whole domain starting from these data. This task is here performed using an interpolation based on ordinary kriging [28]. Our procedure assumes a smooth spatial variation of the sedimentary unit thicknesses. While this hypothesis can be restrictive in practical cases (e.g., in the presence of fault zones), our procedure is still able to handle geological settings of interest such as the occurrence of pinch out layers.

### *3.2. Extraction of the Basin Tilting from the Isostasy Model*

The isostatic movement of the basin is taken into account as a rigid motion. In particular, we locate the computational cell that embed the fine scale geometrical model from the isostatic adjustment simulation, as shown in Figure 3. From the displacement field of the selected cell we evaluate the deformation gradient and the vertical rigid motion. We then isolate the rotational part of the deformation gradient that is uniquely defined by the polar decomposition [29], *F* = *R* · *U*, where *F* is the deformation gradient, *U* is the symmetric stretch tensor and *R*, such that det(*R*) = 1, is the rotation matrix that we want to determine. In particular, from the deformation gradient we evaluate the Right Cauchy-Green Deformation Tensor namely matrix *G* := *F<sup>T</sup>* · *F*. Using the definition of *F*, we obtain that *G* = (*RU*)*<sup>T</sup>* · (*RU*)=( *U<sup>T</sup>* · *R<sup>T</sup>*) · (*R* · *<sup>U</sup>*), and using the orthogonality of the rotation matrix, namely *R<sup>T</sup>* = *<sup>R</sup>*−1, we obtain *G* = *U<sup>T</sup>* · *U*. We recall that *U* is a diagonalizable matrix so can be represented by *U* = *Q*−<sup>1</sup> · Λ · *Q*, where *Q* is the square matrix in which the *i*th column is the eigenvector *qi* of *U* and Λ is the diagonal matrix composed by the corresponding eigenvalues, namely Λ*ii* = *ζi*. We recall that *U* is a symmetric matrix so *U<sup>T</sup>* · *U* = *U* · *U*, so using the spectral decomposition we obtain that *G* = ( *Q*−<sup>1</sup> · Λ<sup>2</sup> · *Q*) where Λ<sup>2</sup> is a diagonal matrix defined by Λ2*ii* = *ζ*2 *i* . From the previous considerations, it follows that by computing the eigenvalues and the eigenvectors of *G*, which is a known matrix, we calculate the matrix *Q*, Λ2, and, as a consequence, Λ. Using *U* = *Q*−<sup>1</sup> · Λ · *Q* we compute the matrix *U* so that we finally retrieve the rotation matrix *R* using *R* = *F* · *U*−1.

The collection of the rotation matrices together with the axial displacement evaluated at every time step of the isostatic adjustment simulation, defines the rigid motion that we consider in fine scale model. In particular we use a Lagrangian and an Eulearian approach for the mechanical and the thermal problems, respectively. More precisely, in the mechanical problem the imposed rigid motion does not causes any additional stress so that the calculation of stress and the strain can be performed in the reference (static) configuration, while the rotation matrix is used to change the orientation of the gravity vector. In this framework, we let pressure, flow and displacement evolve under the action of the weight of the basin and of the ice. Concerning the thermal evolution of the basin we account of the variation of the heat flux with the isostatic movement of the sedimentary basin by means of the Eulerian approach, that is we actually move the domain in the computational model. In this way, both the vertical displacement and the rotation matrix are taken into account in the evaluation of the thermal source, as it will be discussed later.

### *3.3. A Poromechanical Approach to Coupled Hydro-Mechanical Effects*

For the mechanical evolution of the basin we rely on the theory of poroelasticity introduced by Biot in [30] under the quasi-static assumption for modeling a linearly elastic fully saturated porous medium. Such approach is widely used in literature, see [1,16,31,32] for a non exhaustive list of examples. Given a domain Ω ∈ <sup>R</sup>*d*, we consider for simplicity an isotropic material (named the *skeleton*) filled with an isothermal single-phase fluid. A sketch of a typical domain together with the frame of reference is shown in Figure A1. In this framework the momentum equation reads

$$-\nabla \cdot \sigma(\mathbf{u}) + \mathfrak{a} \nabla p = \mathbf{f} \tag{5}$$

$$
\partial\_t \left( \frac{p}{M} + a \nabla \cdot \mathbf{u} \right) - \nabla \cdot K \nabla p = 0 \tag{6}
$$

where (with little abuse of notation with respect to the isostatic adjustment model) here **u** denotes the solid matrix displacement vector and *p* is the variation of pore pressure from the hydrostatic load. We notice that *∂t* denotes the standard partial derivative with respect to time in the Eulerian framework. The parameters *α*, *M* and *K* are the Biot coefficient, the Biot modulus and the hydraulic conductivity, respectively. We recall that the hydraulic conductivity is defined as the ratio between the permeability *ks* and the dynamic viscosity of the fluid *μf* , namely *K* = *ks*/*μf* . Finally **f** is the gravity load of the porous material evaluated as **f** = (*ρs* − *<sup>ρ</sup>f*)**g**, where *ρ<sup>s</sup>*, *ρf* and **g** are the fluid density, the solid density and the gravity vector, respectively. We remark that in the mechanical model the isostasy movement coming from the isostatic adjustment simulation is taken into account by means of a Lagrangian approach. As a consequence, the gravity vector **g** varies during the simulation according to a prescribed profile. Such profile is defined by the angles evaluated from the polar decomposition of the deformation gradient of the computational cell, of the large scale simulation, that contains the portion of the sedimentary basin we consider. To complete the definition of the problem we recall the linear elasticity behavior for the skeleton. This implies that the stress tensor *σ*, appearing in (5), is defined by *σ*(**u**) := <sup>2</sup>*με*(**u**) + *λ*∇ · **u** , where *μ* and *λ* are the Lamé coefficients and *ε*(**u**) is the symmetric gradient of the skeleton displacement. For further details on poromechanical modeling, the interested reader is referred to e.g., [32,33].

For a well-posed problem we must complement the previous governing equations with appropriate boundary and initial conditions. Concerning the initial condition, the following constraints **u** = 0 and *p* = 0 are considered at the initial time *t* = *t*0. Let us label with Γ the top surface of the basin, while *∂*Ω*b* and *∂*Ω*l* are the bottom and the lateral boundary of the domain Ω, as show in Figure A1. According to this notation, we consider the following boundary conditions, *p* = 0 , *σ*(**u**) · **n** = *<sup>σ</sup>ice*(*t*) on <sup>Γ</sup>(*t*), **u** = 0 , ∇*p* · **n** = 0 on *∂*Ω*b*(*t*), **u** · **n** = 0 , ∇*p* · **n** = 0 on *∂*Ω*l*(*t*), where **n** is the unit outward normal to the boundary and *<sup>σ</sup>ice*(*t*) is the load resulting from the ice sheet on top of the basin. We further assume that the load relative to the solid component of the ice sheet transfers only to the solid matrix of the basin, while the fluid phase at the top surface is subject to the hydrostatic load.
