Towards the application of Stokes flow equations to structural restoration simulations

Structural restoration is commonly used to assess the deformation of geological structures and to reconstruct past basin geometries. For this, geomechanical restoration considers faults as frictionless contact surfaces. To bring more phys5 ical behavior and better handle large deformations, we build on a reverse time Stokes-based method, previously applied to restore salt structures with negative time step advection. We test the applicability of the method to structures including sediments of variable viscosity, faults and non-flat topogra10 phy. We present a simulation code that uses a combination of Arbitrary Lagrangian-Eulerian methods and Particle-In-Cells methods, and is coupled with adaptive mesh refinement. It is used to apply the reverse time Stokes-based method on simple two-dimensional geological cross-sections and shows 15 that reasonable restored geometries can be obtained.


Introduction
The Earth's subsurface is the result of millions of years of deformation. Determining the deformation history from present-day structures has been a concern for geoscien- 20 tists who try to understand and quantify basin evolution. Restoration is an ensemble of methods which allow such quantification, by reversing processes that led to the current geometry of a geological region (e.g., Chamberlin, 1910;Dahlstrom, 1969). It covers a number of different processes 25 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, several methods have been developed to take into account the effects of other important parameters, like erosion and 30 deposition of sediments (e.g., Dimakis et al., 1998), 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 35 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 evaluate the consistency of a model and test the hypotheses which lead to its construction, in order to generate paleo-basin geometries consistent with present-day 40 observations for use in more elaborate hydro-mechanical forward models (e.g., Bouziat et al., 2019). In this article, we focus on the structural restoration based on unfolding and unfaulting. 45 Since the beginning of the last century, unfolding and unfaulting has been mostly done with geometric and kinematic rules (e.g., Chamberlin, 1910;Dahlstrom, 1969;Gratier, 1988;Rouby, 1994;Groshong, 2006;Lovely et al., 2018;Fossen, 2016). The first implementations in two 50 dimensions (2D) used balanced restoration, which relies on the conservation of layer bed area and thickness (e.g., Chamberlin, 1910;Dahlstrom, 1969;Groshong, 2006). Map restoration was then developed to study deformations which are mainly horizontal; it can be qualified as a 55 2.5D method (e.g., Cobbold and Percevault, 1983;Rouby, 1994;Ramón et al., 2016). Later, three dimensional (3D) geometrical methods have been proposed (Massot, 2002;Muron, 2005;Lovely et al., 2018), allowing the tracking of internal volumetric deformation. Such methods are all 60 based on the minimization of horizon deformation and on volume conservation, and therefore considerably simplify rock deformation mechanisms, ignore mechanical layering effects and are limited when considering salt basins. In this light, numerous authors have stressed out the necessity of 65 incorporating more physical principles into the restoration of geological models (Fletcher and Pollard, 1999;Ismail-Zadeh et al., 2001;Muron, 2005;Maerten and Maerten, 2006;Moretti, 2008;Guzofski et al., 2009;Al-Fahmi et al., 2016).
Several solutions have emerged for volumetric mechanicsbased restoration, that differ in terms of computational techniques and scale of the area of interest. These solutions can be divided in two main approaches, that have been developped to adress two different problematics of restoration. 10 They differ both in the mechanical laws used to compute the motion of rock layers, and in how these mechanical laws are applied to restore geological models.
The first approach considers the restoration of sediment layers assumed to deform elastically between frictionless fault surfaces. It 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 et al., 2010, 2013aTang 20 et al., 2016;Chauvin et al., 2018). In this approach, internal deformation is not known a priori, and the strain is computed from the mechanical behavior of rocks and the applied boundary conditions. The model is parameterized with elastic properties to mimic the response of rocks to mechani- 25 cal stresses and the restoration displacement is computed by solving the equation of motion, in which the Cauchy stress tensor is defined by Hooke's law. The restoration itself is performed by applying specific boundary conditions to constrain the model. These conditions, usually imposed on the 30 displacement, rely on the following assumptions: the uppermost horizon was flat and horizontal at deposition time, and it was not faulted. Other conditions can be introduced as complementary geological knowledge, such as direction and scale of deformation, or amount of lateral displacement 35 (Chauvin et al., 2018). Although these methods offer significant advances in the structural restoration of geological models, they still present many limitations. First, the boundary conditions set to unfold and unfault the medium are unphysical as the imposed depth of the free surface is the main 40 driver of the deformation (Lovely et al., 2012;Chauvin et al., 2018). These conditions are convenient hypotheses which do not necessarily reflect the paleo-stress state, hence they can be questionned (Durand-Riard et al., 2010;Lovely et al., 2012;Durand-Riard et al., 2013a). Secondly, geomechani-45 cal restoration so far only considers elastic rock properties, neglecting other possible behaviors, such as viscous, viscoelastic or plastic deformation (Gerbault et al., 1998). Transverse isotropic behavior also affects strain localization during restoration (Durand-Riard et al., 2013a), but such a behav- 50 ior is rarely applied in practice. These physical issues raise the question of the capability of this restoration approach to properly recover paleo-deformation. As a consequence, there are no clear guidelines on which method to choose between geometric and kinematic restoration and geomechani- 55 cal restoration, despite the more physical approach of the second one (Maerten and Maerten, 2006;Guzofski et al., 2009). Moreover, in spite of its name, geomechanical restoration is extensively controlled by geometric considerations: flattening of the top layer and a geometric unfaulting based on fric-60 tionless contact conditions to stitch the horizon cutoff lines accross each fault. Another practical issue is the need for a valid volumetric mesh of the structural model, including a boundary representation of the geological domain with the horizons and faults as boundaries (e.g., Muron, 2005), even if 65 the use of implicit horizons relaxes this constraint (Durand-Riard et al., 2010). Such a mesh is difficult to generate, as shown for example by Pellerin et al. (2014), Zehner et al. (2016 and Anquez et al. (2019). Since restoration deals with large deformations, the model evolves and may need to be 70 remeshed. The remeshing algorithms, however, are limited because key structural elements like faults and horizons must be preserved for geomechanical restoration to be used as an interpretation validation tool. To sum up, this restoration approach has overcome some limitations of the "classical" ge-75 ometric restoration process, by taking some of the internal movement of the layers into account for example, but it still needs to be improved to better account for different rheologies, larger deformations, faults, salt tectonics, and boundary conditions. 80 The second approach was introduced in 1999 as a way to improve the restoration of salt structures (Kaus and Podladchikov, 2001;Ismail-Zadeh et al., 2001, 2004Ismail-Zadeh and Tackley, 2010). It relies on considering the rocks as viscous fluids to compute the motion, and applying negative 85 timesteps. It is motivated both by the fact that rock salt and some sediment overburdens behave as viscous fluids over time-scales of millions of years, and by the reversibility of the Stokes equations, which allow the backward timestepping. The first implementations used a linear viscous (New-90 tonian) rheology, and proved to be able to restore 2D seismic cross-sections of salt diapirs (Ismail-Zadeh et al., 2001), and 3D Rayleigh-Taylor instabilities (Kaus and Podladchikov, 2001;Ismail-Zadeh et al., 2004). Since then, the method has been used for 3D unfolding in the absence of gravity (e.g., 95 Schmalholz, 2008), extended to non-linear (power-law) viscous behavior (e.g., Lechmann et al., 2010;Fernandez Terrones), or used to study the reverse modelling of flanking structures (e.g., Kocher and Mancktelow, 2005). Overall, this approach has proven to allow the unfolding of sediment lay-100 ers and the restoration of salt structures, both in 2D and in 3D. In the various previous applications, however, faults are either not present or not taken into account in the restoration process. Also, the top surface in contact with air stays flat during the restoration process as the sedimentation and ero-105 sion process are mostly considered fast enough to flatten the arising topography.
In this paper, we investigate a way of addressing some of the challenges raised by the first approach. We show that it is possible to push the second approach further and apply it 110 to models with faults and a non-flat free surface. For simplicity, we neglect the influence of temperature and consider the rocks as having a (linear) Newtonian rheology. While this considerably simplifies any non-Newtonian or visco-elastoplastic behavior in rocks, we show in the manuscript that this consideration is sufficient in simple setups. In the case of more complex overburdens, the method proved in the litterature to be able to restore various structures with powerlaw viscosities (e.g., Lechmann et al., 2010;Fernandez Terrones). We introduce a numerical scheme combining features 10 of the Arbitrary Lagrangian-Eulerian (ALE) and Particle-In-Cell (PIC) approaches, and using adaptive grid refinement. This specific implementation is motivated by the need for a moving topography, as well as the high accuracy needed for the computation of motion around the faults (). We show that 15 this scheme is accurate enough to consistently restore various geological setups, including faults.
The outline of this paper is as follows: we first present the concepts of Stokes flow-based restoration and its physical underpinnings. In a second part, we introduce the numer-20 ical code we developed for this application. Finally, we show the results that were obtained on an upscaled version of the model presented by van Keken et al. (1997), on a model with no prior knowledge on the material properties and boundary conditions to apply, and on a model with faults and a non-flat 25 free top surface.
and the mass conservation equation for incompressible fluids (continuity equation) where ∇ is the del operator, σ is the stress tensor, f is the specific body force (usually the volumetric weight ρg), and v is the velocity. The stress consists of a deviatoric part τ and an isotropic pressure p: where I is the identity tensor. In the viscous flow assumption, the deviatoric part of the stress is with η the dynamic viscosity and D the infinitesimal strain 45 rate tensor defined by Assembling Eq. (1), (3), (4) and (5), the momentum conservation equation can be written Here, we deal with materials that are highly viscous (with a viscosity η over 10 17 Pa.s), over time scales of thousands to millions of years, so these equations neglect the inertial part of the Navier-Stokes equations (Massimi et al., 2006). As such, they describe a steady-state flow and their resolu-55 tion provides the velocity of a fluid at a specific position and time. When different fluids are present, the conditions that are applied at their boundaries, as well as their differences in density, can create instabilities such as Rayleigh-Taylor instabilities. These instabilities make the flow non-stationary 60 as they advect the viscosity and density fields in time.

