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 physical 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 topography. We present a simulation code that uses a combination of arbitrary Lagrangian–Eulerian methods and particle-in-cell 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 that reasonable restored geometries can be obtained.

1 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 geoscientists 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 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 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 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 observations for use in more elaborate hydromechanical forward models (e.g., Bouziat et al., 2019). In this article, we focus on the structural restoration based on unfolding and unfaulting.

Since the beginning of the last century, unfolding and unfaulting have 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 dimensions (2-D) 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 2.5-D method (e.g., Cobbold and Percevault, 1983; Rouby, 1994; Ramón et al., 2016). Later, three-dimensional (3-D) geometrical methods have been proposed (Massot, 2002; Muron, 2005; Lovely et al., 2018), allowing the tracking of internal volumetric deformation. Such methods are all 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 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 mechanics-based 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 developed to address two different problematics of restoration. 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, 2013a, b; Tang 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 mechanical 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 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 (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 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 questioned (Durand-Riard et al., 2010; Lovely et al., 2012; Durand-Riard et al., 2013a). Secondly, geomechanical restoration so far only considers elastic rock properties, neglecting other possible behaviors, such as viscous, visco-elastic or plastic deformation (Gerbault et al., 1998). Transverse isotropic behavior also affects strain localization during restoration (Durand-Riard et al., 2013a), but such behavior 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 geomechanical 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 frictionless contact conditions to stitch the horizon cutoff lines across 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 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 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” geometric 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.

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, 2004; Ismail-Zadeh and Tackley, 2010). It relies on considering the rocks as viscous fluids to compute the motion and applying negative time steps. It is motivated both by the fact that rock salt and some sediment overburdens behave as viscous fluids over timescales of millions of years, and by the reversibility of the Stokes equations which allows the backward time stepping. The first implementations used a linear viscous (Newtonian) rheology and proved to be able to restore 2-D seismic cross-sections of salt diapirs (Ismail-Zadeh et al., 2001) and 3-D Rayleigh–Taylor instabilities (Kaus and Podladchikov, 2001; Ismail-Zadeh et al., 2004). Since then, the method has been used for 3-D unfolding in the absence of gravity (e.g., Schmalholz, 2008), extended to non-linear (power-law) viscous behavior (e.g., Lechmann et al., 2010), or used to study the reverse modeling of flanking structures (e.g., Kocher and Mancktelow, 2005). Overall, this approach has proven to allow the unfolding of sediment layers and the restoration of salt structures, both in 2-D and in 3-D. 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 erosion processes 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 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-elasto-plastic behavior in rocks, we show in the paper that this consideration is sufficient in simple setups. In the case of more complex overburdens, the method proved in the literature to be able to restore various structures with power-law viscosities (e.g., Lechmann et al., 2010). We introduce a numerical scheme combining features 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 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 numerical 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 free top surface.

2 Using creeping flow equations for geomechanical restoration

2.1 Creeping flow equations

The standard equations for creeping flows are the Stokes equations, consisting of the momentum conservation equation,

$$\begin{array}{}\text{(1)}& \mathbf{\nabla }\cdot \mathit{\sigma }+\mathit{f}=\mathbf{0},\end{array}$$

and the mass conservation equation for incompressible fluids (continuity equation),

$$\begin{array}{}\text{(2)}& \mathbf{\nabla }\cdot \mathit{v}=\mathrm{0},\end{array}$$

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:

$$\begin{array}{}\text{(3)}& \mathit{\sigma }=\mathit{\tau }-p\mathbf{I},\end{array}$$

where I is the identity tensor. In the viscous flow assumption, the deviatoric part of the stress is

$$\begin{array}{}\text{(4)}& \mathit{\tau }=\mathrm{2}\mathit{\eta }\mathbf{D},\end{array}$$

with η the dynamic viscosity and D the infinitesimal strain rate tensor defined by

$$\begin{array}{}\text{(5)}& \mathbf{D}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathbf{\nabla }\mathit{v}+(\mathbf{\nabla }\mathit{v}{)}^T\right].\end{array}$$

