## Abstract

Accurately positioning faults in a geological model is a major concern because they are responsible for offsets of geological sequences. In the tetrahedral models studied in this paper, faults are discontinuities: faces of tetrahedra on either side of a fault are disconnected. Building tetrahedral models can require a large amount of time, especially when there are many faults. We present a tool for making small, real-time, modifications of faults in tetrahedral models arising from geometrical changes required either by new subsurface data or by new interpretations of existing subsurface data. Fault editing is achieved by moving control points on the fault in the tetrahedral grid and by computing a distortion property over an editable volume relative to the control point and spreading this distortion to neighbouring points using the Discrete Smooth Interpolation technique. The editable volume in which tetrahedron vertices are allowed to move is defined by a given distance to the fault. This approach provides a means of editing faults and fault-related features, such as branch-lines.

Grid building remains one of the main challenges of reservoir modelling. Geometric accuracy must sometimes be overlooked in favour of simple grids, which require fewer computing resources. Generally, reservoirs are modelled as curvilinear grids – Cartesian grids in which the layers of cells are deformed to match horizon and fault geometry. These grids offer a fair compromise between precise geometry and easy computation for property simulation, fluid flow and velocity modelling. However, in some instances, reservoirs are so complex that a curvilinear grid does not reflect the geometry of layers and faults accurately, and sometimes cannot even be constructed without making unacceptable approximations. One option for these complex geometries is to use unstructured and irregular meshes, such as polygonal or tetrahedral grids (Lepage 2002; Prévost *et al.* 2004). The resolution of the mesh can be adapted to the structures to be represented: the mesh is coarse in simple regions and denser where the structures or heterogeneity are more complex.

Geological properties are also affected by faults. Some properties, usually related to the geometry of the model, for example distance from a well, may be continuous from one side of the fault to the other. In this case, values on both sides of the fault are identical. Other properties, generally related to rock type, such as porosity, may be continuous before displacement by the fault, but may not be continuous after fault displacement. A method involving vectorial links (Moyen 2005) restores the continuity of the property using fault throw (Fig. 1). Points that were at the same location before faults offset the geological sequences are linked so that property values are equal on both sides of faults. In the first part of this paper, tetrahedral grid building and volume distortion methods are reviewed. Based on the concepts defined in these methods, a fault editing algorithm is presented in the second part. Finally, examples of consistent geological model editing are explained and illustrated.

## Building and editing tetrahedral grids

### Geological applications based on tetrahedral grids

Because their geometry is flexible, tetrahedral meshes have several applications in geology. Tetrahedra can be made to fit the geometry of geological structures and as they are the simplest possible simplexes in three dimensions, tetrahedra are interesting to perform computations on.

#### Velocity modelling

Velocity models play a key role in processing seismic data, especially in pre- or post-stack time and depth migration, as well as in time-to-depth conversion (Yilmaz 1987). Quantitative interpretation of seismic data uses elastic properties derived from velocities. Comparing the results of velocity modelling with real seismic data for a study area can highlight defects in the velocity model or in the layout and geometry of geological structures. Macy & Smith (1998) use a tetrahedral model in which several regions with different geological attributes are separated by velocity discontinuities. The ray path is calculated in each intersected tetrahedron in turn, using the gradient of the velocity function. Velten (1998) defines tetrahedron columns inside the volume model and computes ray paths using several options for the velocity law inside columns.

#### Three-dimensional restoration

Three-dimensional structural restoration provides strain fields that help validate structural interpretations. Muron developed a restoration method for tetrahedral models based on volume conservation and strain minimization, taking into account mechanical rock properties (Muron 2005). Most other volume restoration methods are not fully three-dimensional: stacks of horizons linked by geological constraints are restored simultaneously (Griffiths *et al.* 2002).

The strain fields generated by three-dimensional restoration can be used as input for fracture simulation in tetrahedral models (Macé *et al.* 2005, 2004). These fractures are modelled using a combination of stochastic methods and geological rules (Cacas *et al.* 1990; Josnin *et al*. 2002).

#### Fluid flow modelling

