GB2587909A - Method of three dimensional (3D) modelling a reservoir system using dip-domain boundaries - Google Patents

Abstract

Defining an extended surface boundary description (SBD) comprises obtaining an initial SBD and extension distance. The SBD is transformed to a transformed SBD oriented such that its long-axis is parallel to a first direction and an average normal vector is parallel to a second direction, the directions being mutually orthogonal and defining a first plane, and wherein the transformed SBD comprises source vertices (V1, V2…) separated by boundary segments (S1, S2…). A bisector (i.e. initial extension vector, F1, F2…) is generated from each source vertex, projecting orthogonally to the first plane and bisecting an angle subtended by the adjacent boundary segments. A new bisector (F6 & F7) is generated (repeatedly if required) at each intersection point (IP12, IP14) of a pair of bisectors, when intersections occur within the extension distance. An extended boundary vertex (V1’ in Fig. 8) is generated on each bisector at the extension distance (d in Fig. 8), the extended SBD being generated by connecting the extended boundary vertices. An extended boundary region defined between the extended SBD and the transformed SBD is tessellated, and the region defined by the extended and initial SBDs is transformed such that the initial SBD is in its original location.

Description

METHOD OF THREE DIMENSIONAL (3D) MODELLING A RESERVOIR SYSTEM USING DIP-DOMAIN BOUNDARIES The present invention relates to systems and methods for 3D modelling a reservoir system using dip-domain boundaries, and in particular for conditioning a 3D reservoir model to avoid common mesh artefacts.

BACKGROUND

Extraction of economically valuable resources requires an understanding of the shapes and locations of structures in the Earth's subsurface. These structures primarily consist of the contacts between geological formations, herein referred to as horizons, and discontinuities where contacts are offset by breaking, herein referred to as faults. Horizons and faults can be represented in two dimensions (2D) as lines and in three dimensions (3D) as curvilinear surfaces. Note that, in geology, a horizon may refer to either a bedding surface where there is marked change in the lithology within a sequence of sedimentary or volcanic rocks, or a distinctive layer or thin bed with a characteristic lithology or fossil content within a sequence; whereas a fault may refer to a planar fracture or discontinuity in a volume of rock, across which there has been significant displacement as a result of rock-mass movement, In the upstream oil and gas industry, reservoir modelling is widely adopted to, for example, estimate oil/gas reserves, forecast future production, make decisions on the development of the field, and choose suitable reservoir management. Typically, reservoir modelling involves a process of constructing a computer model of a petroleum reservoir. There are two main categories of reservoir models: geological models and reservoir simulation models. Geological models are computer models that use geophysical and geological information gathered on and below the Earth surface to create representations of certain portions of the Earth's crust (e.g., reservoirs). 'Whereas, reservoir simulation models are computer models that are created to simulate the flow of fluids (e.g., oil, water and gas) within a reservoir. Reservoir simulation models are often implemented with finite difference methods, which use both structured and unstructured grids to accurately represent the geometry of a reservoir. The results of geological models are often used as input for reservoir simulation models.

A typical geological modelling process may comprise, for example, following three main steps: interpretation, conditioning or quality control, and validation. A new horizon (or surface) may be first interpreted based on recognition of geological spatial relationships in raw data, e.g., such as for example seismic data and well data. Then, the interpreted horizon data (e.g. regular (rectangular) grids or triangulated irregular networks (TINs)) may be conditioned using editing techniques to tie an interpretation as closely as is reasonable to raw data and remove artefacts. Finally, the conditioned data may be validated via application of more advanced scientifically based methods to ensure that an interpretation conforms to known geological and physical principles.

When performing geological modelling as has been described above, for example to understand the location or form of a hydrocarbon reservoir in the upstream oil and gas industry, it is often desirable to include a 3D model building phase. During the 3D model building phase, the interpreted horizon data is conditioned or edited in order to generate a more valid and complete geological interpretation which subsequently leads to a more reliable and accurate 3D geological model. Conditioning of a geological model typically requires extension of existing 3D interpreted horizons, e.g., extension of TINs.

However, the quality of horizon extension is strongly dependent on whether the horizon is in the vicinity of a fault. This is because faults can cause breaks in seismic images of the horizons of interest, resulting in locally reduced data quality and thus increased interpretational uncertainty. Consequently, horizon interpretations in the vicinity of faults are likely to be incomplete or unreliable either because part of the horizon data may be missing or because the data may contain a large amount of noise which leads to data artefacts. If such an unconditioned geological model, e.g., having data artefacts or data missing, continues through downstream processes in the workflow, it tends to become increasingly difficult to work with. For example, the unconditioned model can cause digital analysis workflows to fail or introduce unnecessary additional uncertainty and risk, resulting in physically unreasonable or meaningless results. This could introduce unrecognised errors or cause analysis problems for, e.g. other petroleum specialists working in downstream processes, where the geological model may be used as the basis of subsequent processes such as reservoir modelling.