Assembling Eqs. (1), (3), (4) and (5), the momentum conservation equation can be written

$$\begin{array}{}\text{(6)}& \mathbf{\nabla }\cdot \left[\mathit{\eta }(\mathbf{\nabla }\mathit{v}+(\mathbf{\nabla }\mathit{v}{)}^T)\right]-\mathbf{\nabla }p=-\mathit{\rho }\mathit{g}.\end{array}$$

Here, we deal with materials that are highly viscous (with a viscosity η over 10¹⁷ Pa s), over timescales 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 resolution 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, as they advect the viscosity and density fields in time.

2.2 Restoration idea

In forward simulation schemes, the Stokes equations (Eqs. 6 and 2) are solved for pressure and velocity, and the material representation of the geological model is advected from the velocity 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

$$\begin{array}{}\text{(7)}& \mathit{x}(t+\mathrm{\Delta }t)=\mathit{x}\left(t\right)+\mathit{v}\left(t\right)\cdot \mathrm{\Delta }t,\end{array}$$

with x(t) and v(t) the position and the computed velocity of 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 approximation 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 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

$$\begin{array}{}\text{(8)}& \mathit{x}(t-\mathrm{\Delta }t)=\mathit{x}\left(t\right)-\mathit{v}\left(t\right)\cdot \mathrm{\Delta }t\end{array}$$

for the advection of the points of the material model, at each time step, like in Fig. 1.

