A Novel Method for Locally Updating an Earth Model While Geosteering

Author(s)
Erich Suter^{1,2},
Eric Cayeux^{1},
Helmer A. Friis^{1},
Terje Kårstad^{2},
Alejandro Escalona^{2},
Erlend Vefring^{1}

Affiliation(s)

^{1}
IRIS (International Research Institute of Stavanger)/DrillWell (The Drilling and Well Centre for Improved Recovery), Stavanger, Norway.

^{2}
University of Stavanger, Stavanger, Norway.

ABSTRACT

Earth models are important tools for support of decision making processes for optimal exploitation of subsurface resources. For geosteering and other real-time processes where time is a major constraint, effective model management is decisive for optimal decision support. During drilling, subsurface information is received which should optimally be used to modify the 3D earth model. Today this model is typically not altered during the operation. We discuss the principles of a novel method that enables a populated earth model grid to be locally modified when the topology (connectivity) of the geological structure is locally altered. The method also allows local updates of the grid resolution. The modelled volume is split into closed regions by the structural model. Each region is individually discretized and obtains its own subgrid. Properties are stored in separate functions, e.g. for each layer, and transferred into each subgrid via a mapping. A local update of the geological structure implies that only subgrids in regions that are directly affected by the updated structure must be discarded and rebuilt, and the rest of the populated earth model grid is retained. Our focus is on decision support for optimal well placement while geosteering. The proposed method aims to manage multiple model realizations that are never fixed and always locally updated with the most recent measurements and interpretations in real-time, and where each realization is always kept at an optimal resolution.

Earth models are important tools for support of decision making processes for optimal exploitation of subsurface resources. For geosteering and other real-time processes where time is a major constraint, effective model management is decisive for optimal decision support. During drilling, subsurface information is received which should optimally be used to modify the 3D earth model. Today this model is typically not altered during the operation. We discuss the principles of a novel method that enables a populated earth model grid to be locally modified when the topology (connectivity) of the geological structure is locally altered. The method also allows local updates of the grid resolution. The modelled volume is split into closed regions by the structural model. Each region is individually discretized and obtains its own subgrid. Properties are stored in separate functions, e.g. for each layer, and transferred into each subgrid via a mapping. A local update of the geological structure implies that only subgrids in regions that are directly affected by the updated structure must be discarded and rebuilt, and the rest of the populated earth model grid is retained. Our focus is on decision support for optimal well placement while geosteering. The proposed method aims to manage multiple model realizations that are never fixed and always locally updated with the most recent measurements and interpretations in real-time, and where each realization is always kept at an optimal resolution.

KEYWORDS

Geosteering, Optimal Well Placement, Earth Modelling, Local Update, Geological Structure, Connectivity, Uncertainty

Geosteering, Optimal Well Placement, Earth Modelling, Local Update, Geological Structure, Connectivity, Uncertainty

1. Introduction

The planning of a new well is based on information known prior to drilling such as surface seismic and offset wells. The interpretation of the information is captured in an earth model that is used to support decision processes while drilling. Due to incomplete knowledge of the subsurface, the interpretation is always burdened with uncertainty. During the geosteering process, the trajectory of the planned well is adjusted based on measurements obtained from logging tools during the ongoing drilling operation. The new information reduces uncertainty and allows revisions of the geological interpretations made prior to the drilling operation. This requires effective interpretation, integration and utilisation of the new information within the timeframe set by the on-going drilling operation.

However, commercially available three-dimensional earth modelling tools have limited capabilities for real-time updates. Model modifications are complex and labour intensive (see Section A.3 in the Appendix for details), and the time needed for updating the model exceeds the time available during drilling operations. This results in sub-optimal utilization of the measurements obtained while drilling, and is a large drawback for decision making processes that require the most current and precise information.

In Section A.1 in the Appendix, two examples of how a more effective decision loop can contribute to safer and more effective drilling and geosteering are discussed. In Section A.2 in the Appendix, a potential future workflow for effective decision support is indicated. It aims to support decisions for safer and faster drilling, increased future production and reduced drilling costs. The aim of the method suggested in this paper is to enable effective local updates of the earth model grid and local scale uncertainty management around and ahead of the bit, including uncertainties in the topology (connectivity) of the geological structure. Effective earth model management is essential to shorten the geosteering decision loop.

1.1. Current Methods for Geosteering

The recently developed ultra-deep directional electro-magnetic (deep EM) tech- nology is sensitive to resistivity contrasts up to tens of meters around the wellbore (see e.g. [1] [2] [3] ). Deep EM measurements is a technological step change that gives more information about geological structures, fluid contacts and reservoir properties deeper into the formation than traditional Logging While Drilling (LWD) measurements, even ahead of the bit [4] . The technology contributes to closing the gap between seismic scale measurements and well scale measurements. Next, we briefly review two current workflows for updating an earth model in real-time, and how the model is used for decision support for optimal well placement.

In [2] , it is explained that the full-field 3D model typically ignores or blends small-scale features such as faults and facies changes that can impact drilling. For drilling support, the full-field model can be refined around the well or a separate sector model on a fine scale around the well can be constructed. The detailed model is used to extract a 2D section along the planned well path, containing geological structure and properties, which is used for drilling support. During the drilling process, four different interpretation methods are used for local scale structural interpretation, such as dip interpretation and remote boundary detection. Moreover, the results of the inversion of the deep EM measurements acquired while drilling are visualized together with and compared with the seismic. After drilling is completed the 2D well-scale structural model and petrophysical model are modified, as basis for the following update of the 3D full-field model.

In the geosteering workflow discussed in [3] , a 3D sector model around the planned well is generated prior to drilling. The volumetrically smaller model allows an increase in the horizontal and vertical resolution, and the high-resolu- tion model is populated based on full-field petrophysical properties. It is emphasized that 3D modelling enables geologists to envision the result also of lateral changes in the well trajectory, as opposed to the simple vertical changes that 2D models allow, and thus to mitigate misinterpretations in the imaging of facies models. During drilling, the geometries of conformable layers can be updated to locally adjust their depths. Also properties are updated. Geosteering decisions are supported by sharing the continuously updated 3D view among the members of the decision making tea m.

In the mind of the interpreters, geological modelling is a process that takes place in three-dimensional space. If a 2D model is used as the main tool for supporting decisions, important information may be ignored. 3D sector models based on the standard tools used in the industry today do not allow effective updates of more complex geological structures within the short time available during drilling (see discussion in Section A.3 in the Appendix). In particular when drilling in complex geology, effective handling of complex structural uncertainties can be decisive for the outcome of the geosteering operation.

1.2. Main Contributions of This Paper

In the tools and methods for earth modelling that are today standard in the industry, a global grid is used to capture both properties and structure. This technique implies that the entire grid is invalidated even by small changes in the geological structure (see Section A.3 in the Appendix).

In the suggested approach, we use the geological structure to split the modelled volume into a set of regions. Each region is individually handled and discretized without reference to other regions. Moreover, properties are not handled in a single grid but in individual property functions that each represents only a small part of the subsurface. Each property function is handled separately without reference to other property functions. A region and a property function are linked via a mapping. All mappings are also handled separately, without reference to other mappings. As a result, local updates of the geometries and topology of the geological structures that are captured in the earth model grid, as well as local updates of the grid resolution, can be performed locally within a time frame independent of the size of the grid.

When using conventional earth modelling tools, local updates of the structural topology, e.g. insertion or removal of a subseismic fault or layer, require the grid to be globally regenerated and populated. Enabling local updates of the grid when the structural topology is locally modified aims to drastically shorten the time required for such updates. Local control with grid resolution (for each individual rock property) aims at improving the control with the trade-off between numerical accuracy and the time required for handling the grid(s). Moreover, local updates aim to enable effective, local scale uncertainty handling. When uncertainties are handled locally, less time is required for updating the grid. Furthermore, separate handling of each region opens up for parallel computer implementations. Such developments are important to speed up and streamline the earth modelling process when targeting real-time workflows. Rudimentary discussions of the strategy were presented in [5] [6] .

It is emphasized that the proposed method is at an early stage of development, yet too immature to handle complex geological configurations and uncertainties. Therefore, our aim in this article is not to describe a complete method for geosteering support, but to discuss principles for how a populated earth model grid can be effectively modified when the topology (connectivity) of the geological structure is locally altered. Structural modelling can be a complex task, and could e.g. be handled in combination with the methodology described in [7] (see Section A.2 in the Appendix). The principles of the suggested method are demonstrated in 2D using a simplified geological structure.

We start by providing an overview of the proposed method in Section 2. In Section 3, we show how the strategy is applied in a geological setting containing layers and faults. Next, the required input for the grid construction process is discussed in Section 4. In Section 5, the construction of a populated grid is described, including pseudo code. In Section 6, the procedure for locally updating the grid is discussed, and in Section 7 examples are provided. In Section 8, the mapping is discussed in more detail. Section 9 summarizes the paper. In the Appendix, two examples that highlight the need and potential benefits of improved workflows for drilling are presented in Section A.1. In Section A.2, we indicate a potential future geosteering workflow. The conventional 3D earth modelling methodology used in the industry today is reviewed in Section A.3. In Section A.4, more recent technological developments within 3D earth modelling are discussed. Finally, in Section A.5, various strategies for gridding of the earth model and their implications for model management are discussed.

2. Overview of the Proposed Method

