Guillaume Caumon

Pierre Muron

Surface restoration as a

come

Introduction

The geometry of geological models is a significant source of uncertainty in petrophysical modeling and flow simulation. For instance, a small variation in the geometry of a horizon may significantly change the modeling output (volume in place, production forecast, drilling strategy, etc).

Part of this geometrical uncertainty can be assessed from a geophysical standpoint, for instance by considering the velocity model used during migration and time-to-depth conversion, the resolution of seismic data, etc. [START_REF] Thore | Structural uncertainties: detramination, management and applications[END_REF]. These geophysical sources of uncertainty affect significantly the geometry (position and shape) of the geological interfaces, hence must be accounted throughout the exploration / production process.

In this paper, we focus on another type of unknown: for a given geometry of the geological surfaces, it is difficult to estimate the transverse slip of a fault, unless markers are available on both sides of this fault (Figure 1. Then, several possible fault displacement fields, can be obtained by varying the transverse fault slip By default, 3D modelers assume, at least locally, a dip slip displacement, or align the displacement field on the domain boundary (Figure 2). This arbitrary decision is generally not correct from a structural standpoint, hence should be reconsidered.

Choosing the wrong transverse slip may be consequential for other 3D subsurface modeling tasks:

• Modeling of sedimentary facies and properties is typically done in depositional space, using a stratigraphic transform such as Geochron [START_REF] Mallet | Space-time mathematical framework for sedimentary geology[END_REF]. In this context, the fault throw impacts the way distances are computed in stratigraphic space, hence the petrophysical models output by geostatistical algorithms. Typically, errors in the transverse slip can lead to errors in the fault transmissibilities. A 2 Although the transverse throw is not known, it is generally assumed to be locally collinear to the domain boundary, which is not necessarily equivalent to a dip slip displacement (B).

3 F 1 F 2 F 1 F 2 1 2 3 B F 1 2 A B dip throw throw
• Geomechanical approaches to assess the density of fractures from strain [START_REF] Mace | Fast generation of 3d discrete fravture networks in respect of geology[END_REF] rely on the definition of compatible deformations. Errors in the structural model may impede the predictive power of such fracture models.

We propose to use surface restoration to stochastically compute several alternative models of fault displacement on a faulted horizon with a known geometry. These models are obtained by a Monte Carlo Markov Chain sampling method [START_REF] Tarantola | Inverse Problem Theory[END_REF]. This method calls for a representation of fault throw that can be perturbed (Section 1), and a means to assess the compatibility of the displacement field for each sample model (Section 2). With this, the sampling algorithm can be applied to select one or several acceptable models (Section 3). 1 Fault slip representation

E 1 E 2 E 3 E 1 E 2 E 3

Surface borders

A 3D structural model can be defined by a set of surfaces representing geological interfaces. These surfaces consist for instance of faults, stratigraphic surfaces, unconformities, salt boundaries, etc. In this paper, we consider one faulted horizon H of some structural model. The boundary of this horizon is subdivided into k borders {∂H 1 , . . . , ∂H k }, which identify the relations of H with its surroundings. For instance, a horizon completely cut by a fault has two borders on both sides of that fault, plus external borders along the domain of interest. A border is defined by two border extremities. To modify the extent of a border, it is sufficient to displace one of its extremities along its two incident borders.

Representation of fault slip

Typically, the slip of a fault F is defined for a given horizon H by linking both borders of H located on the hanging wall and on the foot-wall of F . For convenience, this matching of borders often relies on geometrical considerations. For instance, border extremities are set on the horizons at the corners of fault blocks, and the throw is distributed evenly between these borders. This is achieved by associating locations having the same curvilinear abscissa on both hanging wall and foot wall borders (Fig. 3-A).

Such an automatic layout of border extremities may not always be realistic, especially when two points are matched only because they lie on the model boundary. In our example, one could use metric curvilinear coordinates from the fault tip and obtain a different throw direction (Fig. 3-B). Instead of curvilinear coordinates, a search for common horizon features on both sides of the fault could also be used.

Fault slip perturbation

In complex configurations with several fault compartments on the horizon, the method suggested in Figure (Fig. 3-B) is not applicable. Therefore, finding an acceptable displacement calls for perturbing the transverse fault slip.

In the general case, the fault slip configuration on some faulted horizon H is defined by:

• a set E 1 , E 2 , . . . , E k of border extremities on H
• a set of throw links associating pairs of borders across each fault [Mallet, 2002, p. 272] For all three types of corner, the transverse slip is made variable by moving one border extremity E i , within its two adjacent borders. The fault throw is then updated using the normalized curvilinear coordinate method depicted in Figure 3. The position of the movable border extremity can be defined relatively to some reference position by a curvilinear coordinate s. The orientation of the horizon can be used to define positive or negative curvilinear coordinates around this reference position. For type 2 corners, positive values of S i will concern say the border extremity on the foot-wall, while negative values will affect the position of the border extremity on the hanging-wall (Fig. 4-B).

For perturbation purposes, we propose to model the position of the movable border extremity E i by a Random Variable S i . For each free border extremity E i , the prior probability distribution P 0 (S i = s) of the curvilinear abscissa is provided as input (Fig. 4). The prior distributions P 0 (S 1 ), . . . , P 0 (S n ) of all n corners can be sampled by a Monte Carlo Method to generate a large number of fault displacement fields. However, all these models are not equally likely to occur. To assess the likelihood of each model, we propose using surface restoration techniques.

Compatibility criterion 2.1 Restoration of geological surfaces

Classically, balanced restoration of a faulted and folded horizon transforms the horizon into a plane subdomain. This can be achieved by computing a 2D curvilinear coordinate system (or parameterization) on the horizon (Fig. 5). This operation basically cancels out fault displacements (unfaulting, or sewing operation), and deforms the original surface to obtain a planar surface (unfolding, or ironing operation). Some types of surfaces are developable, meaning that the unfolding operation can be achieved without distortions. Developability is a good

I R 3 (x,y,z) x y z D (u,v) Φ(x,y,z) x(u,v)
x -1 (x,y,z)

Figure 5: The restoration of a horizon is a parameterization that maps the horizon onto a plane. Note that the fault throw is cancelled out during this operation.

The unfolding problem has been solved by three main types of approaches:

• Packing methods [START_REF] Gratier | Restoration and balance of a folded and faulted surface by best-fitting of finite elements: principles and applications[END_REF][START_REF] Williams | Restoration and balance of complex folded and faulted rock volumes: flexural flattening, jigsaw fitting and decompaction in three dimensions[END_REF][START_REF] Rouby | 3-D restoration of complexly folded and faulted surfaces using multiple unfolding mechanisms[END_REF] independently project or rotate all surface triangles to a plane, according to some deformation mechanism (e.g., flexural slip or inclined shear), then minimize gaps and overlaps between triangles by rotating and translating, then shearing triangles.

• Global geometric optimization methods [START_REF] Lévy | Non-distorted texture mapping for sheared triangulated meshes[END_REF][START_REF] Massot | Implémentation de méthodes de restauration équilibrée 3D[END_REF][START_REF] Lévy | Least square conformal maps for automatic texture generation[END_REF][START_REF] Thibert | A direct method for modeling and unfolding developable surfaces and its application to the ventura basin (california)[END_REF][START_REF] Sheffer | Abf++: fast and robust angle based flattening[END_REF]] compute a 2D curvilinear coordinate system (u, v) on the 3D surface. Three numerical criteria have been formulated to compute this coordinate system so as to minimize surface distortions between present (x, y, z) and restored (u, v) states:

-Isometric mapping [Lévy andMallet, 1998, Massot, 2002] defines two constraints :

||∇u|| = ||∇v|| = 1 (conservation of distances) ∇u • ∇v = 0 (conservation of angles)
-Conformal mapping [START_REF] Lévy | Least square conformal maps for automatic texture generation[END_REF] aims at preserving angles using complex number formulation; a consequence of this conservation is that ∇u and ∇v should have the same norm, which is not necessarily equal to 1.

-Angle-based flattening [START_REF] Sheffer | Abf++: fast and robust angle based flattening[END_REF][START_REF] Thibert | A direct method for modeling and unfolding developable surfaces and its application to the ventura basin (california)[END_REF] minimizes variations of angles directly in angle space, which produces less distortion than previous approaches, but requires solving a non-linear problem.

• Geomechanical methods [START_REF] Dunbar | Palinspastic reconstruction of structure maps: an automated finite element approach with heterogeneous strain[END_REF]] model the geological horizon as an elastic membrane and use finite element methods to restore the surface with a minimum strain. This global approach is an elegant solution, but it is difficult to assess whether the physical model is appropriate: a geological surface is but an interface between two deformed layers. Therefore, we are more confident in geomechanical approach applied on geological volumes [START_REF] Muron | Handling faults in 3d structural restoration[END_REF].