Flow simulation is a key step in understanding the dynamic of an oil or water reservoir. Verma (1996) discusses flow simulation in reservoirs using several types of grids, including Voronoi grids. A Voronoi grid is the dual mesh of a tetrahedral grid. It is obtained by replacing tetrahedra with their barycentres, and linking these points together (Okabe *et al.* 2000). Voronoi grids are used to model fluid flow in reservoirs using the finite elements or finite volume methods (Palagi *et al.* 1993). Two and three-dimensional unstructured cells can be adapted to reservoir heterogeneity (Prévost *et al.* 2004). Streamline-based simulations are an alternative to finite elements and finite volume methods (Matringe *et al.* 2006) and help assess the accuracy of upscaling in these grids and provide an insight into flow behaviour (Baker *et al.* 2001; Blunt *et al.* 1996).

#### Geochron framework

Geometry of geological structures and property modelling can be separated using the Geochron framework (Mallet 2002*b*, 2004; Mallet *et al.* 2004). A three-dimensional parameterization with geological considerations such as fault–fault and fault–horizon contacts and sedimentary discontinuity surfaces is computed on a tetrahedral model of the study area. This parameterization associates each node in the tetrahedral grid to its position in the depositional space. If no fault was active during deposition the layers are perfectly flat in the depositional space, so geological properties can be modelled in this space using a fine Cartesian grid. Computing petrophysical properties such as porosity in the depositional space reduces errors introduced by gridding limitations. The parameterization enables mapping high resolution geological properties in the depositional space onto the coarser tetrahedral volume.

### Creating tetrahedral grids

There are three main groups of tetrahedral grid creation methods (Owen 1998): octree-based, Delaunay-based and advancing front methods. The most popular meshing methods for geological applications are Delaunay-based.

First, interpreters create fault and horizon surfaces from seismic and well data and from geological and regional knowledge (Fig. 2a). Faults represent the structural information in the study area, and partition the subsurface in fault blocks (Fig. 2b). If required, horizon surfaces can be integrated in the model as implicit surfaces. Then, a tetrahedral mesh of the model may be generated using a Delaunay criterion (Lepage 2003). As shown in Figure 2c, faults do not have to extend through the whole tetrahedral volume. Because they are unstructured, tetrahedral volumes can represent complex geometries incorporating, for example, discontinuous faults and fault branching. Rock properties can be stored and visualized directly on the tetrahedral volume or using the Geochron framework (Mallet 2002*b*, 2004; Mallet *et al.* 2004), a three-dimensional texture mapping enables visualization of properties stored in a regular grid (Fig. 2d).

Tetrahedral models in geology can be divided into different fault blocks. Inside a fault block, tetrahedra are all connected, and each tetrahedron shares vertices and faces with its neighbours. A fault is modelled as a topological discontinuity in the volume (Fig. 3a). Tetrahedra on both sides of a fault do not share faces and vertices are duplicated, although the geometrical position of the nodes is the same on both sides. To make property modelling easier on a tetrahedral volume, faces of tetrahedra that are on faults have the same geometry as the triangulated surfaces from which the model was created. In other words, on a fault, three entities have exactly the same shape and location: the triangle of the fault surface and the faces of tetrahedra on both sides of the fault (Fig. 3b).

### Model updating

New information obtained for a geological volume can make a model obsolete. For instance, a new interpretation of seismic data can show that the geometry of a surface is not accurate, or a new well can provide another data point for the location of a surface. Faults are discontinuities in the grid, hence fault editing requires modification of the geometry and possibly the topology of the tetrahedral mesh. The current way of dealing with this problem is to edit the structural model and build a new tetrahedral volume. However, obtaining a high-quality mesh that conforms to geological structures requires time and effort. Local editing techniques on geological models provide a way of making small adjustments to models without having to go through the whole meshing process.

Volume distortion is a very active research field, especially in computer graphics and animation. Free-form deformation (FFD) incorporates a number of techniques for modifying the shapes of objects by the manipulation of control nodes or curves. Barr (1984) introduced a combination of simple object-editing transformations. FFD was made easier to use by Sederberg & Parry (1986), Coquillart (1990) and Chang & Rockwood (1994), who used different kinds of control lattices or curves in order to distort objects. But these methods did not take into account the nature of objects being edited. Hirota *et al.* (1999) added physical laws to FFD in order to preserve the total volume of objects during geometrical distortion.