Restoration idea
In forward simulation schemes, the Stokes equations (6) and (2) are solved for pressure and velocity, and the material representation of the geological model is advected from the ve-65 locity at each time step. The simplest way to do it is by using an Euler scheme, the position x(t + ∆t) of each point of the material model after one time step being computed as with x(t) the position and v(t) the computed velocity of 70 the point at time t, and ∆t the time step (while higher-order methods exist (e.g., Ismail-Zadeh and Tackley, 2010), particularly to stabilize the advection scheme in the case of large time steps, we choose to present the restoration idea with this one for simplicity). This Finite-Difference approxima-75 tion relies on the idea that if the chosen time step ∆t is small enough, we can approximate the velocity of a particle as a constant over this time step (∆t is usually calculated using a Courant-Friedrichs-Lewy (CFL) condition (Courant et al., 1928) to ensure it). Since the Stokes equations are linear and 80 do not depend on previous time steps for the computation of the velocity, we can extend this approximation to backwards simulations. This is the basis of backward time stepping restoration schemes: instead of applying Eq.(7), we apply for the advection of the points of the material model, at each time step, like in Fig. 1. In this light, using viscous fluid properties instead of elastic properties to represent the mechanical behavior of geo-90 logical materials holds several advantages, such as the use of boundary conditions that are closer to reality, like a free surface on top, or the account of other rheologies like a salt layer. 3 Implementation in a specific code