A main component of our strategy is to separate the modelled volume into disjoint regions that are individually handled. The regions are defined using the geological interfaces in the structural model (see Section 4.1). Each region R gets its own grid, or even a set of grids if different properties should be discretized at different resolutions. A grid for a particular region is called a subgrid G_{Z}, and it should be generated at the resolution and quality decided by the application of the model. The subscript Z indicates the properties of each specific subgrid, e.g. its resolution. The resolution of each subgrid is in general independent of the other subgrids for this region and of subgrids for neighbouring regions. However, if two subgrids in neighbouring regions are connected by the sharing of faces along their common boundary, dependencies between the two neighbouring subgrids are introduced. Typically, the subgrid is an unstructured grid.

The properties are not represented in a single global grid, but in a set of individual property functions (see Section 4.2). Each property function represents a specific property, e.g. porosity, for a geologically defined small rock volume, e.g. a depositional layer. A property, say porosity, can be represented at multiple resolutions by application of multiple functions (or a multi-resolution function). Control with resolution of both the subgrid and the property function provides a large degree of flexibility for locally controlling the earth model resolution. A given property could be represented at varying resolution depending on the location in the model, and different properties can be represented at different resolutions within the same region.

A mapping links a region with the corresponding part of a property function, and allows population of the subgrid with properties from the property function. For example, if a layer is split by faults into multiple regions, a set of mappings link each region with the corresponding part of a property function. (Note that a “function” and a “mapping” is the same in the mathematical terminology. In this article we let a “function” represent values of a rock property. A “mapping” is used to link a subgrid with a property function.)

When the geological structure is modified, only the regions and subgrids that are in direct contact with the modified parts of the structure must be discarded and rebuilt (see Section 6). The rest of the existing subgrids are retained. Moreover, a local update of the grid resolution implies that only the subgrids in the affected regions are discarded and that new subgrids are established for these regions. Again, the rest of the existing subgrids are kept. An update of a property function only implies that a small set of the existing subgrids must be repopulated. The method is independent of the particular strategy used for structural modelling as long as a structural model (see Section 4.1) can be extracted.

To update the populated grid as a fully local operation requires that all involved data structures can be locally updated. In our strategy, the grid, the mappings and the property functions are data structures that represent information for only a small part of the subsurface. Local updates of the populated grid are thus achieved by avoiding the use of any globally defined data structure that cannot be locally updated. In this paper, simple examples of local updates of the fault network and the stratigraphy are shown.

In the proposed method, each individual region is assigned its own mapping that links a subgrid for the region with a property function for a stratigraphic layer. First, this individual handling of each layer aims to enable more flexibility for controlling and locally altering the interpretation of the stratigraphic record. Second, a region has a relatively simple shape and does not contain any geological interfaces in its interior. As a result, the mapping can be kept simple without the complexity of taking interior geometries into consideration. A particularly attractive choice of mapping is therefore one that only requires knowledge of the geometric boundaries of the region and the geometric boundaries of the corresponding parameter domain of the property function (see Section 5.2) for its construction. Such mappings have been developed, verified, optimized and documented elsewhere in the scientific community, see Section 8. The use of general purpose mappings also allows us to capitalize from future developments within the field. Moreover, also gridding strategies for the discretization of each region may benefit when regions do not contain any interior geometries that the grid must honour. Gridding strategies are discussed in more detail in Section A.5 in the Appendix.

3. Application of Method to Layered Media with Faults

Next, we exemplify the approach for a geological setting comprising layered depositions with faults. A region is then assumed to be a part of a layer within a fault block (note, however, that a region could also be defined differently within the proposed framework). In this setting, a layer therefore consists of a set of regions separated by faults. Each property is handled in a separate property function Φ that represents a property for a layer L. But more flexible designs are also possible using the proposed framework. For example, a single property function could cover a set of layers, a specific fault block, or some other geologically defined volume of rock.

When a structural element such as a fault or another geological interface is updated, i.e. geometrically perturbed, removed, inserted, or its connectivity with other structural elements is modified, only the regions whose boundaries are affected by the structural update are involved in the update of the grid. The principle is indicated in Figure 1, where the insertion of faults in the interior of the fault block denoted “A” only affects the subgrids in this particular fault block.

In Figure 2, the population of an earth model grid is explained. The top layer L^{1} is split into two regions and by the fault, and similarly for the bottom layer L^{2}. Each of the four regions is individually discretized and we obtain four subgrids. Two properties, represented by the two property functions Φ^{1} and Φ^{2} for the two layers L^{1} and L^{2} respectively, are linked with each region using four corresponding mappings. Each mapping links a region with the corresponding part of the parameter domain of the property function. For example, the subgrid in is linked with the part of the parameter domain of Φ^{1} via. Informally, we can say that the mapping deforms the interior of into the shape of. This allows the subgrid to be populated by sampling values from Φ^{1}. This deformation is further explained in Section 5. The effect of deforming e.g. into the shape of is that it is elongated in horizontal direction and squeezed together in vertical direction. The white near-vertical

Figure 1. At the top is a model representing a faulted reservoir with alternating sands and shales. At the bottom the initial model is locally updated by inserting three new vertical faults in fault block A.

Figure 2. A structural model with two depositional layers (L^{1} and L^{2}) split by a fault contains four regions R that are individually gridded. Property values are interpolated from the property functions Φ^{1} and Φ^{2} to the right and transferred into the grid in each region by the use of mappings (). The result is a populated grid as shown at the bottom. The colors in the figure indicate property values.

line that splits Φ^{1} into and is the image of a part of the fault, correspondingly for Φ^{2}. Note that the direction of the arrows in the figure indicates how property values evaluated from the property functions are transferred to the subgrids. The mappings are in fact directed in the opposite direction (see Section 5.2). Also note that we have zoomed in on the fault so that only parts of the layers are shown, they extend further both to the right and the left. As a result, only the corresponding parts of the properties around the fault are shown in the populated earth model at the bottom of the figure.

In the simple example in Figure 2, created using a rudimentary software prototype, the property functions used to populate each layer are identical. Clearly, for realistic modelling, separate property functions will exist for each layer.

4. Input for Earth Model Construction

The required input for generating a populated earth model grid is a structural model and associated property functions.

4.1. The Structural Model

The structural model captures an interpretation of the geological structure at a specific resolution. The geometries in the structural model create a partition of space into disjoint polytopic regions (polygons in 2D). The boundary Ω_{R} of each region is represented by a set of geometric patches, where each patch is a part of a geological interface. For example, a region can be a layer in a fault block and be bounded by parts of fault geometries and stratigraphic interfaces. Each region is a continuous closed volume, and it is the smallest volumetric object in the partition of space as it cannot be subdivided by any other geometric patch in. Geometric patches do not cross, and each patch stop into another patch (including patches representing the model boundary). A layer can e.g. be locally split in two as a result of hiatus or erosion so that it has zero thickness. Then the two parts are handled as two separate regions. The part of the layer with zero thickness is not represented by a region and is therefore not part of the earth model grid. Neither faults nor layers are required to cross the entire model. The described structural model is similar to the sealed structural model described in 3D in e.g. [8] [9] [10] [11] .

In the initial strategy described in this paper there are no geological interfaces in the interior of a region. This implies that faults must terminate into other geological interfaces. To enable faults to terminate in the interior of a layer, adjustments are required when mapping properties to subgrids. Furthermore, when discretizing a region, this part of the fault must be taken into consideration.

An important future target is to allow integration of the proposed method with conventional tools for 3D earth modelling. The rules for the described structural model allow us to import from and possibly integrate with the standard tools used in the industry today. Potentially, the rules could also be adapted to other earth modelling approaches.

4.2. Property Functions

For each layer L we have at least one property function Φ. The function represents a physical property, say porosity, density, velocity or saturation. A bivariate scalar property function is defined by over a parameter domain D^{Φ}, for example the unit square.

Each property function is independent of the other property functions, and can be managed at its own resolution. For example, porosity could be represented at a finer resolution than saturation. The same property can be represented by property functions at different resolutions. Property functions can be defined over parameter domains D^{Φ} of different shapes. This may be useful for more optimal handling of facies and property distributions of more complex geological shapes. Different functions, typically over the same parameter domain and with the same resolution, can be used to represent multiple realizations of the same property.

A set of functions can be constructed from the property representation of existing earth models. Each function is then handled separately, which provides flexibility e.g. for locally updating the stratigraphic record and for handling each property at its own resolution. However, to modify existing properties in a geologically reasonable manner to match an updated interpretation of stratigraphic interfaces is not straight-forward. Property functions are mathematical constructions without the burdening requirement of carrying direct geological meaning. Geological meaning is only assumed after the functions have been used to populate the grid. This provides an extra degree of freedom such that the functions can be set up in ways that are mathematically convenient. On the other hand, it requires that the mapping of properties from the functions to the grid ensures that the geological meaning is restored when the grid is populated from the property functions. This issue is discussed in Section 5.3.

The property functions can be set up in ways that are mathematically convenient. Φ could in principle be any type of function, as long as it can be evaluated everywhere within D^{Φ}. Potential strategies include that the function is represented over a grid, or is a uniform value as in [12] . The function could also be a complex analytical function, or a multi-resolution function as suggested in [5] [13] [14] . A property function can be constructed from grid-based property distributions that are imported from external tools. It could also be derived from the interpretation of well logs in real-time.

In the example in Figure 2, the geological structure is imported from Petrel and the property functions are generated from an imported Petrel grid. Currently, the construction of the property functions take place in a simple manner and is only intended for demonstrating principles; by identifying a layer in an imported corner-point grid, the value in each cell in the grid is simply transferred into a regular grid with the same topology as the corner-point grid (i.e. the same number of cells in all directions). The populated regular grid then constitutes a property function. But in this manner, e.g. collapsed cells are not properly handled. An improved procedure would sample the properties from the geological space where the corner-point grid has been adapted to the geological structure.