While unfolding is always treated as an optimization problem, unfaulting is addressed in various manners. [START_REF] Thibert | A direct method for modeling and unfolding developable surfaces and its application to the ventura basin (california)[END_REF] manually rotate and translate rigid fault blocks using a graphical program. [START_REF] Meeting | Least-squares palinspastic restoration of regions of normal faulting = application to the Campos basin (Brazil)[END_REF] propose a packing method to minimize gaps and overlaps between rigid, planar fault blocks. This method is similar in spirit to the triangle packing method of [START_REF] Gratier | Restoration and balance of a folded and faulted surface by best-fitting of finite elements: principles and applications[END_REF], but is adapted to handle complex polygonal fault blocks. Both the automatic and manual packing methods are interesting since they clearly separate the unfolding and unfaulting steps [START_REF] Rouby | 3-D restoration of complexly folded and faulted surfaces using multiple unfolding mechanisms[END_REF]. However, rigid fault blocks approximations are not exact, since gaps and overlaps may remain after restoration; packing methods also require faults terminating within the horizon to be extrapolated to define separate fault blocks.

In geometric optimization methods [Lévy andMallet, 1998, Massot, 2002], the fault displacement field is provided as input before the unfolding. Indeed, optimizing fault throw direction within a global geometric optimization is difficult, because of non linear effects on the objective function. In Section 3, we will see how inversion theory can be used to address this problem and capture some uncertainty about the fault displacement field.

Thanks to their finite element geomechanical approach, [START_REF] Dunbar | Palinspastic reconstruction of structure maps: an automated finite element approach with heterogeneous strain[END_REF] allow the horizon borders on the hanging wall and foot wall of a fault to slide against each other. They also define a set of cinematic rules for transverse slipping along fault contacts, but rely on interpretive input to define fault slip when a fault cut the horizon from end to end as F 2 in Figure 2. From these rules, the optimization code restores the horizon with minimum strain.

Strain calculation

For all restoration methods, the dilatation θ can be computed everywhere on the horizon H by considering the change of area of a unit disc dS (0) in the restored state when deformed into the ellipsoid dS at present: θ = dS -ds (0) dS The (orthogonal) directions of the strain ellipsoid dS are the principal strain vectors, and the radii of the ellipsoid give the principal strain λ 1 and λ 2 .

Restoring the same surface using several transverse fault slip configurations yields several maps of dilatation and strain. Several global metrics can then be proposed to assess the likelihood of a fault slip configuration. For these metrics, we propose considering some measure of spread (standard deviation, interquartile range) of one or several quantities over the horizon H, for instance:

• The dilatation θ around 0.

• The ratio of the principal elongations λ 1 /λ 2 around 1, since the strain tensor should be isotropic due to the geometric constraints.

