A Using a Stokes flow code to simulate other types of rheologies

Get Complete Project Material File(s) Now! »

Geomechanical restoration of sedimentary basins

In sedimentary basins, structural models consist mainly of sedimental rock layers which have been folded and faulted. Other deformation mechanisms also include erosion and deposition of sediments, and compaction of the layers. In these setups, geomechan-ical restoration was introduced as a mechanical simulation in which the rock layers are modeled with elastic properties and the faults by free-slip surfaces (e.g., Maerten and Maerten 2001; Muron 2005; Moretti et al. 2006; Durand-Riard et al. 2013b). The restoration itself is performed with specific boundary conditions, such as the flattening of the top surface or the removal of fault throw by binding the hanging-wall and footwall cutoﬀ lines. While these methods introduce more physical strains than those of purely geometrical schemes, they are still limited in several ways. Firstly, the deformation drive for the simulations is still based on geometric unfolding and unfaulting. These simplifying conditions are unphysical and do not necessarily respect the paleo-stress state. Secondly, rock behavior is only modeled with elastic properties. While these are a good assessment of rock deformation at timescales up to a few thousand years, they neglect the important viscous and plastic behavior that can appear at geologic timescales. Moreover, they force the use of overly simplifying hypotheses, such as con-sidering the interface with a salt layer as a free surface (Chauvin et al. 2018). Lastly, current geomechanical restoration methods rely on a surface or volume mesh for the computation of the model strain and deformation. These meshes are used because they can precisely account for the diﬀerent scales and shapes of geological objects. Their construction, however, can be complicated, particularly in 3D. Moreover, the precision on computations using them relies on their validity. In restoration, the model under-goes large deformations, distorting the mesh and invalidating it. Consecutively, the model needs to be edited and remeshed, which t akes time and limits the application of geomechanical restoration as a validation tool.

Geometric and kinematic restoration

Since the beginning of the last century, structural restoration methods have been mostly based on geometric and kinematic assumptions (e.g., Chamberlin 1910; Dahlstrom 1969; Gratier 1988; Rouby 1994; Groshong 2006; Lovely et al. 2018; Fossen 2016).
Their first implementation was the balancing restoration, which relies mostly on the conservation of layer bed area and thickness to restore cross-sections (e.g., Chamber-lin 1910; Dahlstrom 1969; Groshong 2006) (Fig. 1.1). Its major drawback, as every other 2D method, is to consider that deformation can only happen in the plane of the cross-section. This is possible only if the out-of-plane deformation can be consid-ered negligible, which precludes applications in settings such as strike-slip systems for example.
Map restoration has been developed to study deformations that are mainly horizon-tal (e.g., Cobbold and Percevault 1983; Gratier et al. 1991; Rouby 1994; Ramón et al. 2016). It considers structural models as a stack of layers modeled by surfaces that are folded and faulted, and aims at flattening these surfaces. This method is based on ge-ometrical conservation and recovery of fault oﬀsets, and restoration of horizon surfaces underlying the top layer is performed using kinematic laws such as flexural slip or sim-ple shear. The resulting computed displacement can be used to analyze the strain and stress along the horizon surface. Although all the horizons are defined in a 3D space, this restoration process considers only surfaces and not the entire volume of the model; it is therefore qualified as 2.5D.
Finally, fully 3D geometrical methods have also been developed (Massot 2002; Muron 2005; Medwedeﬀ et al. 2016), allowing the tracking of internal volumetric deformation. They are all based on the minimisation of deformation and on volume conservation. Among them, Massot (2002) relies on unfolding and unfaulting the uppermost hori-zon and then interpolating the deformation in the rest of the model on a regular grid representing the geological volume. This interpolation is made using Discrete Smooth Interpolation (DSI) (Mallet 1989, 1992, 1997), and was later extended to tetrahedral meshes. Another method of volumic restoration by Medwedeﬀ et al. (2016) and Lovely et al. (2018) is based on depositional coordinates (GeoChron model) computed using two types of deformations: minimal deformation and flexural slip (Mallet 2004; Moyen 2005). As 3D methods only work in specific setups, however, most applications of geometric and kinematic restoration remain 2D (cross-section and map restoration).
The main issue with all these methods is that they considerably simplify rock de-formation mechanisms. Indeed, they rely on a global simple shear (or flexural slip), assumption which generally ignores mechanical layering eﬀects. Moreover, the principle of masse conservation is reduced to the conservation of the surface or volume of the layers. In this light, numerous authors have stated the necessity of incorporating me-chanics into the restoration of geological models (e.g., Fletcher and Pollard 1999; Muron 2005; Maerten and Maerten 2006; Guzofski et al. 2009; Al-Fahmi et al. 2016).

