Validating novel boundary conditions for three-dimensional mechanics-based restoration: An extensional sandbox model example

Geomechanical restoration methods are dependent on boundary conditions to ensure geological consistency of the restored model in terms of geometry and strain. Classical restoration boundary conditions, such as attening a datum horizon, may lead to inconsistent displacement and strain elds. We restore a laboratory structural sandbox model with known deformation history in order to develop guidelines for de nition of boundary conditions that produce improved results from geomechanical restorations. The sandbox model has a basal silicone layer, includes syn-kinematic deposition, and is characterized by structures analogous to those found in supra-salt extensional environments. The deformed geometry is interpreted from 3D tomography imaging, and a time-series of cross-section tomography images provides a benchmark to quantify restoration error and inform boundary conditions. We con rm that imposing a lateral displacement equal and opposite to fareld tectonic shortening or extension provides a more accurate restoration. However, the amount of displacement may not be known in real cases. We therefore test several established methods, using only the unrestored geometries, to assess the amount of shortening that should be used to guide geomechanical restorations. An accurate estimation is provided by the area-depth method and potentially by a dilatation analysis. Additionally, novel fault compliance boundary conditions produce improved results in the vicinity of crossing and branching faults. Application of similar methods should produce improved restoration of natural geologic structures.

This uncertainty in appropriate boundary conditions is particularly problematic if one intends to analyze the corresponding stress or strain for fracture analysis or other purposes (Maerten and Maerten, 2006;Mejía-Herrera et al., 2014;Stockmeyer et al., 2017). The problems illustrated by Lovely et al. (2012) are: (1) there may be instances when the classical boundary conditions, as dened above, may be unphysical and (2) there is no specic guideline to choose appropriate boundary conditions in restoration. Durand-Riard (2010); Lovely et al. (2012) and Durand-Riard et al. (2013b) suggest that these classical boundary conditions may be insucient to restore geologically consistent and physical strain. They show on synthetic models that a lateral displacement boundary condition along a boundary wall is necessary to recover the expected strain in compressive, extensional, and strike-slip and oblique-slip contexts. In addition, they show that the restoration displacement eld is only consistent with the forward displacement eld when the amount of the displacement condition on a wall is equal to the amount of forward displacement. These works highlight the need for additional constraints derived from geologic or tectonic insights for mechanics-based restoration. The main challenge for dening these additional constraints is that they require knowledge of the deformation history, which is rarely accessible and, ideally, should be an output of the mechanics-based restoration. Moreover, these studies of boundary conditions were applied to numerical or synthetic models, which are typically idealizations of natural geologic structures, and present additional uncertainties and assumptions (structural interpretation, deformation path, etc.).
Models from laboratories are often used to validate restoration methods (e.g., Schultz-Ela, 1992;Yamada and McClay, 2003;Maerten and Maerten, 2006;Groshong et al., 2012;Moretti and Callot, 2012). As they are laboratory experiments, the forward boundary conditions and the mechanical behaviors are known, and the kinematic evolution of structures may be recorded. Moreover, the interpretation uncertainties of these deformed forward models are generally small, such that the applied boundary conditions can be considered the primary source of uncertainty in restoration attempts. This is a signicant benet compared to restoration attempts of natural structures, where restoration uncertainties result from the interplay of boundary condition uncertainties and of structural interpretation uncertainties (Gratier et al., 1991;Schultz-Ela, 1992).
In this study, we performed a sequential restoration on an analog model deformed in laboratory. Computed tomography (CT) images capture the sequential development of this analog model. These images constrain the forward deformation path of each structure, capturing paleo-geometries through time on one edge of the analog model, and, thus, provide a reasonable benchmark for restoration quality and boundary condition testing. In the following, we describe the structural sandbox model and our tests of boundary conditions with the goal of restoring deformed geometries and related fault slip that are consistent with the reference model paleogeometries. New fault compliance boundary conditions are proposed to handle the complex fault network identied on the CT images. In addition, we propose methods to dene lateral displacement boundary conditions without detailed knowledge of the forward deformation path, improving the viability of the 3D geomechanical restoration method for use with natural geologic structures. Geologists typically have access only to the present-day state of deformed strata, often times informed by sparse and uncertain data. Such data may be consistent with multiple interpretations that may vary signicantly (e.g., Frodeman, 1995;Bond et al., 2007;Wellmann et al., 2010;Bond, 2015;Cherpeau and Caumon, 2015). Thus, the analysis of rock deformation through time is made dicult by the lack of direct information on paleo-structures and the limitations of the available data. To overcome some of these concerns, laboratory analog models are widely used to model viable deformation styles and paths of natural geologic structures (e.g., McClay, 1990).
They provide a way to follow the evolution of well-known geometries under known physical mechanisms through time. X-ray tomography is a common, non-destructive method to image the 3D structures of a deformed structural sandbox (e.g., Colletta et al., 1991;. X-ray tomography resolution on an analog model is generally sucient to study deformed structures with minimal geometric uncertainties. Distinct horizons and faults can be observed due to density contrasts in the model's stratigraphy. Thus, as pointed by Colletta et al. (1991), X-ray tomography is a valuable tool to analyze the temporal evolution of laboratory models. Moreover, an analog model must be properly dened to reproduce the behavior of geological structures. Scale, mechanical materials, physical processes, and timing are examples of parameters to consider for the purpose of assessing the degree to which analog models represent natural structures.
In this paper, we restore a laboratory model analogous to supra-salt extensional structures observed in salt basins around the world, such as the Gulf of Mexico, Angola and Morocco. It is well established that dry sand (no cohesion) is a viable material for modeling brittle and ductile rock deformation in sedimentary systems (e.g., Panien et al., 2006;Victor and Moretti, 2006;Dooley et al., 2007;Moretti and Callot, 2012;Darnault et al., 2016). Moreover, Weijermars et al. (1993); Victor and Moretti (2006) and Moretti and Callot (2012), among others, have shown that an analog composed by a stack of sand above silicone can produce structures representative of natural salt basins. Silicone has a very weak rheology relative to sand. Thus, sand layers deform and may penetrate into the silicone. This eectively reproduces the subsurface at the interface between a viscous salt layer and overlying brittle rocks (Weijermars et al., 1993).
Our work is based on a deformed structural sandbox done in laboratory by IFPEN (http://www.ifpe nergiesnouvelles.fr) and C&C Reservoirs, 2016, DAKS TM -Digital Analogs Knowledge System (http: //www.ccreservoirs.com) ( Figure 1) to reproduce extensional salt structures. The model box was initially composed of two horizontal layers composing a pregrowth stratigraphy: one of silicone at the bottom with a thickness of 1.8 cm (0.71 in), and one of sand above with a thickness of 4 mm (0.2 in). The initial thickness (along Z direction) of the pregrowth strata is 2.2 cm (0.87 in). Along the Y axis the structural sandbox length is 10 cm (3.9 in), and 18 cm (7.1 in) along the X axis. The model box was inclined by 1.5 • (Weijermars et al., 1993;Victor and Moretti, 2006). Deformation was induced initially by gravity sliding along this tilt (toward the eastern side on Figure 1). On the down-dip end of the model there was no wall to restrain the motion of the materials. As the model deformed, alternating layers of pyrex or sand were deposited (one layer every 16 min in mean), further driving deformation by a combination of gravity spreading and gravity gliding (Victor and Moretti, 2006). This deposition of successive stratigraphic horizons during deformation represents syn-tectonic strata (i.e., growth strata). As pyrex and sand strata are deposited above the silicone, this experiment describes supra-salt structures. At each depositional time step, the newly deposited sediments lled the available model space. In total, 12 layers were deposited ( Figure 1) over the course of the forward analog model. The total duration of the experiment is 4h16min. The properties of the silicone, pyrex and sand are provided in Table 1.
The scaling from the analog model scale to real eld scale is approximately 1 cm (0.4 in) for 1 km (0.6 mi), consistently with similar analog models (Ellis and McClay, 1988;Dooley et al., 2007;Wu et al., 2009;Hidayah, 2010;Darnault et al., 2016). See Hubbert (1937) and Ramberg (1981) for more details about the methods used to dene this scaling. It is possible to distinguish silicone, sand, pyrex and the fault osets by tomographic imaging due to their density contrasts. Indeed, faults are visible in the analog model due to sand and pyrex dilatation (areas of lower density) (Colletta et al., 1991;Cobbold and Castro, 1999;Le Guerroué and Cobbold, 2006;Groshong et al., 2012). Moreover, sand and pyrex have a sucient Hounseld density contrast (Table 1) to distinguish them in X-ray tomography (Panien et al., 2006;Darnault et al., 2016), allowing the analysis of fault oset. In addition, we assume that friction is negligible on the edges of the structural sandbox model (Souloumiac et al., 2012).

Interpretation of the structural sandbox model
A 3D X-ray tomography volume of the nal state of the deformed box was produced. The volume is dened by an X-ray tomography section every 3 mm (0.1 in) along the Y axis, and an X-ray tomography section every 0.6 mm (0.02 in) along the X and Z axes, producing a tomography volume that can be interpreted with similar methods as a 3D seismic reection survey ( Figure 2A). Unfortunately, we had access to only a part of the structural sandbox volume. Indeed CT imaging only recorded through time a specic interval of the structural sandbox. Beyond this interval, down-dip, the analog model continued to deform but was out of the scope of the CT imaging. Although we analyzed the majority of the volume, a part on the eastern side could not be considered in our study. From our interpretation, we built a 3D explicit numerical structural model (boundary representation, Figure 2B) using SKUA-GOCAD (Paradigm, 2015), in which horizon and fault surfaces are conformal (Caumon et al., 2004). The geological model is composed of three primary structures: two grabens in the western and central regions of the model, and a series of west-dipping half-grabens in the eastern region of the model. The layer of silicone representing autochthonous salt is not explicitly represented in our numerical model. Within the analog model, we identied 52 faults. As the purpose of the modeling is the restoration and the strain analysis, faults with small osets were ignored for simplicity ( Figure 3B). Additionally, a few faults that dene narrow fault blocks were removed and the horizons were made continuous ( Figure 3C

Structural uncertainties
Although the analog structures are well imaged, uncertainties in our structural interpretations exist. This is largely the result of approximating diuse horizons and faults in the analog model with discrete surfaces in our structural representation. Thus, quantifying our interpretation precision is necessary in order to properly evaluate the quality of subsequent restoration results. Figure 4 illustrates four examples of interpretation uncertainties of the deformed analog model. Boundaries between white layers (sand) and black layers (pyrex or silicone) are typically blurred gray ( Figure 4). The thickness of these gray transition zones provides an estimate of the uncertainty associated with an interpreted horizon between two strata intervals. Although this thickness may vary laterally ( Figure 4), we estimate an error of 0.1 cm (0.04 in) is a representative uncertainty for all of our interpretations. We will be mindful of this precision as we analyze our restoration results. We also note that other restoration uncertainties, such as nite element approximation or mechanical simplications, although present, are not considered in our uncertainty analysis.
2 Restoration settings 2.1 Physical volumetric model We created a 3D mesh ( Figure 5A) from the structural model ( Figure 2B) using the Geogram (Lévy, 2015), RINGMesh (Botella et al., 2016;Pellerin et al., 2017), VorteXLib (Botella, 2016a,b) and TetGen (Si, 2015a,b) libraries. It is composed of 647,558 tetrahedra and the average tetrahedron length is 0.25 cm (0.098 in). We made some eorts to reduce the number of tetrahedra, which impacts the restoration computational time, and to avoid imprecision due to a coarse mesh. The VorteXLib library enabled us to develop a 3D mesh maximizing the quality of the tetrahedra (equilaterality) to avoid numerical issues during restorations (Parthasarathy et al., 1994;Shewchuk, 2002;Munson, 2007). The silicone layer was not represented in the model used for the restoration for two primary reasons. First, its rheology is far weaker than the sand and the pyrex, and thus is not considered to contribute any signicant resistance. Second, the viscous behavior of the silicone interval cannot be properly represented by the elastic constitutive law invoked in our restoration method. Thus, we focused on restoration of the sand and pyrex layers that overlie the silicone, with the base of the model being a free surface that represents the top of the silicone (Stockmeyer and Guzofski, 2014). As sand and pyrex are rheologically similar (Panien et al., 2006), we applied homogeneous elastic properties for the entire model: Young's modulus was set to 70 GPa (10 7 psi) and Poisson's ratio to 0.2 (Holtzman et al., 2009).

Classical boundary conditions
A video of one edge of the analog model was recorded by X-ray computed tomography, allowing us to visualize the deformation and model geometries through time. The deformation front was located on the eastern side of the model. In contrast, the western side was only weakly deformed ( Figure 2A). Thus, we xed the western wall in the X and Y directions during the restoration, allowing it to only move vertically ( Figure 5B). As the experiment is inside a box, no ow occurred in the north-south direction (Y axis) through the northern and southern walls. Therefore, during the restoration, we xed the northern and the southern walls in Y ( Figure 5C), as recommended by Durand-Riard (2010) and Durand-Riard et al. (2013b) in other deformation contexts. At each restoration step, we set a datum boundary condition for the uppermost stratigraphic surface because we know the original depositional gradient ( Figure 5D). As the model had a tilt of 1.5 • toward east, we rotated the entire model before each restoration. This allowed us to set our datuming boundary condition to a constant Z value ( Figure 5D). We rotated our restored models to their proper geometries for proper comparisons between the numerical models and the CT images. The basal horizon, which denes the interface between sand and silicone, i.e., H1, was dened as a free surface (e.g., Stockmeyer and Guzofski, 2014).
We also dened fault contact conditions to tie the hanging wall and footwall cuto lines of the uppermost (attened) horizon ( Figure 5E) and to avoid any gap or penetration along fault surfaces ( Figure 5F). We ensured fault compliance by contact mechanics (Wriggers and Laursen, 2006) which is a master-slave approach adapted to restoration purposes by Muron (2005), and Maerten and Maerten (2006). This method enables us to tie fault blocks without any friction along fault planes (Muron, 2005;Maerten and Maerten, 2006;Wriggers and Laursen, 2006). The slave surface cannot penetrate nor have a gap with the master surface but the contrary is possible when faults are curved, owing to limited mesh resolution. No relative displacement constraint between the master and the slave is dened: motion is bilateral and is a consequence of both energy minimization and the constraint for the two sides of the fault to be in contact. In the case of the contact of the fault cuto lines of the restored (uppermost) horizon ( Figure 5E), the throw is already dened by the datuming condition applied on the Z-component of the uppermost horizon ( Figure 5D). The heave is dened by contact mechanics as explained above.

2.3
Non-classical boundary condition: imposed shortening condition As extension clearly occurred during forward deformation, we test an optional lateral shortening condition applied to the down-dip model boundary during restoration. We consider this boundary condition analogous to those suggested by Durand-Riard (2010); Lovely et al. (2012) and Durand-Riard et al. (2013b). This optional shortening condition is limited to motion along the X axis (i.e., in the west-east direction) and is applied to the eastern wall ( Figure 5G). In this paper, when we refer to a no shortening condition, we refer to a restoration scenario without this shortening condition set to the down-dip wall. In these scenarios, the eastern wall is free to move along any direction and the resultant shortening is the output of the restoration, ultimately controlled by the datum and fault slip conditions.

Non-classical boundary conditions: contacts between faults
The complexity of a model increases substantially with the numbers of faults due to the increasing number of the interactions between them and with the horizons (i.e., cuto relationships, (Pellerin et al., 2015)). The of such a complex fault network is a dicult task during structural modeling, but also during each step of a sequential restoration. To accomplish this task, we present in this paper two additional fault contact boundary conditions that we applied in our restorations. These conditions use contact mechanics, as classical fault contact boundary conditions that we previously dened (Wriggers and Laursen, 2006).

Handling branching faults
In our structural modeling procedure, a fault is represented by two surfaces, one for the hanging wall and one for the footwall. This and the fault contact conditions enable the sliding of the fault blocks along the faults ( Figure 6). As a result, a branch line between two faults is represented on the main fault by two surface internal borders (black dots in (B1) and (C1) in Figure 7 with F1the main fault). In the case of branching contacts between faults, discontinuities may occur in the restored state if care is not taken ( Figure 7B). Therefore, we set contact conditions to tie the internal surface borders and thus avoid internal gaps or overlaps in the restored state ( Figure 7C).

Handling oset fault surfaces
There were several situations where a fault surface was oset by a dierent fault. To properly characterize these faulted faults, we split each oset fault into two or more distinct fault surfaces. For example, Figure 8 shows that fault F1 is cut and displaced by fault F2. In this case, F1 was represented by two independent faults: F1-hw and F1-fw where the labels -hw and -fw respectively refer to the hanging wall and footwall sides of F2. As a result, F1-hw and F1-fw are able to move independently. However, as F1 was originally a single fault surface, and since we know all of the faults are normal faults at all times during the forward model, we know that the slip between F1-hw and F1-fw along F2 should decrease through the restoration until it becomes null. In other words, the distance between F1-hw and F1-fw should decrease along F2 until they merge. Upon removing of all the fault osets, F1-hw and F1-fw should no longer behave independently, but form a single continuous fault surface. We ensure this condition by a set of contact conditions that aims to tie the dierent connected components of an oset fault (Figure 8). In our case, we were able to quickly determine which restoration steps to apply this contact condition for a particular oset fault (presenting apparent continuity) by investigating the CT images that recorded the forward deformation process.
3 Results: restoration of the analog model 3.1 Sequential restoration We performed a partial sequential restoration using RINGMecha , a mechanicsbased restoration library based on the work of Muron (2005) and Durand-Riard (2010). We used a timeindependent nite element solver to perform the restoration (e.g., Zienkiewicz and Taylor, 2000a,b;Belytschko et al., 2013) with a small deformation assumption. After each restoration step, we removed the uppermost, restored layer before performing the subsequent restoration step. Using the classical and newly dened boundary conditions described above, we performed four steps of sequential restoration for our model (Figures 9-12), yielding a restoration of more than half of the growth strata interval. Restorations with a shortening boundary condition are in Figures 9E, 10E, 11E, and 12E. As we had the CT images of the paleo-states of the northern wall, Qualitatively, the general consistency on the northern edge between the restorations with prescribed shortening and the reference CT images is quite good, indicating a robust and accurate restoration. Restorations without a shortening boundary condition are shown in Figures 9F, 10F, 11F, and 12F. In these models, we only avoided the shortening boundary condition for the last restoration step. For example, the result shown in Figure 9F was not used as the starting model for the restoration in Figure 10F. The starting model for Figure 10F was generated by removing the restored, uppermost layer from the model shown in Figure 9E. In this way, we attempt to avoid propagating errors. For each restoration step that does not include the shortening boundary condition ( Figures 9F, 10F, 11F, and 12F), it is clear that there was not enough extension restored to be considered an acceptable restoration result. In each case, the restored faults are too far down-dip relative to the reference position obtained from the CT tomography video. In contrast, the restorations that included the shortening boundary condition provide a better qualitative match between the restored models and the reference CT images Step 1 Step 2 Step 3 Step 4

Validation: quantitative comparison with a reference solution
The visual comparison of the restored models with the references provides valuable insight on the quality of our restorations. However, quantitative analysis is necessary to rigorously and objectively assess the restoration quality Vidal-Royo, 2015, 2016) and uncertainties. To quantify the dierence between the restored models and the reference paleo-geometries, we measured the magnitude of residual dip slip along faults after a restoration step. Dip slip provides a quantitative measure of the recovered strain along each fault. For each fault that crosses the northern edge, we measured the amount of oset for each horizon on the CT pictures. In addition, we computed the dip slip values on each fault in the numerical models at the considered time steps.
The dierence of dip slip values between the restored model and the reference, that we call delta, is calculated with DS res and DS ref respectively the dip slip in the restored state and the dip slip on the reference CT image. The corresponding distributions for the rst four restoration steps (with shortening boundary conditions) are shown in Figure 13. To avoid bias in this analysis, these distributions do not include the dip slip measures at the uppermost horizon, as these dip slip values are dened by input boundary conditions ( Figure 5E). Table 2 presents the mean and median values for each dip slip delta distribution shown in Figure 13, as well as the percentage of dip slip deltas within the picking uncertainty range estimated to be between -0.2 cm (-0.08 in) and +0.2 cm (+0.08 in). This uncertainty value originates from the picking uncertainty (0.1 cm) applied on the footwall and the hanging wall. It also considers similar uncertainties in our interpretations of paleo-geometries and dip slip magnitudes on the CT images. The distributions show maxima near zero delta. In addition, the majority of the residual dip slip measurements are within the uncertainty range considered (Table 2). Therefore, we suggest these restorations, which each included the applied shortening boundary condition to the down-dip model wall, are valid. Nevertheless, some slip measurements from restoration models dier signicantly from the reference solution. These are clearly not a common result, except in the fourth restoration step (white intervals in Figure 13, Table 2, see discussions).

Estimation of shortening
As shown previously, the amount of extension that is restored without the applied boundary condition to the down-dip model wall consistently underestimates the actual amount of extension that occurred in the forward model. For natural structures, the total amount of extension (or shortening) that occurred to yield the presentday geometry is generally unknown. In these cases, an estimation of the amount of displacement can be attempted using 2D kinematic restoration approaches (e.g., Chamberlin, 1910;Dahlstrom, 1969) or 2D areadepth analysis (Epard and Groshong, 1993;Groshong et al., 2003;Groshong, 2006;Groshong et al., 2012).  Table 2: For each restoration step, the mean and the median of the dip slip delta distribution ( Figure 13) Table 3 presents the incremental shortening evaluated by dierent methods, in particular fault heave and bed length conservation. The former corresponds to the required horizontal displacement to tie the uppermost horizon parts as a pure rigid motion of the fault blocks. The latter, in addition to joining the uppermost horizon parts, assumes that this horizon conserves its bed length and is restored to horizontal. In this case, the horizontal displacement is equal to the horizontal extension of the analog model in the X direction (i.e., down-dip direction) before restoration minus the sum of the lengths of the uppermost horizon parts. For each of these two methods, we used the unrestored model (with applied shortening boundary condition) geometry at each restoration step. Both methods provide displacement estimates that are signicantly less than the expected values (Table 3). In other words, rigid motion along faults is not an accurate measure of total tectonic displacement for our model and bed lengths did not remain constant through deformation. This latter conclusion is a known expectation for extensional structures (e.g., Xiao and Suppe, 1992). There is internal deformation accommodated by structures below image resolution or by deformation of a more continuous nature.

Area-depth method
We applied the area-depth method (Epard and Groshong, 1993;Groshong et al., 2003;Groshong, 2006;Groshong et al., 2012) to estimate the total forward extension without the need of a reference paleo-geometry. The areadepth method may be used to calculate the magnitude of shortening or extension of a system above a basal detachment. A benet of this method is that it accounts for the displacement due to faults omitted from the interpretation or tectonic strain that is below our imaging resolution, which has been found to accommodate up to 60% of total extension within a given system (e.g., Kautz and Sclater, 1988;Marrett and Allmendinger, 1992;Baxter, 1998;Groshong et al., 2003). The area-depth method is independent of the mechanical processes and is based on assumptions of area conservation and plane-strain, given that a thin detachment level exists. The area-depth method denes for each horizon a regional depth of detachment and a lost area inside the graben (below the regional datum and above the horizon), as shown in Figure 14. This lost area is equal to the product of the displacement that produced the graben and the depth to the detachment level. It follows that the total extension is given by the lost area divided by the depth to the detachment. We did not plot an area-depth H8 regional  graph of the entire growth sequence, which would integrate each lost area and each distance from the regional to a reference level, since by denition the layers did not undergo the same magnitude of extension (Groshong et al., 2003). Indeed, such a plot enables to evaluate the common displacement and the depth to detachment only for pregrowth strata or for no-growth sequences of growth strata (Groshong, 2015). Thus, we assume that the depth to the detachment is known. In our analog model, the denition of the detachment level is not straightforward, as the silicone layer is thick and may act as a distributed detachment zone. Assuming that no slip occurs along silicone boundaries (Weijermars et al., 1993), we approximated the detachment level to be at the middle of the silicone layer ( Figure 14). The regional level of each horizon is dened by a straight line dipping 1.5 • (parallel to the detachment) and starting from the intersection between the horizon and the most western fault (Figure 14). Our calculations only use the north CT image of the analog nal deformation stage.
The estimates of the shortening magnitude increments for the rst four horizons using the area-depth method and the CT images are given in Table 3. The amounts of displacement predicted by the area-depth method are within 15% of the shortening magnitudes provided by the CT images. We consider this a valid estimate given the structural uncertainties.

3D dilatation analysis
We propose a complementary approach to evaluate the model forward extension from calculations of dilatation, with V r and V u respectively the restored volume and the unrestored volume. As previously discussed, horizontal dilatation is expected during the experiment. Due to the small duration and the scale of the experiment, the strata used in our model are not expected to undergo signicant vertical compaction (Schultz-Ela, 1992).
Thus, we expect the volume to increase during forward deformation experiment. A consequence is that the volume should decrease during restoration, resulting in negative dilatation calculations from Equation (2). We ran a large number of geomechanical restorations varying the magnitude of shortening imposed as a boundary condition. Figure 15 represents the proportion of tetrahedra with a positive dilatation from the unrestored state to the restored state according to magnitudes of the imposed shortening conditions. In this way, we attempt to estimate the magnitude of shortening required to minimize the number of tetrahedra with positive dilatation.
Similarly, Durand-Riard (2010), with a contractional model, used a lateral (elongation) displacement to reduce the number of tetrahedra with a negative dilatation. The shortenings in Table 3 are displayed in the dierent graphs of Figure 15. As expected, the number of tetrahedra with a positive dilatation decreases when the magnitude of the applied shortening increases. However, in each scenario (Figure 15), a plateau of diminishing returns develops with additional applied shortening. The beginning of each plateau, as well as the area-depth estimates, provides a much improved estimate of the shortening magnitude than the other methods investigated above (Table 3). While this conclusion is still empirical, we suggest that dilatation may be an eective tool to estimate the magnitude of the lateral displacement boundary condition and to evaluate the validity of the restored state.  Figure 15: Proportion of tetrahedra with a positive dilatation after each restoration step (steps 1 to 4) according to dierent shortenings (in centimeters). Black dots are data points (restoration simulations). For each restoration, the unrestored state is the restored state at the previous restoration step with the imposed shortening from CT image (and without the uppermost restored layer). Shortenings in Table 3    There are several potential explanations for the requirement of an imposed shortening boundary condition. A rst reason is the granular nature of the growth strata. According to Groshong et al. (2003); Yamada and McClay (2003); Le Guerroué and Cobbold (2006), and Moretti and Callot (2012), dilatation is likely to occur in structural sandbox models when granular materials undergo shear. This eect enables faults to develop in unconsolidated materials (e.g., Colletta et al., 1991;Cobbold and Castro, 1999;Le Guerroué and Cobbold, 2006;Groshong et al., 2012). As deformation progresses, additional shear occurs and more voids develop in the system (Groshong et al., 2003;Le Guerroué and Cobbold, 2006). This disorder is at the origin of an increase of the global volume, and thus must be countered by applied shortening in the restoration. Table 4 provides the forward dilatation v (volumetric strain) of the sand and pyrex strata measured on the CT images for each with A u i and A d i respectively the area of the sand and pyrex strata on the CT image at deposition time of the layer i and just before the deposition of the layer just above the layer i. At each step, the forward dilatation is positive, which means the analog model area on the northern edge increased through time. Nevertheless, we have just such an evidence of dilatation on the northern edge thanks to the CT images; this dilatation may or may not be compensated elsewhere within the volume. However, due to the style of deformation of the sandbox model, it is very probable that the forward dilatation observed on the CT images is representative of the entire volume. Rock dilatation may exist in real extensional elds in which sediments contain uids (e.g., Boerner and Sclater, 1992). Another important reason which explains the need for the shortening boundary condition is that a part of the fault displacement may not be taken into account. Indeed, all the observed faults are not represented, and there may be faults below tomography resolution. Even if their osets are small, accumulated fault heaves may represent signicant forward extension. This is analogous to nature with faults below seismic resolution, which are not imaged and, thus, unable to be represented at the macro-scale (e.g., Groshong et al., 2003).

Residual amounts of fault dip slip values
Although the distributions in Figure 13 are encouraging, the number of inconsistent fault dip slip deltas increases with each successive restoration step. A possible explanation for this is that each residual dip slip on these faults is not corrected between restoration steps to t the dip slip observed on CT images, leading to the accumulation of errors. Another possible explanation is that some faults are kept within the volumetric model whereas they were not present in the analog model at the time corresponding to the restoration step. As faults behave as sliding surfaces in our volumetric mesh, small articial slip may be present on these surfaces, leading to local inconsistent shear strain. We kept these faults to avoid rebuilding a new structural model, which can be quite time-intensive for such complex fault networks (Zehner et al., 2015). Another observation of these distributions suggests that the restorations seemed to have recovered too much dip slip (numerous negative deltas). Since the forward deformation path involves friction on faults, this result may be due to the frictionless contacts in our mechanics-based restoration method (Wriggers and Laursen, 2006).

5.3
Mismatches with the area-depth method As mentioned previously, the area-depth method provides a reasonable estimate of the incremental extension that occurred during the forward model (Table 3). However, the estimates are not perfect, in particular for the restoration step 4 (Figures 12 and 15). Several factors may explain these errors. First, material dilatation observed on the CT images and attested by several authors (Yamada and McClay, 2003;Le Guerroué and Cobbold, 2006) is inconsistent with the area conservation hypothesis underlying the area-depth method. As the upper horizons have accumulated less dilatation than the bottom horizons, the constant area hypothesis deteriorates with each successive restoration step. Second, the silicone layer could migrate laterally and blend with the sand and pyrex, leading to area changes. In the analog model, from the CT image of the restoration step 4 to the CT image representing the rst unrestored state (northern edge), we calculated a forward dilatation of the silicone to ∼6.8%. Third, as a part of the analog model on the eastern side was not available for analysis in our CT tomography images, we could not integrate this data in our area-depth computations. Fourth, the denition of the detachment level, even based on several reasonable assumptions, is uncertain. Fifth, as the units of sand and pyrex can penetrate into the silicone, the resulting subsidence modies the denition of the regional levels. This eect is equivalent to the oating regional mentioned by Groshong (2015) for a buckle-style fold above a thick salt unit.

Boundary conditions
This study suggests that for extensional systems, the combination of classical boundary conditions and a new lateral displacement boundary condition along the dominant transport direction ( Figure 5) may yield consistent restored geometries. In this paper, we also propose the use of novel contact conditions to ensure consistent restoration of complex branching and crossing fault geometries. These new constraints enable eective sequential restoration of four steps of the analog model. Without them, only two steps could have been performed, and quality of these restorations would have been reduced. We believe that the boundary conditions presented in Figure 5 can be applied in compressive contexts with an elongation displacement condition instead of the shortening condition. Indeed, Durand-Riard (2010) shows that an elongation condition is necessary to properly restore a fault-bend fold model. An estimation of the elongation may be done using the area-depth method (Groshong et al., 2012). In case of strike-slip faults, displacement conditions parallel to the strike direction should also be considered, as shown by Durand-Riard et al. (2013b).
In the literature and in this paper, all the boundary conditions correspond to displacement conditions except for the mechanical contact conditions which are a mix between displacement and traction conditions (Muron, 2005;Wriggers and Laursen, 2006;Maerten and Maerten, 2006). Such displacement conditions may lead to unphysical strain elds (Lovely et al., 2012). In reality, rock deformation is a consequence of force constraints. Maerten and Maerten (2006) suggest the possibility of employing mechanical boundary conditions that incorporate the far eld stress as an additional boundary condition. The main diculty of this technique would be to know the intensity of the forces to apply (Muron, 2005). A rst start could be to use the determined displacement condition for a model (e.g., Figure 15) and convert it to a force: dilatation multiplied by Young's modulus in linear elasticity. In addition, the overburden force is not incorporated in our geomechanical restoration method. As our experiment was gravity-driven, it would be interesting to add an overburden body force to the nite element procedure to analyze its impact on the restored geometries.

Conclusions
The restoration of an analog model, in which the structural uncertainties are limited and paleo-geometry is well known, enabled us to dene eective boundary conditions that yield optimal restored models using mechanicsbased restoration. For extensional structures, a shortening boundary condition was applied to obtain a good t with reference paleo-geometries. Such a condition may be estimated by the area-depth method. Our experiments suggest that an analysis of the volumetric dilatation can complement the estimate of the shortening boundary condition magnitude. Moreover, to handle complex fault networks, we propose the application of contact conditions on internal fault borders and between fault connected components. Ultimately, the methods developed in this paper, in particular the lateral displacement boundary condition, should lead to improved results if applied to geomechanical restorations of natural structures.
work. Frantz Maerten, Robert Worthington and Frank Zwaan provided helpful comments that improved the quality of this manuscript.