Multiscale modeling of glacial loading by a 3D Thermo-Hydro-Mechanical approach including erosion and isostasy

We present a computational framework that allows investigating the Thermo-Hydro-Mechanical response of a representative part of a sedimentary basin during a glaciation cycle. We tackle the complexity of the problem, arising by the mutual interaction among several phenomena, by means of a multi-physics, multi-scale model with respect to both space and time. Our contribution addresses both the generation of the computational grid and the algorithm for the numerical solution of the problem. In particular we present a multi-scale approach accounting for the global deformation field of the lithosphere coupled with the Thermo-Hydro-Mechanical feedback of the ice load on a representative part of the domain at a finer scale. In the fine scale model we also include the erosion possibly caused by the ice melting. This methodology allows investigating the evolution of the sedimentary basin as a response to glaciation cycle at a fine scale, taking also into account the large spatial scale movement of the lithosphere due to isostasy. The numerical experiments are based on the analysis of simple scenario, and show the emergence of effects due to the multi-physics nature of the problem that are barely captured by simpler approaches.


Introduction
Reconstructing the stress and deformation history of a sedimentary basin is a challenging and important problem in the geosciences and a variety of applications [47]. The mechanical response of a sedimentary basin is the consequence of complex multi-physics processes involving mechanical, geochemical, geophysical, geological and thermal aspects [49]. The strongly coupled nature of the deformation problem may be understood in terms of the feedbacks underlying crustal dynamics. The pore fluid pressure affects stress, stress changes can lead to fracturing, and fracturing can affect pore fluid pressure [4,46].
Basin scale compaction processes involve mechanical and chemically induced transformations that take place during the accumulation of sediments [49]. In this context a number of approaches have considered the geochemical and mechanical compaction problem from a one-dimensional perspective, i.e. by considering mass, momentum and energy balances along the vertical direction, applied to fluid and solid phases [4,18,21,46]. These simplified one-dimensional approaches may be effective in interpreting qualitatively well data (e.g. [15]), however, they cannot capture inherently three dimensional processes that may arise due to the coupling of mechanical deformations and fluid mechanics in geological bodies that play an important role in the presence of glaciations [32].
Hydro-mechanical effects of continental ice sheets are widely recognized to cause movements and stresses of overridden terrains by ice load. The effect of the ice load on top of the sedimentary basin can be represented by the combination of two effects. The first is a large scale effect where we consider the action of the ice load on the entire lithosphere. The second is a fine scale analysis where we take into account the thermo-hydro-mechano-chemical (THMC) effects of the ice load into a small portion of the crust, such as a sedimentary basin [32].
In the global large scale framework, the interaction between the lithosphere and the glaciation cycle is modeled by means of a viscoelastic model. This choice is based on significant previous efforts devoted to define a proper mathematical model for the description of glacial isostatic adjustment. Initially this problem has been considered by Rayleigh [38] which studied the problem of a pre-stressed elastic compressible layer as an approximation of a "flat" planet. After Rayleigh's work other authors enriched his theory, including many other details, like the effect of viscosity or the stratified structure of the Earth. First of all, Love [28] gave a more detailed theory and defined the basic concepts which are included in more recent works. Peltier and his coauthors in a series of articles [35,36,52] gave a detailed description of a more realistic viscoelastic model of stratified Earth. The Peltier's model is essentially an extension of Love's model, where a viscoelastic rheology is used instead of an elastic one. All the mathematical details of this theory are contained in the works of Biot and, more recently, Ogden [6,33].
In this work we apply these models to describe the global deformation field of the lithosphere and to extract from it the information of the movement of a selected part of the sedimentary basin. Since the spatial scale of such region is very small compared to the global scale, we describe it as a rigid motion. More precisely, the rigid motion of the fine scale basin model is extrapolated from the lithosphere displacement fields and used at run time to move the computational grid in the simulations. In what follows we describe in more details this the workflow of this multiscale approach.
A number of numerical simulation tools have been presented to model THMC processes [31,42,44]. While considering a similar mathematical approach, our work introduces the following new features with respect to previous studies: i) it combines the global deformation of the lithosphere with the local simulation of pressure and temperature fields; ii) it is built on available geophysical and geological information on the whole sedimentary system and relying on information available at selected wells; iii) it allows considering the effects of erosion induced by glaciation, which is often neglected in previous studies.
Our work focuses on the integration of THM simulation of a single glaciation cycle with larger scale information available on basin scale compaction and lithosphere dynamics, thus the proposed THM simulation of glaciations can be cast within a multi-scale geological simulation framework. A visual sketch of the multiscale model outline is provided in Figure 1.

