**4. Results and Discussion**

### *4.1. 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 [11–13]. 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<sup>3</sup> years, according to the variable thickness profile shown in Figure 2. More precisely, we assume that the glaciation phase starts 26 × 10<sup>3</sup> years ago and ends at present. We split this time window in three uniform intervals of 8.6 × 10<sup>3</sup> 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 agreemen<sup>t</sup> with the results reported in [24,25].

**Figure 2.** The evolution of the height of the ice sheet is shown on the left and the prescribed temperature field on top on the basin is reported on the right.

**Figure 3.** The simulation of isostatic adjustment with visualization of the selection of the location where the deformation tensor is extracted.

**Figure 4.** A zoom of the computational cells on the crust layer of the glacial rebound simulations. The magnitude of the vertical displacement is shown (meters). The black edges highlight the computational cell used to calculate the deformation gradient *F*.

**Figure 5.** We show the average displacement field over time in the *x*, *y* and *z* direction at the nodal points of of the selected cell. Precisely, the light regions mark the minimum and the maximum displacement sampled at the nodes of the selected computational cell.

### *4.2. Results of the Kriging Algorithm for the Geological Basin Model Setup*

An example of a simple reconstructed basin geometry is represented in Figure 6. These results are obtained on a synthetic example which considers a domain of 50 by 20 km in planar dimension. The one-dimensional simulations consider sediment deposition for a total period of time of 40 Ma. The sedimentation velocity is assumed to be fixed in time and is assigned at each *x*, *y* location between a maximum value 90 m/Ma (at the domain center) and and a minimum of 40 m/Ma (at the domain boundaries, i.e., at *x* = 0 and *x* = 50 km). Interfaces between different layers are approximated in such a domain by applying ordinary kriging starting from interface depths calculated at 20 locations randomly placed in the computational domain. Four interfaces are considered in total, including the top and bottom surfaces. 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 [6,34,35], are not explicitly included. However, the effects of chemical processes are considered in the one-dimensional model employed to approximate the interfaces [5,36] 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.

**Figure 6.** Kriging-based basin-scale reconstruction of layers interfaces.

### *4.3. 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 ky with a time step of 0.14 ky chosen to follow with enough detail glaciation evolution.


**Figure 7.** On the left we show a sketch of the physical domain, ot top of which is superposed a layer visualizing the ice sheet. The labels 1, 2, 3 mark the different materials. On the right the list of materials properties is reported.

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*,

$$\text{ls}(\mathbf{x}, z, t) = -\left(-0.285\mathbf{x} + 1.6z + 12353\right) + 400(1 + a) \,\tag{8}$$

$$a = 0 \qquad \text{if } t < -8.7 \text{Ky} \,, \tag{9}$$

$$a = \frac{t + 8.7\text{Ky}}{8.7\text{ky}}\qquad\text{if }t \ge -8.7\text{Ky}\,.\tag{10}$$

Finally, in the simulations presented in this section, we consider that the basin is exposed to a spatially dependent heat source *qr*, 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 *qr*(**x**) = *qr*(*z*) = *q*0 *f*(*z*). Following [37], 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* + *qr*,<sup>0</sup> where *φ* is the surface heat flux from the bottom surface, *S*, *V* is the basin volume and *qr*,<sup>0</sup> is the baseline radiogenic source. In this way we obtain *q*0 = 0.1 μW/m<sup>3</sup> .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* = *<sup>q</sup>*0*e*<sup>−</sup>*<sup>z</sup>*/*<sup>D</sup>* where the parameter *D* = 5 km.

**Figure 8.** The evolution of the top surface of the basin due to erosion form configuration *S*1 to *S*2 is shown during the last part of the simulation, from −8.7 ky to the present.

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 [38] that produces the 3D a labeled computational grid shown in Figure 9, right panel.

**Figure 9.** A sketch of the pipeline to build the geological model and the mesh. On the left we show the point cloud and corresponding reconstruction of the horizons. In the middle we show the lateral surfaces together with the internal horizon. On the right we report the final Computational grid of the geological model.

For the analysis and the interpretation of the results we subdivide the simulation in three different phases: phase A (formation of the ice sheet) from −26 to −17.4 ky where the ice is growing, phase B (isostatic adjustment) from −17.4 to −8.7 ky where the ice load is steady and most of the isostatic motion takes place and the phase C (erosion) from −8.7 to 0 ky where the ice on top of the basin vanishes and the erosion takes place. A schematic of this temporal subdivision is reported in Figure 10. The current configuration of the domain <sup>Ω</sup>(*t*), rotated according to the rigid motion coming from the isostatic adjustment simulation at different time steps, is shown in Figure 11. The evolution of displacement and pressure along the phases A, B and C of the simulation is shown in Figure 12. 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 12 and 13, where we show at *tB* = −13 ky 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 13, 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.

**Figure 10.** Timeline of the simulation. On the top we show the height of the ice, on the bottom temperature profile (red) and the thickness of the eroded material (green).

**Figure 11.** Time variation of the relative total thermal source (*qr*) during the evolution of the domain because of isostasy and erosion.

**Figure 12.** 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 *tA* = −21.7*ky tB* = −13*Ky* and *tC* = <sup>−</sup>4.3*Ky*, from left to right.

**Figure 13.** 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.

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 scale determined by the different permeability values.

The comparative study of the two scenarios described above is reported in Figure 14. 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 12, 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 14). Such peak is not present in the second scenario, where we neglect the isostatic movement of the basin, as we can see from Figure 14 (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 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.

**Figure 14.** Pressure field on the lateral (on top) and front (on the bottom) lines *r*1 and *r*2, respectively. Solid (*Piso*) and dashed line (*Pno*−*iso*) mark the results obtained considering and neglecting the effect of erosion and of the isostasy movements, respectively.

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 12, 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 15 that a difference of 60 deg*C* 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 agreemen<sup>t</sup> with the thermal gradients expected in literature, see for example [37,39] for thermal properties and, in particular [37] for expected temperature profiles. For a more detailed analysis of the thermal variations from the initial equilibrium state, we report, in Figure 16, difference of the temperature field at time *tA* = −21.7*ky tB* = −13*Ky* from the one at *t*0 = <sup>−</sup>26*Ky*. We notice that at *tA* and *tB*, the main temperature variations are driven by the increase of surface temperature from −10 deg*C* to 0 deg*C* 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 *tA* and *tB*, we observe that at time *tB* the temperature in the central part of the basin has slightly increased with respect to the one at *tA*, by the effect of the thermal source . The profile at *tC* differs significantly form the others, because of the erosion, active in the time interval from *tB* to *tC*. During erosion, the reference temperature of 0 deg*C* enforced at the contact interface with the ice, progressively shifts downwards, swiping a region of the basin that was previously hotter than 0 deg*C*. For this reason, the temperature in the basin after erosion significantly decreases with respect to the initial time. Finally, in Figure 11 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ω*/ <sup>Ω</sup>(*<sup>t</sup>*0) *q*(**x**)*d<sup>ω</sup>*. These data sugges<sup>t</sup> 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.

**Figure 15.** Comparison of the initial (*t* = −26 ky on the left) and final (*t* = 0 ky on the right) temperature fields. It is observed, at the depth of 4 km, temperature gradient of approximately 60◦C as qualitatively expected in [37].

**Figure 16.** The difference of the temperature field at time *tA* = −21.7 ky, *tB* = −13 ky and *tC* = −4.3 ky from that at *t*0= −26 ky is shown from left to right.

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 [25]. 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.