Geomechanical restoration

Mechanics-based restoration has been developed since the 2000s as a geomechanical simulation with specific boundary values (Maerten and Maerten 2001; De Santi et al. 2002; Muron 2005; Moretti et al. 2006; Maerten and Maerten 2006; Guzofski et al. 2009; Durand-Riard 2010; Durand-Riard et al. 2013b; Tang et al. 2016). In this approach, internal deformation is not known a priori, and the resulting strain is computed from the mechanical behavior of rocks and the applied boundary conditions. Tracking the strain and stress through time gives valuable information on the deformation mechanisms, the potential fracture areas and the possible paleo-geometry. The model is parameterised by elastic properties to mimic the response of rocks to mechanical stresses. The restoration displacement is computed using fundamental laws of motion such as mass conservation and linear momentum conservation, instead of geometric or kinematic rules. In practice, it is done by solving the equation of motion r = fext (1.1).
where r is the del operator, fext is the sum of external forces, and is the Cauchy stress tensor. Under the assumption that rocks behave as elastic materials, is defined by Hooke’s law for elastic motion: = C : ; (1.2).
with C the elastic fourth-order tensor and the (small-)strain tensor associated to the displacement u by = 1 r(u) + r(u)T : (1.3).
The restoration itself is performed by applying specific boundary conditions to con-strain the model (Fig. 1.2). These conditions, usually imposed on the displacement, rely on the following asumptions: the uppermost horizon was flat and horizontal at its deposition time, and it was not faulted. Other conditions can be introduced as comple-mentary geological knowledge, such as direction and scale of deformation, or amount of shortening (Chauvin et al. 2018). Compared to previous methods, it allows the con-sideration of more physical deformation processes for the sedimental layers, leading to more physical strains inside the restored model. It can also be used to study the defor-mation history of the model where it is least known, particularly in regions away from the boundaries where the displacement is imposed.
Although these methods oﬀer significant advances in the structural restoration of geological models, they still present many limitations. First, the true mechanical de-formation of rocks is highly irreversible and the boundary conditions set to unfold and unfault the medium are often unphysical (Lovely et al. 2012). Moreover, these conditions are shortcuts to the paleo-stress state, hence can be easily questioned (Durand-Riard et al. 2010; Lovely et al. 2012; Durand-Riard et al. 2013a). For example, Durand-Riard et al. (2013b) claims that classical boundary conditions are unable to account for strike-slip and highly oblique-slip faults, and oﬀers two alternatives for a better simulation of faults. The first is to constrain fault motion by using known fault-slip data as pierc-ing points. This, however, implies having suﬃcient data, which is rarely the case, and can result in unrealistic strain around the constrained points on the fault. The second method is similar to that of Lovely et al. (2012): it constrains fault displacement using the opposite of the forward fault displacement. This was proved to provide an accurate recovery of faults and may yield a strain field more consistent with the forward strain patterns, which cannot be achieved with classical boundary conditions. Nevertheless, geomechanical restoration only considers elastic rock properties so far, neglecting other possible behaviours, such as visco-elasticity or plasticity (Gerbault et al. 1998). Trans-verse isotropic behaviour also aﬀects results (Durand-Riard et al. 2013a). Moreover, the elastic behavior can only be applied to a model where the internal deformations remain small, which is usually not the case in restoration. These physical issues raise the question of the capability of geomechanical restoration to properly recover paleo-deformation. There is, therefore, no clear guidelines on which method to choose between geometric and kinematic restoration and geomechanical restoration, despite the benefits of the second one (Maerten and Maerten 2006; Guzofski et al. 2009). Moreover, while it is called geomechanical restoration, it is extensively controlled by geometric consid-erations: flattening of the top layer and a geometric unfaulting based on stitching the horizon cutoﬀ lines across each fault.
Another issue is the need for a valid volumetric mesh of the structural model, in-cluding a boundary representation of the geological domain with the horizons and faults as boundaries (e.g., Muron 2005; Durand-Riard et al. 2010). Such a mesh is diﬃcult to generate, as shown by Pellerin et al. (2014) and Anquez et al. (2016). Since restoration deals with large deformations, this mesh evolves and has to be edited a lot or even rebuilt (Fig. 1.3), limiting the applicability of the geomechanical restoration to be used as a validation tool.

The Lagrangian and Eulerian frames