5. Construction of a Populated Grid

The construction of a populated grid takes place by first identifying the closed regions in to populate. Each region R is handled independently of the other regions. Each region can contain multiple subgrids at different resolutions and with different numerical qualities, allowing each property within each region to be handled separately at its own resolution. There are three main steps required for populating a region R with a property, namely 1) construction of a grid in the region with a resolution and quality adapted to the property in question, 2) construction of a mapping to link the region with the corresponding domain of a property function, and 3) populating the grid by transferring values interpolated from the property function using the mapping. Once a boundary polygon for R has been created, step 1) and 2) are independent and could be completed in parallel.

Next, a procedure for populating a subgrid with a property is described. The procedure is repeated for all subgrids within all required regions and for all required properties.

5.1. Construction of Boundary Polygon and Subgrid

A polygon P_{R} representing the boundary Ω_{R} of the region R is required both for gridding and for construction of the mapping. P_{R} is constructed by extracting and joining the geometric patches in that together constitute the boundary of R. Details for its construction are discussed in Section 5.2. The resolution of P_{R} controls the resolution of the subgrid as well as the computational efficiency of the evaluation of the mapping. Once P_{R} is constructed, a subgrid G_{Z} at an appropriate resolution and quality to discretize the interior of R can be constructed using a suitable grid generator. If such a subgrid already exists, it may be reused.

5.2. Construction of the Mapping and Population of the Subgrid

For linking the property function Φ for a layer with the region R, we apply a mapping f. First, let us assume that we have an unfaulted layer L, represented by a single region R. The exact shape of L is a matter of geological interpretation, so we assume that we have arrived at an alternative we call L^{1}. An informal and intuitive way to understand the mapping of properties from a property function to a layer is to imagine that the shape of Φ is deformed into the shape of L^{1} as shown in Figure 3. Then the mapping ensures that the interior follows. If the geological structure is modified from L^{1} to L^{2}, a new mapping f^{2} is generated to deform Φ to fit within L^{2} as shown in the figure. The nodes of the polygonal boundaries of D^{Φ}, L^{1} and L^{2} are also indicated. Note that the direction of the arrows also in this figure indicates how property values are transferred from the property function to the layer.

Next, we describe a preliminary strategy for handling faults. The approach demonstrates how property functions are used to populate faulted layers, but further developments are required for proper handling of updates in the geological structure when addressing realistic model alterations while drilling (see discussion in Section 5.3). Φ represents the property function for the complete layer L, and we assume that n − 1 faults separate L into n regions, i = 1, ∙∙∙, n.

Figure 3. Two mappings f^{1} and f^{2} can informally be said to deform the property function Φ to the right to fit within any of the two alternative interpretations L^{1} and L^{2} of the shape of a sedimentary layer L. P_{R}_{,t} is the geometric patch representing the top stratigraphic boundary of L^{1}, whereas P_{D}_{,t} is the geometric patch representing the top boundary of the parameter domain D^{Φ} of Φ. The other parts of the boundaries of L^{1} and D^{Φ} are correspondingly marked.

A fault may displace multiple layers, but here we only consider the part of a fault that affects the layer that is currently being populated. The parameter domain D^{Φ} of Φ must then be correspondingly split into n subdomains, for i = 1, ∙∙∙, n. Now each in L can be associated with each corresponding. As an example, consider Figure 2 where L^{1} is split into two regions and by the fault. The parameter domain of the property function Φ^{1} is then split into two subdomains and by a curve which is the image in D^{Φ} of the part of the fault that splits the layer. In the example the curve is an almost vertical line, but it can also have a more complex geometry.

The mapping is on the form from a source polygon Ω_{S} to a target polygon Ω_{T}. In the described strategy, Ω_{S} is the polygonal boundary P_{R} of R and Ω_{T} is the polygonal boundary P_{D} of D, for any and corresponding. We call the mapping, and it links any pair R and D so that properties represented by property functions can be used to populate subgrids in R.

The mapping we use (see Section 8) requires that the source polygon P_{R} and the target polygon P_{D} are topologically equivalent, so that P_{R} can be deformed into P_{D}. Here, topological equivalence means that the two polygons have the same number of nodes and that the nodes have the same ordering. Next, a recipe to construct P_{R} and P_{D} is discussed.

In the same manner as P_{R} was constructed (see Section 5.1), P_{D} is constructed by joining the geometric patches that constitute the boundary of D. We let the patches constituting P_{R} and P_{D} be pairwise associated; the top boundary of R is associated with the top boundary of D, the right hand side part of the boundary of R is associated with the right hand side part of the boundary of D, and so on. Each geometric patch in the pairwise association must be topologically equivalent. This is obtained by inserting new nodes into either of the two patches, but typically into the patch being part of the boundary of P_{D}. This is because P_{D} generally has the lowest geometric resolution, as a result of the simple quadratic shape of D^{Φ}. As an example, consider Figure 3 where R is the entire L^{1} and correspondingly, D is the domain D^{Φ} of the entire property function. For this example, L^{1} is thus considered to be the part of a layer in the interior of a fault block. P_{R}_{,t} is the top stratigraphic interface of R, P_{R}_{,l} is the left hand side boundary of R, P_{R}_{,b} is its bottom stratigraphic interface, whereas P_{R}_{,r} is the right hand side boundary of R. Correspondingly, P_{D}_{,t} is the top boundary of D, P_{D}_{,l} is its left hand side boundary, P_{D}_{,b} is its bottom boundary, whereas P_{D}_{,r} is the right hand side boundary of D. P_{D}_{,t} is associated with P_{R}_{,t}, P_{D}_{,l} is associated with P_{R}_{,l}, P_{D}_{,b} is associated with P_{R}_{,b}, and P_{D}_{,r} is associated with P_{R}_{,r}. Note that P_{R}_{,t} and P_{R}_{,b} can have a different number of nodes, and there are no constraints regarding the placement of these nodes. Thus, the top and bottom boundaries of a layer can have different geometric resolutions.

In Figure 3, we see that P_{R}_{,t} contains five nodes. D is a square, so P_{D}_{,t} is a straight line represented only by its two end points at the upper left corner and upper right corner of the square. To ensure that P_{R}_{,t} and P_{D}_{,t} are topologically equivalent, P_{D}_{,t} is refined by inserting new nodes (see Figure 3). We let P_{R}_{,t} and P_{D}_{,t} be parameterized by normalized arc length. The three new nodes in P_{D}_{,t} are inserted at the same parameter values s_{i}, for i = 1, 2, 3, as where P_{R}_{,t} has interior nodes. A similar refinement is applied to P_{D}_{,b} with respect to its counterpart P_{R}_{,b}. Neither P_{R}_{,l} nor P_{R}_{,r} have interior nodes, so neither P_{D}_{,l} nor P_{D}_{,r} need further refinement. Now that all pairwise associated geometric patches are topologically equivalent, we join {P_{D}_{,t}, P_{D}_{,l}, P_{D}_{,b}, P_{D}_{,r}} to form P_{D}. Each of the patches must be oriented so that a valid polygon is formed, and P_{D} must be have the same orientation (clockwise or anticlockwise) as P_{R}. Now P_{R} and P_{D} are topologically equivalent.

can be constructed from the source and target polygons P_{R} and P_{D}, respectively (see Section 8, where a particular type of mapping is discussed). The mapping allows any point to be mapped to its corresponding point, as. A property value can then be evaluated by interpolation of Φ, and the value can be used to populate the subgrid G_{Z} that covers R. Typically, all nodes or cell centres in G_{Z} are given values in this manner. The procedure is called “backward mapping” and is well known from image warping, see e.g. [15] . Note that for handling of more complex fault networks, e.g. where a fault terminates into another fault, the strategy must be generalized to handle more complex topological relationships in the geological structure.

Pseudo code for grid construction is provided in Algorithm 1. The input is a list of regions in where subgrids should be constructed, together with associated property functions. The output is a populated subgrid for each region. The algorithm handles each region independently of the other regions.

5.3. Handling of Properties When the Structure Is Locally Modified

Existing properties are the results of modelling prior to drilling, or of predictions made earlier during the drilling process. But property predictions are burdened by uncertainties, and predictions of the geological structures and properties that are constrained by the most recent measurements and interpretations obtained during the ongoing drilling process are important to support the geosteering decision process.

Algorithm 1. Part 1: Generation of subgrids for a set of regions.

Algorithm 1. Part 2.

Algorithm 1. Part 3.

More complex property distributions are time consuming to generate for the entire grid. Therefore, to reduce the time for locally modifying the structural model, existing properties represented in property functions should be reused if possible. In a property distribution captured in a grid, e.g. facies objects are implicitly represented. The sizes of the objects and the distances between the objects carry geological meaning, and the property distribution is matched to the shape of the layer it originally populated. In the parts of the model where there is no alteration of the geological structure during a local model update, which is typically a very large portion of the model, existing properties could be reused because their geological meaning is not altered. But when a property function is used to populate a subgrid that is adapted to a locally updated geological structure, care must be taken to ensure that the original geological meaning is (approximately) restored so that the updated model is geologically reasonable. This depends on the complexity of the property function, as well as on the complexity of the structural update, e.g. the amount of deformation of an existing layer. A major revision of the structure at local scale may render existing complex property representations locally inapplicable, so that new properties that respect the new measurements should be generated at the local scale. This depends on the application of the model. Such issues have not yet been properly addressed. More basic property representations, such as constant-valued representations, are easier to handle. In the method described in [7] , the aim is to allow automatic modification of both properties and the geological structure.