Figure 1 Example of the restoration scheme for a simple setup (a): as the arrows in panel (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 panel (c).

In this light, using viscous fluid properties instead of elastic properties to represent the mechanical behavior of geological 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

3.1 Presentation

The restoration scheme presented in Sect. 2 has been implemented in the FAIStokes (Finite element Arbitrary Eulerian-Lagrangian Implementation of Stokes) code. It relies almost entirely on the deal.II library (Bangerth et al., 2007; Arndt et al., 2019, 2020) for all finite-element-related algorithms. The material tracking is based on the PIC method (e.g., Asgari and Moresi, 2012; Thielmann et al., 2014; Gassmöller et al., 2018, 2019; Trim et al., 2019). The general workflow of the code is shown in Fig. 2 and details of implementation are discussed in the following subsections. Five benchmarks have been carried out to test the computation parts of the code and are presented in Appendix Sects. A, B, C, D and E.

Figure 2 Schematic workflow of the FAIStokes code structure. The pre-refinement step occurs at the beginning of the simulation (or during a re-initialization of the grid) to ensure that the velocity used for the advection step is computed using the adaptively refined grid.

3.2 Finite element discretization

The finite element method (FEM) was introduced in the late 1950s (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 domain is discretized on a set of quadrilateral elements, on 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₂×Q₁ elements that satisfy the 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 at which structural 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. We use Dirichlet and Neumann boundary conditions that we adapt (e.g., rigidity, 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.

3.3 Material discretization

The geomechanical simulation of a specific domain requires to choose an appropriate kinematic description to follow the displacement inside the geological layers. Continuum mechanics 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 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. It has various formulations and implementations, both in 2-D (e.g., Willett et al., 1993; Poliakov et al., 1996; Massimi et al., 2006, 2007; Fillon et al., 2013; Rose et al., 2017) and, more recently, in 3-D (e.g., Braun, 2003; Thieulot, 2011; Thieulot et al., 2014). Most of these methods 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 a possible free surface. It is particularly useful in geomechanics, where the vertical deformation is generally 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 grid has an ALE part, as it can adapt to follow the movement of the free surface.

3.4 The PIC method

During mechanical simulations, the material properties inside the model are tracked using particles; each of these particles discretizes 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 to solve the Stokes equations on the grid. Following this, the particles are advected using the solution on the grid.

At the beginning 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 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 does not 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 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 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 finite element grid to reasonable accuracy.

3.5 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 ¹ and step-32 ². The grid is created first as a quadrilateral from the coordinates of the bottom left and top right corners of the domain. This quadrilateral is then split in order to get cells 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 preconditioner 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 term on the right-hand side, the norm of the gravity vector g of Eq. (6) is always 9.81 m s⁻² in our simulations, and its direction is always downwards. The matrix system is solved using an iterative Flexible Generalized Minimal Residual (FGMRES) solver preconditioned by a block matrix involving the Schur complement (Kronbichler et al., 2012). This solution is then used to refine and coarsen the grid adaptively using deal.II's features, based on a gradient estimator in order to minimize the local error. 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 and E show the results of benchmarks that tested the computation of the velocity on different setups.

3.6 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 grid cell with a Q₂ interpolation scheme. Depending on whether the simulation is forward or backward, the displacement of each particle for a time step Δt is computed using Eqs. (7) or (8). The value of Δt is computed from the CFL condition. The default value for the CFL number is 0.085, 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 second-order Runge–Kutta scheme in space: at each time step, the particles 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 again by half of this new displacement. This scheme reduces the error in the advection process without need for simulation time 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 and E show the results of benchmarks that tested the interpolation of the velocity in time-dependant problems.

3.7 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 2-D cases). It is advected at 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 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. Figure 3 illustrates the whole process.

Figure 3 Process for the update of the free surface (the motion is exaggerated for the sake of the explanation and is less extreme in reality). (a) Initial state where the velocity is computed on the grid. (b) The point swarm tracking the free surface is advected according to the computed velocity. (c) The grid nodes at the top of the free surface are moved vertically to match the point swarm. (d) The deformation of the grid is diffused to the rest of the nodes.

Since our models are isothermal, no special processing is required to correct the temperature field during 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 (referred to as FSSA in the rest of the paper) developed by Kaus et al. (2010) and showcased in Quinquis et al. (2011) has been implemented in FAIStokes; we benchmark it in Appendix E.

4 Results

In addition to the benchmarks presented in the Appendix, which mainly check the algorithms of the code, we tested our 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 displacements 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.

4.1 Diapiric growth model

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.

Figure 4 Setup of the model scaled up from van Keken et al. (1997).

The upper layer represents sediments that are denser than the lower layer which contains salt (ρ_o=2600 kg m⁻³ for the sediment layer and ρ_s=2150 kg m⁻³ for the salt layer). A sinusoidal instability initiates the movement at the beginning of the simulation. The model is limited to a 10 km × 9.142 km domain (the 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² 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 5 times more particles around the interface between the two layers 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 (${\mathit{\eta }}_{\mathrm{o}}={\mathit{\eta }}_{\mathrm{s}}={\mathrm{10}}^{\mathrm{19}}$ Pa s), and the second one with material properties closer to reality with a lower viscosity for salt (${\mathit{\eta }}_{\mathrm{o}}=\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{19}}$ Pa s for the sediment layer and ${\mathit{\eta }}_{\mathrm{s}}=\mathrm{1.4}\times {\mathrm{10}}^{\mathrm{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 obtained at the end of the simulation. The state obtained after 6×10⁶ years for the first test and 1.5×10⁶ years for the second test, as well as the restored models, are shown in Fig. 5.

Figure 5 Particle swarms for the two synthetic diapiric growth experiments. The darker grey and brown parts on the swarms are due to the higher density of particles around the interfaces. The particles have the same initial position (a) in the two experiments, with different material properties. The result of the forward simulation after 6×10⁶ years for the first experiment is shown in panel (b). Panel (c) shows the result of the forward simulation for the second experiment after 1.5×10⁶ years. The results for the restoration simulations are shown in panels (d) and (e) for the first and second experiments, respectively.

We can see that while the isoviscous experiment has a rather smooth forward result, the second experiment with a less viscous 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 on the left-hand side of the model.

To check the quality of the restoration in the two experiments, 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. The mean value for this distance is 14 m (0.1 % error) for the isoviscous case and 201 m (2 % error) for the variable viscosity case, and the maximum value is 143 m (1.5 % error) for the isoviscous case and 4947 m (49 % error) for the variable viscosity case. While these results are quite good for the isoviscous case, we could think that the variable viscosity case restoration is too inaccurate. Histograms for the errors in the two experiments are given in Figs. 6 and 7 and help explain this phenomenon.

Figure 6 Error logarithmic distribution for the first experiment (isoviscosity) on the diapiric growth model.

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

Figure 7 Error logarithmic distribution for the second experiment (variable viscosity) on the diapiric growth model.

4.2 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 mimics passive diapirism structures created by syn-deformation 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 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 whether the boundary conditions 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 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 results of the restoration, we tested different possibilities. For the density, the value for salt rock is ρ_salt=2160 kg m⁻³, while the value for sediments can vary depending on the type and origin of deposition mechanisms; we considered here a value ${\mathit{\rho }}_{\mathrm{o}}\in \left[\mathrm{2600};\mathrm{3300}\right]$ kg m⁻³. For simplicity, we set the viscosity of the salt layer at η_salt=10¹⁷ Pa s and only vary the viscosity of the sediments ${\mathit{\eta }}_{\mathrm{o}}\in \left[{\mathrm{10}}^{\mathrm{19}};{\mathrm{10}}^{\mathrm{21}}\right]$ Pa s in order to test the effect of the contrast.

Figure 8 Setup of the simulation for the model generated with the method proposed by Clausolles et al. (2019). The initial model is sampled on a regularly spaced particle swarm.

We did five experiments with different values of ρ_o and η_o:

• Exp. 1: ρ_o=2600 kg m⁻³, η_o=10¹⁹ Pa s

• Exp. 2: ρ_o=3300 kg m⁻³, η_o=10¹⁹ Pa s

• Exp. 3: ρ_o=2950 kg m⁻³, η_o=10²⁰ Pa s

• Exp. 4: ρ_o=2600 kg m⁻³, η_o=10²¹ Pa s

• Exp. 5: ρ_o=3300 kg m⁻³, η_o=10²¹ Pa s

As this is a simplification of a real case application, and there is no information on the type of sediments, in each experiment the density and viscosity are homogeneous in the sediment and salt layers.

Figure 9 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.

The results for the five experiment simulations are given in Fig 10. Depending on the experiment, we choose to stop the 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.

Figure 10 Results of the five restoration experiments done on the salt model setup of Fig. 8, after different time spans t_end.

Overall, the restoration process removes the diapir and leaves a salt scar, while the sediment layers remain globally flat. Since this setup is generated by a method for syn-deformation 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. In this case, the stress state inside the model being incorrect, the sediment and salt layers could not 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 the results are still quite convincing despite the high level of simplification, this shows that the salt diapir is a result of upbuilding and not downbuilding. The analysis of the five experiments shows that in this setup, the viscosity contrast between salt and sediment and the density of the sediments 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 4 orders 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 viscosity 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 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 analogue models (Weijermars et al., 1993; Dooley et al., 2005).

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.

4.3 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, as shown in Fig. 11.

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.

It consists of a layered overburden underlain by salt and cut by two 60^∘ faults. As the intent of the paper is to focus on the restoration of structural models, we do not consider the formation of the faults with plasticity but rather 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 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. 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 6.25 m × 6.25 m cells near the faults), as shown in Fig. 12 for the first time step.

Figure 13 Results of the simulations for the simplified graben model at (a) the end of the forward simulation and (b) the end of the restoration simulation.

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 different density and viscosity. Material properties inside the model are as follows:

• Overburden type 1 layer: ${\mathit{\eta }}_{\mathrm{o}\mathrm{1}}=\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{19}}$ Pa s, ρ_o1=2550 kg m⁻³.

• Overburden type 2 layer: ${\mathit{\eta }}_{\mathrm{o}\mathrm{2}}=\mathrm{5.0}\times {\mathrm{10}}^{\mathrm{19}}$ Pa s, ρ_o2=2600 kg m⁻³.

• Salt layer: ${\mathit{\eta }}_s=\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{17}}$ Pa s, ρ_s=2160 kg m⁻³.

• Faults: ${\mathit{\eta }}_{\mathrm{f}}=\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{16}}$ Pa s, ρ_f=2200 kg m⁻³.

Like in the diapiric growth model, we first do a forward simulation and then apply the restoration scheme on the model obtained. The lateral flow is set to 10 mm yr⁻¹ 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 8 times more 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 inward 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 simulations is shown in Fig. 13.

Figure 14 Error on the restoration of the interfaces in the model. The lower the number of the interface, the lower the interface it corresponds to. For example, interface 0 is the interface between the salt and the lowest sediment layer.

Figure 14 shows the difference between the position of the interfaces before the forward simulation and at the end of the backward simulation. The numbering of the interfaces follows their position in the model, with the interface 0 being the lowest salt–sediment interface.

The mean and maximum errors for each interface are given in Table 1.

Table 1 Mean and maximum errors for the restoration of each interface.

Overall, the results for the restoration simulation are quite good, with mean errors around 1 % of the forward deformation for the layer interfaces. Figure 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 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 occurring on the top of the faults, where high viscosity contrasts occur. The error resulting from these instabilities, however, is quite small compared to the amount of deformation in the model (around 200 m of slip on the faults).

5 Discussion

While the results of the three test models in the previous section are promising, their purpose is not to correctly compute 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 unphysical 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 et al., 2001; Lechmann et al., 2010). As such, restoration using 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 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 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, 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 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 timescales are too short to reflect the slow movement occurring at geological timescales. The values we took are inspired from numerical simulations, but they have a large uncertainty (at least 1 order of magnitude) (e.g., Massimi et al., 2006; Kronbichler et al., 2012), as they are calibrated using postglacial rebound data, for example. Working on analogue 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 phenomena occur could guide the choice of viscosity values in subsurface models.