Presentation
The restoration scheme presented in Sect. 2 has been implemented in the FAIStokes 1 code. It relies almost entirely on the deal.II library (Bangerth et al., 2007;Arndt et al., 2019, 5 2020) for all Finite Element related algorithms. The material tracking is based on the Particle-In-Cell (PIC) method (e.g., Asgari and Moresi, 2012;Thielmann et al., 2014;Gassmöller et al., 2016Gassmöller et al., , 2018Gassmöller et al., , 2019Trim et al., 2019). The general workflow of the code is shown in Fig. 2 and details of implementa-10 tion are discussed in the following sub-sections. Five benchmarks have been carried out to test the computation parts of the code and are presented in Appendices A, B, C, D and E.

Finite Element discretization
The Finite Element Method (FEM) was introduced in the 15 late 1950's (Hughes, 2012). Since then, it has emerged as one of the most powerful methods for solving Partial Differential Equations (PDEs) numerically. In FAIStokes, the FEM algorithms are based on the deal.II library. The do-1 Finite element Arbitrary Eulerian-Lagrangian Implementation of Stokes main is discretized on a set of quadrilateral elements, on 20 which Finite Element (FE) basis functions are defined. The aim of this paper is not to do a thorough review of the FEM, so only the specifications of the FAIStokes code will be presented here. For solving the Stokes equations, we use quadrilateral Taylor-Hood Q 2 × Q 1 elements that satisfy the 25 Ladyzhenskaya-Babuška-Brezzi (LBB) condition for stability (Donea et al., 2004). Contrarily to many creeping flow codes that are used to study the subsurface, we do not solve the heat transport equation, both for simplicity and because it is likely to have only a small effect on the strain at the scale 30 at which structural restoration is generally applied (i.e. basinscale, close to the surface). Moreover, there may be important temperature diffusion at geological time scales, particularly in salt layers, and it is not reversible. We use Dirichlet and Neumann boundary conditions that we adapt (e.g. rigidity, 35 free-slip, free surface, specific traction or velocity) for each boundary to the different problems at hand. Appendices A, B, D and E showcase results of the FE benchmarking.

Material discretization
The geomechanical simulation of a specific domain requires 40 to choose an appropriate kinematic description to follow the displacement inside the geological layers. Continuum me- . Schematic workflow of the FAIStokes code structure. The pre-refinement step occurs at the beginning of the simulation (or during a reinitialization of the grid) to ensure that the velocity used for the advection step is computed using the adaptively refined grid. chanics first distinguished 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, but neither of them is specif-5 ically 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 Arbitrary Lagrangian Eulerian (ALE) formulation (Fullsack, 1995;Donea et al., 2004), which inherits features from both methodologies, was developed. It has various formulations and implementations, both in 2D (e.g., Willett et al., 1993;Poliakov et al., 1996;Massimi et al., 2006Massimi et al., , 2007Fillon et al., 2013;Rose et al., 2017) and, more recently, in 3D (e.g., Braun, 2003;Thieulot, 2011;Thieulot et al., 2014). Most of these methods rely on keeping 15 track of the material properties in a Lagrangian way, while computing the displacement on a grid that can only deform vertically to account for a possible free surface. It is particularly useful in geomechanics, where the vertical deformation is generally small compared to the horizontal deformation, 20 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 grid has an ALE part as it can adapt to follow the movement of the free surface. 25 During mechanical simulations, the material properties inside the model are tracked using particles; each of these particles discretize the small part of the model around them and its properties. At each time step, the material properties of the particles are projected onto the grid. They are then used 30 to solve the Stokes equations on the grid. Following this, the particles are advected using the solution on the grid.

The PIC method
At the begining of the simulations, FAIStokes either creates a model from a function giving the distribution of the material parameters or loads a particle swarm from a file. In 35 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 given function is then used to associate the material properties to the particles depending on their position. Since the particle swarm 40 doesn't directly track the interfaces, it has to be dense enough to recover accurately the material properties of the model; depending on the simulation, some parts of the model can therefore be densified to keep the appropriate accuracy. At each time step, the material properties are interpolated from 45 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 distance depending on the smallest element 50 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 effect of very high viscosity differences (possibly of several orders of magnitude) on the solver and is more computationally efficient despite the higher grid refinement needed (Deubelbeiss and Kaus, 2008;Thielmann et al., 2014;Heister et al., 2017). In the simulations we present hereafter, we were able to verify that this averaging verifies the conservation of the volume and mass in the model. Appendices A, B, D and E test the interpolation of the material properties from the particle swarm to the 10 finite element grid to reasonable accuracy.

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 2 and step − 32 3 . The grid is created first as a quadrilateral 15 from the coordinates of the bottom left and top right corners of the domain. This quadrilateral is then split in order to get cells closest to a square (depending on the model bounding box size) and refined and coarsened adaptively several times to construct the initial grid. The FE matrix, its precon-20 ditioner and the right-hand side force-vector are constructed using the material properties interpolated from the particle swarm as described in the previous subsection. In the righthand side, the norm of the gravity vector g of Eq. (6) is always 9.81 m.s −2 in our simulations, and its direction is al- 25 ways downwards. 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 30 in order to minimize the local error. 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. Appendices A, B, D, E show the results of benchmarks that tested the compu-35 tation of the velocity on different setups.