While drilling, when time is limited and an updated model is urgently needed for decision support, methods that provide approximate solutions within short time are typically preferred to methods that provide more exact solutions long after a decision has been made and the drilling operation has progressed. Therefore, careful consideration with respect to the modelling requirements set by the decision process is required. The time needed to generate a model-based decision recommendation must be weighed against the desired quality of the recommendation.

Finally, it is important to note that the current focus of the proposed strategy is on local updates to support decisions for the well being drilled, not on updating the model globally during the drilling process. In [16] , a novel strategy for multi-resolution earth model gridding is described. It is based on the methodology described in this paper, and one of its objectives is to aid in the localization of the updates by arranging the regions in ^{ }in a hierarchical manner. Outside a region at any level of detail in the hierarchy, the model will not be modified. Moreover, the approach in [16] aims to increase the modelling efficiency while drilling by locally controlling the resolution of the geological structure and the grid.

6. Local Updates of the Earth Model Grid

With Algorithm 1 in mind, we explain the general method to locally update the geological structure in a populated grid. Let be the pre-update structural model, while the corresponding post-update structural model is denoted. is obtained by locally updating the structural connectivity and/or geometry in. is the set of regions in that are affected by the local structural update, whereas is the set of regions in that are established or deformed in the update.

Thus, a local update of implies that one or more regions in are affected; they are either deformed or new regions have taken their place. In both cases, is the set of regions in that must be attended. The rest of the regions in already exist in. The general method is to first remove the subgrids for the regions in. To reestablish the global earth model grid, Algorithm 1 is run to generate populated subgrids only for the regions in. If a region is only slightly deformed during the structural update, it could be possible to deform its existing populated subgrids.

Structural modelling can be challenging, and the particular technique to locally updating is not considered here. For example, such updates could be controlled in combination with an external process where multiple model realizations are constrained by new measurements obtained while drilling as discussed in [7] . In this article we focus on how the populated global grid can be locally (rather than globally) modified when a local update of the structural model is performed.

7. Examples of Local Updates of an Earth Model Grid

Next, some basic examples are shown that demonstrate the principles of how the grid is locally modified when the structural topology is locally altered. Such local updates cannot be performed using the methodologies on which the commercially available earth modelling tools are based. It requires further developments to handle more realistic geological configurations and uncertainties than those shown in the examples.

In Figure 4, a new layer is inserted as a local model update. The at the top has two layers L_{1} and L_{2}, and a fault splits the layers into totally four regions

Figure 4. An example of a local stratigraphic update. Top: an initial model with two layers L_{1} and L_{2} split by a fault. Bottom: a pinch-out in yellow marked L_{3} is inserted. The update affects only the two regions that together constitute L_{1}, the subgrids in the two regions for L_{2} are retained. The colors in the figure indicate property values.

that are separately gridded and populated with properties. In the updated at the bottom of the figure, the volume formerly covered by the top layer L_{1} is now occupied by L_{1} and the pinch-out L_{3} in yellow. consists of the two regions representing L_{1}. The subgrids for were discarded while the subgrids for L_{2} were kept. consists of the two updated regions for L_{1} as well as the two new regions for L_{3}. A new property function for L_{3} was generated before was sent for gridding using Algorithm 1.

In the example in Figure 4 the property functions used to populate each layer are identical just as in Figure 2, so that the effect of a local structural update can be examined. The procedure for populating the subgrid in any R using properties represented in the corresponding D is based on deformation of the boundary polygon of R into the shape of the boundary polygon of D (see Section 5.2). In Figure 4, the thickness of the layer L_{1} is locally decreased in the local update. The figure shows that the polygon deformation implies that the property representation in the interior of L_{1} is correspondingly deformed, so that the vertical distances between the objects that are implicitly represented in the grid are decreased. See Section 5.3 for further discussion.

In Figure 5, an example is shown where the grid resolution and the fault network are locally modified within two separate fault blocks. The synthetic model consists of alternating sands (in orange) and shales (in gray), and all faults are vertical. The five layers in the middle of the model contain the faulted reservoir rocks. The two layers at the top and the two layers at the bottom of the model do not contain faults. Let us assume that these four boundary layers are of less interest for the modelling purpose at hand. Their boundary geometries were originally at the same resolution as the geometries closer to the reservoir, but they were coarsened to allow coarser subgrids. A coarser subgrid requires less computational time for its generation and population, but comes at the expense that fewer details can be captured in the grid. The grid resolution is generally finer within the sand layers than within the shale layers. This is useful if variations in the rock heterogeneity should be captured at different levels of detail for the two facies.

At the top is the initial model, it has finer grid resolution within the fault block denoted “A”. The middle and the bottom models in the figure were obtained by locally updating the model at the top. First, the grid resolution in the interior of fault block “A” was coarsened as a local model update. Then, for each of the two models, fault block “B” was modified in a local update by inserting new vertical faults with small displacements (a local update of the structural topology). The subgrids for all new regions within the fault block were constructed at a fine resolution. In fault block “B”, the number of new faults and their respective displacements are different for each of the two updated models. The assumption is that the two updated models both represent possible realizations of the fault network, but with significant differences in the reservoir connectivity and with corresponding consequences for the resulting flow patterns. Within the fault block, parameters such as number of faults, fault location and fault

Figure 5. Top: an initial model with faulted reservoir rocks. The grid resolution varies across the model. Middle and bottom: the fault network and the grid resolution were locally updated within fault blocks A and B. Gray indicates shale, orange indicates sand rocks.

displacement were used to automatically update the fault network in the interior of the fault block. The vertical faults were inserted using a simple fault operator that locally moved the stratigraphic interfaces vertically up or down. All model updates were accomplished in a fully automatic fashion that also supports uncertainty modelling.

8. The Selected Mapping

Mapping between polytopes (in 2D they are planar polygonal domains) is a well known problem in computer graphics and geometric modelling. Numerous mappings with various numerical properties exist in the literature. One group of mappings is based on barycentric coordinates. Barycentric coordinates are frequently used to represent a point in the interior of a polygon as an affine combination of the nodes of the polygon. The coordinates are unique for triangles and tetrahedra, but for arbitrary simple polygons there are many generalizations that each has a different set of numerical properties. The examples shown in this paper are generated using a mapping based on mean value coordinates (MVC). It was first described in [17] , while pseudo-code can be found in [15] . It has also been extended to 3D, see [18] . In [19] , a 2D mesh that conforms to the geological interfaces is used for seismic restoration. When the structure is deformed in the restoration process, the MVC-based mapping ensures that the positions of the nodes in the interior of the triangulated mesh follow the restored interfaces. Then the properties stored in the grid are always available during restoration.

The general procedure to populate a subgrid with values from a property function was described in Section 5.2. As shown there, the mapping from a source polygon to a target polygon takes the form. When using the MVC- based mapping, barycentric coordinates for a given source point with respect to P_{R} are calculated (see [15] for details). Then the coordinates are kept fixed while P_{R} is deformed into P_{D}. The target point can then be evaluated by applying the barycentric coordinates with respect to P_{D}.

The MVC-based mapping has many favourable properties. One of the most prominent is its computational efficiency; it has a closed form and it can easily be parallellized, allowing multiple property values to be interpolated simultaneously. It is not based on a grid and thus independent of the resolution of such a grid. However, the mapping discussed in [15] [17] is not bijective. A bijective mapping that also enable the application of source and target polygons of any shape may allow e.g. to address facies distributions of complex shapes in an easier manner. In [20] , smooth and bijective mappings that also extend to 3D are discussed.

9. Conclusions

Decision making to optimize the exploitation of subsurface resources is challenging in particular when targeting more complex fields and reservoirs. Three- dimensional grid based earth models are routinely used for decision support in workflows where time is not a major constraint. In geosteering operations, new measurements received while drilling should be used to modify the pre-drill interpretation captured in the earth model and support right-time decision making based on the most recent measurements and interpretations. But today’s methods fall short in the attempt to update the model in a timely manner.

We have described a novel method that aims to enable real-time local updates of the topology of the geological structure and the grid resolution in a 3D earth model while drilling. The main principles for locally updating the populated grid when the geological structure is modified have been discussed, but the developments have not yet come far enough to address more realistic geological problems. Examples have been shown for basic cases. If the method can be further developed to handle realistic geology, it would offer a number of advantages for increased grid handling efficiency:

・ The grid can be locally modified when the topology (connectivity) or geometry of the structural model is updated.

・ The grid resolution can be locally updated.

・ Each property can be handled at its own resolution.

・ Grid handling can be parallellized to further reduce the computational time required for managing the model.

・ Uncertainties in the structure and properties can be handled at a local scale while the rest of the model is kept unaltered.

The method is described in 2D. The mapping described in Section 8 has a three-dimensional counterpart. Therefore the basic principles for locally updating the grid as discussed in this article, applied to a simplified geological structure, should be possible to extend to 3D. Future work should first and foremost focus on management of more complex and realistic geology. The potential for improved modelling efficiency provided by local model updates, control with model resolution and parallel processing should be further explored. Moreover, modelling for real-time applications could be achieved by a high degree of automation so that the need for time-consuming manual work is minimized. Automation should be addressed by algorithms that update the model, both structure and properties, in a geologically realistic manner. Such algorithms should be controlled by geological parameters that are also intuitive to geoscience experts, to allow capturing of the geological reasoning behind the interpretation directly in the model.