Borrel & Rappoport (1994) introduced *simple constrained deformation* (SCODEF). Constraint points possessing a radius of influence and desired distortion determine local B-spline basis functions. The entire space, along with the objects in it, is then distorted according to a linear combination of these functions, creating bumps in space. Real-time geometrical distortion is possible with this method. Grosse (2002) adapted FFD and SCODEF to geological problems and described a number of solutions for interactive editing of geological models. Caumon *et al.* (2004) developed tools for editing a geological model in real-time whilst retaining other essential constraints, such as contacts between horizons and fault surfaces.

One of the problems faced when modifying the geometry of a set of vertices in a tetrahedral volume is that tetrahedra may flip over. When that occurs, the mesh is no longer valid, because the sign of the algebraic volume of inverted tetrahedra has changed, and because some edges may intersect faces of these tetrahedra. Figure 4 shows an example of tetrahedron inversion where one point is moved through the opposing face of the tetrahedron. Shontz & Vavasis (2003) had to tackle this problem when they modelled the beating of a canine heart using a tetrahedral mesh. They applied geometrical distortion to the triangulated surface bounding the tetrahedral mesh. Then the points inside the volume were moved using linear weighted Laplacian smoothing. A constraint on the function used to deform the boundary of the volume ensured that the mesh did not flip over. They obtained satisfactory results for the canine heart, but the linear weighted Laplacian smoothing could not prevent inversion when the geometrical distortion was large.

In conclusion, some of these volume editing methods can be adapted with varying degrees of success to geological models. In FFD methods, control curves are distorted and the object being edited is updated accordingly. SCODEF is an interesting method but it creates bumps centred on constrained points, whilst editing fault geometry requires smooth geometrical distortion in the control volume. Our method tries to solve these problems and to provide an intuitive, real-time distortion method adapted to geological problems.

## Fault-editing method

The fault-editing tool presented here permits small adjustments to fault geometry in real-time. An editable volume around the fault is defined and it is distorted by picking and dragging arbitrary control points. The geometrical distortion of the tetrahedral mesh is instantly visible and editing can be stopped when the new geometry is acceptable. Although the fault-editing tool was developed as a plug-in to the gocad software, similar approaches could, in principle, underpin methods implemented in other software systems.

### Fault editing

Before modifying fault geometries in a tetrahedral model, user-defined control points must be created on the fault which is to be edited and the volume in which the tetrahedron nodes will be allowed to move must be selected. In order to define this volume, a property is computed from the distance from the fault and stored at the nodes of the tetrahedral mesh. This distance property can be adapted to provide different shapes of control volumes. If only the distance from the fault is used, the control volume has boundaries that are parallel to the fault. Using the distance from a point at the centre of the fault will provide a spherical control volume. Any combination of these distance functions can be computed, thus varying the shape of the control volume. A user-defined threshold value for the distance property defines the editable volume: tetrahedra with property values below this threshold are inside and tetrahedra with property values above the threshold are outside the editable volume and are therefore fixed. A low value enables small changes on a part of the fault, whereas a higher value enables general reshaping of the fault. If the threshold is too low, there will not be enough tetrahedra inside the editable volume to enable fault editing. In order to make the volume selection step easier, visualization tools can also be used to display the distance property on surfaces within a model (e.g. Frank 2006).

As a control point is edited, by user-defined dragging inside the control volume, a distortion vector is computed between the starting position of the control point and its current position. This vector **V** provides distortion direction and magnitude, with the direction *W*=(1/|**V**|). **V** stored for later use and the magnitude computed as a one-dimensional property on the control volume. This distortion property is interpolated from the value set by the control point to zero on the borders of the editable volume using the Discrete Smooth Interpolation technique (see below; Mallet 1992, 2002*a*). This means the distortion varies smoothly from a maximum value at the control point to zero on the borders of the editable volume. Tetrahedra outside the editable volume are not distorted, whereas the vertices of the tetrahedral mesh that are inside the editable volume are moved along **W**, using the distortion property computed for each vertex.