Most geomechanical restoration schemes consider basin rocks to have an elastic behavior, whereas the rheology used herein does not display Poisson effects. Incorporating elasticity in viscous flow has been done 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), for example. The problem is that these schemes, like every implementation of elasticity, use values of the 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. For example, the incompressibility constraint, which implies a Poisson ratio of 0.5, can be relaxed, which could be used to account for lesser values of the Poisson ratio. Regarding the rheology of faults, we cannot directly use their usual forward modeling implementation considering the rock as having a plastic behavior. Indeed, the previous stress history is needed to simulate such 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 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 geomechanical restoration schemes that consider faults as free-slip surfaces.

A significant issue with the boundary conditions in geomechanical simulations is the difficulty to estimate the paleo-forces at play several kilometers underground. We therefore 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 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 (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 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., 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 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 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 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 models shown in Sect. 4.1 and 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 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 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, 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 it performs similarly to the free surface on the benchmarks presented in Appendices D and E. In restoration simulations, it can delay the instabilities that appear in models driven exclusively by gravity but does not remove them.

6 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 allowed various tests on its implementation. 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 including 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 results 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 viscosity and density of geological layers and to upgrade the specific implementation of faults.

Appendix A: Taking into account small scales inside a model: the Rayleigh–Taylor instability benchmark

This benchmark is based on the analytical solution of a 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 boundary condition, while the sides have a free-slip boundary condition.

The first layer, made of fluid 1 with properties (ρ₁,η₁), overlays the second layer, made of fluid 2 (ρ₂,η₂). An initial sinusoidal disturbance of the interface between the two layers is introduced, characterized by an amplitude Δ and a wavelength λ, as shown in Fig. A1.

Figure A1 Rayleigh–Taylor instability benchmark initial setup.

Under these conditions, the velocity of the diapiric growth v is given by Ramberg (1968):

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{v}{\mathrm{\Delta }}}=-K{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{1}}-{\mathit{\rho }}_{\mathrm{2}}}{\mathrm{2}{\mathit{\eta }}_{\mathrm{2}}}}{h}_{\mathrm{2}}g,\end{array}$$