• The scalar product of the curvilinear frame vectors (∇u•∇v), since the angles should be preserved (per the conformal constraint.

Since it is difficult to assess which maximum values are acceptable, these measures of likelihood are considered by the sampling process not in absolute values, but relatively from one realization to the next.

Sampling algorithm 3.1 A Monte-Carlo approach

We propose to sample possible transverse fault slip configurations using a Markov Chain Monte-Carlo (MCMC) sampling algorithm [START_REF] Tarantola | Inverse Problem Theory[END_REF]. The algorithm starts with a set of prior un-certainty distributions {P 0 (S 1 ), . . . , P 0 (S n )} (Figure 4), and a vector of curvilinear coordinate values {s 10 , . . . , s n0 }. The algorithm then proceeds as follows: (f) If the current realization j is retained, m prev ← m cur ; j ← j + 1

1. j ← 1; m cur ← ∞; m prev ← ∞;
This MCMC sampling technique considers the displacement of all faults at once, which makes full sense in a global restoration framework. No bias is introduced thanks to the use of a random path. Performance may be an issue with large surfaces; however, geological expertise can drastically help increasing performance by freezing some fault throws deemed well defined, by setting the prior uncertainty of some corners to 0.

Although MCMC has a reputation of slowness, the performance of the algorithm is very reasonable: hundreds of configurations are sampled in a few minutes, due to the speed of the restoration algorithm.

Dependency among model parameters

Classically, Monte Carlo sampling consists of independent drawings of model parameters. However, drawing all values s i(j-1) independently from each other may lead to flipping of borders, if two border extremities cross each other. This type of inconsistency could be left for the restoration algorithm to be detected.

Instead, because such border flipping is obviously wrong and should not impede the method performance, one could update the prior distribution P 0 (S i ) for the corner c i by the values sampled for the adjacent corners c i -and c i + , and then sample from P (S i |S i -, S i + ). A more convenient solution is to use normalized curvilinear coordinates between -1 at c i -and 1 at c i + along the two borders incident to c i .

Implementation details

When geological interfaces are represented as triangulated surfaces, the displacement of a border extremity may slightly change the geometry of the horizon H. Therefore, the reference horizon geometry must be cached by the sampling system. Reverting to the geometry of the reference surface between two perturbations is needed not to introduce any bias in the sampling.

Results

For validation purposes, the method has been applied on a synthetic faulted horizon shown in Figure 6. Each fault block is planar, so the deformation comes only from the geometry and extension of borders. 300 iterations were run with a prior uncertainty of 20% of the total border length. The perturbation from one realization to the next one was limited to 5% of the prior uncertainty range. The restoration quality metric m cur was the interquartile range of the dilatation. Although the sampling algorithm is successful in producing several realistic fault throw models, some of the "best" realizations are only optimal with regard with the proposed criterion. In particular, some of the models in Figure 6 have large dilatation locally, which is not seen using the interquartile range. Visual inspection is then needed to select acceptable models, before better criteria are defined.

The method can naturally be applied to more complex models. In Figure 7, the throw on a horizon built from true seismic data is considered (reference model courtesy of Total). The workflow parameters are the same as in the previous case, except for the quality metric, which is the p 99 -p 1 of the dilataion to avoid local extrema. Borders have been split inside the faults to introduce more flexibility. Results obtained show that all models are very close to each other, which suggests that the configuration corresponds to a quite stable minimum of the objective function.

Suggested improvements

• Both geometric [Lévy andMallet, 1998, Massot, 2002] and geomechanical [START_REF] Dunbar | Palinspastic reconstruction of structure maps: an automated finite element approach with heterogeneous strain[END_REF] global approaches do not explicitly separate the horizon strain due to the unfaulting from the strain due to unfolding. For this, we suggest first unfolding each fault block separately without extrapolating dying faults. Then, reapplication of the optimization method would force a perfect continuity along faults, and allow checking the influence of transverse fault slip.

• The suggested Monte Carlo sampling technique is reputably slow to explore the space of uncertainty. While it works in a reasonable amount of time for horizons affected by a few faults only (about 10 minutes for the case study depicted in Figure 7), the computational cost to obtain significant results would be much higher in heavily faulted domains. In this case, a divide-andconquer approach could be used to optimize the sampling process: instead of randomizing the all fault displacements between each realization, the optimization could be done fault by fault.

Discussion and Conclusion

Subsurface uncertainty modeling has long been limited to petrophysical uncertainty assessment through geostatistical stochastic modeling. However, it is clear that a significant part of subsurface uncertainty stems from the lack of structural and sedimentological information and from model simplifications.

Proper uncertainty assessment then relies on randomizing structures [START_REF] Lecour | Modeling of stochastic faults and fault networks in a structural uncertainty study[END_REF][START_REF] Thore | Structural uncertainties: detramination, management and applications[END_REF], geological scenarios [START_REF] Massonnat | Can we sample the complete geological uncertainty space in reservoir-modeling uncertainty estimates[END_REF][START_REF] Caumon | Assessment of global uncertainty for early appraisal of hydrocarbon fields[END_REF], and seemingly departing from the principle of parsimony by considering parameter-rich models.

This last point is worth a remark. The principle of parsimony states that one should always choose the simplest explanation (model with as few parameters as possible) for a phenomenon. This principle is often wrongly simplified by taking the simplest explanation for observations, which forces considering a new model with more parameters whenever a new observation is made. The distinction does not make much of a difference in deterministic modeling, which uses optimization theory. It does in stochastic modeling: the challenge of all good uncertainty models is to choose from the start enough parameters to possibly include all past and future observations. The model should then be suited to the phenomenon / geological concept, and not to some partial and noisy data.

In this paper, we have proposed a method to sample, from a given geometry, several ways of associating fault gaps on a horizon by stretching or shrinking the faulted borders. Acceptable throw models are selected according to the Metropolis rule on the distortion observed after 2D restoration. The limiting assumptions of the method are therefore intrinsically connected with those of the 2D restoration: the deformation and faulting should be connected to one major tectonic event, and the initial deposition geometry should be relatively planar. This method also uses one single scalar criterion to screen surfaces; using a combination of several restoration quality metrics would definitely improve the results obtained.

This work is just one step towards a better characterization of structural uncertainties. Future works include sensitivity analysis to assess the impact of fault throw uncertainty as compared to the other sources of uncertainty. A possible -but time-consuming-improvement for this work would be randomizing the throw direction jointly on several horizons, which would amount to use 3D restoration as a model selection criterion.

Figure 1 :

 1 Figure 1: Anatomy of a fault. The net slip can be decomposed into dip slip and strike slip in the fault plane, or into horizontal (strike slip + heave) and vertical throws. In this paper, we consider the transverse slip as seen from a faulted horizon.

Figure 2 :

 2 Figure 2: Typical default definition of fault displacement in Gocad (left: global view of a three-block synthetic model; right: details of fault displacement vectors): along the contact A, the fault slip of F 1 must be locally aligned on the fault contact F 1 ∩ F 2 for compatibility.Although the transverse throw is not known, it is generally assumed to be locally collinear to the domain boundary, which is not necessarily equivalent to a dip slip displacement (B).

Figure 3 :

 3 Figure 3: Definition of the fault slip on a faulted horizon border. By default (A), the displacement vectors are defined using normalized curvilinear coordinates. It would make more sense to use metric coordinates from the fault tip (B).

Figure 4 :

 4 Figure 4: Sampling the position of a border extremity in two current corner configurations: a branching fault contact (left) and the contact of a fault on the model boundary (right). Borders across a fault gap are associated by dashed lines. The probability distribution for each curvilinear abscissa can be randomly sampled from the input cdf.

  Define a random path visiting all n corners (b) For each corner c i , draw randomly the abscissa s ij correlated with the previous abscissa s i(j-1) , and update the corner geometry (c) Unfold the surface using the appropriate surface restoration approach (simple shear or global optimization technique) (d) Assign the global surface distortion to m cur (e) Compare m cur and m prev , and apply Metropolis' rule to decide whether to keep or discard the current realization {s 1j , . . . , s nj }.

Figure 6 :

 6 Figure 6: Simple case of a three fault block flat horizon. The reference model and the 8 "best" realizations with regard to the considered criterion (interquartile range of the dilatation), out of 300 sample models.

Figure 7 :

 7 Figure 7: The reference model and the 5 "best" realizations (out of 300) obtained on an actual horizon.

Acknowledgements

The authors want to thank the members of the Gocad Consortium for their support. We also thank EarthDecision, Bruno Lévy, Jéroôme Massot and Julien Massenet for providing part of the software used in this work, and Albert Tarantola for his faith in MCMC techniques.