If required, **W** can be modified so that each vertex is moved in a geologically consistent way. This is of importance when a Geochron parameterization is available: horizons in the tetrahedral mesh are represented by a pseudo-time function. If this function is distorted, so is the geometry of the horizons. In order to solve this problem, when a Geochron parameterization is provided, the component of vector **W** that is normal to the horizons is computed for each vertex and removed. This ensures that the vertices are moved along constant values of the pseudo-time function, thus preserving the geometry of the horizons.

Vertices must also be moved in a geometrically consistent way, otherwise tetrahedra could invert. In order to prevent such inversions in the tetrahedral mesh, checks are performed at each step of the distortion. If there are no mesh inversions, editing can proceed and the geometry is updated. If the mesh does invert, the distortion is computed from the distortion vector stored in the previous step and the tool is terminated. Editing is carried out by moving control points around in the selected control volume. When a control point is moved, the editing tool enters a loop in which the distortion property is computed over the control volume and the geometry is updated: these editing steps are shown in Figure 5.

### Background on Discrete Smooth Interpolation

The Discrete Smooth Interpolation (DSI) method, described in detail in Mallet (1992, 2002*a*), is used as a linear solver to compute the distortion function. The tetrahedral mesh provides a finite set of interconnected nodes constituting a linear model. The DSI method minimizes the degree of violation of soft constraints set on the linear model. A constraint is a set of equations linking the nodes in the linear model and translating an idea on how the interpolated property should vary between the nodes that are involved. For instance, one of the basic DSI constraints is the roughness constraint, which states that the interpolated function should vary smoothly from one node of the model to the other. This roughness constraint is translated as a set of coefficients linking neighbouring nodes in the linear model.

The DSI problem can be expressed as solving the set of linear equations
1
with **A**_{ci} a vector containing the coefficients of constraint c_{i} at each node of the linear model, here the tetrahedral mesh. *φ* is the unknown function at the nodes. *k* is the number of constraints on function * φ*. The goal is to find, for each constraint c

_{i}, the value of

*φ*for which the linear combination

**A**

_{ck}

^{t}·

*φ*is closest to target value b

_{ci}. The system of linear equations is solved using a least squares method which minimizes an error criterion computed from the values of the unknown function at the nodes of the tetrahedral mesh.

This method has two main features that make it particularly useful for interpolating the distortion function. First, it has proved its robustness in many different problems because it is implemented in the Gocad software. Second, designing new constraints for the DSI method is only a matter of translating mathematical formulae into C++ code that can be managed by the software. The mathematical framework of DSI makes it easy to develop new geological constraints as well as geometrical constraints.

In addition to the soft constraints described above, which are honoured in a least squares sense, the DSI method can also honour hard constraints such as: 2

## Interpolating distortion function

The most important step in the editing loop illustrated in Figure 5 is interpolation of the distortion function. This section explains which DSI constraints are set on the distortion function for the interpolation to yield the expected result.

Three constraints are set on distortion property *φ*, defined for each node *α* of the tetrahedral volume Ω. For a constraint set on the distortion function, Equation 1 can be written as
3

A distance constraint transfers the distortion of the control points to the volume. Assuming that control point *P* is located inside tetrahedron 𝒯 with nodes {*α*_{0}, *α*_{1}, *α*_{2}, *α*_{3}}, that {*g*_{0}, g_{1}, g_{2}, g_{3}}, are the barycentric coordinates of *P* inside 𝒯 and that the distortion property value at control point *P* is *φ*_{p}, the DSI constraint associated to control point *P* can be written as:
4

A smoothness constraint ensures the gradient of *φ* is constant over the volume, thus propagating the distortion induced by the control points to the editable volume. On the border of the editable volume, a hard constraint specifies that the distortion is zero. The distortion property is constrained to vary smoothly from the value given by the control point to zero on the border of the editable volume by the constant gradient constraint. This constraint was described in (Mallet 2003; Moyen 2005), and it is set on each pair of tetrahedra inside the editable volume having one common face.

As a tetrahedral model is edited, tetrahedra must not be flattened or inverted. If that were to happen, tetrahedron edges would cross and algebraic tetrahedron volumes would change signs. A non-inversion constraint for the tetrahedra ensures that they do not invert when pushed against the border of the editable volume or against a tetrahedron pinned by another control point. The distorted tetrahedron 𝒯 has vertices *α*_{0}, *α*_{1}, *α*_{2}, *α*_{3} located at positions **x**(*α*_{0}), **x**(*α*_{1}), **x**(*α*_{2}), **x**(*α*_{3}). (vectors **x** are three-dimensional coordinate vectors, *φ* is the intensity of distortion in direction **W**. 𝒯* is the tetrahedron after distortion). The goal of the non-inversion constraint is to prevent any point from getting on the other side of the opposite face. On a tetrahedron, this amounts to preventing the sign of the algebraic volume from inverting.