with K, the dimensionless growth factor, given by

$$\begin{array}{}\text{(A2)}& K={\displaystyle \frac{-d_{\mathrm{12}}}{c_{\mathrm{11}}{j}_{\mathrm{22}}-d_{\mathrm{12}}{i}_{\mathrm{21}}}},\end{array}$$

which involves the following factors:

$$\begin{array}{rl}& {\mathit{\varphi }}_{\mathrm{1}}={\displaystyle \frac{\mathrm{2}\mathit{\pi }{h}_{\mathrm{1}}}{\mathit{\lambda }}}\\ & {\mathit{\varphi }}_{\mathrm{2}}={\displaystyle \frac{\mathrm{2}\mathit{\pi }{h}_{\mathrm{2}}}{\mathit{\lambda }}}\\ & c_{\mathrm{11}}={\displaystyle \frac{\mathrm{2}{\mathit{\eta }}_{\mathrm{1}}{\mathit{\varphi }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\eta }}_{\mathrm{2}}\left(\cosh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}\right)-\mathrm{1}-\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}^{\mathrm{2}}\right)}}-{\displaystyle \frac{\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}^{\mathrm{2}}}{\cosh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}\right)-\mathrm{1}-\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}^{\mathrm{2}}}}\\ & d_{\mathrm{12}}={\displaystyle \frac{{\mathit{\eta }}_{\mathrm{1}}\left(\sinh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}\right)-\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}\right)}{{\mathit{\eta }}_{\mathrm{2}}\left(\cosh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}\right)-\mathrm{1}-\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}^{\mathrm{2}}\right)}}+{\displaystyle \frac{\sinh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}\right)-\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}}{\cosh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}\right)-\mathrm{1}-\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}^{\mathrm{2}}}}\end{array}$$

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& i_{\mathrm{21}}={\displaystyle \frac{{\mathit{\eta }}_{\mathrm{1}}{\mathit{\varphi }}_{\mathrm{2}}\left(\sinh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}\right)+\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}\right)}{{\mathit{\eta }}_{\mathrm{2}}\left(\cosh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}\right)-\mathrm{1}-\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}^{\mathrm{2}}\right)}}+{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{2}}\left(\sinh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}\right)+\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}\right)}{\cosh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}\right)-\mathrm{1}-\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}^{\mathrm{2}}}}\\ & j_{\mathrm{22}}={\displaystyle \frac{\mathrm{2}{\mathit{\eta }}_{\mathrm{1}}{\mathit{\varphi }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{2}}}{{\mathit{\eta }}_{\mathrm{2}}\left(\cosh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}\right)-\mathrm{1}-\mathrm{2}{\mathit{\varphi }}_{\mathrm{1}}^{\mathrm{2}}\right)}}-{\displaystyle \frac{\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}^{\mathrm{3}}}{\cosh \left(\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}\right)-\mathrm{1}-\mathrm{2}{\mathit{\varphi }}_{\mathrm{2}}^{\mathrm{2}}}}.\end{array}\end{array}$$