There are two main frames for describing a physical phenomenon: the Eulerian frame of reference, also known as the spatial description, and the Lagrangian frame of reference, also known as the material description (e.g., Cornet 2015).
Let us consider a body B and a particle p inside it. We consider their respective position B0 and x0 at time t = 0, and B and x at time t. The Lagrangian representation of the body B would be a function that gives the position x of any particle at time t as a function of its position at time t = 0 : x = (x0; t). This representation is generally used in solid mechanics, where the paths of the particles are of interest. In the Eulerian description, the movement of the particles is evaluated at several fixed positions xi. As time passes, diﬀerent particles occupy points xi, and the objective is to describe the motion, or the temperature, or any feature of the particles at points xi. The initial position of the particle that occupies point x at time t is then x0 = 1(x; t) (cf Fig. 1.7). This description is commonly used in fluid dynamics.
The Eulerian methodology is based on a fixed mesh on the whole domain in study, so it can implement accurate and eﬃcient computation of displacements in a geological model (e.g., Deubelbeiss and Kaus 2008; Braun et al. 2008). However, Eulerian methods leave no trace of the movement of particles, so discontinuities such as faults, horizons or other moving surfaces cannot be tracked through time, which is critical when dealing with large deformations. The material properties are also not tracked directly, so their value at the nodes of the mesh can become unaccurate and the computation of the velocity from the Stokes equations loses accuracy after a few time steps.
The Lagrangian approach, on the other hand, incorporates the interfaces, material properties and geometric description by tracking the reference configuration at each step (e.g., Podladchikov et al. 1993; Mello and Henderson 1997; Ismail-Zadeh et al. 2001; Massimi et al. 2006). In this approach, a trace of the path followed by the points of the model can be kept. However, repeated application of mesh movements when the model evolves and particles are displaced may lead to bad quality or even invalid mesh. Since the error in any computation performed in a mesh is linked to its quality, in order to keep the error as low as possible, one then needs remeshing at some time steps, which requires significant eﬀorts. In large deformations, remeshing may even be needed at each time step. To relax such an eﬀort, some methods have been developed to mitigate the loss of numerical precision due to mesh distortion in large deformation simulations (Braun and Sambridge 1994, 1995), but they are limited to 2D simulations where the mesh is a Delaunay triangulation.
When dealing with large deformations of geological models, neither of these two representations is well adapted, since one needs both the Lagrangian ability to keep track of the material properties and interfaces, and the Eulerian precision and absence of remeshing in case of large displacement. This is why methods mixing the two approaches have been developed when dealing with creeping flow models involving large strains.

READ Simulation Environment

The Arbitrary Lagrangian Eulerian representation

A notable mix of the two previous frames is the Arbitrary Lagrangian Eulerian (ALE) formulation (Hughes et al. 1981; Fullsack 1995). This method has been developed to overcome the limitations of the Eulerian and Lagrangian formulations by inheriting features from both methodologies. It has often been used in geomechanical simulations involving creeping flow phenomena such as mantle flow simulations. It has had various formulations and implementations, in 2D [Willett et al. (1993); Fullsack (1995); Poliakov et al. (1996); Gemmer et al. (2004); Massimi et al. (2007); Fillon et al. (2013)] and more recently in 3D [Braun (2003); Braun et al. (2008); Braun and Yamato (2010); Ismail-Zadeh et al. (2004c); Thieulot et al. (2014)].
Most of these methods rely on keeping track of the material model in a Lagrangian way, while computing the displacement or the velocity on a Eulerian mesh. Contrarily to the Eulerian frame, however, in the ALE formulation, the Eulerian mesh has an arbitrary velocity and can be deformed. This is particularly interesting to follow specific surfaces inside the computational domain and better account for the material model at interfaces with high contrasts of properties. In particular, allowing for moving interfaces enables the use of free surfaces, moving boundaries or internal interfaces. Since the computational Eulerian mesh is independent of the Lagrangian material, it can be done while keeping deformation inside the computational domain to a minimum. Erosion and deposition of mass is also simplified. Indeed, the Eulerian mesh can easily be deformed on its outer shape in the ALE formulation (e.g., Braun 2003; Thieulot 2011), whereas in the Lagrangian approach it would imply removing or creating mesh elements at each time step.
Depending on the mechanical simulation, the implementation varies as the quan-tities that need to be tracked are not the same. Lagrangian nodes can be used to track interfaces through time, by advecting them at each time step from the velocity computed on the Eulerian mesh. These Lagrangian interfaces can then be used for the interpolation of material properties on the Eulerian mesh. The Eulerian mesh can also be modified to better fit some of them, for example by deforming the Eulerian mesh so its nodes fit the tracked interfaces (e.g., Fullsack 1995; Gemmer et al. 2004), or by using adaptive mesh refinement (e.g., Braun et al. 2008; Bangerth et al. 2019) (Fig. 1.8). In some other implementations, no internal tracking is done and the Lagrangian part of the code is there only to deform the Eulerian mesh boundaries, to account for a free surface at the top and possibly a movement condition on the sides of the model (Fig. 1.9) (e.g., Gemmer et al. 2004; Thieulot et al. 2014). This is particularly useful in mechanical simulations of the mantle, where the viscosity and density of the materials don’t need to be tracked as they depend mainly on the temperature field, which can be computed at each time step by solving the heat equation on the Eulerian mesh.