Modelling efficiency is important to support optimal well placement while geosteering. The suggested approach could potentially form an essential part of a methodology for effective uncertainty modelling where an ensemble of three- dimensional earth model realizations is always kept up-to-date with the most recent measurements, interpretations and uncertainty estimates, and is always at an optimal resolution, even during real-time workflows.

Acknowledgements

The authors acknowledge the Research Council of Norway, ConocoPhillips, AkerBP, Lundin Norway, Statoil and Wintershall for financing the work through the research centre DrillWell―Drilling and Well Centre for Improved Recovery, a research cooperation between IRIS, NTNU, SINTEF and UiS. Moreover, experts from the industry have guided in the identification of the research path. For generating subgrids in our examples, we have used Triangle [21] . We also thank Schlumberger for providing access to their Petrel software for generating the pre-update models in our examples.

Appendix

A.1. Two Geosteering/Operational Challenges

Example 1: Effective updates of the geological uncertainty while drilling aim to offer support for improved geosteering workflows. Today, drilling speed (ROP, rate of penetration) is often reduced to allow time for geological interpretation and decision making while geosteering (see [22] ), with the consequence that the drilling performance decreases. But with longer drilling durations there is an increased probability for hole collapse before the production liner is set, in particular when drilling long horizontal sections. With high speed bidirectional telemetry like wired drill-pipe it is possible to quickly reprogram the rotary steerable system to new settings, and therefore it is possible to shorten the directional update loop considerably. This allows higher ROP so that the probability for drilling problems is reduced, but implies that less time is available for geological interpretation and steering decisions. A more effective decision loop will thus contribute to safe and effective drilling, and minimize the time spent in the open hole section.

Example 2: Geological interpretation is uncertain. Assume that measurements indicate that the drill bit may have penetrated a fault or the roof of the reservoir and drilled into a formation associated with high risk for drilling problems. The longer the drilling continues in this formation, the higher is the risk of encountering problems. However, it is important to stay in the pay zone to maximize future production. The decision to be taken is if 1) a side-track should be drilled to reduce the risk of serious drilling problems, 2) retain the strategy and continue along the currently planned trajectory, or 3) continue drilling to collect more information and continuously re-evaluate the decision to side-track. In the latter case, when more information is available it may already be too late to avoid drilling problems. Also the uncertainty in the well trajectory should be considered. To support decision making, the geosteering decision support tool indicated in Section A.2 could be applied to continuously monitor the operation and provide real-time model-based recommendations.

A.2. A Potential Future Geosteering Workflow

In [7] , a novel workflow is discussed where a set of earth model realizations are used to capture uncertainties in the locations and displacements of faults, in the shapes of stratigraphic interfaces and in water saturation. In the workflow the realizations are automatically conditioned to deep electromagnetic (EM) measurements while drilling. To minimize the time spent for managing multiple realizations representing complex geology in real-time, the ability to effectively update both structure and properties is crucial.

A decision analytic framework for geosteering is discussed in [23] [24] [25] . It aims to provide unbiased and consistent decision support under uncertainty by optimization with respect to multiple weighed geosteering and drilling objectives. Such objectives may include e.g. to minimize the probability of drilling into formations associated with drilling problems, maximize reservoir exposure and future production, minimize dogleg severity, minimize the cost of drilling, etc. Local earth model updates in real-time is a key element in the approach.

A possible highly automated future workflow for real-time geosteering and drilling decision support is to 1) if required, generate new realizations of the geology around and ahead of bit by locally updating the earth model, 2) locally condition all realizations by the recently received measurements as described in [7] , and 3) employ decision analytics methods to provide decision support under uncertainty. This process should be run in a continuous loop to assess the current risks and aid the optimization of the drilling operation. For support of the workflow excellent control with the geological structure and grid is paramount, which is the theme addressed in this paper. The faster the situation can be analyzed, the uncertainties and probabilities can be calculated, and a decision recommendation can be produced, the faster the modelling results can be applied by the drilling/geosteering team. This will contribute to safer and faster drilling, increased future production and reduced drilling costs.

A.3. Current Methods for Earth Modelling

The interpretation of the geological structure is frequently burdened by first order uncertainties, see e.g. [10] [13] [26] [27] [28] [29] . Poor assessment of structural uncertainties can thus have dramatic effects on the decisions to be taken, in particular when drilling in more complex geology. Uncertainty in the topology (connectivity) of the geological structure includes how geological interfaces, e.g. stratigraphic and fault surfaces, are connected. Topological uncertainty also includes if particular faults and layers exist or not, the lateral magnitude and depth of an erosional event (which layers that are eroded), complex fault patterns around a salt dome, if fault segments are linked or not, stratigraphic correlation between wells (e.g. if layers pinch out or not), if there is communication between layers across a fault or not, and so on. For example, when geosteering in seismically obscured areas, such as below salt or gas, interpretation uncertainties are often higher and the measurements and interpretations obtained while drilling become more important to guide the steering of the well.

Numerous numerical methods require a grid for their discretization algorithms. A grid is by nature a rigid numerical construction where the topological relationship between its cells cannot be modified in a flexible manner, and grid construction is an area of active research.

In today’s 3D earth modelling methodologies, implemented in software tools such as Petrel and IRAP RMS, a globally defined corner-point grid at a specific resolution is constructed early in the modelling process [12] [27] [30] . The grid construction is based on a deterministic representation of the geological structure, often denoted a base case or reference model. Rock properties are then distributed in this grid. The grid thus represents both structure and properties. The strategy implies that modifications in the topology of the geological structure cannot be effectively transferred to the earth model without invalidating the existing grid [12] [27] [30] . The updated structure is incompatible with the connectivity between the cells in the grid, and the cell connectivity cannot be modified in a general manner. Therefore, for each such update, a new grid based on the updated structure must be constructed. Moreover, all properties must be recomputed over the new grid. The reconstruction of the grid and distribution of properties may require much computational time, depending on the size of the model. In addition, much manual work is typically required to handle more complex updates of the geological structure.

As a result of the slow and complicated management, a crucial class of geological uncertainties may be underestimated or overlooked. Today, structural uncertainties are normally addressed by perturbing the geometry of a grid while the topology remains unaltered (see e.g. [31] ).

In [30] , a framework for modelling uncertainties in fault location and fault geometry in a structural model is presented. Here it is explained that if the base case structural model is updated, there are two possible procedures to construct a grid that match the new structure. In general, the grid must be entirely rebuilt and the new grid must be populated with properties as explained above. The second alternative is to deform the grid so that it matches the updated structure. This is an attractive option because the grid is not invalidated and the properties stored in the grid remain intact. In [32] [33] it is demonstrated how alterations in the displacement of a fault can be accommodated by grid deformation. But this only works for grids with simple fault geometries and when there are no changes to the topology of the geological structure [30] [32] .

Recently, a system for closed-loop reservoir modelling has been developed [33] [34] [35] [36] . In the history-matching workflow a set of model realizations that are used for capturing geological uncertainty, including geometric uncertainties in the geological structures, are updated. A main advantage of the workflow is that multiple model realizations can be automatically generated in batch and in parallel, in a fully reproducible manner. In [34] it is discussed how the geometry and depth of a stack of stratigraphic interfaces can be modified in a geologically realistic manner. In [33] it is explained that for each realization of the earth model grid, all individual modelling steps are still performed as in the conventional modelling process. Whenever the structural model is updated, each realization of the grid is constructed from scratch.

A.4. Recent Earth Modelling Methods

In the approach described in [13] [27] , the structure is split from the property representation and all properties are stored in a globally defined rectilinear grid. A single globally defined mapping, called the uvt-transformation, links the property grid (in a parametric uvt-space) with a geological grid (in the geological xyz-space) that conforms to the geological structure. The mapping enables population of the geological grid with values from the property grid. However, structural uncertainty modelling is performed by geometric perturbation of the base case structural model [13] [37] . Geometric perturbation does not include modifications in the topology of the fault network or layering.

In [14] [38] , a strategy for seismic interpretation based on capturing the geological evolution in the model is presented. The evolution is described as a sequence of geological processes that take place through geological time. For each step in geological time, a structural model can be constructed over a computational grid. The computational grid is a regular grid with the same resolution everywhere. In [14] it is explained that the properties are handled in a parameter space, separate from the structure. A globally defined bijective mapping links the existing properties with the restored structure, and the properties can be mapped into the computational grid where the structure resides. The strategy permits local updates of faulting by the application of a fault operator that is used to insert and remove faults by locally updating the mapping.

Earth modelling strategies where the structure and the properties are separated and connected via a mapping introduce a new level of flexibility for updating the earth model. When the structure is modified, the existing properties can be reused without the need for a full reconstruction of the properties. However, the numerical characteristics of the mapping determine e.g. its computational efficiency and its ability to handle local updates in the structure. In the two methods explained in [13] [14] , properties are represented in a globally defined grid and linked to the structure via a single mapping. Our approach is similar in that the structure and the properties are handled separately. However, we do not apply a global strategy.