Isostatic glacial rebound model
From the mechanical point of view the interior of the Earth can be considered as composed of four main layers: the inner and outer core, the mantle and the lithosphere [48]. During the growth of a continental ice sheet, the lithosphere under the ice load is deformed into the mantle and the removal of the ice load during deglaciation initiates a rebound process. The uplift is well known in formerly glaciated areas, e.g. North America and Scandinavia, and in currently deglaciating areas, e.g. Alaska, Antarctica, and Greenland. Compared to water, the mantle viscosity is 10 22 − 10 25 times higher, therefore the uplifting will be slowed down and continue long time after the ice has gone. The entire process of subsidence during the glacial growth, followed by uplift during and after deglaciation, is referred to as glacial isostatic adjustment. The glacial isostatic adjustment process is dependent on the viscosity structure of the mantle, as well as the elastic thickness of the lithosphere. Observations of this process can therefore be used to gain insight into these properties of the Earth and this is important for an understanding of the dynamics of the Earth's interior.
A well established assumption for the computation of the solid Earth response to surface ice loads over glacial timescales is that the Earth can be considered as a viscoelastic body ( [37]). In particular the lithosphere can be assumed to be elastic and the solid mantle beneath behaves as a viscous fluid. A complete review of the state of the art concerning the modeling and simulation of the glacial rebound can be found in [50,51] whereas the importance of this phenomenon in the context of basin simulation is discussed in [17,26,53].

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 3 since the dynamics due to the glacial isostatic adjustment does not involve the core of the planet. Following the approach presented in [43] we assume that the viscoelastic stress tensor is given σ = σ e − q, where σ e 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 σ e is given by the sum of the deviatoric and the volumetric parts σ e = 2µ e(u) − pI , where I denotes the identity matrix. The deviatoric part is the product between the shear modulus µ and the deviatoric strain e(u) defined as The volumetric part depends on the pressure p which is defined by the equation where ν is the Poisson ratio. The internal variable q is defined by the evolution equation whereq 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: By means of the previous expressions the viscoelastic stress tensor σ is evaluated through a convolution integral defined as 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: in Ω, 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 where G is the universal gravitational constant and the boundary conditions are g = −g surface n 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 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 ρ(∇ x g)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 σ 0 (which is not important for our purposes) and a dynamic term σ d which is the solution of Problem (3) with f = ρ(∇ x g)u and the boundary conditions.

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.

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 [15] to interpret vertical distributions of system properties with a one-dimensional model, which was developed in [18] starting from classical approaches to vertical compaction modeling (e.g., [49]).
At each location such one-dimensional model provides an approximation of layers interface locations, whose characterization under uncertainty is investigated in detail in [14]. 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 [39]. 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.

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 [22], 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 T · F . Using the definition of F , we obtain that G = (RU ) T · (RU ) = (U T · R T ) · (R · U ), and using the orthogonality of the rotation matrix, namely R T = R −1 , we obtain G = U T · U . We recall that U is a diagonalizable matrix so can be represented by U = Q −1 · Λ · Q, where Q is the square matrix in which the i th column is the eigenvector q i 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 T · U = U · U , so using the spectral decomposition we obtain that G = (Q −1 ·Λ 2 ·Q) where Λ 2 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 −1 · Λ · 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.

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 [5] under the quasi-static assumption for modeling a linearly elastic fully saturated porous medium. Such approach is widely used in literature, see [7,16,30,47] for a non exhaustive list of examples. Given a domain Ω ∈ R 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 17. In this framework the momentum equation reads 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 k s and the dynamic viscosity of the fluid µ f , namely K = k s /µ f . Finally f is the gravity load of the porous material evaluated as f = (ρ s − ρ f )g, where ρ s , ρ 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) := 2µε(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. [13,16].
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 17. According to this notation, we consider the following boundary conditions, p = 0 , σ(u) · n = σ ice (t) on Γ(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 σ 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.

Thermal effects
The heat transfer is modeled through the following advection diffusion equation where T is the temperature field, ρ f and c f are the fluid density and specific heat, respectively; ρ b and c b are the bulk density and the bulk specific heat defined as , ρ s and c s being the porosity, the solid density and specific heat capacity; v D is the Darcy velocity defined as v D = −K∇p, b is the bulk thermal conductivity which is an average of the conductivity of the solid and the fluid phase and finally q r is the heat source. Concerning the initial condition, we consider an homogeneous temperature field, namely T = 0 ∀x ∈ Ω(t = t 0 ). According to the nomenclature shown in Figure 17, we consider the following boundary conditions, where n is the unit outward normal to the boundary and T ice (t) is the basin top temperature. Following a widely used approach in literature, see for example [17,30], the presence of the ice on top of the basin is taken into account by means of a variation of the top temperature of the basin (T ice (t)) during the glaciation cycle.
4 Results and discussion

Results of the glacial rebound model
The glacial rebound simulation addresses the deformation of the lithosphere of the whole Earth, based on a viscoelastic model. The forcing term of such simulation is the load of the ice sheet under a glaciation cycle.
The approach presented in this work is different from the one presented in the previously cited articles [35,36,52]. In these articles the authors describe a model for the deglaciation, the initial condition of their model is a glaciated Earth, the authors take into account the presence of this load at the initial condition defining a prestressed configuration. We rather simulate a full cycle of glaciation-deglaciation process to avoid the definition of a prestressed configuration, in our model the initial condition considered is a fully relaxed configuration, without enforced load.
In this case we consider a benchmark problem, where a circular sheet of radius R = 1111 Km (equivalent to an angular sector of 10 • ) centered at the North pole is formed and melt over a time window of 26 10 3 years, according to the variable thickness profile shown in Figure 2. More precisely, we assume that the glaciation phase starts -26 10 3 years ago and ends at present. We split this time window in three uniform intervals of 8.6 10 3 years each. In the first one we assume formation of the ice sheet up to a maximal thickness of 4 Km. In the central phase, we assume that the ice is static, while in the last term we model ice melting with a linearly decreasing profile of the ice thickness. The results of such simulation are reported in Figure 3. In Figure 4 we show a zoom of the computational cells on the crust layer of the glacial rebound simulations. The zoom is taken in correspondence of the computational cell used to calculate the deformation gradient F . More precisely, the numerical calculation of the deformation gradient involves the displacement field in the x, y and z direction at the nodal points of of the selected cell. The time history of the displacement at these time points is shown in Figure 5. Even if the presented test is just a synthetic example, the order of magnitude of the obtained vertical displacement is reasonable and in good agreement with the results reported in [17,53].  We assume that two interface collapse on each other, i.e. the second and third interface correspond for x > 35 km. This means that one of the layers is not found at locations x > 35 km, and that a layer pinch-out occurs. Note that in the simulation approach presented in Section 3 the effects of chemical compaction processes, such as quartz precipitation and smectite to illite transformation [19,21,45], are not explicitly included. However, the effects of chemical processes are considered in the one-dimensional model employed to approximate the interfaces [18,41] at selected locations, and therefore are implicitly embedded in the system geometry. The geometrical reconstruction of the sedimentary system obtained in this way can then be further enriched by mapping mineral compositions, porosity, permeability or other properties which may be available at selected location, e.g. through compositional kriging. These data are not considered in the following for simplicity.

Results of the thermo-hydro-mechanical effects of glaciation
In this section we consider the application of the solver described in Section 3. As introduced in Section 3.3, the mechanical effect of the glaciation is taken into account by means of a variable load. This load is related to the height of the ice accumulated on top of the basin and follows the curve shown in the left part of Figure 2. Concerning the thermal problem, adiabatic conditions are considered at the bottom and the lateral surfaces of the physical domain while a time-dependent temperature profile is imposed in the top boundary of the basin. This profile models the thermal effect of the ice and it is shown in the right panel of Figure 2.
The configuration of the basin simulated in this section and the physical parameters used to initialize the THM model are reported in Figure 7. More precisely, to perform the three-dimensional THM simulation, we consider a portion of the basin of 4 × 4 Km size with an average depth of of 2.8 Km located at the x, y coordinates (32.5, 36.5)× (8,12) km of the basin shown in Figure 6. The domain is split into three layers (numbered as 1,2 and 3 from bottom to top, as shown in the left part of Figure 7) and one of them is not continuous across the whole planar domain size resulting in a geological model with a pinch-out, consistent with the data in Figure 6. The material properties of all the different components of the basin are summarized in the right panel of Figure 7. We assume that the intermediate layer has a permeability that is significantly higher than the other ones and we consider the evolution of the basin during a time window of 26.000 years (26Ky) with a time step of 0.14Ky chosen to follow with enough detail glaciation evolution.
In the example of this section, the surface erosion on top of the basin is active during the last 8.7 Ky of the simulation. It is modeled as the prescribed evolution of the upper part of the sedimentary basin. More precisely, we assume that during the last part of the simulation a certain amount of material is removed from the upper part of the basin. As a consequence the top surface of the physical domain evolves from S 1 to S 2 , as shown in Figure 8. From the modeling standpoint, the motion of the top surface is described by means of a the level set function, that is defined in terms of spatial coordinates x, y and time t, ls(x, z, t) = −(−0.285x + 1.6z + 12353) + 400(1 + a) , Finally, in the simulations presented in this section, we consider that the basin is exposed to a spatially  x (k m ) y ( k m ) z (km) Figure 6: Kriging-based basin-scale reconstruction of layers interfaces dependent heat source q r , determined by the heat flux from the mantle and by the internal radiogenic thermal source. The combination of these factors is accounted as a volumetric thermal source of this form q r (x) = q r (z) = q 0 f (z). Following [11], the baseline source q 0 is determined by means of a heat balance equation that distributes over the basin volume the heat flux coming form the bottom surface (i.e. the one closer to the mantle) and the radiogenic source. Precisely, we set q 0 = φS/V + q r,0 where φ is the surface heat flux from the bottom surface, S, V is the basin volume and q r,0 is the baseline radiogenic source. In this way we obtain q 0 = 0.1µW/m 3 . To account for the exponentially decreasing radioactivity with depth, this source is evaluated as a function of the actual position of the domain according to the following empirical formula q = q 0 e −z/D where the parameter D = 5Km.
The kriging-based horizon reconstruction addressed before approximates the interface locations on a user-defined spatially uniform grid (with uniform size equal to 0.5 km in our example), so that each interface is represented by a point cloud. Such point clouds are used to generate the internal surfaces of the basin as shown in the left panel of Figure 9. Then the boundaries of the horizons are used for the creation of the lateral surfaces of the geological model as shown in the middle panel. These operations result in the definition of a water tight geological model and are performed using the platform GOCAD. The geological model is then used as the input of the mesh generator RINGMesh [34] that produces the 3D a labeled computational grid shown in Figure 9, right panel.
For the analysis and the interpretation of the results we subdivide the simulation in three different Young Modulus E 1 10 11 P a E 2 2 10 10 P a E 3 10 10 P a Rock density ρ s 2.2 10 3 Kg/m 3 Rock permeability K 1 10 −18 m 2 K 2 10 −12 m 2 K 3 10 −16 m 2 Thermal diffusivity b 1 10 −6 m 2 /s b 2 2 10 −6 m 2 /s b 3 10 −6 m 2 /s Water density ρ l 10 3 Kg/m 3 Water viscosity µ l 10 −3 P a s Dimension l 4 10 3 m Radiogenic source q 10 −7 W/m 3 K  Figure 10. The current configuration of the domain Ω(t), rotated according to the rigid motion coming from the isostatic adjustment simulation at different time steps, is shown in Figure 16. The evolution of displacement and pressure along the phases A, B and C of the simulation is shown in Figure 11. In the first phase (A), the action of the ice load generates the mechanical compaction of the basin. According to the Biot model the compaction leads to an increase in the pore pressure in the different layers of the basin. This effect is more evident in the top layers, 2 and 3, while layer 1 is almost unperturbed, because it is the most impermeable and it is subject to zero displacement conditions on the bottom surface. Moreover, we notice that the pressure field in the pinch out layer (layer 2) is almost uniform and equal to the value at the interface with the upper layer (layer 3). This effect is due to the fact that layer 2 is the most permeable. According to the Darcy law, this is the region where the smallest pressure gradients are observed.
In the second phase (B) we can observe the effect of the isostatic adjustment. The interaction between the isostatic motion end the poromechanical problem is illustrated in Figures 11 and 12, where we show at t B = −13Ky the pressure field in an internal slice of the domain, together with the displacement direction inside layer 3. We observe that the imposed rotation generates a tangential component of the load applied on the top of layer 3, namely σ ice n · t = 0. This effect, combined with the boundary conditions enforcing zero normal displacement on the lateral surfaces of the basin, generates a mechanical compression effect along the tangential direction of the top and bottom surface planes of layer 3. Because of the poromechanic coupling, this tangent stress induces the pressure peak observed at the left corner of the basin. Furthermore, in the right part of Figure 12, an internal slice along the xz plane is shown. From this view, the transition from layer 1 to layer 2 (at the pinch out) is visible. In this visualization, the effect of the permeability jump can be appreciated. More precisely, the pressure increases only along the interface between layers 1 and 3, while the pressure gradients are significantly smeared out along the contact line with layer 2 as a consequence of the high permeability of it. Finally, we notice that these effects are difficult to appreciate in the bottom layer, where the high young modulus limits the displacement field and the low permeability almost blocks the diffusion of the pressure peak occurring in the top layer.
To appreciate the impact of taking into account the isostatic response, we compare the pressure field in the top layer with the one obtained through a different numerical experiment in which we neglect the isostatic adjustment and the erosion. In the second scenario the basin model is completely static. In such case, because the load on top of the basin is constant, the general temporal trend of the pressure is to progressively reach a steady state in all the layers of the basin according to the characteristic temporal Rigid motion Displacement Pressure t A t B t C Figure 11: Evolution of the isostatic adjustment, the displacement and the pressure fields along the different phases of the simulation. In particular, we show on the top we show an the imposed rigid motion. The displacement (middle row) and the pressure fields (the variation from the hydrostatic pressure profile is shown, bottom row) are reported at times t A = −21.7ky t B = −13Ky and t C = −4.3Ky, from left to right.
scale determined by the different permeability values. The comparative study of the two scenarios described above is reported in Figure 13. In this Figure  we compare the pressure field along two transversal lines r 1 and r 2 , obtained switching on (first scenario, solid line) and off (second scenario, dashed line) the isostatic motion of the basin. As previously observed in the central panel of Figure 11, we notice that in the first scenario the pressure reaches a peak in the left corner of the layer 3 (as also illustrated in Figure 13). Such peak is not present in the second scenario, where we neglect the isostatic movement of the basin, as we can see from Figure 13 (dotted line). From that comparison we notice that limited spatial variations of the pressure are generate without taking into account the rigid motion of the basin, Conversely, the introduction of the isostatic rebound significantly Figure 12: Front (on the left) and lateral (on the right) view of the basin. The color marks the pressure field, the arrows shows the direction of the solid displacement (in the interior and on the boundary of the domain, that is u · n = 0) and the black lines show the direction of gravity. The configuration of the basin layers is also shown, to facilitate the interpretation of the results.
increases the pressure variability, generating local peaks that depend on the morphology of the basin. Then, we conclude that the motion due to the isostatic adjustment simulation is the primary reason of the overpressure generation.
We now focus on the final phase (namely phase C) when the ice sheet is disappearing from the top of the basin and the erosion takes place. In this phase the mechanical load due to the ice weight goes to zero. The reduction of mechanical compaction is also augmented by the erosion of the top part of the basin. The combination of these effects leads to a general decrease of the pore pressure in all layers, as shown in the right panel of Figure 11, right column. More precisely, the lowest values of pressure (in terms of variations from the hydrostatic load) are located in the softest layer.
Concerning the temperature field, the model is initiated with a temperature field at thermal equilibrium and it evolves according to the (time dependent) heat equation under the action of a space dependent heat source. To analyze and validate the simulation of the thermal field we address two indicators: the overall thermal gradient from the top to the bottom of the basin and the temporal variations of the temperature from the steady state. Concerning the vertical temperature gradient, we observe from Figure 14 that a difference of 60 degC along the vertical axis is established as a consequence of thermal balance between the thermal source and the imposed temperature at the top of the basin. This results is in good agreement with the thermal gradients expected in literature, see for example [11,20] for thermal properties and, in particular [11] for expected temperature profiles. For a more detailed analysis of the thermal variations from the initial equilibrium state, we report, in Figure 15, difference of the temperature field at time t A = −21.7ky t B = −13Ky from the one at t 0 = −26Ky. We notice that at t A and t B , the main temperature variations are driven by the increase of surface temperature from −10 degC to 0 degC at the top of the basin, due to the protective effect of the ice cap, that is progressively forming in this phase. Comparing more in detail the profiles at t A and t B , we observe that at time t B the temperature in the central part of the basin has slightly increased with respect to the one at t A , by the effect of the thermal source . The profile at t C differs significantly form the others, because of the erosion, active in the time interval from t B to t C . During erosion, the reference temperature of 0 degC enforced at the contact interface with the ice, progressively shifts downwards, swiping a region of the basin that was previously hotter than 0 degC. For this reason, the temperature in the basin after erosion significantly decreases with respect to the initial time. Finally, in Figure 16 we report the variation of the total thermal source term, relative to the initial state, that is the spatial integral during the evolution of the domain Ω(t) of the heat source Ω(t) q(x)dω/ Ω(t 0 ) q(x)dω. These data suggest that variation of the basin volume due to erosion decreases the thermal source more significantly than the progressive increase of basin depth due to isostasy.
The analysis of the thermal field confirms that, in the short time window of one cycle of glaciation, the surface temperature, possibly modulated by erosion, is the main diving factor that determines temperature variations from the equilibrium state. This finding is in accordance with the detailed analysis of the effects of glaciations on sedimentary basins provided in [17]. It should also be observed that using a time dependent heat equation correctly models the natural inertia of the basin to change its equilibrium state. As a result, it is expected that the repetition of many glaciation cycles is required to significantly cool down the entire basin.

Conclusions
We have developed a fully coupled multiphysics and multiscale description of the evolution of a sedimentary basin under the effect of a glaciation cycle. To the best of our knowledge, thermo-hydro-mechanical effects are combined with isostatic adjustment and erosion within a fully time dependent three-dimensional simulation for the first time. Although the geological model that has been considered does not represent yet the complexity of a real sedimentary basin, we have described a pipeline of steps that could handle the most complex cases as well, thanks to a multiscale approach that decouples phenomena occurring at very different space and time scales. For example, our methodology establishes a quantitative framework to transfer information from the definition of a large scale geological architecture to local fluid displacement and deformation dynamics. Preliminary numerical results suggest that the combination of all these phenomena reveals the emergence of effects that were not expected or predictable using simpler approaches. In the considered synthetic test case we have quantified the effects of large scale isostatic displacement on local pressure and displacement field, in the presence of layer a pinch-out. Moreover, our results demonstrate Figure 14: Comparison of the initial (t = −26Ky on the left) and final (t = 0Ky on the right) temperature fields. It is observed, at the depth of 4 Km, temperature gradient of approximately 60 degC as qualitatively expected in [11]. 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 [1,3]. 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 h be a face shared by two elements K + and K − , define the unit normal vector n + and n − on the e pointing exterior to K + and K − , respectively. With φ ± = φ| K ± we set the average operators as together with the jump operators:  Equation (4) is discretized using the standard SIPG method [2]: 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 [25]. The Equations (3) are discretized using a standard second order accurate, one-step and unconditionally stable scheme [43] 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: 2τ e(u n − u n−1 ) − p n I + e − ∆tn τ h n−1 + ρ(∇ x g)u n = 0 in Ω, 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 − ∆tn 2τ and all the known terms in the first equation are written as a unique term F n : F n = ∇ · 2μ n e(u n−1 ) − e − ∆tn τ h n−1 .
These equations can be rewritten in a full discrete form using the standard SIPG-method for the linear elasticity [40]: Where the bilinear forms are defined by 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 This system is solved using the GMRES method and it is preconditioned with the block-triangular preconditioner WhereÂ is the matrix associated to the Discontinuous Galerkin approximation of the vector-valued Laplace operator in V k+1 h and M is the mass matrix in Q k h . The implementation of P −1 requires the computation of the inverse matrix forÂ and M . The mass matrix M can be inverted easily, because a proper choice of the base functions for the space Q k h and the quadrature rule to evaluate the integral leads to a diagonal mass matrix. The inverse of matrixÂ is replaced with a geometric multigrid preconditioner [10].
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.

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 [9]. 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 Ω(t), into a larger computational domain Ω(t) ∪ Ω out as shown in Figure 17. 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, Ω(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.
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   where P 1 denotes the space of scalar piecewise linear polynomials on T h . We introduce the following bilinear forms:a(u, v) := 2 Ω(t) µε(u) : ε(v)dΩ + Ω(t) λ(∇ · u)(∇ · v)dΩ and c(p, q, K) := Ω in K∇p · ∇qdΩ. Let us also introduce the operators D p (p, q, u) = 1/M (p, q) + α(∇ · u, q) , D T (T, w, v) = ρ b c b (T, w) + ρ f c f (v · ∇T, w) and Θ p, q K Γ = − Γ K∇p · nqdΓ + γh −1 Γ pqdΓ − Γ K∇q · npdΓ, where (·, ·) is the standard inner product in the space L 2 (Ω(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 [8,23,24,27]. 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 n , n = 1, 2, ..., N can be written as follows: find (u h , p h , T h ) ∈ V h × Q h × W h such that: where Λ(p, q, k) = c(p, q, K) + Θ p, q K Γ and τ is the computational time step. The solution of equations (11) is not trivial. We decouple the solution of poromechanical problem (i.e. the first two equations of (11)) from the solution of the thermal problem (last equation of (11)). 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 [7,29]. 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,k h , p n,k h ) ∈ V r h × Q s h we evaluate p n,k+1 h ∈ Q s h where · k is the index of the fixed stress iteration. In the second step, given the new pressure where η p and η u are the desired tolerances (here chose equal to 10 −7 ). A schematic of the complete algorithm for the solution of the thermal poro mechanical problem in a compact form is shown in Figure  18. For further details about the performance of this method we remand to [12].