Project Topics
Get Complete Project Material File(s) Now! »
Structural restoration
Geological model restoration covers a number of different processes and methodologies. The classical techniques are unfolding and unfaulting using length/area preservation in order to remove the effects of tectonic forces. In addition to this, a large number of methods have been developed to take into account the effects of other important param-eters, like erosion and deposition of sediments (e.g., Dimakis et al. 1998; Chauvin 2017), isostasy compensation (e.g., Allen and Allen 2013), thermal subsidence due to mantle thermal effect (Royden and Keen 1980; Allen and Allen 2013), rock decompaction due to a change of load (e.g., Athy 1930; Durand-Riard et al. 2011; Allen and Allen 2013), or, at a smaller scale, the reverse migration of channelized systems (e.g., Parquer et al. 2017). These methods allow us to assess the temporal evolution of a structural model, in order to evaluate its consistency and test the hypotheses which led to its construction. It can also serve to quantify the shortening of an area, or analyze the strain in a current model. Paleo-basin geometries consistent with present-day observations can then be generated, for use in more elaborate hydro-mechanical forward models (e.g., Bouziat et al. 2019). In this manuscript, the focus is put on structural restoration based on unfolding and unfaulting. The following section is largely based on the bibliographic review on the subject by Chauvin (2017).
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).
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 effects. 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 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:
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 offer 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 offers 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 sufficient 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 affects 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 cutoff 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 difficult 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.
Figure 1.3: 2D synthetic example of how deforming a mesh through time can make it invalid for various computations, implying the need for a remeshing algorithm to be applied. (a) and (b) are the states of a mesh during displacement computation and after, (c) and (d) are the same at the next displacement iteration. (e) shows the problems occurring due to the mesh deformation: invalid crossing triangles and triangles that are too thin for the computation to be accurate on them.
To sum up, geomechanical restoration has overcome some limitations of the “clas-sical” geometric restoration process, but there is substantial need for improvement to better account for physical processes such as different rheologies, larger deformations, faults, and salt tectonics.
Stokes flow
While most rocks act as elastic materials when submitted to small deformations and for timescales of a few milliseconds to hundreds of years, they may exhibit a viscous fluid-type behavior when dealing with large deformations that occur at timescales of thousands, millions or more years (Moretti and Froidevaux 1986; Massimi et al. 2006; May et al. 2014; Cornet 2015; May et al. 2015). In particular, rock salt reacts as a viscous fluid when submitted to high stress (Nalpas and Brun 1993; Weijermars et al. 1993; Hudec and Jackson 2007). This has motivated the use of viscous flow simulations to study the evolution of salt structures, as well as other geological processes, both forward and backward in time. Moreover, we show in Appendix 1.A that a creeping flow simulation code can also be used to simulate other rheologies with only a few additional computations.
Application to restoration
The idea of a restoration method based on creeping flow was introduced at the EGS General Assembly in Hague in 1999 and later published in Ismail-Zadeh et al. (2004a). It was introduced as a way to improve the restoration of salt structures. Indeed, salt rock is less dense than most sedimental rocks, and behaves as a viscous fluid over geological timescales (Nettleton 1934; Woidt 1978; Poliakov et al. 1996). As such, salt layers form Rayleigh-Taylor instabilities with their overburdens, which can lead to salt structures taking a wide variety of shapes depending on the environment (Jackson and Hudec 2017). In particular, factors such as the type of overburden, the position and rate of erosion and deposition of sediments over it, the tilt of the depositional surface or the tectonic setting can have a large impact on salt flow and the resulting structures (e.g., Zaleski and Julien 1992; Podladchikov et al. 1993; Poliakov et al. 1996). This rheology and wide variety of complicated shapes make the conservation of layer area and thickness unable to perform restoration on models where salt has a large role. In this context, considering rocks as highly vicous fluids in a restoration scheme can, not only hold more physical results, but allow the restoration of some models to make more geological sense. Moreover, as said in the previous subsection, Stokes equations describe a steady-state flow where time appears only in the application of the velocity to the density and viscosity fields. As the computation of the velocity does not depend on previous times (contrarily to elasticity constraints which depend on previous strain), the advection can be applied with a negative time. This is the basis of backward time stepping restoration schemes: during the numerical advection, instead of applying Eq.
Figure 1.4: Example of the restoration scheme for a simple setup (a): as the arrows in (b) represent the velocity computed at a specific time step for a forward scheme, the advection of the material model in a restoration scheme is done with the opposite of the computed velocity, shown in (c).
In this light, representing the mechanical behavior of geological materials with vis-cous fluids holds several advantages, such as the use of ‘natural’ boundary conditions, like a free surface on top, and the account of rheologies like salt layers.
The first creeping flow restoration implementations used a linear viscous (Newto-nian) rheology, and proved to be able to restore 2D seismic cross-sections of salt di-apirs (along with a backstripping method for the sedimentation processes, see Fig. 1.5) (Ismail-Zadeh et al. 2001) and 3D Rayleigh-Taylor instabilities (Kaus and Podladchikov 2001; Ismail-Zadeh et al. 2004c). They were also used in the restoration of mantle ther-mal behaviors (e.g., Ismail-Zadeh et al. 2004b; Korotkii and Tsepelev 2007; Ismail-Zadeh et al. 2016). Since then, the method has been used for 3D unfolding in the absence of gravity (e.g., Schmalholz 2008), extended to non-linear (power-law) viscous behavior (e.g., Lechmann et al. 2010; Fernandez 2014) (Fig. 1.6), or used to study the reverse modeling of flanking structures (e.g., Kocher and Mancktelow 2005). Overall, this ap-proach has proven to allow the unfolding of sediment layers and the restoration of salt structures, both in 2D and in 3D. In the various previous applications, however, faults are either not present, not taken into account in the restoration process, or simply considered as frictionless surfaces. Also, the top surface in contact with air stays flat during the restoration process as the sedimentation and erosion processes are mostly considered fast enough to flatten the arising topography.
Numerical methods
In order to perform creeping flow restoration on geological models, one needs to solve the Stokes equations Eq. (1.11) in models where the velocity field evolves in time. However, there is usually no analytical solution for these equations in the studied models, due to their complexity. The solution then has to be approximated. In order to have enough precision and to deal with very large models, this is usually done numerically. The model is then discretized on a set of points with a prescribed connectivity, and the solution of the Stokes equations is approximated using methods such as the FEM or FDM.
Table of contents :
Introduction
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
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
4.A Impact of the Poisson coefficient
General conclusions
Bibliography