In [12] [39] , a surface-based method for adaptive gridding during fluid flow simulation is presented. Here the modelled volume is split into separate rock volumes by surfaces that represent geological interfaces to capture e.g. stratigraphic and diagenetic heterogeneities. Aiming to avoid upscaling, each property within each volume is uniform. Each rock volume is separately gridded, and the grid resolution can be locally updated within each volume. A characteristic aspect of the approach is that it avoids the complexity of handling property representations captured in a global grid. Numerically, the strategy suggested in this paper enables the same functionalities. But it also offers an additional level of flexibility as it allows capturing non-uniform property distributions within each rock volume. This provides an alternative when capturing e.g. gradational changes in different directions, or more complex trends and distributions.

A.5. Gridding Strategy

Geological heterogeneities have complex geometries. It is emphasized e.g. in [12] that using strictly rectangular (Cartesian) grids, approximately rectangular (corner-point) grids or PEBI grids of a given spatial resolution often provide a poor representation of geological heterogeneity. Local updates and uncertainty handling of complex geological structures and other heterogeneities at well scale are main motivations behind the suggested strategy. A simplex-based subgrid (triangulation in 2D, tetrahedralization in 3D) typically offers better flexibility to adapt to complex structural geometries, and enables local control with grid resolution and quality. To optimize the grid to its application, it should be possible to use different constraints in the gridding algorithm to obtain subgrids of different qualities. Different quality parameters such as shape and orientation of grid cells, grid resolution and how well the grid can be adapted to complex geological structures are important for many numerical schemes.

The geometry and topology of the structural model has severe consequences for the construction of the grid. In [12] [27] it is discussed how a corner-point grid may fail to capture complex structural architectures. A complex structural model may result in a grid of too poor quality to support various simulations, or even inhibit the generation of a grid. Moreover, in the trade-off between numerical accuracy and the time spent for computations, grid resolution is a main factor. But the structural resolution dictates the resolution of the grid because the size of the grid cells cannot be coarser than the distance between individual elements in the geological structure. The requirements to the grid size may therefore dictate a coarse structural resolution that fails to capture important structural features. Furthermore, in [12] , it is discussed how the resolution of the grid cannot be modified in a flexible manner. In [27] , it is emphasized that a normal procedure to reduce the complexity and size of the grid is to simplify or leave out known structural elements such as faults (as exemplified in [2] ). For drilling in complex geology, where structural accuracy in the model around and ahead of bit is of particular importance, such limitations and practices are far from optimal. Optimally, the model capturing the geological interpretation should not be obstructed by limitations in the modelling method.

Gridding of complex geological structures is well known to be problematic. In the proposed strategy each region is individually discretized, potentially with different subgrids for each property. The subgrids in Figure 5 are generally not matching across their shared boundaries. For many applications, this is not a problem. For visualization, even small gaps between regions are acceptable. However, many numerical schemes require that the faces of neighbouring subgrids match across their shared boundaries. In [12] [39] it is discussed how grid resolution is independently controlled for each rock volume, and how grids for separate regions are linked together. When more constraints are used in the gridding process to obtain high quality grids, more computation time is generally required. For optimal performance in a real-time environment, it could be important to carefully tune the input parameters of the grid generator to the application of the grid. Moreover, similar to the method in [39] , our method allows populated subgrids to be generated in parallel. It may thus benefit from approaches for domain partitioning, see e.g. [40] , to further streamline the gridding process.

Flow modelling is not our primary concern. Yet, we believe that also such applications could potentially benefit from the proposed strategy when effective assessment of more complex structural uncertainties is required. In [39] , the generation of unstructured grids for use in the next generation of unstructured- mesh fluid flow simulators is discussed. Unstructured grids allow the capabilities of such simulators to be fully utilized in the modelling of complex reservoir architectures. In [41] , different workflows for construction of tetrahedral grids for simulations on complex structures are evaluated.

Submit or recommend next manuscript to SCIRP and we will provide best service for you:

Accepting pre-submission inquiries through Email, Facebook, LinkedIn, Twitter, etc.

A wide selection of journals (inclusive of 9 subjects, more than 200 journals)

Providing 24-hour high-quality service

User-friendly online submission system

Fair and swift peer-review system

Efficient typesetting and proofreading procedure

Display of the result of downloads and visits, as well as the number of cited articles

Maximum dissemination of your research work

Submit your manuscript at: http://papersubmission.scirp.org/

Or contact ijg@scirp.org

Cite this paper

Suter, E. , Cayeux, E. , Friis, H. , Kårstad, T. , Escalona, A. and Vefring, E. (2017) A Novel Method for Locally Updating an Earth Model While Geosteering.*International Journal of Geosciences*, **8**, 237-264. doi: 10.4236/ijg.2017.82010.

Suter, E. , Cayeux, E. , Friis, H. , Kårstad, T. , Escalona, A. and Vefring, E. (2017) A Novel Method for Locally Updating an Earth Model While Geosteering.

References

[1] Bittar, M. and Aki, A. (2015) Advancement and Economic Benefit of Geosteering and Well-Placement Technology. The Leading Edge, 34, 524-528.

https://dx.doi.org/10.1190/tle34050524.1

[2] Arata, F., Gangemi, G., Mele, M., Tagliamonte, R.L., Tarchiani, C., Chinellato, F.J., Denichou, M. and Maggs, D. (2016) High-Resolution Reservoir Mapping: From Ultradeep Geosteering Tools to Real-Time Updating of Reservoir Models. Paper SPE-183133-MS Presented at the Abu Dhabi International Petroleum Exhibition & Conference, 7-10 November, Abu Dhabi, UAE.

https://dx.doi.org/10.2118/183133-MS

[3] Bashir, F., Al-Hawi, M., Bhuana, I., Abbas, M. and Ghamdi, Y. (2016) Real Time 3D Modeling to Optimize Geosteering in Clastic Reservoir—Case Study. Paper OTC-26635-MS Presented at the Offshore Technology Conference Asia, Kuala Lumpur, Malaysia, 22-25 March 2016.

https://dx.doi.org/10.4043/26635-MS

[4] Constable, M.V., Antonsen, F., Stalheim, S.O., Olsen, P.A., Fjell, O.Z., Dray, N., Eikenes, S., Aarflot, H., Haldorsen, K., Digranes, G., Seydoux, J., Omeragic, D., Thiel, M., Davydychev, A., Denichou, J.-M., Salim, D., Frey, M., Homan, D. and Tan, S. (2016) Looking Ahead of the Bit While Drilling: From Vision to Reality. SPWLA 57th Annual Logging Symposium, Reykjavik, 25-29 June 2016.

[5] Suter, E., Cayeux, E., Vefring, E., aesheim, L., Friis, H., Escalona, A. and Karstad, T. (2010) An Efficient Approach for Earth Model Updates. Paper SPE-136319-MS Presented at the SPE Russian Oil and Gas Conference and Exhibition, Moscow, Russia, 26-28 October 2010.

https://dx.doi.org/10.2118/136319-MS

[6] Suter, E., Cayeux, E., Escalona, A., Karstad, T. and Vefring, E. (2012) A Strategy for Effective Local Updates of the Geological Structure in an Earth Model during Drilling. Extended Abstract Presented at the 74th EAGE Conference and Exhibition Incorporating EUROPEC, Copenhagen, Denmark, 4 June 2012.

https://dx.doi.org/10.3997/2214-4609.20148222

[7] Luo, X., Eliasson, P., Alyaev, S., Romdhane, A., Suter, E., Querendez, E. and Vefring, E. (2015) An Ensemble Based Framework for Proactive Geosteering. SPWLA 56th Annual Logging Symposium, Long Beach, California, USA, 18-22 July 2015.

[8] Caumon, G., Lepage, F., Sword, C.H. and Mallet, J.-L. (2004) Building and Editing a Sealed Geological Model. Mathematical Geology, 36, 405-424.

https://dx.doi.org/10.1023/B:MATG.0000029297.18098.8a

[9] Caumon, G., Collon-Drouaillet, P., Le Carlier de Veslud, C., Viseur, S. and Sausse, J. (2009) Surface-Based 3d Modeling of Geological Structures. Mathematical Geosciences, 41, 927-945.

https://doi.org/10.1007/s11004-009-9244-2

[10] Caumon, G. (2010) Towards Stochastic Time-Varying Geological Modeling. Mathematical Geosciences, 42, 555-569.

https://doi.org/10.1007/s11004-010-9280-y

[11] Pellerin, J., Caumon, G., Julio, C., Mejia-Herrera, P. and Botella, A. (2015) Elements for Measuring the Complexity of 3d Structural Models: Connectivity and Geometry. Computers & Geosciences, 76, 130-140.

https://doi.org/10.1016/j.cageo.2015.01.002

[12] Jackson, M.D., Hampson, G.J., Saunders, J.H., Elsheikh, A., Graham, G.H. and Massart, B.Y.G. (2013) Surface Based Reservoir Modelling for Flow Simulation. Special Publications, 387, 271-292.

https://doi.org/10.1144/SP387.2

[13] Mallet, J.-L. (2014) Elements of Mathematical Sedimentary Geology: The GeoChron Model. EAGE Publications.

[14] Hjelle, O. (2014) A Hamilton-Jacobi Framework for Modeling Geological Folding and Deformation. PhD Dissertation, Department of Imaging Physics, Faculty of Applied Sciences, Delft University of Technology.

https://dx.doi.org/10.4233/uuid:e5932fea-ce78-4aa7-8224-7eaad640b843

[15] Hormann, K. and Floater, M.S. (2006) Mean Value Coordinates for Arbitrary Planar Polygons. ACM Transactions on Graphics, 25, 1424-1441.

https://doi.org/10.1145/1183287.1183295