Velocity interpolation
Once the grid refinement has been completed, the particle swarm is advected by the obtained solution. In FAIStokes, the interpolation of the velocity is done separately in each 40 grid cell with a Q 2 interpolation scheme. Depending on whether the simulation is forward or backward, the displacement of each particle for a time step ∆t is computed using Eq. (7) or (8). The value of ∆t is computed from the CFL condition. The default value for the CFL number is 0.085, 45 but it can be reduced depending on the simulation (for example, the results shown in the next section use a CFL number of 0.0085, while the benchmarks in the Appendix use a CFL number of 0.042). The advection is done with a 2 nd -order 2 https://dealii.org/9.0.0/doxygen/deal.II/step_31.html 3 https://dealii.org/9.0.0/doxygen/deal.II/step_32.html Runge-Kutta scheme in space: at each time step, the parti-50 cles are first advected by half the computed displacement; the velocity is then interpolated on their new position to update the displacement, and particles are advected them again by half of this new displacement. This scheme reduces the error in the advection process without need for simulation time 55 step refinements. It is computationally efficient because the computation of the displacement on the particle swarm is inexpensive as compared to solving the FE matrix system. Appendices C, D, E show the results of benchmarks that tested the interpolation of the velocity in time-dependant problems. 60

Free surface implementation
In the case of a free surface on the top of the model, the top surface is tracked by a separate point swarm. This point swarm is denser than the material particle swarm and is one dimension lower (i.e. a line in our 2D cases). It is advected at 65 each time step the same way as the particle swarm that represents the geological model. After its displacement or during the setup of the grid, the free surface point swarm is used as a reference to move vertically the nodes of the grid at the top of the model, so that they match the free surface. This 70 vertical displacement is then propagated to the rest of the grid so that the grid cells stay as close to squares as possible, while not affecting the other boundaries. Fig. 3 illustrates the whole process. Since our models are isothermal, no special processing is required to correct the temperature field during 75 this process. Appendix D shows the results of a benchmark that tests the free surface implementation along with other computational parts of the code. The free surface stabilization algorithm (refered to as FSSA in the rest of the paper) developed by Kaus et al. (2010) and showcased in Quinquis 80 et al. (2011) has been implemented in FAIStokes; we benchmark it in Appendix E.

Results
In addition to the benchmarks presented in the Appendices, which mainly check the algorithms of the code, we tested 85 our specific restoration scheme on three simple models. In those experiments, the boundary conditions are simplified and quite unrealistic, but the goal here is to check the behavior of the reverse-time modeling in simple settings. In particular, we choose to neglect basal and lateral displace-90 ments in the first two models, which are known to play a role in salt tectonics (Koyi, 1996;Ismail-Zadeh et al., 2004) but would require a calibration and would increase the degrees of freedom of the problem. 95 The first model is scaled-up from van Keken et al. (1997). The setup consists of a simple two-layered system driven by gravity, as shown in Fig. 4. The upper layer represents sediments that are denser than the lower layer which contains salt (ρ o = 2600 kg.m −3 for the sediment layer and ρ s = 2150 kg.m −3 for the salt layer). A sinusoidal instability initiates the movement at the begining of the simulation.

Diapiric growth model
The model is limited to a 10 km × 9.142 km domain (the 5 width value is given by van Keken et al. (1997) to yield the largest growth rate for the diapir) with free slip boundary conditions on the sides and no slip boundary conditions on the top and bottom sides. The grid has 32 2 initial elements and two levels of additional adaptive refinement. The particle swarm has a heterogeneous particle density: it is first sampled regularly in the model and then densified to five times more particles around the interface between the two layers 5 to facilitate the tracking of material properties. The average distance between two particles near the interface is 14.3 m. The total number of particles is 64, 000. Two experiments were performed in this model: the first one as a test with isoviscous materials (η o = η s = 10 19 Pa.s), the second one 10 with material properties closer to reality with a lower viscosity for salt (η o = 2.8 × 10 19 Pa.s for the sediment layer and η s = 1.4 × 10 17 Pa.s for the salt layer).
For each experiment, we first did a forward simulation, and then we applied the restoration scheme to the results ob-15 tained at the end of the simulation. The state obtained after 6×10 6 years for the first test and 1.5×10 6 years for the second test, as well as the restored models, are shown in Fig. 5. We can see that while the isoviscous experiment has a rather smooth forward result, the second experiment with a less vis-20 cous salt leads to the creation of a salt weld (surface where the salt layer thickness has reached or almost reached zero, the salt having creeped away) at the bottom and left-hand side of the model.
To check the quality of the restoration in the two exper- 25 iments, we compute for each particle the distance between its original position before the forward simulation and its position at the end of the restoration process. given in Fig. 6 and Fig. 7, and help explain this phenomenon. The high error values in the variable viscosity case are due to the creation of a basal weld, which mixes the particles at the bottom of the model. Some of these particles are not well restored and stay at the bottom of the model, creating very 40 large errors (hence the error bars of 1 to 20 particles with an error higher than 500 m in Fig. 7). The basal weld in itself creates large distortions which explain the overall large errors at the interface. However, if we look at the model at the end of the experiments in a global way, not taking into account 45 small irregularities, and study only the boundary between the two layers, the maximum distance between the initial model and the restored model is only 50 m (0.5% error) for the first experiment and 125 m (1.25% error) for the second, which is acceptable considering the large amount of total deforma-50 tion.

Stochastically generated salt diapir model
This model was generated with the method proposed by Clausolles et al. (2019). It consists of a salt diapir that mimicks passive diapirism structures created by syndeformation 55 differential sediment loading. The input for the salt diapir is a seismic image interpreted to segment it in three regions: salt, sediment, and uncertain. The salt-sediment interface is then generated in the uncertain zone, from available data, geological knowledge and a random scalar field that takes into 60 account the uncertainties. The setup is quite simple but interesting for two reasons. First, this model was not created by a forward viscous simulation, and the rheology of the salt and sediments is not known. Second, this model has a high uncertainty and it is uncertain wether the boundary condi-65 tions we apply can restore it or not. Therefore, this test case can be assimilated to the simplification of a real case application. The initial particle swarm contains 102, 510 particles regularly sampling the model, and we apply free-slip boundary conditions on the top and side model boundaries, and a 70 no slip boundary condition on the bottom. Figure 8 shows the initial state of the model. The grid has 48×80 initial elements and three levels of additional adaptive refinement; its state at the beginning of the simulation is shown in Fig. 9. In order to assess the influence of the value of the parameters on the 75 results of the restoration, we tested different possibilities. As this is a simplification of a real case application, and there 90 is no information on the type of sediments, in each experiment the density and viscosity are homogeneous in the sediment and salt layers. The results for the 5 experiment simulations are given in Fig 10. Depending on the experiment, we choose to stop the 95 restoration process after different durations t end . Indeed, as the viscosity and density vary from one experiment to the other, so does the model relaxation time.
Overall, the restoration process removes the diapir and leaves a salt scar, while the sediment layers remain glob-100 ally flat. Since this setup is generated by a method for syndeformation diapirs, a full restoration of the model should have taken into account the deposition of the sediments at the same time as the formation of the diapir, by removing the sediment layers one by one. For simplification purposes and in order to test the process with simple boundary conditions, such sedimentation processes were not implemented.