We set ρ₁=3300 kg m⁻³, ρ₂=3000 kg m⁻³, η₁=10²¹ Pa s, $L_x=h_{\mathrm{1}}+h_{\mathrm{2}}=\mathrm{512}$ km, and Δ=3 km. We make η₂ vary between 1.25×10²⁰ and 2.5×10²³ Pa s, while λ takes three values: $L_x/\mathrm{2},L_x/\mathrm{4}$ and L_x∕8 (Fig. A2).

Figure A2 Initial setup of the Rayleigh–Taylor instability benchmark with three different wavelengths: (a) $\mathit{\lambda }=L_x/\mathrm{2}$, (b) $\mathit{\lambda }=L_x/\mathrm{4}$ and (c) $\mathit{\lambda }=L_x/\mathrm{8}$.

A first run is done, where the FEM grid is fixed to 80×80 elements, each containing 10² regularly spaced particles. In 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² regularly spaced particles, and then the swarm is densified to 64 times more particles around the interface between the two fluids. The results are shown along with the analytical ones in Fig. A3.

Figure A3 Comparison between numerical and analytical results for the Rayleigh–Taylor instability benchmark. The numerical results are computed for a 80×80 element grid and for the same grid with three levels of adaptive refinement.

Overall, results show good agreement between the computed solution and the reference, especially in the case of adaptive refinement, where the relative error falls beneath 2.5% for all the curves. Since ϕ₁ 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 in viscosity.