[16] Suter, E., Friis, H.A., Escalona, A., Karstad, T. and Vefring, E.H. (2017) A Novel Method for Multi-Resolution Earth Model Gridding. Paper SPE-182687-MS Presented at the SPE Reservoir Simulation Conference, Montgomery, Texas, USA, 20-22 February 2017.

[17] Floater, M.S. (2003) Mean Value Coordinates. Comput. Aided Geom. Des., 20, 19-27.

https://doi.org/10.1016/S0167-8396(03)00002-5

[18] Ju, T., Schaefer, S. and Warren, J. (2005) Mean Value Coordinates for Closed Triangular Meshes. ACM Transactions on Graphics, 24, 561-566.

https://doi.org/10.1145/1073204.1073229

[19] Gilardet, M., Guillon, S., Jobard, B. and Komatitsch, D. (2013) Seismic Image Restoration Using Nonlinear Least Squares Shape Optimization. Procedia Computer Science, 18, 732-741.

https://doi.org/10.1016/j.procs.2013.05.237

[20] Schneider, T. and Hormann, K. (2015) Smooth Bijective Maps between Arbitrary Planar Polygons. Computer Aided Geometric Design, 35-36, 243-254.

https://doi.org/10.1016/j.cagd.2015.03.010

[21] Shewchuk, J.R. (1996) Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator. In: Lin, M.C. and Manocha, D., Eds., Applied Computational Geometry: Towards Geometric Engineering, Vol. 1148 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, 203-222.

[22] Antonsen, F., Barbosa, J.E.P., Morani, B., Klein, K., Kjolleberg, M., McCann, A., Olsen, P.A., Constable, M.V., Eidem, M., Gjengedal, J.A., Antonov, Y., Hartmann, A., Larsen, D., Skillings, J. and Tisley-Baker, R. (2015) Well Placement and Reservoir Characterization Challenges on Peregrino Mitigated by Development of Extra-Deep Azimuthal EM Technology. SPWLA 56th Annual Logging Symposium, Long Beach, California, USA, 18-22 July 2015.

[23] Kullawan, K., Bratvold, R. and Bickel, J.E. (2014) A Decision Analytic Approach to Geosteering Operations. SPE Drilling and Completion, 29, 36-46.

https://doi.org/10.2118/167433-PA

[24] Kullawan, K., Bratvold, R. and Bickel, J.E. (2016) Value Creation with Multi-Criteria Decision Making in Geosteering Operations. International Journal of Petroleum Technology, 3, 15-31.

https://dx.doi.org/10.15377/2409-787X.2016.03.01.2

[25] Kullawan, K., Bratvold, R. and Nieto, C. (2016) Decision-Oriented Geosteering and the Value of Look-Ahead Information: A Case-Based Study. SPE Journal.

http://dx.doi.org/10.2118/184392-PA

[26] Bond, C., Gibbs, A., Shipton, Z. and Jones, S. (2007) What Do You Think This Is? “Conceptual uncertainty” in Geoscience Interpretation. GSA Today, 17, 4-10.

https://doi.org/10.1130/GSAT01711A.1

[27] Mallet, J.-L. (2008) Numerical Earth Models. EAGE Publications.

[28] Caers, J. (2011) Modeling Uncertainty in the Earth Sciences. Wiley.

https://doi.org/10.1002/9781119995920

[29] Bond, C.E. (2015) Uncertainty in Structural Interpretation: Lessons to Be Learnt. Journal of Structural Geology, 74, 185-200.

https://doi.org/10.1016/j.jsg.2015.03.003

[30] Roe, P., Georgsen, F. and Abrahamsen, P. (2014) An Uncertainty Model for Fault Shape and Location. Mathematical Geosciences, 46, 957-969.

https://doi.org/10.1007/s11004-014-9536-z

[31] Caumon, G. (2014) Thoughts on Geological Uncertainty Assessment in Integrated Reservoir Modeling. Extended Abstract Presented at the Second EAGE Integrated Reservoir Modelling Conference, Dubai, UAE, November 2014.

https://doi.org/10.3997/2214-4609.20147463

[32] Seiler, A., Aanonsen, S.I., Evensen, G. and Lia, O. (2010) An Elastic Grid Approach for Fault Uncertainty Modelling and Updating Using the Ensemble Kalman Filter. Paper SPE-130422-MS Presented at the SPE EUROPEC/EAGE Annual Conference and Exhibition, Barcelona, Spain, 14-17 June 2010.

https://dx.doi.org/10.2118/130422-MS

[33] Pettan, C. and Stromsvik, J.F. (2013) The Peregrino Challenge: How to Keep Reliable Models While Drilling Eight Wells per Year. Offshore Technology Conference, Rio de Janeiro, Brazil, 29-31 October 2013, OTC-24522.

[34] Stenerud, V.R., Kallekleiv, H., Abrahamsen, P., Dahle, P., Skorstad, A. and Viken, M.H.A. (2012) Added Value by Fast and Robust Conditioning of Structural Surfaces to Horizontal Wells for Real-World Reservoir Models. Paper SPE-159746-MS Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, 8-10 October 2012.

https://dx.doi.org/10.2118/159746-MS

[35] Hanea, R., Evensen, G., Hustoft, L., Ek, T., Chitu, A. and Wilschut, F. (2015) Reservoir Management under Geological Uncertainty Using Fast Model Update. Paper SPE-173305-MS Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, USA, 23-25 February 2015.

https://doi.org/10.2118/173305-ms

[36] Aarnes, I., Midtveit, K. and Skorstad, A. (2015) Evergreen Workflows That Capture Uncertainty—The Benefits of an Unlocked Structure. First Break, 33, 89-92.

[37] Caers, J. (2015) Is the Geochron Model the Way Forward for Mathematical Sedimentary Geology? First Break, 33, 95-96.

[38] Hjelle, O., Petersen, S.A. and Bruaset, A.M. (2013) A Numerical Framework for Modeling Folds in Structural Geology. Mathematical Geosciences, 45, 255-276.

https://doi.org/10.1007/s11004-013-9452-7

[39] Jackson, M., Percival, J., Mostaghimi, P., Tollit, B., Pavlidis, D., Pain, C., Gomes, J., Elsheikh, A., Salinas, P., Muggeridge, A. and Blunt, M. (2015) Reservoir Modeling for Flow Simulation by Use of Surfaces, Adaptive Unstructured Meshes, and an Overlapping-Control-Volume Finite-Element Method. SPE Reservoir Evaluation and Engineering, 18, 115-132.

https://doi.org/10.2118/163633-PA

[40] Loseille, A., Alauzet, F. and Menier, V. (2017) Unique Cavity-Based Operator and Hierarchical Domain Partitioning for Fast Parallel Generation of Anisotropic Meshes. Computer-Aided Design, 85, 53-67,

https://dx.doi.org/10.1016/j.cad.2016.09.008

[41] Zehner, B., Borner, J.H., Gorz, I. and Spitzer, K. (2015) Workflows for Generating Tetrahedral Meshes for Finite Element Simulations on Complex Geological Structures. Computers & Geosciences, 79, 105-117.

https://doi.org/10.1016/j.cageo.2015.02.009

[1] Bittar, M. and Aki, A. (2015) Advancement and Economic Benefit of Geosteering and Well-Placement Technology. The Leading Edge, 34, 524-528.

https://dx.doi.org/10.1190/tle34050524.1

[2] Arata, F., Gangemi, G., Mele, M., Tagliamonte, R.L., Tarchiani, C., Chinellato, F.J., Denichou, M. and Maggs, D. (2016) High-Resolution Reservoir Mapping: From Ultradeep Geosteering Tools to Real-Time Updating of Reservoir Models. Paper SPE-183133-MS Presented at the Abu Dhabi International Petroleum Exhibition & Conference, 7-10 November, Abu Dhabi, UAE.

https://dx.doi.org/10.2118/183133-MS

[3] Bashir, F., Al-Hawi, M., Bhuana, I., Abbas, M. and Ghamdi, Y. (2016) Real Time 3D Modeling to Optimize Geosteering in Clastic Reservoir—Case Study. Paper OTC-26635-MS Presented at the Offshore Technology Conference Asia, Kuala Lumpur, Malaysia, 22-25 March 2016.

https://dx.doi.org/10.4043/26635-MS

[4] Constable, M.V., Antonsen, F., Stalheim, S.O., Olsen, P.A., Fjell, O.Z., Dray, N., Eikenes, S., Aarflot, H., Haldorsen, K., Digranes, G., Seydoux, J., Omeragic, D., Thiel, M., Davydychev, A., Denichou, J.-M., Salim, D., Frey, M., Homan, D. and Tan, S. (2016) Looking Ahead of the Bit While Drilling: From Vision to Reality. SPWLA 57th Annual Logging Symposium, Reykjavik, 25-29 June 2016.

[5] Suter, E., Cayeux, E., Vefring, E., aesheim, L., Friis, H., Escalona, A. and Karstad, T. (2010) An Efficient Approach for Earth Model Updates. Paper SPE-136319-MS Presented at the SPE Russian Oil and Gas Conference and Exhibition, Moscow, Russia, 26-28 October 2010.

https://dx.doi.org/10.2118/136319-MS

[6] Suter, E., Cayeux, E., Escalona, A., Karstad, T. and Vefring, E. (2012) A Strategy for Effective Local Updates of the Geological Structure in an Earth Model during Drilling. Extended Abstract Presented at the 74th EAGE Conference and Exhibition Incorporating EUROPEC, Copenhagen, Denmark, 4 June 2012.

https://dx.doi.org/10.3997/2214-4609.20148222