5
In this case, the stress state inside the model being incorrect, the sediment and salt layers couldn't be restored to a completely flat state. For example, the shallow sediments should have been removed early in the restoration process, and as such were deformed beyond their sedimentation point. While do not have a big impact on the shape of the model after the restoration process. Only the shape of the sediments at the base of the diapir is slightly different from one experiment to the other. Experiments 4 and 5 have serrated shapes that are not geologically probable, probably because of the four or-20 ders of magnitude of viscosity contrast between the salt and sediments. The main difference between the experiments is the relaxation time for each restoration process. If the duration of the formation of the diapir was known, it could then be used to reduce the uncertainty on which density and vis-25 cosity to use. It also seems that the curvature of restored layers changes. This could provide another criterion to further evaluate the results (but would call for variable sediment viscosity testing). However, this is a difficult path forward, because sediment deposition clearly plays a major role during 30 salt displacement (e.g., Giles and Lawton, 1999;Hudec and Jackson, 2007;Rowan et al., 2012). Moreover, recovering the full deformation path during sediment deposition would  also call for further studies, possibly on laboratory analogs (Weijermars et al., 1993;Dooley et al., 2005).

Simplified graben model
The last model is a simplification of the creation of a graben in sediments submitted to lateral flow in an extensive context, 5 as shown in Fig. 11.
It consists of a layered overburden underlain by salt and cut by two 60°faults. As the intent of the manuscript is to focus on the restoration of structural models, we do not consider the formation of the faults with plasticity, but rather 10 start with two faults already present. The domain size is 6 km horizontally and 2 km vertically, and the right boundary is subjected to lateral flow. This is modeled by a Dirichlet condition applying a specific value for the velocity in the horizontal direction (the vertical value for the velocity is still 15 free, as free slip boundary conditions). The left and bottom boundaries are set to free slip, and the top boundary is considered a free surface. For the model to evolve without interference with the lateral flow on the right boundary, the faults are positioned at one third of the model width from the left. 20 In order to capture the geometry of the faults, which is essential in this setup, the adaptive refinement of the grid is an important feature of the proposed implementation. As such, the grid is refined specifically near the faults, with 60 × 20 initial elements and four levels of adaptive refinement (up to 25  . Adaptively refined grid for the first time step of the simulation. We can see that the grid is refined to a high level at the interface between the salt and the sediment overburden, where the highest velocity gradients appear. On the contrary, it is coarsened where the velocity has small gradients, particularly in the upper right and upper left corners. 6.25 m × 6.25 m cells near the faults), as shown in Fig. 12 for the first time step. The faults are considered as shear bands with a lower viscosity and density than the rest of the overburden. The overburden is layered with two types of rocks with slightly dif- Like in the diapiric growth model, we first do a forward simulation, and then apply the restoration scheme on the model 15 obtained. The lateral flow is set to 10 mm/year outwards and the model is observed for 35, 000 years, both in forward and backward simulations, to have sufficient deformation. The particle swarm has a heterogeneous particle density, with a regular sampling in most of the model and eight times more 20 particles near the faults and around the interface between the salt and the overburden. The average distance between two particles in the densified zone is 2.5 m. The total number of particles at the first step of the simulation is about 330, 000. This number decreases during the forward simulation, as the particles are removed once they flow outside of the model. During the restoration simulation, new particles are added near the right boundary when the particle swarm flows in-5 ward due to the negative time stepping. Their material properties are determined from the particles moving inwards and their motion. The model at the end of the forward and backward simulation is shown in Fig. 13. Figure 14 shows the difference between the position of the interfaces before the 10 forward simulation and at the end of the backward simulation. The numbering of the interfaces follows their position in the model, interface 0 being the lowest salt-sediment in-terface. The mean and maximum error for each interface are given in Table 1.

15
Overall, the results for the restoration simulation are quite good, with mean errors around 1% of the forward deformation for the layer interfaces. Fig. 14 shows that the graben part of the model (between the two faults) is approximatively 7 m lower than it should be at the end of the restoration. This 20 is due to the model topography being slightly tilted from the right-hand side of the model towards the faults at the end of the restoration. The largest errors are located on the two faults. In particular, the free surface behaves well during the simulation, except for some small instabilities occuring on 25 the top of the faults, where high viscosity contrasts occur. The error resulting from these instabilities, however, is quite Faults (η f , ρ f )

o o
Salt layer (η s , ρ s ) Figure 11. Setup of the simulation for the simplified graben model. The initial particle swarm is densified near the faults and at the interface between salt and sediment. Figure 12. Adaptively refined grid for the first time step of the simulation. We can see that the grid is refined to a high level near the faults, where the velocity gradient is high. On the contrary, it is coarsened where the velocity has small gradients, particularly in the lower right and left corners. small compared to the amount of deformation in the model (around 200 m of slip on the faults).