After distortion, the algebraic volume 𝒱 (𝒯*) of tetrahedron 𝒯* is: 5

If we define three vectors such that: 6 equation (5) simplifies to: 7

Noting that **d**_{1} × **d**_{2} · **d**_{3} = 6𝒱(𝒯) and keeping only first order terms yields:
8

Our goal is to ensure that the sign of 𝒱 (𝒯*) is the same as the sign of 𝒱 (𝒯). One way of achieving this is to constrain *φ* so that:
9
Using equation (8):
10
In terms of a DSI constraint on *φ*:
11

As second-order terms are neglected, b_{c} is chosen as 6𝒱^{2} (𝒯) (*ε*−1) with *ε* ≈ 0.1. The value of *ε* controls how close a tetrahedron can go towards inverting. A low value for *ε*, such as 0.05 for example, will result in nearly flat tetrahedra. A higher value, such as 0.2, will ensure all tetrahedra keep a certain volume, but will stop editing closer to the starting point. This constraint is installed on all tetrahedra inside the editable volume.

The non-inversion constraint has dramatic effects on tetrahedral volumes. In two dimensions, the user can evaluate when triangles will invert and stop before that point. But in three dimensions, it is much more difficult to evaluate the instant when editing goes too far and tetrahedra begin to flip over. The non-inversion constraint stops editing when flipping begins and goes back one step in the editing to a state where the model is valid.

## Consistent editing of geological models

In this part, three fault editing examples are shown. The interpreted fault in the first example is not consistent with seismic data for the study area. Its geometry is therefore edited using the fault-editing tool and the result is a better match to the seismic data. In the second example, a single fault is edited while a porosity simulation is mapped on the tetrahedral model. Steps must be taken so that the final model is still geologically valid, so that it honours both the seismic data and the porosity simulation data. In the last example, a fault branch line is edited and it is shown how editing of a fault branch line in a geologically meaningful way is more complicated than editing a single fault.

### Adjusting a fault to seismic information

After the tetrahedral model was created from interpreted surfaces, fault surfaces within the tetrahedral mesh can be displayed along with seismic data. This can show discrepancies between the seismics and the geometry of the tetrahedral model (Fig. 6b). If this happens, the horizon and fault surfaces used for the creation of the tetrahedral mesh can be edited, and the model can be tessellated again. An alternative time-saving process is to use the fault-editing tool to adjust the geometry of the fault to seismic data.

As it is difficult to show a tetrahedral mesh being distorted, only a slice of the model is displayed. Figure 6a shows a global view of the sliced tetrahedral model and the point of view the other images were taken from. Displaying a section in the seismic cube provides a convenient basis for editing, as the user can move control points along a seismic section.

Once a control point is set on the fault to be edited (Fig. 7a) it can be dragged until the geometry of the fault matches seismic information (Fig. 7b). The control point is constrained to move along the normal to the fault being edited. The vertices of the tetrahedral mesh are moved along this direction minus the normal to the horizons at each vertex, which is computed from the Geochron parameterization pseudo-time function. This ensures that the pseudo-time function is not distorted, thus preserving horizon geometry. Figure 8 shows superimposed slices of the initial tetrahedral model and of the tetrahedral model after the fault was distorted. The fault geometry is much more consistent with seismic information after editing.

### Three-dimensional property mapping

The three-dimensional mapping in the Geochron model introduces another challenge. In the Geochron framework, geological data are mapped from an image of the depositional space, where horizons are flat, onto a faulted and folded geometry. If a property simulation is mapped to the tetrahedral volume, this simulation can have constraints such as values known from well data. These values must remain constant throughout editing.