This benchmark ensures the validity of the code in the presence of large viscosity contrasts, even if those contrasts are located on deformations that are small compared to the size of the model. It also validates the averaging of the density and viscosity from the particles to the finite element grid.

Appendix B: Taking into account viscosity changes: the falling block benchmark

This benchmark appears in Gerya (2019) and is presented in Thieulot (2011). It consists in modeling the fall of a block of fluid of properties (ρ₁,η₁) inside another fluid of properties (ρ₂,η₂), with ρ₁>ρ₂. The domain is a square of size $L_x=L_y=\mathrm{500}$ km, and the block (a square in 2-D) of size 100×100 km is initially centered at point (x=250 km, y=400 km), as shown in Fig. B1.

Figure B1 Falling block benchmark initial setup.

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 5² particles in each initial element and then densifying it up to 64 times more particles around the zone of interest (i.e., the falling block). Free-slip boundary conditions are imposed on all sides of the domain. We carry out five experiments:

• Exp. 1: η₂=10²⁰ Pa s, ρ₁=3220 kg m⁻³;

• Exp. 2: η₂=10²¹ Pa s, ρ₁=3300 kg m⁻³;

• Exp. 3: η₂=10²² Pa s, ρ₁=6600 kg m⁻³;

• Exp. 4: η₂=10²³ Pa s, ρ₁=3300 kg m⁻³;

• Exp. 5: η₂=10²⁴ Pa s, ρ₁=9900 kg m⁻³.

In all the experiments, the density of the surrounding fluid is ρ₂=3200 kg m⁻³ and the viscosity of the block is varied between 10¹⁹ and 5×10²⁷ Pa s. The velocity of the falling block is measured at its center at t=0 for all experiments. Following physical intuition, one expects the velocity of the block to act as follows: (1) decrease when the viscosity of the surrounding fluid η₂ increases (i.e., when going from Exp. 1 to Exp. 5) and (2) increase with the density contrast (ρ₁−ρ₂) in each experiment. To check this behavior, we measure $v{\mathit{\eta }}_{\mathrm{2}}/({\mathit{\rho }}_{\mathrm{1}}-{\mathit{\rho }}_{\mathrm{2}})$ as a function of the viscosity contrast log ₁₀(η₂∕η₁). The results of the benchmark are plotted in Fig. B2.

Figure B2 Velocity measurements as a function of the viscosity contrast between surrounding medium and block for the experiments of the falling block benchmark.

We can see that the experimental points line up on a single curve, which shows that FAIStokes can deal with gravity-driven simulations where 0.6$\mathit{\%}\le ({\mathit{\rho }}_{\mathrm{1}}-{\mathit{\rho }}_{\mathrm{2}})/{\mathit{\rho }}_{\mathrm{2}}\le \mathrm{210}$ %, and the viscosity contrasts are as strong as ${\mathrm{10}}^{-\mathrm{6}}\le {\mathit{\eta }}_{\mathrm{2}}/{\mathit{\eta }}_{\mathrm{1}}\le {\mathrm{10}}^{\mathrm{5}}$ in a consistent manner.

Appendix C: Advecting particles: the rotation benchmark

The last benchmark aims at assessing the error in the advection part only.

Figure C1 Setup for the rotation benchmark, assessing errors on the advection of particles.