Discussion
While the results of the three test models in the previous section are promising, their purpose is not to correctly compute 5 the deformation of the subsurface in a forward mechanical simulation, but rather to assess the validity of the proposed restoration scheme and the underlying concepts in various geological cases.
In this paper, we have made some links between two different types of structural restoration approaches. On one side, geomechanical restoration methods have been relying on considering rocks as elastic material and flattening the top surface of the horizons (e.g., Guzofski et al., 2009;Lovely et al., 2012;Chauvin et al., 2018), and may lead to unphysi-15 cal strains. On the other side, dynamic restoration methods have considered viscous fluid rheologies for the rock layers, the deformation being driven by density contrasts and boundary conditions, and applied with a backward advection scheme (e.g., Kaus and Podladchikov, 2001;Ismail-Zadeh 20 et al., 2001;Lechmann et al., 2010). As such, restoration us-ing Stokes equations is expected to provide more physical strains, given that the boundary conditions and the rheology inside the model are close enough to reality. This method is not new (e.g., Kaus and Podladchikov, 2001;Ismail-Zadeh 25 et al., 2001), but has been restricted to the restoration of salt structures and small-scale folds, in environments where the topography remains flat. Here, we are interested in applying it to environments with faults and large displacements of the topography. In this scope, we introduce a numerical 30 scheme combining features of the ALE and PIC approaches, and using adaptive grid refinement. We show that it is accurate enough to produce consistent results on the restoration of models with a viscous backward advection approach. When applying a reverse-time Stokes restoration scheme, 35 two important questions appear: what are the material properties of the geological objects inside the model, and what type and intensity of boundary conditions should be applied to these geological objects? Regarding the material properties, the diapir test model of Sect. 4.2 gave a first idea of 40 how to choose them, and previous articles have considered the question on specific setups (e.g., Lechmann et al., 2010). The density of the subsurface depends on the type of rocks that are present, and its estimation is relatively easy. The viscosity, however, is not trivial, as laboratory observation time 45 scales are too short to reflect the slow movement occuring at geological time scales. The values we took are inspired from numerical simulations, but they have a large uncertainty (at least one order of magnitude) (e.g., Massimi et al., 2006;Kronbichler et al., 2012), as they are calibrated using postglacial rebound data for example. Works on analog sandbox experiments and further experiments on models with more geological knowledge should prove to be useful in estimating a proper viscosity for the restoration of different rock rheologies. In particular, the duration over which the geologic phe- Most geomechanical restoration schemes consider basin rocks to have an elastic behaviour, whereas the rheology used herein does not display Poisson effects. Incorporating elastic-15 ity in viscous flow has been done, for example by using an effective viscosity to account for the elastic part of the material while minimizing the modifications to the viscous flow code (e.g., Moresi et al., 2003). The problem is that these schemes, like every implementation of elasticity, use values of the 20 stress and strain at previous time steps. The elastic behavior then depends completely on the stress state at the beginning of the simulation, which is not available in restoration schemes. However, specific material properties could still be taken into account in other ways in Stokes-based restoration. 25 For example, the incompressibility constraint, which implies a Poisson's ratio of 0.5, can be relaxed, which could be used to account for lesser values of the Poisson's ratio. Regarding the rheology of faults, we cannot directly use their usual forward modeling implementation considering the rock as hav-5 ing a plastic behavior. Indeed, the previous stress history is needed to simulate such a behavior, and it is not available in restoration, which studies backward movement. Here, we used a specific viscosity for the implementation of faults in restoration, which holds two advantages. First, since all the 10 faults are already identified at the beginning of the restoration process, we do not need to allow the creation of faults in backward simulations. Second, using an effective viscosity for the faults allows for a more realistic simulation of shear band and damage zone behavior, compared to previous ge-15 omechanical restoration schemes that consider faults as freeslip surfaces. A significant issue with the boundary conditions in geomechanical simulations is the difficulty to estimate the paleoforces at play several kilometers underground. We therefore 20 need to choose Dirichlet and Neumann boundary conditions that best fit the tectonic knowledge about the region of study. For example, deformation is generally strongly influenced by the horizontal stress state, implying compressive or extensive structures and the need for corresponding conditions 25 on the side boundaries (Chauvin et al., 2018). Another example is the top surface of the models, which can be considered as being on ground level, and is therefore in contact with air. This interface is complicated to handle due to the several orders of magnitude in the material property contrast 30 (very high density and viscosity for rocks versus almost null density and viscosity for the air). In geomechanical simulations, several approaches exist to model its behavior. The simplest topographic surface solution is to set a free-slip condition which removes the normal component of the velocity 35 at the boundary. This simplification is mostly used in cases where the movement of the top surface is negligible compared to the rest of the model. In order to do more realistic simulations, two main approaches are available: the implementation of a free surface, or the "sticky air" method (e.g., 40 Crameri et al., 2012, for a benchmark and a comparison of the two methods). The sticky air method considers a layer of material with a low viscosity and zero density, the difficulty being that this viscosity needs to be sufficiently low to be negligible compared to the rest of the model, but high 45 enough for the solvers to converge. The free surface method considers that no force is applied on the surface of the computational mesh. While this is theoretically simple, it is numerically complicated to implement, as it also means that the computational mesh needs to honor the movement of the free 50 surface. In FAIStokes, the free surface method is applied by tracking the movement of the top surface and allowing the grid nodes to move vertically (Sect. 3.7). In order to stabilize its movement and avoid some of the instabilities that can appear, the free surface stabilization algorithm presented in 55 Kaus et al. (2010) has been implemented. The free surface implementation has been benchmarked and performs well in forward simulations (Appendices D and E). In restoration simulations, the results are more mitigated. Indeed, in models where the only drive is a density contrast (such as the 60 models shown in Sect. 4.1 and Sect. 4.2), the free surface shows instabilities. This appears particularly when working with models that have a near-horizontal or initially horizontal top surface. In those setups, any small computational error in the computation of the vertical part of the velocity can lead to 65 instabilities that increase exponentially in reverse time. Several approaches involving specific tractions on the top surface have been tested to remove or correct this instability, but we have not yet devised any efficient means to prevent it in such setups. In particular, the FSSA delays this phenomenon, but 70 does not suppress it altogether. In models where other drives for the deformation occur, such as lateral flow, the results are more promising, as shown in Sect. 4.3. Indeed, in such setups, the boundary conditions introduce larger strains that dampen the instabilities. For test and comparison purposes, 75 the sticky air method has also been implemented and coupled with the FSSA. It uses the moving grid feature of the free surface so that the cells containing air and rock layers are distinct, and performs similarly to the free surface on the benchmarks presented in Appendices D and E. In restoration 80 simulations, it can delay the instabilities that appear in models driven exclusively by gravity, but doesn't remove it.