The texture properties which define the three-dimensional mapping on the tetrahedral volume are attached to the vertices of the tetrahedral mesh. As their geometry is independent of the attached texture, the texture properties must be updated in order to avoid distorting the mapping along with the geometry. To achieve this, the texture properties are copied to a background volume before editing begins. As the vertices in the editable volume move, their texture properties are retrieved from the background volume. The vertices move, but the property mapping is not distorted because their texture properties are kept up-to-date. Figure 9a shows a slice in a tetrahedral model with a mapped porosity simulation in which a fault intersects a channel. When the fault is edited, the texture properties are updated so that the channel is not distorted. Not updating the properties distorts the porosity mapping and the channel is stretched on one side of the fault and compressed on the other (Fig. 9b). When the texture properties are updated, the channel keeps its original shape and is transferred in part from one side of the fault to the other (Fig. 9c). If other faults are detected inside the editable volume, the value of the distortion property on their nodes is fixed to zero so that their geometry is not modified.

### Editing fault branch lines

Editing faults in a geological model is more demanding than just ensuring that the tetrahedra do not flip over. If the object to be edited is not a single fault but a fault branch line, it must be edited in such a way that the main fault does not change shape, so that a secondary fault glides along the main fault without deforming it. When a control point is close to a fault branch line, or even on the branch line, additional constraints ensure that the shape of the main fault is preserved during editing. Branch lines are detected automatically in the initialization phase. If such a line is found inside the editable volume, the normal to the main fault is interpolated as a three-dimensional vector field on the editable volume. If other faults are inside the editable volume, their geometry is fixed. When a control point moves inside the editable volume, the component of the distortion vector **W** which is colinear with the normal to the main fault is removed for each vertex (Fig. 10). This ensures vertices only move along the main fault, preserving the general shape of the main fault while editing the secondary fault. Figure 11 shows an example of branch line editing, in which the branch line is straightened out whilst retaining the basic shape of the main fault.

## Discussion and conclusions

When modelling faulted study areas with tetrahedral meshes, a small geometrical correction to a fault may require a complete rebuild of the model and can be very time-consuming. The fault-editing tool presented in this paper enables direct geometrical adjustments in real-time. The user can control the aspect of the mesh precisely while editing, whereas existing FFD techniques usually involve control lattices or control curves. The magnitude of the editable volume and additional control points can be used to specify precisely which regions of the model can change, and which regions should not be modified.

The fault-editing tool also features safety precautions which prevent the introduction of inconsistencies in the model. Tetrahedra are prevented from flipping over and making the mesh invalid. Vertices are moved along the pseudo-time function of the Geochron parameterization, thus preserving horizon geometry on both sides of edited faults. The texture properties for the Geochron parameterization are updated as the model is edited so that texture-mapped properties are not distorted. And finally, fault branch lines are edited in such a way that the secondary fault glides along the main fault, leaving the general shape of the main fault undisturbed.

The tool provides graphic feedback fast enough to enable comfortable fault editing. Table 1 gives benchmarks of the tool. The time given for each set of nodes inside the editable volume is the average time spent in one tool loop, from the update of the coordinates of the control point to the graphic update (Fig. 5). The frame rate decreases when the editable volume expands, but a volume with about 700 nodes (about 4200 tetrahedra) can be edited in real-time on standard hardware. The frame rate only depends on the number of nodes inside the editable volume: the actual tetrahedral mesh usually is much larger than the sub-volume used for editing (about 50000 tetrahedra for the model tested here).

In some cases, larger modifications have to be made to tetrahedral models. A fault may have to be translated over such a distance that keeping existing tetrahedra would lead to tetrahedra of very poor quality. Faults may also have to be added to or removed from the mesh. These changes cannot be made while keeping the same tessellation: if a fault is added to the mesh, a topological discontinuity has to be introduced in the model, and if a fault is removed, this discontinuity has to be sealed. The next step for fast model editing is to develop tools that enable topological changes in tetrahedral models without having to recompute the whole tessellation.

## Acknowledgments

This work is part of a Ph.D. sponsored by the Association Scientifique pour la Géologie et ses Applications through the gOcad consortium. Consortium members are hereby acknowledged. Thanks to EarthDecision for providing the gOcad development environment. Thanks to T. Frank for visualization tools, R. Moyen for the Geochron parameterization and B. Leflon for the porosity simulation.

- © The Geological Society of London 2007