[7] Luo, X., Eliasson, P., Alyaev, S., Romdhane, A., Suter, E., Querendez, E. and Vefring, E. (2015) An Ensemble Based Framework for Proactive Geosteering. SPWLA 56th Annual Logging Symposium, Long Beach, California, USA, 18-22 July 2015.

[8] Caumon, G., Lepage, F., Sword, C.H. and Mallet, J.-L. (2004) Building and Editing a Sealed Geological Model. Mathematical Geology, 36, 405-424.

https://dx.doi.org/10.1023/B:MATG.0000029297.18098.8a

[9] Caumon, G., Collon-Drouaillet, P., Le Carlier de Veslud, C., Viseur, S. and Sausse, J. (2009) Surface-Based 3d Modeling of Geological Structures. Mathematical Geosciences, 41, 927-945.

https://doi.org/10.1007/s11004-009-9244-2

[10] Caumon, G. (2010) Towards Stochastic Time-Varying Geological Modeling. Mathematical Geosciences, 42, 555-569.

https://doi.org/10.1007/s11004-010-9280-y

[11] Pellerin, J., Caumon, G., Julio, C., Mejia-Herrera, P. and Botella, A. (2015) Elements for Measuring the Complexity of 3d Structural Models: Connectivity and Geometry. Computers & Geosciences, 76, 130-140.

https://doi.org/10.1016/j.cageo.2015.01.002

[12] Jackson, M.D., Hampson, G.J., Saunders, J.H., Elsheikh, A., Graham, G.H. and Massart, B.Y.G. (2013) Surface Based Reservoir Modelling for Flow Simulation. Special Publications, 387, 271-292.

https://doi.org/10.1144/SP387.2

[13] Mallet, J.-L. (2014) Elements of Mathematical Sedimentary Geology: The GeoChron Model. EAGE Publications.

[14] Hjelle, O. (2014) A Hamilton-Jacobi Framework for Modeling Geological Folding and Deformation. PhD Dissertation, Department of Imaging Physics, Faculty of Applied Sciences, Delft University of Technology.

https://dx.doi.org/10.4233/uuid:e5932fea-ce78-4aa7-8224-7eaad640b843

[15] Hormann, K. and Floater, M.S. (2006) Mean Value Coordinates for Arbitrary Planar Polygons. ACM Transactions on Graphics, 25, 1424-1441.

https://doi.org/10.1145/1183287.1183295

[16] Suter, E., Friis, H.A., Escalona, A., Karstad, T. and Vefring, E.H. (2017) A Novel Method for Multi-Resolution Earth Model Gridding. Paper SPE-182687-MS Presented at the SPE Reservoir Simulation Conference, Montgomery, Texas, USA, 20-22 February 2017.

[17] Floater, M.S. (2003) Mean Value Coordinates. Comput. Aided Geom. Des., 20, 19-27.

https://doi.org/10.1016/S0167-8396(03)00002-5

[18] Ju, T., Schaefer, S. and Warren, J. (2005) Mean Value Coordinates for Closed Triangular Meshes. ACM Transactions on Graphics, 24, 561-566.

https://doi.org/10.1145/1073204.1073229

[19] Gilardet, M., Guillon, S., Jobard, B. and Komatitsch, D. (2013) Seismic Image Restoration Using Nonlinear Least Squares Shape Optimization. Procedia Computer Science, 18, 732-741.

https://doi.org/10.1016/j.procs.2013.05.237

[20] Schneider, T. and Hormann, K. (2015) Smooth Bijective Maps between Arbitrary Planar Polygons. Computer Aided Geometric Design, 35-36, 243-254.

https://doi.org/10.1016/j.cagd.2015.03.010

[21] Shewchuk, J.R. (1996) Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator. In: Lin, M.C. and Manocha, D., Eds., Applied Computational Geometry: Towards Geometric Engineering, Vol. 1148 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, 203-222.

[22] Antonsen, F., Barbosa, J.E.P., Morani, B., Klein, K., Kjolleberg, M., McCann, A., Olsen, P.A., Constable, M.V., Eidem, M., Gjengedal, J.A., Antonov, Y., Hartmann, A., Larsen, D., Skillings, J. and Tisley-Baker, R. (2015) Well Placement and Reservoir Characterization Challenges on Peregrino Mitigated by Development of Extra-Deep Azimuthal EM Technology. SPWLA 56th Annual Logging Symposium, Long Beach, California, USA, 18-22 July 2015.

[23] Kullawan, K., Bratvold, R. and Bickel, J.E. (2014) A Decision Analytic Approach to Geosteering Operations. SPE Drilling and Completion, 29, 36-46.

https://doi.org/10.2118/167433-PA

[24] Kullawan, K., Bratvold, R. and Bickel, J.E. (2016) Value Creation with Multi-Criteria Decision Making in Geosteering Operations. International Journal of Petroleum Technology, 3, 15-31.

https://dx.doi.org/10.15377/2409-787X.2016.03.01.2

[25] Kullawan, K., Bratvold, R. and Nieto, C. (2016) Decision-Oriented Geosteering and the Value of Look-Ahead Information: A Case-Based Study. SPE Journal.

http://dx.doi.org/10.2118/184392-PA

[26] Bond, C., Gibbs, A., Shipton, Z. and Jones, S. (2007) What Do You Think This Is? “Conceptual uncertainty” in Geoscience Interpretation. GSA Today, 17, 4-10.

https://doi.org/10.1130/GSAT01711A.1

[27] Mallet, J.-L. (2008) Numerical Earth Models. EAGE Publications.

[28] Caers, J. (2011) Modeling Uncertainty in the Earth Sciences. Wiley.

https://doi.org/10.1002/9781119995920

[29] Bond, C.E. (2015) Uncertainty in Structural Interpretation: Lessons to Be Learnt. Journal of Structural Geology, 74, 185-200.

https://doi.org/10.1016/j.jsg.2015.03.003

[30] Roe, P., Georgsen, F. and Abrahamsen, P. (2014) An Uncertainty Model for Fault Shape and Location. Mathematical Geosciences, 46, 957-969.

https://doi.org/10.1007/s11004-014-9536-z

[31] Caumon, G. (2014) Thoughts on Geological Uncertainty Assessment in Integrated Reservoir Modeling. Extended Abstract Presented at the Second EAGE Integrated Reservoir Modelling Conference, Dubai, UAE, November 2014.

https://doi.org/10.3997/2214-4609.20147463

[32] Seiler, A., Aanonsen, S.I., Evensen, G. and Lia, O. (2010) An Elastic Grid Approach for Fault Uncertainty Modelling and Updating Using the Ensemble Kalman Filter. Paper SPE-130422-MS Presented at the SPE EUROPEC/EAGE Annual Conference and Exhibition, Barcelona, Spain, 14-17 June 2010.

https://dx.doi.org/10.2118/130422-MS

[33] Pettan, C. and Stromsvik, J.F. (2013) The Peregrino Challenge: How to Keep Reliable Models While Drilling Eight Wells per Year. Offshore Technology Conference, Rio de Janeiro, Brazil, 29-31 October 2013, OTC-24522.

[34] Stenerud, V.R., Kallekleiv, H., Abrahamsen, P., Dahle, P., Skorstad, A. and Viken, M.H.A. (2012) Added Value by Fast and Robust Conditioning of Structural Surfaces to Horizontal Wells for Real-World Reservoir Models. Paper SPE-159746-MS Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, 8-10 October 2012.

https://dx.doi.org/10.2118/159746-MS

[35] Hanea, R., Evensen, G., Hustoft, L., Ek, T., Chitu, A. and Wilschut, F. (2015) Reservoir Management under Geological Uncertainty Using Fast Model Update. Paper SPE-173305-MS Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, USA, 23-25 February 2015.

https://doi.org/10.2118/173305-ms

[36] Aarnes, I., Midtveit, K. and Skorstad, A. (2015) Evergreen Workflows That Capture Uncertainty—The Benefits of an Unlocked Structure. First Break, 33, 89-92.

[37] Caers, J. (2015) Is the Geochron Model the Way Forward for Mathematical Sedimentary Geology? First Break, 33, 95-96.

[38] Hjelle, O., Petersen, S.A. and Bruaset, A.M. (2013) A Numerical Framework for Modeling Folds in Structural Geology. Mathematical Geosciences, 45, 255-276.

https://doi.org/10.1007/s11004-013-9452-7

[39] Jackson, M., Percival, J., Mostaghimi, P., Tollit, B., Pavlidis, D., Pain, C., Gomes, J., Elsheikh, A., Salinas, P., Muggeridge, A. and Blunt, M. (2015) Reservoir Modeling for Flow Simulation by Use of Surfaces, Adaptive Unstructured Meshes, and an Overlapping-Control-Volume Finite-Element Method. SPE Reservoir Evaluation and Engineering, 18, 115-132.

https://doi.org/10.2118/163633-PA

[40] Loseille, A., Alauzet, F. and Menier, V. (2017) Unique Cavity-Based Operator and Hierarchical Domain Partitioning for Fast Parallel Generation of Anisotropic Meshes. Computer-Aided Design, 85, 53-67,

https://dx.doi.org/10.1016/j.cad.2016.09.008

[41] Zehner, B., Borner, J.H., Gorz, I. and Spitzer, K. (2015) Workflows for Generating Tetrahedral Meshes for Finite Element Simulations on Complex Geological Structures. Computers & Geosciences, 79, 105-117.

https://doi.org/10.1016/j.cageo.2015.02.009