Conclusions
We have presented a scheme that exploits the reversibility of Stokes equations to perform structural restoration on different geological setups. While the principle of the method is not new, we have shown that it may be applied on models with faults and a non-flat topography. As such, it may improve some of the issues with the current geomechanical restoration implementations that are used for such environments. The FAIStokes code was developed to apply this restoration scheme and allow various tests on its implemen-10 tation. Among those tests, we presented three simple models and the results we obtained with them. Those results are encouraging, although the numerical method has difficulties dealing with the restoration of salt in the presence of welds. The free surface is well managed in our experiment includ-15 ing lateral flow, but also leads to instabilities in the restoration process when the flow inside the models is driven exclusively by density. Overall, we still show that combining adaptive grid refinement with the PIC and ALE approaches gives enough accuracy to produce consistent restoration re-20 sults on different model setups.
We intend to follow this work by applying the method to more complex models, starting with the restoration of sandbox experiments (e.g., Colletta et al., 1991). This will allow us to do more precise tests on the value to choose for the vis-25 cosity and density of geological layers, and to upgrade the specific implementation of faults.
Code availability. The code corresponding to this paper is available to members of the RING consortium in the FAIStokes software. The Rayleigh-Taylor instability by Ramberg (1968) and was carried out in various numerical studies (Deubelbeiss and Kaus (2008); Thieulot (2011)). It consists of a two-layer system driven by gravity, the density of the bottom layer being smaller. The bottom and top boundaries have a no slip bound-45 ary condition, while the sides have a free slip boundary condition.
The first layer, made of fluid 1 with properties (ρ 1 , η 1 ), overlays the second layer, made of fluid 2 (ρ 2 , η 2 ). An initial sinusoidal disturbance of the interface between the two 50 layers is introduced, characterized by an amplitude ∆ and a wavelength λ, as shown in Fig. A1.
Under these conditions, the velocity of the diapiric growth v is given by Ramberg (1968): with K the dimensionless growth factor given by which involves the following factors: We set ρ 1 = 3300 kg.m −3 , ρ 2 = 3000 kg.m −3 , η 1 = 60 10 21 Pa.s, L x = h 1 +h 2 = 512 km, and ∆ = 3 km. We make η 2 vary between 1.25 × 10 20 and 2.5 × 10 23 Pa.s, while λ takes three values: L x /2, L x /4, L x /8 (Fig. A2). A first run is done, where the FEM grid is fixed to 80 × 80 elements, each containing 10 2 regularly spaced particles. In 65 order to test the influence of adaptive refinement, we conduct a second run with a grid starting at 80 × 80 elements and three levels of adaptive refinement. We also refine the particle swarm adaptively: each initial cell is first filled with 5 2 regularly spaced particles, and then the swarm is densified to 70 64 times more particles around the interface between the two fluids. The results are shown along with the analytical ones in Fig. A3.
Overall, results show a good agreement between the computed solution and the reference, especially in the case of 75 adaptive refinement, where the relative error falls beneath 2.5% for all the curves. Since φ 1 is inversely proportional to the wavelength λ, it means that the code can account well for small disturbances, especially with the use of adaptive refinement on the parts with higher velocity and high contrasts 80 in viscosity.
This benchmark ensures the validity of the code in the presence of large viscosity constrasts, even if those constrasts  This benchmark appears in Gerya (2019) and is presented in Thieulot (2011). It consists in modelling the fall of a block of fluid of properties (ρ 1 , η 1 ) inside another fluid of properties (ρ 2 , η 2 ), with ρ 1 > ρ 2 . The domain is a square of size L x = L y = 500 km, and the block (a square in 2D) of size 10 100 × 100 km is initially centered at point (x = 250 km, y = 400 km), as shown in Fig. B1.
The simulation is carried out on a 50 × 50 element grid that is adaptively refined three times. Like in the previous benchmark, the particle swarm is created by first introducing In all the experiments, the density of the surrounding fluid 25 is ρ 2 = 3200 kg.m −3 and the viscosity of the block is varied between 10 19 and 5 × 10 27 Pa.s. The velocity of the falling block is measured in its centre at t = 0 for all experiments. Following physical intuition, one expects the velocity of the block to act as follows: (a) decrease when the vis-30 cosity of the surrounding fluid η 2 increases (i.e. when going  from Exp.1 to Exp.5), and (b) increase with the density contrast (ρ 1 − ρ 2 ) in each experiment. To check this behavior, we measure vη 2 /(ρ 1 − ρ 2 ) as a function of the viscosity contrast log 10 (η 2 /η 1 ). The results of the benchmark are plotted in Fig. B2. 5 We can see that the experimental points line up on a single curve, which shows that FAIStokes can deal with gravitydriven simulations where 0.6% ≤ (ρ 1 − ρ 2 )/ρ 2 ≤ 210% and the viscosity contrasts are as strong as 10 −6 ≤ η 2 /η 1 ≤ 10 5 in a consistent manner. 10 Appendix C: Advecting particles : the rotation benchmark The last benchmark aims at assessing the error in the advection part only. The setup of the model is a square of size 10 × 10 km, where we study the advection of a single par-15 ticle, starting at coordinates (8 km, 5 km) and doing a 2π rotation around the center point (5 km, 5 km) (Fig. C1). A velocity field is prescribed in the domain and discretized on the grid: on each grid point, the velocity has a constant norm The grid is not adaptively refined here, and is composed of 16×16 elements. In order to have scales that are geologically 5 relevant, we choose v 0 = 3 cm.year −1 and vary the time step ∆t between 500 and 2000 years (in this setup, the CFL numbers chosen for our simulations would give a timestep be-tween 175 years for the lowest CFL number and 1753 years for the highest CFL number). The second order Runge-Kutta 10 scheme presented in Sect. 3.6 is used at all time steps. We then evaluate the distance ∆r = |r(θ = 0)−r(θ = 2π)|. This distance gives us a measure of the error made in the computation of the particle advection, and allow us to compare different advection schemes. Figure C2 shows the results obtained 15 for a 2π rotation of the particle with different interpolation schemes. We can see that reducing the timestep linearly reduces the error on the radius r(θ). In this setup, the type of  Figure C2. Results for the rotation benchmark obtained with different time steps and advection schemes.
interpolation mostly impacts the stability of the interpolation, and not the accuracy.
Appendix D: Taking into account the top surface in contact with air : the free surface benchmark This benchmark is presented in Crameri et al. (2012), where 5 it is applied on several numerical codes to compare their implementation of the free surface, and evaluate the use of the 'sticky air' method. It will be used here to evaluate the quality of our approximation and interpolation of the free surface. It consists on a cosine-shaped layer of homogeneous litho-10 sphere overlaying a homogeneous layer of mantle. For this type of model, Ramberg (1981) gives an analytical solution for the maximal height of the topography at each time t : where γ is the relaxation rate and h initial is the value of h 15 at the beginning of the simulation. The model setup for the benchmark is shown in Fig. D1. The bounding box of the model spans 2800 km by 707 km. The underlaying mantle layer is 600 km thick, while the lithosphere has a thickness between 93 and 107 km. The 20 lithosphere's top surface is cosine-shaped with an amplitude of 7 km and a wavelength of the size of the domain. The mantle and lithosphere have a density of ρ M = ρ L = 3300 kg.m −3 and a viscosity of η M = 10 21 Pa.s and η L = 10 23 Pa.s, respectively. We set free slip boundary condi- 25 tions for the sides and a no slip condition on the bottom of the model. The initial grid is made of 16 × 64 elements and is adaptively refined 3 times. The particle swarm contains 484, 160 particles; it is constructed by first sampling regularly the domain, and then adaptively densifying the swarm 30 to 64 times more particles in the lithosphere and upper part of the mantle. In this setup, Crameri et al. (2012) gives a characteristic relaxation rate γ = 0.2139 × 10 −11 s −1 and a characteristic relaxation time t relax = 14.825×10 3 year. The results obtained with FAIStokes are given in Fig. D2. 35 The numerical results are close to the analytical ones, with only a 1.3% error at the characteristic relaxation time. This shows the capacity of FAIStokes to compute the solution of Stokes equations with a free surface for small vertical deformation, and to advect the particles inside the model. 40 It also gives another evaluation of the handling of gravitydriven flow, this time with the addition of evolution inside the model. This benchmark is presented in Kaus et al. (2010), where it is used to assess the results of the free surface stabilization algorithm (FSSA) presented in the same article. It is used here to verify the implementation of the same algorithm in our code, as well as check the behavior of the free surface 50 in another setup. The benchmark model is another Rayleigh-Taylor instability with a dense, more viscous layer sinking into a less dense fluid (Fig. E1).
The model span is 500 km × 500 km; the side boundaries have a free-slip condition, the lower boundary is no-55 slip, and the top boundary is a free surface. The initial perturbation between the two layers is sinusoidal with an amplitude of 5 km. The computation is carried out on a grid with 25 × 25 initial elements and three adaptive refinement steps. The particle swarm counts 25, 000 particles; it is constructed 60 by first sampling regularly the model and then densifying it to 64 times more particles around the interface. The specificity of this benchmark is the apparition of a sloshing instability (also refered to as the "drunken sailor" instability) if the simulation time step is too large. Specifically here, without the FSSA, the forward simulation is stable with a time 5 step ∆t of 2500 years, but with ∆t = 5000 years, an instability emerges as the velocity pattern changes direction from one time step to the other (Fig E2). In order to follow the evolution of the free surface, we keep trace of the altitude of the most top-left point over time. Re-10 sults of a 0.5 Myr simulation, for different time steps ∆t, with and without the FSSA, are shown in Fig. E3. We can see that the implementation of the FSSA stabilizes the sloshing behavior of the free surface that appeared whith a time step ∆t of 5000 years, and keeps the free surface stable 15 even with higher time steps. Moreover, the results show great similarities to those that can be found in Kaus et al. (2010) and Thieulot (2019). This validates the implementation of the free surface stabilization algorithm. It also gives another evaluation of the handling of gravity-driven flow with a free 20 surface, this time with the additional resolution of an instability that can occur with free surfaces.
Author contributions. The writing of this paper, as well as the code of the software and the results, was done by M S-S. CT helped with useful discussions and new points of view on the subject of the pa-25 per, as well as help on the code structure, benchmarks, and methods. GC and PC provided useful discussions on the subject, guidance on solving the different problems, and funding for the PhD thesis that allowed this research through the RING-Gocad consortium. All authors contributed to the conceptualization, the design of 30 experiments and the analysis of results. CT, GC and PC also helped on the manuscript with useful feedback and corrections.
Competing interests. The authors declare that they have no conflict of interests.