A Novel Method for Locally Updating an Earth Model While Geosteering

Earth models are important tools for support of decision making processes for optimal exploitation of subsurface resources. For geosteering and other realtime 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.


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.

Current Methods for Geosteering
The recently developed ultra-deep directional electro-magnetic (deep EM) technology 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-resolution 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.

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.

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.

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 1 1 R and 1 2 R 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 { } , , , f f f f .Each mapping links a region with the corresponding part of the parameter domain of the property function.For example, the subgrid in 1  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.

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

The Structural Model
The structural model S  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 ∈  R 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 S  .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.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.

Property Functions
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.

Construction of a Populated Grid
The construction of a populated grid takes place by first identifying the closed regions  in S  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.

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 S  that together constitute the boun- dary 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.

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  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 R in L can be associated with each corresponding i D .As an example, consider Figure 2 where L 1 is split into two regions R 1 1 and R 1  2 by the fault.The parameter domain of the property function Φ 1 is then split into two subdomains D 1  1 and D 1 2 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 : f Ω Ω from a source polygon Ω S to a tar- get 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 f , , 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.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 S  where subgrids should be constructed, together with as- sociated property functions.The output is a populated subgrid for each region.
The algorithm handles each region independently of the other regions.

Handling of Properties When the Structure Is Locally Modified
Existing properties are the results of modelling prior to drilling, or of predictions 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  S in a hierarchical manner.Out- side 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.

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 Structural modelling can be challenging, and the particular technique to locally updating S   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.that are separately gridded and populated with properties.In the updated S   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 re- gions 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.

Examples of Local Updates of an Earth Model Grid
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

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 fre-quently 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 R D P P f , .When using the MVCbased mapping, barycentric coordinates for a given source point ∈ R x 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 ′ ∈ D x 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.

Conclusions
Decision making to optimize the exploitation of subsurface resources is challenging in particular when targeting more complex fields and reservoirs.Threedimensional 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 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]  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 ex-isting 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 res-olution 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 unstructuredmesh fluid flow simulators is discussed.Unstructured grids allow the capabilities

1 R is linked with the part 1 1 D 1 D 1 D 1 RFigure 1 .
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 . 1 D and 1 2 D
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 ( f ).The result is a populated grid as shown at the bottom.The colors in the figure indicate property values.

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 ter domain D Φ , for example the unit square.
. 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 R , i = 1, •••, n.

Figure 3 .
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.
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 ∈ Rx to be mapped to its corresponding point ′ 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.


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.be a structural model Let  be the regions from S  to be gridded and populated Let Φ be the property functions Φ for each L containing an R in  Let  be a set of grids G Z ← ∅  Let  be a set of mappings f ← ∅  Let  be a set of arrays a for storing grid-based property values parallel for 12: return a 13: end function 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.


S be the pre-update structural model, while the corresponding post-update structural model is denoted S is obtained by locally updating the structural connectivity and/or geometry in S   .  is the set of regions in S that are affected by the local structural update, whereas   is the set of regions in S   that are established or deformed in the update.Thus, a local update of S   implies that one or more regions   in S   are affected; they are either deformed or new regions have taken their place.In both cases,   is the set of regions in S   that must be attended.The rest of the regions in S   already exist in S  .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.

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

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

A. 1 .
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 Appendix 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 [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 in-to 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.