The setup of the model is a square of size 10×10 km, where we study the advection of a single particle, starting at coordinates (8, 5 km) and doing a 2π rotation around the center point (5, 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 and is always normal to the line connecting the point to the model center:

$$\begin{array}{}\text{(C1)}& \mathit{v}=v_{\mathrm{0}}\cdot {\mathit{e}}_{\mathit{\theta }}=\left(\begin{array}{c}{v}_{\mathrm{0}}\cdot \sin \mathit{\theta }\\ v_{\mathrm{0}}\cdot \cos \mathit{\theta }\end{array}\right).\end{array}$$

The grid is not adaptively refined here and is composed of 16×16 elements. In order to have scales that are geologically relevant, we choose v₀=3 cm yr⁻¹ and vary the time step Δt between 500 and 2000 years (in this setup, the CFL numbers chosen for our simulations would give a time step between 175 years for the lowest CFL number and 1753 years for the highest CFL number). The second-order Runge–Kutta scheme presented in Sect. 3.6 is used at all time steps. We then evaluate the distance $\mathrm{\Delta }r=\left|r\right(\mathit{\theta }=\mathrm{0})-r(\mathit{\theta }=\mathrm{2}\mathit{\pi }\left)\right|$. This distance gives us a measure of the error made in the computation of the particle advection and allows us to compare different advection schemes.

Figure C2 Results for the rotation benchmark obtained with different time steps and advection schemes.

Figure C2 shows the results obtained for a 2π rotation of the particle with different interpolation schemes. We can see that reducing the time step linearly reduces the error on the radius r(θ). In this setup, the type of 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 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 lithosphere 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:

$$\begin{array}{}\text{(D1)}& h_{\mathrm{analytical}}\left(t\right)=h_{\mathrm{initial}}\exp (-\mathit{\gamma }t),\end{array}$$

where γ is the relaxation rate and h_initial is the value of h at the beginning of the simulation. The model setup for the benchmark is shown in Fig. D1.

Figure D1 Model setup for the 2-D free surface benchmark.

The bounding box of the model spans 2800 × 707 km. The underlaying mantle layer is 600 km thick, while the lithosphere has a thickness between 93 and 107 km. The 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 ${\mathit{\rho }}_{\mathrm{M}}={\mathit{\rho }}_{\mathrm{L}}=\mathrm{3300}$ kg m⁻³ and a viscosity of η_M=10²¹ Pa s and η_L=10²³ Pa s, respectively. We set free-slip boundary conditions 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 three times. The particle swarm contains 484 160 particles; it is constructed by first sampling regularly the domain and then adaptively densifying the swarm 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 ($\mathit{\gamma }=\mathrm{0.2139}\times {\mathrm{10}}^{-\mathrm{11}}$ s⁻¹) and a characteristic relaxation time ($t_{\text{relax}}=\mathrm{14.825}\times {\mathrm{10}}^{\mathrm{3}}$ years). The results obtained with FAIStokes are given in Fig. D2.

Figure D2 Comparison between the analytical and numerical results of the maximum topography over time for the free surface benchmark.

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. It also gives another evaluation of the handling of gravity-driven flow, this time with the addition of evolution inside the model.

Appendix E: Upgrading the free surface movement: the sloshing benchmark

This benchmark is presented in Kaus et al. (2010), where it is used to assess the results of the 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 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).

Figure E1 Sloshing free surface benchmark initial setup.

The model span is 500 km × 500 km; the side boundaries have a free-slip condition, the lower boundary has a no-slip condition, 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 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 referred 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 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).

Figure E2 Simulation evolution for Δt=5000 years, showing the sloshing instability: the velocity pattern changes from one time step to the other, with the velocity norm increasing each time.

In order to follow the evolution of the free surface, we keep trace of the altitude of the most top-left point over time. Results of a 0.5 Myr simulation, for different time steps Δt, with and without the FSSA, are shown in Fig. E3.

Figure E3 Altitude of the most top-left point of the grid over time, for the sloshing free surface benchmark, for different time steps with and without the FSSA.

We can see that the implementation of the FSSA stabilizes the sloshing behavior of the free surface that appeared with a time step Δt of 5000 years and keeps the free surface stable 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 surface, this time with the additional resolution of an instability that can occur with free surfaces.