Effective execution of model analysis, or acceptable quality for visualization of the 3D model, often requires that horizon surfaces fit precisely to faults, meaning that there should be no locations with space or "gaps" between the horizon and fault. The aforementioned conditioning step is able to refill the missing data or remove the data artefacts in the fault zone and thereby allow the horizon of interest (e.g., a source TIN) to be extended towards the fault plane (e.g., a target TIN) such that the extended or new horizon is attached to the fault surface. Alternatively or additionally, horizons may be desired to be extended into free space such that surface meshes, for example, can be grown into areas where there is either no horizon data or horizon data is missing. Extension of horizons into free space is particularly useful when two or more separate areas of horizon data are to be joined or models with horizons of equal extents occupying a volume are to be created.

In a conventional horizon or surface extension method, a new horizon boundary line, extended from an original horizon boundary line, is first generated. Subsequently, the area between the new and original horizon boundary lines is tessellated. The new horizon boundary line comprises a plurality of new boundary points which are extended along the outward direction from their corresponding original boundary points by a predetermined extension distance. The extension vector, having information of the extension direction and distance, is determined to be the cross product of the normal of a triangle mesh and the outward edge boundary. However, with this method, extension vectors are not "aware" of other vectors and can cross and overlap when the edge geometries of the source horizon are irregular or when the extension distance is increased. The crossing of extension vectors results in die generation of undesired features during the horizon extension such as overlapping mesh areas with overturned elements and/or empty mesh areas. Therefore, the conventional or existing method is limited to small extension distances when used to extend horizons with irregular edge geometries.

The aforementioned problem can be mitigated with further processing of horizon data using manual editing techniques, such as smoothing out, resampling or cutting out die irregular edge geometries of the source surface before applying surface extension. However, such data processing is a manual process and thus requires a significant amount of time to complete. Moreover, the fact that each separate surface requires individual processing further compounds the issue.

Therefore, it is desirable to have a method that is able to overcome problems of existing horizon extension methods while not introducing new issues that would compromise the performance of geological modelling.

SUMMARY OF INVENTION

In a first aspect of the invention, there is provided a computer implemented method for defining a 3D extended surface boundary description of a subsurface volume comprising: a) obtaining an initial surface boundary description and an extension distance describing the distance that the initial surface boundary description is to be extended; b) transforming the initial surface boundary description to obtain a transformed initial surface boundary description which is oriented such that a long-axis of the transformed initial surface boundary description is substantially parallel to a first direction and an average normal vector of die transformed initial surface boundary description is substantially parallel to a second direction, the first direction and the second direction being mutually orthogonal and defining a first plane, and wherein the transformed initial surface boundary description comprises a plurality of source vertices with a boundary segment between each adjacent pair of source vertices; c) generating a bisector from each of said source vertices, each of these bisectors being source bisectors that are projectablc to a third direction orthogonal to the first plane and bisecting an angle subtended by said boundary segments immediately adjacent its respective source vertex; d) extending each bisector until said bisectors are extended the extension distance, while generating a new bisector at each intersection point where pairs of adjacent said bisectors intersect within said extension distance, in place of the intersecting pair of bisectors; e) generating an extended boundary vertex on each bisector at the extension distance, said bisectors at the extension distance being final bisectors; 0 generating said extended surface boundary description by connecting said extended boundary vertices; and g) tessellating an extended boundary region defined between said extended surface boundary description and said transformed initial surface boundary description; h) transforming a region defined by said extended surface boundary description and said initial surface boundary description such that said initial surface boundary description is in its original location.

BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments of the invention will now be described, by way of example only, by reference to the accompanying drawings, in which: Figure 1 schematically depicts a conventional method for determining a new boundary vertex extended from an existing boundary vertex on a to-be-extended surface boundary; Figure 2 shows (a) a plan view and (b) a 3D view of an example surface extension of an irregular horizon boundary using the method of Figure 1; Figure 3 shows an extreme example of the artefact problem where the method of Figure 1 is used to extend a concave or embayed surface boundary; Figure 4 shows a flowchart of the horizon or surface extension method in accordance with an embodiment; Figure 5 (a) schematically illustrates transformation of a selected boundary region; and (b) the corresponding transformed boundary region in accordance with an embodiment; Figure 6 shows a flowchart of an implementation of step 403 of Figure 4 in accordance with an embodiment; Figure 7 schematically illustrates the process of determination of (a) initial bisectors and (6, c) new bisectors based on an existing surface boundary which is part of step 403 of Figure 4, in accordance with an embodiment; Figure 8 schematically illustrates the process of determination of new boundary vertices based on the bisectors determined in the process of Figure 7, in accordance with an embodiment; Figure 9 schematically illustrates (a) an example extended or new boundary line after following the step 403 of Figure 4 and (b) an example tessellation of the extended surface area defined in Figure 9(a), in accordance with an embodiment; and Figure 10 shows a flowchart of an implementation of step 405 of Figure 4, in accordance with an embodiment.

DETAILED DESCRIPTION OF THE EMBODIMENTS

Figure 1 schematically depicts the aforementioned conventional method for determining a new boundary vertex extended from an existing boundary vertex on a to-be-extended surface boundary.