Material discretization of a geological model

The geomechanical simulation of a specific domain requires to choose an appropriate kinematic description to follow the displacement inside the geological layers. Contin-uum mechanics distinguish two main frames: the Eulerian frame of reference, also known as the spatial description, and the Lagrangian frame of reference, also known as the material description (Cornet 2015). Both methods have their advantages and disadvantages, as shown in Section 1.4.1.1, but neither of them is specifically adapted in the case of large displacements over time, such as those studied here. In order to overcome the limitations of the two approaches, the ALE formulation (Fullsack 1995; Donea et al. 2004), which inherits features from both methodologies, was developed (see Section 1.4.1.2). Most of its implementations rely on keeping track of the material properties in a Lagrangian way, while computing the displacement on a grid that can only deform vertically to account for the deformation of a possible free surface. It is particularly useful in geodynamics, in models where the vertical deformation is small compared to the horizontal deformation, and in the case of highly viscous fluids in the mantle, for which the density and viscosity depend mostly on the temperature and depth. In FAIStokes, the computation of the velocity is done by solving the Stokes equations on a quadrilateral grid. This grid has an ALE part as its nodes can move vertically to follow the movements of the top surface of the model.
During mechanical simulations, the material properties inside the model are tracked using particles; each of these particles discretizes the properties of a small part of the model around it. At the beginning of the simulations, FAIStokes either creates a par-ticle swarm from a distribution of the material parameters or loads it from a file. In the first case, a regularly distributed particle swarm is generated, with a density of particles depending on the size of the smallest element of the computation grid (the two parameters can be tuned in the parameter file, as shown in Appendix 2.A). The given material distribution is then used to associate the material properties to the particles depending on their position. In the second case, one can load a model obtained through previous simulations, or create one. In order to create a model easily, I supervised and worked with two MSc students in a lab internship on the development of a GUI (coded in python and PyQt) for the creation of particle swarms for FAIStokes (Fig. 2.2).
Figure 2.2: Graphical User Interface for the creation of particle swarms to use as input for FAIStokes. In the present example, an image has been loaded, its transparency has been changed, and the user has started digitizing the interfaces and the faults of the cross-section. Since the particle swarm does not directly track the interfaces, it has to be dense enough to recover accurately the material properties of the model; some parts of the model can be manually densified to keep the appropriate accuracy.
At each time step, the material properties are interpolated from the particle swarm to the grid in order to build the FE matrix and its preconditioner. For each element, the density is interpolated on the quadrature points using an arithmetic mean of the densities of the particles around the quadrature points (closer than a fourth of the diagonal of the smallest elements of the domain). The viscosity is recovered for each element using a harmonic mean of the viscosities of the particles inside the element. This reduces the eﬀect of very high viscosity diﬀerences (possibly of several orders of magnitude) on the solver (compared to interpolating it on the quadrature points using the particles nearby), and is more computationally eﬃcient despite the higher grid refinement needed (Deubelbeiss and Kaus 2008; Thielmann et al. 2014; Heister et al. 2017). In the simulations we did with FAIStokes, we have checked that this averaging verifies the conservation of the volume and mass in the model. Subsections 2.3.1, 2.3.2, 2.3.4 and 2.3.5 test the interpolation of the material properties from the particle swarm to the finite element grid and show reasonable accuracy.

FEM discretization

