Dealing with the free surface in a backward creeping flow scheme

Get Complete Project Material File(s) Now! »

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 offsets, 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; Medwedeff 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 Medwedeff 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).

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: (1.2).
with C the elastic fourth-order tensor and the (small-)strain tensor associated to the displacement u by 2= 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.

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, different 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 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.
When using a Lagrangian frame to track the material model, the material param-eters need to be transfered from the Lagrangian frame to the Eulerian mesh nodes at each time step. The Lagrangian mesh can then be heavily distorted, as long as it is valid. The other way around, the velocity needs to be interpolated on the Lagrangian nodes to advect them. As mentionned above, the Eulerian mesh can move as well. Although it implies additional computation, tracking precisely the interfaces while min-imizing the need for a remeshing of the Lagrangian mesh make this method valuable for geomechanical simulations with very large deformations, compared to the Eulerian or Lagrangian methods.

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).

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 effect 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 diffusion 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 different problems at hand. Subsections 2.3.1, 2.3.2, 2.3.4 and 2.3.5 showcase results of the FE benchmarking.

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.
Appendices
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 

GET THE COMPLETE PROJECT