This method generates the extended triangulated irregular network (TIN) triangles by computing an extended or new surface boundary using extension vectors. As illustrated in Figure 1, the extension vectors are calculated in 3D space as the outer pointing cross product of the edge vectors and normal vectors of the associated TIN triangles. Their relationship is expressed by following equations: Ex1 = d(E1 x n1). 1.11 [21 Ex2 = d(E2 x n2).

where d is the extension distance. El and E2 are the unit vectors of the existing boundary edges, n1 and n2 are the unit normal vectors of the TIN triangles T1 and T2, respectively. The directions of the resultant extension vectors E1 and Ex2 point outwards from the existing surface edges. The normal vectors of the TIN triangles (not shown in Figure I) are located at the label positions, e.g., Ti and T2, corresponding to the TIN triangle centroids. At each label position, e.g., position of the label TI, the normal vector, e.g., T. projects perpendicularly from the TIN triangle face, e.g., TI and out of the page.

The calculated extension vectors can then be used to determine intermediate boundary vertices which are subsequently used to determine the new extended boundary vertex. For example, in the case of extending the existing boundary vertex V, positions of the two intermediate vertices V1 and V2 are determined by: V1= V + Ex1, [3] V2 = V + Ex2. [4] The new extended boundary vertex is located at the midpoint V12 between such t' o intermediate vertices VI and V2.

When the above-described surface extension method is implemented in geological models, e.g., to determine the extended surface boundary where the existing boundary is not straight or gently curved, extension vectors can cross and unwanted artefacts, such as overturned or overlapping TIN triangles, are often produced. Such artefacts are more likely to occur either when the existing surface comprises neighbouring boundary triangles, the normal vectors of which vary greatly in direction, or when the required surface extension distance is greater than three times of the average TIN triangle width. Figure 2(a) and 2(b) respectively show a plan view and a 3D view of an example surface extension of an irregular horizon boundary using the foregoing method. Both overturning and overlapping TIN triangles are produced. Depending on how irregular the existing surface boundary is, the above-described method is limited to allow only small extension distances, e.g., smaller than three times of the average TIN triangle width, for those surface boundaries with irregular geometries. In some cases, the extension vectors may cross one another, resulting in overturned sections of the extended surface boundary. This is evident in Figure 2 where many overturned boundary edges are present in the tangled areas of the extended surface.

Figure 3 shows an extreme example of this artefact problem where the foregoing method is used to extend a concave or embayed surface boundary. Figure 3(a) shows the original surface (in dark grey) with an irregular surface boundary. The embayment (in light grey) in the boundary of the surface is to be filled with the extended surface. Figure 3(b) shows the embaymcnt in the surface boundary being covered by a plurality of extension vectors calculated by using equations [1] and [2]. When surface boundaries to be extended comprise partially enclosed or embayed sections with acute and obtuse internal angles, the application of the foregoing method generates extension vectors that cross and result in the extended surface being generated with overlapping and/or overturned triangles. Figure 3(c) shows an extended surface (formed by TIN triangles) determined using the foregoing method. As shown in Figure 3(c), shading is applied to indicate degree of triangle overlapping or extension vector crossing in the extended surface. The empty space or empty TIN zone in the middle of the surface where no TIN triangles are generated results from partial surface extension. As mentioned in the Background Section, further data processing using manual editing techniques can be used to condition the problematic sections of the surface boundary before extending the surface. However, such manual methods are subjective and time-consuming to implement.

In the present disclosure, a method that is based on dip-domain boundaries is proposed to extend boundaries of 3D horizons while preventing crossing of extension vectors. Dip-domain boundaries have proven to be useful for extending boundaries of horizons in geological models. After interpretation, a horizon may be represented by, for example, a vector based TIN which is constmcted by triangulating a set of vertices (points). The vertices (points) are connected to form a network of contiguous and non-overlapping triangular facets. By way of example, extension of a horizon boundary can be obtained by growing a TIN, which is constructed to represent the original horizon, toward a desired/predetermined new surface boundary. The growth of a TIN requires the determination of extension vectors and the subsequent generation of new triangular facets following the determined extension vectors. Described below are embodiments of the proposed method for conditioning horizon data (e.g., extension of existing (or source) 3D horizons) used in 3D geological models. The embodiments disclosed here use dip-domain boundaries to avoid common mesh artefacts caused by the crossing of extension vectors and therefore can accurately, and efficiently determine new extending boundaries. This is achieved by the fact that boundary segments are utilized to determine local extension vectors for all existing boundary vertices that are to be extended.

Figure 4 shows a flowchart of the horizon or surface extension method in accordance with an embodiment. The horizon extension method may comprise for example six main steps, i.e. steps 401 to 406 of Figure 4.

At step 401, a boundary region of the interpreted 3D horizon represented by a vector based TIN (or other grid representation such as a regular grid/rectangular grid) is selected and subsequently transformed to approximately in the X-Y plane for horizon extension. According to an embodiment, the selection of the boundary region is achieved first by determining which grid polygons are boundary polygons (e.g., which triangles are boundary triangles in the case of a TIN). To be a boundary polygon/triangle, a polygon should comprise at least one boundary vertex, where a boundary vertex could be described as any vertex on the surface boundary to be extended. The following description will be explicitly applicable to TINs, but will be generally applicable to other grid representations in which case "triangle" should be read as "polygon".

Once the boundary triangles are all determined, computational methods may be used to search for all triangles that are within a pre-defined distance of the boundary triangles. The pre-defined distance may be equivalent to a total width of e.g., between 5 and 25 triangles, such as 20 triangles, 10 triangles, or 5 triangles. Figure 5(a) schematically illustrates such a selected boundary region in accordance with an embodiment. Note that, the selected boundary region is a 3D surface Once selected, at step 402, the boundary region may be transformed to be substantially in the X-Y plane so that it is ready for implementing boundary extension. In an embodiment, the transformation of the selected boundary region may comprise rotation of the region such that the resulting average normal vector n of the selected region is substantially parallel to the Z axis (e.g., the axis perpendicular to the surface being extended). As such, the selected boundary region substantially overlaps with the X-Y plane.

in a typical embodiment, such as the embodiment shown in Figure 5, the transformation of the selected boundary region may comprise for example two steps: at step 1, the selected boundary region is rotated about the Z-axis so that the average normal vector n of the selected region is aligned with the X-Z plane (e.g., rotation of angle a in Figure 5(a)); at step 2, the selected boundary region is further rotated about the Y-axis so that the avenge normal vector n of the selected region coincides with the Z axis (e.g., rotation of angle 13 in Figure 5(a)).

Note that, the mutually orthogonal axes X. Y, and Z are defined here according to the exemplary coordinate system shown in Figure 5(a); different orientations of the boundary region are also applicable. Note that, the embodiments disclosed herein will work as long as bisectors generated during the process of surface extension are projectable in a third direction orthogonal to a plane oriented in a first and second direction defined respectively by the transformed initial surface boundary description long-axis and the average normal vector of the selected boundary region, i.e. the selected boundary region is extended substantially in 2D or an approximation of a 2D representation.

The average normal vector of a surface region may be calculated by first summing together all the normal vectors of triangular facets and then dividing the sum by the total number of triangular facets. A boundary' segment may comprise a portion or length of boundary with no change of direction (e.g., comprising a ID flue). The normal of an individual triangle with vertices e.g., vO, vl, v2, is given as the cross product of (v2-v1) and (v2-v0). Note that, triangles are arranged such that vertices vO, v 1, v2 are arranged in a consistent order, and consequently the calculated normal vectors point in consistent directions.

Figure 5(b) shows a plan view of the boundary region after the aforementioned two-step transformation. Now, the average normal vector of the region stays in the X-Z plane and is parallel to the Z axis of the coordinate system. In other words, the selected 3D boundary region overlaps substantially with the horizontal plane or X-Y plane and the extension direction is parallel to the X axis.

At step 403, the transformed boundary region may be used as a source surface for surface extension. Note that, the transformed boundary region is still a 3D surface. According to an embodiment, step 403 may be implemented in several sub-steps. Figure 6 shows a flowchart of an implementation of step 403 in accordance with an embodiment. In this embodiment, step 403 may further comprise six sub-steps, i.e. steps 601 to 607.

At step 601, an initial set of extension vectors may be determined and subsequently generated on the boundary vertices of the surface boundary to be extended, one extension vector being generated on each boundary vertex. The direction of each initial extension vector (hereafter referred to as a "bisector") may be the direction that bisects the angle subtended by the surface boundary segments incident on each boundary vertex. Note that, the direction of each initial extension vector may not necessarily be the final extension direction for the corresponding boundary vertex. The final extension directions may be determined in combination of the defined extension distance and whether any adjacent bisectors intersect within the defined extension distance.

Figure 7(a) schematically illustrates the surface boundary segments and the initial bisectors (or extension vectors) that are determined by the surface boundary segments. As shown in Figure 7(a), the surface boundary comprises a plurality of boundary segments, e.g., Si to S6. The boundary segments are connected via a plurality of boundary vertices, e.g., VI to VI Each pair of adjacent boundary segments fonns a folding angle that is used to determine the direction of a corresponding bisector. For example, the boundary segments SI and S2 are connected by the vertex VI and thereby form the folding angle Al2. The folding angle Al2 is used to determine the direction of the bisector Fl. Similar considerations apply to other boundary segments for determining initial bisectors thereof.

At step 602, the plurality of bisectors may he evaluated sequentially or in parallel to see if any two adjacent bisectors intersect along the extension direction. If one or more intersections are confirmed, the flow proceeds to step 603; or else the flow proceeds to step 604.

At step 603, if any two adjacent bisectors intersect along the extension direction and within the target extension distance, they do not cross, but instead a new bisector may be created. By way of an example, Figure 7(b) and 7(c) schematically illustrate two examples of determination of the direction of a new bisector. The direction of the new bisector F6 bisects the angle A13 subtended by the direction of the segment SI on the immediate left of the left bisector F I and the segment S3 on the immediate right of the right bisector F2; i.e., the segments SI, S3 adjacently on either side of the segment S2 between the two boundary vertices VI, V2 at the base of the adjacent bisectors F1, F2.

More specifically, the two initial bisectors F I and F2 intersect before the target extension distance is reached, resulting in generation of an intersecting point 1P12. At the intersecting point 1P12, a new folding angle Al3 is subtended by the direction of the segment Si on the immediate left of the left bisector F I and the direction of the right segment S3 on the immediate right of the right bisector F2. The new folding angle A13 is then used to determine the direction of the new bisector F6 (e.g., in the direction bisecting the folding angle A13). According to an embodiment, at the time a new bisector is generated, step 602 may be repeated to check if the newly generated bisector intersects the adjacent bisector before the target extension distance is reached, if more intersections are found. step 603 is repeated.

Figure 7(c) schematically illustrates an example of such a scenario. As shown in Figure 7(c), before the target extension distance is reached, the new bisector F6 further intersects the adjacent bisector F3 resulting in generation of an intersecting point IP 14. At the intersecting point TP14, a new folding angle Al4 is subtended by the direction of the left segment Si on the immediate left of bisector F6 and the direction of the right segment S4 on the immediate right of bisector F3. The new folding angle A14 is subsequently used to determine the direction of the new bisector F7 (e.g. as bisecting this angle).

Once all the bisectors are generated as a result of mutual intersections within the target extension distance, the extended or new surface boundary may be created at step 604 by calculating segment vertices for the extended surface boundary at the target extension distance.

Figure 8 schematically illustrates the process of determination of new boundary vertices in accordance with an embodiment. When the target extension distance is reached and the initial bisector incident on an existing (or source) vertex results in at least one intersection with another bisector within such distance, the extension direction of the vertex follows the bisector path which comprises the initial bisector as well as newly generated bisectors.

For example, the bisector path of the vertex VI is shown as a bold dashed line in Figure 8. The bisector path comprises the initial bisector Fl and the new bisectors F6 and F7 which are generated as a result of intersections at the intersecting points TP12 and IP14 respectively. As such, the bisectors [F1, F6, F7] or [F2, F6, F7] or [F3, F7] together define a bisector path between the source vertex V1 and new extended vertex V1'.

The method comprises determining each new or extended vertex of the extended boundary at an extension distance d. This extended vertex VI' is defined at a distance along the final bisector F7 (e.g., at extension distance) which is determined as follows. One (or both) of the outermost initial bisectors (i.e., F I or F3 in this example, where the outermost initial bisectors are those drawn from the outermost source vertices V1, V3 contributing to the final bisector) which have contributed to the final bisector (i.e., bisector F7) in the corresponding bisector paths is selected. The selected outermost initial bisector is then extended by the extension distance d. In the example shown, this has been done to initial bisector Fl. The extension distance d is thus the distance between the source vertex VI and a point P on the extended initial bisector. The extended vertex VI' is then defined at the point on the final bisector F7, which when projected to the extended initial bisector (i.e., at a 90 degree angle from the extended initial bisector) corresponds to the point P on the extended initial bisector. Of course, where the final bisector is an initial bisector (e.g., the initial bisector is extended the full extension distance without intersecting another) then this point on the final bisector in a direction perpendicular to a point P will be the same point as point P. In an embodiment, once the extended boundai-y vertex of a particular source vertex is determined, the source boundary segments involved in determination of the extended boundary vertex may be identified. Following that, the next source vertex to be extended may be determined as the one to the right of the rightmost of the identified boundary segments. For example, referring to Figure 8 again, boundani segments Si. S2, S3 and S4 are all involved in determining the new extended boundary vertex V1'. Accordingly, the rightmost boundary segment is S4 and the source vertex on the right of the boundary segment S4 is V4. Therefore, the next source vertex to be extended is the vertex V4. Note that, such a determination process is equally applicable for the next source vertex to the left of the leftmost of the identified boundary segments. The next source vertex to the left of the leftmost of the identified boundary segments may be determined at the same time as the next source vertex to the right of the rightmost of the identified boundary segments.

In the case where the bisector on a particular source vertex does not intersect any adjacent bisectors, e.g., the initial bisector F4 on the source vertex V4 in Figure 8, the bisector path of the source vertex follows the direction of the initial bisector on that source vertex and hence overlaps with tlw reference line extended from the same initial bisector. As such, the extended or new boundary vertex corresponds to a point on the reference line where its distance to the source vertex equals the extension distance Once all the existing or source vertices have been extended, step 606 may be performed such that the resultant extended vertices are systematically joined or connected to form an extended surface boundary line. Between any two adjacent extended boundary vertices, there exists an extended boundary segment. Figure 9(a) shows an example extended boundary line which is extended from the boundary of the transformed region of Figure 5(b) and is generated after following the step 403 of Figure 4 (or steps 601 to 606 of Figure 6).

The extended surface boundary and the original boundary may form a polygonal outline. At step 404, the extended surface area defined by the polygonal outline may be tessellated in order to create a suitable TIN. Hence, vertices of the triangles may be determined so as to create suitable triangles which form the desired TIN. According to an embodiment, a first grid of points may be produced from the polygonal outline. The grid size may be adjustable to cover the bounding rectangle of the polygon (or the area defined by the polygonal outline). The number of grid points may be determined such that the grid spacing is similar to or substantially the same as the average segment size of the polygonal outline.

Once the first grid is established, the grid points may be used as the vertices for subsequent area tessellation. In an embodiment, the extended surface, i.e, the area defined by the polygonal outline, is constructed by creating, for example a constrained 2D Delaunay Triangulation with the segments of the outer polygon being constraints and the grid points inside the bounding polygon being vertices of triangles. Figure 9(b) shows an example tessellation of the extended surface area defined by the polygonal outline of Figure 9(a). After tessellation, the extended surface (light grey area in Figure 9(b)) is represented by a TIN comprising a plurality of triangles. Note that, the extended surface created at the end of step 404 is a 2D surface in e.g., the X-Y plane.

To convert the foregoing 2D extended surface (e.g., in the X-Y plane) to a 3D surface, Z values may be determined for all vertices of surface triangles of the extended surface arca. This may be implemented at step 405 which may further comprise several sub-steps. Figure 10 shows a flowchart of an implementation of step 405 in accordance with an embodiment.

At step 1001, a horizontal (or X-Y plane) regular grid may be generated to sample the Z values of the surface triangles. in an embodiment, grid points of the horizontal regular grid may form a plurality of rectangular cells arranged to cover both the existing boundary region, e.g., dark grey region in Figure 9(b), and the extended surface region, e.g., the light grey region in Figure 9(b).

The number of (e.g., rectangular) cells on the regular grid may be determined such that the grid spacing is similar to or substantially the same as the average edge length of the TIN triangles. In an embodiment, each surface triangle may comprise at least one grid point of the regular grid.

At step 1002; a plurality of vertical rays, propagating along the Z axis, may be projected from the grid points to the existing 3D boundary region, e.g., light grey region in Figure 9(b). In an embodiment, a plurality of intersecting points may be generated on the TIN triangles of the existing boundary region.

At step 1003, the generated intersecting points on the existing surface may be used to determine the Z values of the grid points in the existing boundary region. In an embodiment, the Z values of the grid points may be determined by calculating the distance between the grid points to their corresponding intersecting points on the existing surface. Subsequently, the determined Z values may be assigned to the grid points in the existing boundary region.

At step 1004, the grid points which have been assigned with the Z values at step 1003 may be used to determine the Z values of the grid points in the extended surface. Note that, grid points with the Z values may be described as assigned grid points In an embodiment, step 1004 may be carried out in following sub-steps: Identifying border grid points which should be unassigned grid points around the border of the assigned area (where grid points are assigned with the Z values), and each of which should neighbour with at least one assigned grid point; ii. Identifying all the assigned grid points within an area defined by a search radius for each boundary grid point; Using the identified assigned grid points to create a best fit 3D plane; iv. Determining the Z value for the boundary grid point by projecting a ray from the boundary grid point to the best fit 3D plane such that an intersecting point is generated on the best fit 3D plane; Assigning the determined Z value to the corresponding boundary grid point; vi. Repeating sub-steps i. to v. using the newly assigned grid points undl all die unassigned grid points are assigned with Z values.

At step 1005, by the time all the grid points are assigned with Z values, the Z value of each triangle vertex of the extended surface may be calculated. According to an embodiment, calculation of the Z value of each triangle vertex may comprise hi-linear interpolation of the Z values of the nearest four assigned grid points. At the end of step 1005, a 3D extended surface may be created.

Finally, at step 406 of Figure 4, the whole surface region including the resultant 3D extended surface and the existing boundary region is transformed such that the existing boundary region is returned to its original location and orientation in the 3D geological model. In an embodiment, such transformation may comprise the steps described in step 402 in reverse.

To summarise, the methods disclosed herein significantly improve the accuracy and efficiency of conditioning of 3D horizon data used in 3D geological models. The methods use dip-domain boundaries, which intersect but do not cross, to avoid common mesh artefacts caused by crossing of extension vectors and hence can quickly and accurately determine new extended boundaries for an existing 3D horizon. As such, better physical modelling of a reservoir system is achieved.

One or more steps of the methods and concepts described herein may be embodied in the form of computer readable instructions for running on suitable computer apparatus, or in the form of a computer system comprising at least a storage means for storing program instructions embodying the concepts described herein and a processing unit for performing the instructions. As is conventional, the storage means may comprise a computer memory (of any sort), and/or disk drive, optical drive or similar. Such a computer system may also comprise a display unit and one or more input/output devices The concepts described herein find utility in all aspects of modelling, optimisation and prediction of hydrocarbon reservoir and well systems, and may aid in, and form part of, methods for extracting hydrocarbons from such hydrocarbon reservoir and well systems. For example, the methods herein may be used in the locating of a hydrocarbon reservoir or determining the form of the hydrocarbon reservoir. The determined form of the hydrocarbon reservoir can then be used, for example, in optimizing hydrocarbon recovery, for example, by selecting one or more of said different production strategies which are determined to maximize hydrocarbon production and/or minimize production costs. Different strategies may include optimizing or determining well placement and/or deciding on how the hydrocarbon is best extracted (e.g., optimize parameter values such as wellhead pressure or gas lift rate, types of producing fluid used.) it should be appreciated that the above description is for illustration only and other embodiments and variations may be envisaged without departing from the spirit and scope of the invention.

Claims (19)

1. CLAIMS1. A computer implemented method for defining a 3D extended surface boundary description of a subsurface volume, comprising: a) obtaining an initial surface boundary description and an extension distance describing the distance that the initial surface boundary description is to be extended; b) transforming the initial surface boundary description to obtain a transformed initial surface boundary: description which is oriented such that a long-axis of the transformed initial surface boundary description is substantially parallel to a first direction and an average normal vector of the transformed initial surface boundary description is substantially parallel to a second direction, the first direction and the second direction being mutually orthogonal and defining a first plane, and wherein the transformed initial surface boundary description comprises a plurality of source vertices with a boundary segment between each adjacent pair of source vertices; c) generating a bisector from each of said source vertices, each of these bisectors being source bisectors that are projectable in a third direction orthogonal to the first plane and bisecting an angle subtended by said boundary segments immediately adjacent its respective source vertex; d) extending each bisector until said bisectors are extended the extension distance, while generating a new bisector at each intersection point where pairs of adjacent said bisectors intersect within said extension distance, in place of the intersecting pair of bisectors; c) generating an extended boundary vertex on each bisector at the extension distance, said bisectors at the extension distance being final bisectors; generating said extended surface boundary description by connecting said extended boundary vertices; g) tessellating an extended boundary region defined between said extended surface boundary description and said transformed initial surface boundary description; and h) transforming a region defined by said extended surface boundary description and said initial surface boundary description such that said initial surface boundary description is in its original location.
2. 2. A computer implemented method as claimed in claim 1, wherein each extended boundary vertex is generated at a position on its respective final bisector defined as the point which when projected on an extended outermost source bisector, corresponds to a point at the extension distance along said extended outermost source bisector, wherein said extended outermost source bisector comprises at least one of the outermost pairs of source bisectors which contributes to the generation of said final bisector, having been extended the extension distance.
3. 3 A computer implemented method as claimed in claim 2, wherein, the direction of each bisector generated at each of said intersection points is that which bisects an angle subtended by the directions of the boundary segments on the immediate outside of the outermost pairs of source bisectors which contribute to the generation of the pair of intersecting bisectors defining said intersection poi nt.
4. 4. A method as claimed in any preceding claim wherein said step of obtaining an initial surface boundary description comprises: obtaining an initial polygonal grid representation of the subsurface volume; determining a boundary region of the polygonal grid representation by determining which grid polygons of the grid are boundary polygons which comprise at least one vertex on said surface boundary to be extended; and determining the initial surface boundary description from said boundary region.
5. 5. A method as claimed in claim 4, comprising defining the boundary region as extending a predefmcd distance from the boundary polygons.
6. 6. A method as claimed in claim 5, wherein the predefined distance is determined by a predefined number of polygons.
7. 7. A method as claimed in any preceding claim, wherein step b) comprises: rotating said initial surface boundary description around said second direction such that said average normal vector of said transformed initial surface boundary description is substantially parallel to a second plane defined by said second direction and said third direction; and rotating said initial surface boundary description around said first direction such that said average normal vector of said transformed initial surface boundary description is substantially parallel to said second direction.
8. 8. A computer implemented method as claimed in claim 7, wherein step g) further comprises determining values for said extended boundary region in said second direction.
9. 9. A computer implemented method as claimed in claim 8, wherein said determination of said values for said extended boundary region comprises: generating a second polygonal grid of points in a third plane defined by said first direction and said third direction, said second grid of points being arranged to cover both said boundary region and said extended boundary region; projecting a plurality of rays in said second direction from said second grid of points to said boundary region and subsequently intersecting said boundary region, wherein at least one ray is projected from each grid point of said second arid; determining a value in said second direction for said each grid point of said second grid that intersects with said boundary region; and identifying border grid points of said second grid, wherein each of said border grid points are those which have no value in said second direction but are neighboured with at least one grid point of said second grid having said value in said second direction.
10. 10. A computer implemented method as claimed in claim 9, further comprising: for each of said border grid points, defining a search area with a search radius and searching within said search area for all grid points of said second grid having said values in said second direction; for at each of said border grid points, calculating a best fit 3D-plane for said search area defined by said search radius; from each of said border grid points, projecting a ray in said second direction towards said best fit 3D-plane and subsequently intersecting said best fit 3D-plane, and determining a value in said second direction for each of said border grid points of said second grid that intersects said best fit 3D-plane.
11. 11. A computer implemented method as claimed in claim 10 wherein transforming saidextended surface boundary description comprises:orienting said extended surface boundary description to match orientation of said initial surface boundary description prior to it being transformed.
12. 12. A computer implemented method as claimed in any preceding claim, wherein each boundary segment comprises a substantially straight boundary portion.
13. 13. A computer implemented method as claimed in any preceding claim, wherein said subsurface volume comprises a hydrocarbon reservoir and the method comprises a step of using said extended surface boundary description in optimizing hydrocarbon recovery from said hydrocarbon reservoir.
14. 14. A computer implemented method as claimed in claim 13, comprising locating said hydrocarbon reservoir using said extended surface boundary description.
15. 15. A computer implemented method as claimed in claim 13 or 14, wherein the method comprises using said extended surface boundary description to determine a form of said subsurface 20 volume.
16. 16. A computer implemented method as claimed in any preceding claim herein said step of generating a bisector comprises generating said bisector in said third plane.
17. 17. A computer implemented method as claimed in any preceding claim, comprising using said form of said subsurface volume to optimize a production strategy so as to maximize hydrocarbon production and/or minimize production costs using said determined form of the hydrocarbon reservoir.
18. IS. A computer program comprising computer readable instructions which, when run on suitable computer apparatus, cause the computer apparatus to perform the method of any preceding claim.
19. 19. A computer program product comprising the computer program of claim 18.