In FAIStokes, the FEM (see Section 1.4.2 for a presentation of the method) algorithms are based on the deal.II library. The domain is discretized on a set of quadrilateral elements, on which FE basis functions are defined. As a review of the application of the FEM to the resolution of Stokes equations is provided in Section 1.4.2, only the spec-ifications of the FAIStokes code are presented here. For solving the Stokes equations, we use quadrilateral Taylor-Hood Q2 Q1 elements that satisfy the LBB conditions for stability (Donea et al. 2004). Contrarily to many creeping flow codes that are used to study the subsurface, the heat transport equation is not solved, both for simplicity and because it is likely to have only a small eﬀect on the strain at the scale at which struc-tural restoration is generally applied (i.e., basin-scale, close to the surface). Moreover, there may be important temperature diﬀusion at geological timescales, particularly in salt layers, and it is not reversible. As for the boundary conditions, we can use Dirichlet and Neumann conditions that we tune (e.g. rigidity, free-slip, free surface, specific trac-tion or velocity) for each boundary of the diﬀerent problems at hand. Subsections 2.3.1, 2.3.2, 2.3.4 and 2.3.5 showcase results of the FE benchmarking.

Grid and solvers

The grid and solvers come from the deal.II code, and their use is highly inspired from the deal.II tutorials step 31 1 and step 32 2. The grid is created first as a quadrilateral from the coordinates of the bottom left and top right corners of the domain. This quadrilateral is then split in order to get cells as close to a square as possible, and refined and coarsened adaptively several times to construct the initial grid. One of the advantages of using deal.II is that it ensures the continuity of the solution between the cells at the hanging nodes (nodes that are corners of a refined cell and in the middle of an unrefined cell line) by using specific constraints in the assembly (see lecture 16 from Wolfgang Bangerth’s online lectures for more details on the implementation in deal.II). The FE matrix, its preconditioner and the right-hand side force-vector are constructed using the material properties interpolated from the particle swarm, as described in Subsection 2.2.1. In the right-hand side of Eq. (1.11), the norm of the gravity vector g is always 9:81 m:s 2 in our simulations. Its direction is downwards by default, but it can be rotated with an angle given in the parameter file, in order to account for tilted models. The matrix system is solved using an iterative FGMReS solver preconditioned by a block matrix involving the Schur complement (Kronbichler et al. 2012). This solution is then used to refine and coarsen the grid adaptively using deal.II’s features, based on a gradient estimator in order to minimize the local error. The grid can also be adaptively refined depending on the viscosity changes or the presence of faults in the model. Depending on the input level of refinement, the cycle of building the matrix system, solving it, and adaptively refining and coarsening the grid is repeated several times, as shown in Fig. 2.1. Subsections 2.3.1, 2.3.2, 2.3.4 and 2.3.5 show the results of benchmarks that challenge the computation of the velocity on diﬀerent setups.

Table of contents :

1 State of the art
1.1 Introduction
1.2 Structural restoration
1.2.1 Geometric and kinematic restoration
1.2.2 Geomechanical restoration
1.3 Stokes flow
1.3.1 General equations
1.3.2 Application to restoration
1.4 Numerical methods
1.4.1 Material representations
1.4.2 The FEM method
Appendices
1.A Using a Stokes flow code to simulate other types of rheologies
2 FAIStokes
2.1 Introduction
2.2 Code and specificities
2.2.1 Material discretization of a geological model
2.2.2 FEM discretization
2.2.3 Grid and solvers
2.2.4 Velocity interpolation
2.2.5 Free surface implementation
2.3 Benchmarking a creeping flow code
2.3.1 Taking into account small scales inside a model : the Rayleigh- Taylor instability benchmark
2.3.2 Taking into account viscosity changes : the falling block benchmark
2.3.3 Advecting particles : the rotation benchmark
2.3.4 Taking into account the top surface in contact with air : the free surface benchmark
2.3.5 Upgrading the free surface movement : the sloshing benchmark
2.A FAIStokes parameter file
3 First results and first issues
3.1 Introduction
3.2 Results on simple models
3.2.1 The upscaled Van Keken model
3.2.2 The stratified overburden diapir model
3.3 Introducing more geological conditions
3.3.1 Dealing with the free surface in a backward creeping flow scheme
3.3.2 The simplified graben model
3.4 First assessment of the choice of parameters
3.5 Discussion on the results
4 Restoration of an analogue model
4.1 Introduction
4.2 From the laboratory experiment to the numerical model
4.2.1 Available data
4.2.2 Creation of the numerical model
4.3 Choosing the boundary conditions
4.3.1 Restoration using kinematic boundary conditions
4.3.2 Upgrading the kinematic conditions to natural boundaries
4.4 Choosing the model parameters
4.5 Discussion on the analogue model results
Appendices
4.A Impact of the Poisson coefficient
General conclusions
Bibliography