Patent Description:
Many subsurface formations include some degree of fracturing, i.e., the presence of faults, joints, cracks and other discontinuities that separate rock within the subsurface formation. Fractures generally have greater permeability and porosity than solid rock, so accounting for the effects of fractures is generally desirable for accurate fluid flow simulation. In this regard, a number of different fracture abundance measures have been proposed to represent the relative amount of fracturing within a subsurface formation, including, for example, fracture density, fracture intensity, fracture porosity, etc. Some conventional approaches, for example, calculate a fracture density as a P<NUM> value (number of fractures per unit length along a scanline) from wells. In addition, in some approaches a P<NUM> value (sum of fracture area per unit volume) is inferred from the P<NUM> value by making an assumption that fractures entirely intersect a borehole as well as corrected from borehole deviation and then using a statistical method to interpolate P<NUM> in a three-dimensional (3D) grid as an input for Discrete Fracture Network (DFN) generation.

The P<NUM> value is desirable in many applications because fracture size is accounted for in the value and does not depend on borehole trajectory. However, accurate fracture sizes within a borehole are generally difficult to obtain from borehole images and core logging, and generally result in the calculation of only a "relative" P<NUM> measurement from wells. Furthermore, interpolation of this measurement generally creates large uncertainties within the 3D grid that generally cannot be easily quantified.

Therefore, a need exists in the art for improved evaluation of P<NUM> and other fracture abundance parameters, and in particular, an improved evaluation having greater accuracy and/or greater computational efficiency than convention approaches. <CIT> describes a system, method and computer readable medium capable of extracting, characterizing and modeling fracture networks in a subterranean formation. Seismic data points indicating one or more discontinuities in the subterranean formation may be identified using processed seismic data and arranged according to a tree structure. The connected discontinuity points are utilized to generate discontinuity planes or surfaces by interpolation. A 3D mesh is utilized to visualize the extracted discontinuity planes or surfaces, where horizontal planes intersect the discontinuity planes or surfaces at sample lines, the discontinuity points of adjacent sample lines are linked to form a 3D mesh.

According to the present invention there is provided a method of evaluating fracture abundance in a subsurface formation according to claim <NUM>, an apparatus for evaluating fracture abundance in a subsurface formation according to claim <NUM>, and a program product for evaluating fracture abundance in a subsurface formation according to claim <NUM>. One or more preferred embodiments are defined in the appended dependent claims. Thus, the present invention is defined by the appended claims. Embodiments not falling within the scope of these claims are exemplary only.

For a better understanding of the invention, and of the advantages and objectives attained through its use, reference should be made to the Drawings, and to the accompanying descriptive matter, in which there is described example embodiments of the invention.

The herein-described embodiments utilize a number of techniques to evaluate fracture abundance in a subsurface formation, e.g., a region or volume of the Earth such as a volume potentially incorporating recoverable hydrocarbons. Fracture abundance, in particular, may be evaluated using a three-dimensional fracture network defined using a plurality of geometric primitives such as triangular elements disposed within a three-dimensional volume, and the evaluation of fracture abundance may result in the generation of one or more fracture abundance parameters for the subsurface formation. A fracture abundance parameter, within the context of the invention, may be any parameter that is indicative of the abundance or amount of fracturing in a subsurface volume, e.g., based on fracture intensity, fracture density, fracture porosity, etc. In some embodiments, a fracture abundance parameter may be based, for example, on ratios between different dimensional values, e.g., various PXY values, where x is the dimension of the measured value or feature and y is the dimension of the sampling region. In the embodiments discussed below, for example, a fracture abundance parameter may be a P<NUM> value of fracture density, based on the sum of the areas of fractures in a unit volume such as a grid cell. As other measurements or parameters may be used to represent relative amounts of fracturing in a subsurface volume, however, the invention is not limited to P<NUM> fracture density values.

Other variations and modifications will be apparent to one of ordinary skill in the art.

Turning now to the drawings, wherein like numbers denote like parts throughout the several views, <FIG> illustrates an example data processing system <NUM> in which the various technologies and techniques described herein may be implemented. System <NUM> is illustrated as including one or more computers <NUM>, e.g., client computers, each including a central processing unit (CPU) <NUM> including at least one hardware-based processor or processing core <NUM>. CPU <NUM> is coupled to a memory <NUM>, which may represent the random access memory (RAM) devices comprising the main storage of a computer <NUM>, as well as any supplemental levels of memory, e.g., cache memories, non-volatile or backup memories (e.g., programmable or flash memories), read-only memories, etc. In addition, memory <NUM> may be considered to include memory storage physically located elsewhere in a computer <NUM>, e.g., any cache memory in a microprocessor or processing core, as well as any storage capacity used as a virtual memory, e.g., as stored on a mass storage device <NUM> or on another computer coupled to a computer <NUM>.

Each computer <NUM> also generally receives a number of inputs and outputs for communicating information externally. For interface with a user or operator, a computer <NUM> generally includes a user interface <NUM> incorporating one or more user input/output devices, e.g., a keyboard, a pointing device, a display, a printer, etc. Otherwise, user input may be received, e.g., over a network interface <NUM> coupled to a network <NUM>, from one or more external computers, e.g., one or more servers <NUM> or other computers <NUM>. A computer <NUM> also may be in communication with one or more mass storage devices <NUM>, which may be, for example, internal hard disk storage devices, external hard disk storage devices, storage area network devices, etc..

A computer <NUM> generally operates under the control of an operating system <NUM> and executes or otherwise relies upon various computer software applications, components, programs, objects, modules, data structures, etc. For example, a petro-technical module or component <NUM> executing within an exploration and production (E&P) platform <NUM> may be used to access, process, generate, modify or otherwise utilize petro-technical data, e.g., as stored locally in a database <NUM> and/or accessible remotely from a collaboration platform <NUM>. Collaboration platform <NUM> may be implemented using multiple servers <NUM> in some implementations, and it will be appreciated that each server <NUM> may incorporate a CPU, memory, and other hardware components similar to a computer <NUM>.

In one non-limiting embodiment, for example, E&P platform <NUM> may implemented as the PETREL Exploration & Production (E&P) software platform, while collaboration platform <NUM> may be implemented as the STUDIO E&P KNOWLEDGE ENVIRONMENT platform, both of which are available from Schlumberger Ltd. and its affiliates. It will be appreciated, however, that the techniques discussed herein may be utilized in connection with other platforms and environments, so the invention is not limited to the particular software platforms and environments discussed herein.

In general, the routines executed to implement the embodiments disclosed herein, whether implemented as part of an operating system or a specific application, component, program, object, module or sequence of instructions, or even a subset thereof, will be referred to herein as "computer program code," or simply "program code. " Program code generally comprises one or more instructions that are resident at various times in various memory and storage devices in a computer, and that, when read and executed by one or more hardware-based processing units in a computer (e.g., microprocessors, processing cores, or other hardware-based circuit logic), cause that computer to perform the steps embodying desired functionality. Moreover, while embodiments have and hereinafter will be described in the context of fully functioning computers and computer systems, those skilled in the art will appreciate that the various embodiments are capable of being distributed as a program product in a variety of forms, and that the invention applies equally regardless of the particular type of computer readable media used to actually carry out the distribution.

Such computer readable media may include computer readable storage media and communication media. Computer readable storage media is non-transitory in nature, and may include volatile and non-volatile, and removable and non-removable media implemented in any method or technology for storage of information, such as computer-readable instructions, data structures, program modules or other data. Computer readable storage media may further include RAM, ROM, erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other solid state memory technology, CD-ROM, DVD, or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium that can be used to store the desired information and which can be accessed by computer <NUM>. Communication media may embody computer readable instructions, data structures or other program modules. By way of example, and not limitation, communication media may include wired media such as a wired network or direct-wired connection, and wireless media such as acoustic, RF, infrared and other wireless media. Combinations of any of the above may also be included within the scope of computer readable media.

Various program code described hereinafter may be identified based upon the application within which it is implemented in a specific embodiment of the invention. However, it should be appreciated that any particular program nomenclature that follows is used merely for convenience, and thus the invention should not be limited to use solely in any specific application identified and/or implied by such nomenclature. Furthermore, given the endless number of manners in which computer programs may be organized into routines, procedures, methods, modules, objects, and the like, as well as the various manners in which program functionality may be allocated among various software layers that are resident within a typical computer (e.g., operating systems, libraries, API's, applications, applets, etc.), it should be appreciated that the invention is not limited to the specific organization and allocation of program functionality described herein.

Furthermore, it will be appreciated by those of ordinary skill in the art having the benefit of the instant disclosure that the various operations described herein that may be performed by any program code, or performed in any routines, workflows, or the like, may be combined, split, reordered, omitted, and/or supplemented with other techniques known in the art, and therefore, the invention is not limited to the particular sequences of operations described herein.

Those skilled in the art will recognize that the example environment illustrated in <FIG> is not intended to limit the invention. Indeed, those skilled in the art will recognize that other alternative hardware and/or software environments may be used without departing from the scope of the invention.

<FIG> illustrate simplified, schematic views of an oilfield <NUM> having subterranean formation <NUM> containing reservoir <NUM> therein in accordance with implementations of various technologies and techniques described herein. <FIG> illustrates a survey operation being performed by a survey tool, such as seismic truck <NUM>, to measure properties of the subterranean formation. The survey operation is a seismic survey operation for producing sound vibrations. In <FIG>, one such sound vibration, sound vibration <NUM> generated by source <NUM>, reflects off horizons <NUM> in earth formation <NUM>. A set of sound vibrations is received by sensors, such as geophone-receivers <NUM>, situated on the earth's surface. The data received <NUM> is provided as input data to a computer <NUM> of a seismic truck <NUM>, and responsive to the input data, computer <NUM> generates seismic data output <NUM>. This seismic data output may be stored, transmitted or further processed as desired, for example, by data reduction.

<FIG> illustrates a drilling operation being performed by drilling tools <NUM> suspended by rig <NUM> and advanced into subterranean formations <NUM> to form wellbore <NUM>. Mud pit <NUM> is used to draw drilling mud into the drilling tools via flow line <NUM> for circulating drilling mud down through the drilling tools, then up wellbore <NUM> and back to the surface. The drilling mud may be filtered and returned to the mud pit. A circulating system may be used for storing, controlling, or filtering the flowing drilling muds. The drilling tools are advanced into subterranean formations <NUM> to reach reservoir <NUM>. Each well may target one or more reservoirs. The drilling tools are adapted for measuring downhole properties using logging while drilling tools. The logging while drilling tools may also be adapted for taking core sample <NUM> as shown.

Surface unit <NUM> may also collect data generated during the drilling operation and produces data output <NUM>, which may then be stored or transmitted.

Generally, the wellbore is drilled according to a drilling plan that is established prior to drilling. The drilling plan sets forth equipment, pressures, trajectories and/or other parameters that define the drilling process for the wellsite. The drilling operation may then be performed according to the drilling plan. However, as information is gathered, the drilling operation may need to deviate from the drilling plan. Additionally, as drilling or other operations are performed, the subsurface conditions may change. The earth model may also need adjustment as new information is collected.

Surface unit <NUM> may include transceiver <NUM> to allow communications between surface unit <NUM> and various portions of the oilfield <NUM> or other locations. Surface unit <NUM> may also be provided with or functionally connected to one or more controllers (not shown) for actuating mechanisms at oilfield <NUM>. Surface unit <NUM> may then send command signals to oilfield <NUM> in response to data received. Surface unit <NUM> may receive commands via transceiver <NUM> or may itself execute commands to the controller. A processor may be provided to analyze the data (locally or remotely), make the decisions and/or actuate the controller. In this manner, oilfield <NUM> may be selectively adjusted based on the data collected. This technique may be used to optimize portions of the field operation, such as controlling drilling, weight on bit, pump rates, or other parameters. These adjustments may be made automatically based on computer protocol, and/or manually by an operator. In some cases, well plans may be adjusted to select optimum operating conditions, or to avoid problems.

As shown, the sensor (S) may be positioned in production tool <NUM> or associated equipment, such as christmas tree <NUM>, gathering network <NUM>, surface facility <NUM>, and/or the production facility, to measure fluid parameters, such as fluid composition, flow rates, pressures, temperatures, and/or other parameters of the production operation.

While <FIG> illustrate tools used to measure properties of an oilfield, it will be appreciated that the tools may be used in connection with non-oilfield operations, such as gas fields, mines, aquifers, storage, or other subterranean facilities.

Part, or all, of oilfield <NUM> may be on land, water, and/or sea.

Data plots <NUM>-<NUM> are examples of static data plots that may be generated by data acquisition tools <NUM>-<NUM>, respectively, however, it should be understood that data plots <NUM>-<NUM> may also be data plots that are updated in real time. These measurements may be analyzed to better define the properties of the formation(s) and/or determine the accuracy of the measurements and/or for checking for errors. The plots of each of the respective measurements may be aligned and scaled for comparison and verification of the properties.

Static data plot <NUM> is a seismic two-way response over a period of time. Static plot <NUM> is core sample data measured from a core sample of the formation <NUM>. The core sample may be used to provide data, such as a graph of the density, porosity, permeability, or some other physical property of the core sample over the length of the core. Tests for density and viscosity may be performed on the fluids in the core at varying pressures and temperatures. Static data plot <NUM> is a logging trace that generally provides a resistivity or other measurement of the formation at various depths.

A production decline curve or graph <NUM> is a dynamic data plot of the fluid flow rate over time. The production decline curve generally provides the production rate as a function of time. As the fluid flows through the wellbore, measurements are taken of fluid properties, such as flow rates, pressures, composition, etc..

While a specific subterranean formation with specific geological structures is depicted, it will be appreciated that oilfield <NUM> may contain a variety of geological structures and/or formations, sometimes having extreme complexity. In some locations, generally below the water line, fluid may occupy pore spaces of the formations. Each of the measurement devices may be used to measure properties of the formations and/or its geological features. While each acquisition tool is shown as being in specific locations in oilfield <NUM>, it will be appreciated that one or more types of measurement may be taken at one or more locations across one or more fields or other locations for comparison and/or analysis.

The data collected from various sources, such as the data acquisition tools of <FIG>, may then be processed and/or evaluated. Generally, seismic data displayed in static data plot <NUM> from data acquisition tool <NUM> is used by a geophysicist to determine characteristics of the subterranean formations and features. The core data shown in static plot <NUM> and/or log data from well log <NUM> are generally used by a geologist to determine various characteristics of the subterranean formation. The production data from graph <NUM> is generally used by the reservoir engineer to determine fluid flow reservoir characteristics. The data analyzed by the geologist, geophysicist and the reservoir engineer may be analyzed using modeling techniques.

<FIG> illustrates an oilfield <NUM> for performing production operations in accordance with implementations of various technologies and techniques described herein. As shown, the oilfield has a plurality of wellsites <NUM> operatively connected to central processing facility <NUM>. The oilfield configuration of <FIG> is not intended to limit the scope of the oilfield application system. Part or all of the oilfield may be on land and/or sea. Also, while a single oilfield with a single processing facility and a plurality of wellsites is depicted, any combination of one or more oilfields, one or more processing facilities and one or more wellsites may be present.

As noted above, evaluation of fracture abundance parameters such as P<NUM> fracture density may be limited in conventional approaches due in part to difficulties associated with accurately accounting for fracture sizes from borehole or other formation data. Embodiments consistent with the invention, on the other hand, may utilize a three-dimensional approach incorporating various features that facilitate evaluation of fracture abundance in a subsurface formation in a more computationally efficient and accurate manner than such approaches.

Embodiments consistent with the invention, in particular, are based in part upon a determination of the areas of geometric primitives that are used to represent a fracture network within a three-dimensional volume, e.g., within a three-dimensional model of a subsurface formation. The geometric primitives are implemented as two-dimensional triangles defined by collections of three points in the three-dimensional volume, although other two-dimensional shapes may be used as geometric primitives in other examples not forming part of the present invention. The areas furthermore are determined on a subvolume-by-subvolume basis, e.g., with primitives that fall entirely within a subvolume having areas corresponding to the areas of the entire primitives, and with primitives that fall partially within a subvolume being clipped at the boundaries of the subvolume such that the areas are of the clipped portions of the primitives. In some embodiments, for example, a fracture network may be overlaid into a regular grid of cubic cells, and as such, areas of primitives representing a fracture network are determined on a cell-by-cell basis.

In the illustrated embodiments discussed hereinafter, the areas of each of the primitives (or clipped portions thereof) within each cell are summed, and then a ratio may be taken against the volume of each cell to generate a "real" P<NUM> fracture density for each cell. In contrast with many conventional approaches, the P<NUM> fracture density may be a more directly-calculated or "true" value rather than an inferred and interpolated value, and may be used to better constrain and validate a fracture network and/or to define correction factors to correct any relative P<NUM> values inferred from a borehole and/or interpolated in a 3D grid. As noted above, however, the techniques described herein may be used to calculate other fracture abundance parameters, so the invention is not limited to the particular fracture density calculations discussed herein.

<FIG>, for example, illustrates a sequence of operations <NUM> capable of being implemented in data processing system <NUM> to evaluate fracture abundance within a subsurface formation. In block <NUM>, a fracture network is defined within the cells of a three-dimensional model using a plurality of geometric primitives. As will become more apparent below, in the embodiments, a fracture network is generated based on data, for example in part from outcrop and/or seismic data, which are used to generate two-dimensional (2D) polylines. In addition, in the embodiments, the fracture network is generated by growing or expanding the 2D polylines in one or more directions and representing the resulting shapes using geometric primitives, with all polylines grown in the same direction (e.g., horizontally or vertically), or with different polylines grown in different directions. A fracture network may be generated manually (e.g., through a computer interface) and/or based on collected data (e.g., from geomechanical properties).

A polyline, in this regard, refers to a line comprised of one or more line segments, and a two-dimensional polyline is a polyline comprised of one or more line segments that lie within the same plane. Thus, it will be appreciated that while the polylines are referred to as two-dimensional polylines, such polylines may still be one-dimensional entities in some instances, e.g., where such polylines include only one line segment or where the segments of such polylines extend along the same axis.

Next, in block <NUM>, the combined areas of the geometric primitives within at least a subset of the cells of the 3D model are determined by summing together the areas of individual geometric primitives within each of the cells. Then, in block <NUM>, a fracture abundance parameter is generated for the fracture network. For example, in some embodiments, the fracture abundance parameter may include a fracture density such as a P<NUM> value for one or more cells in the 3D model. In other embodiments, other values indicative of fracture abundance may be generated from the determined combined areas.

The fracture abundance parameter may then be used for various purposes in various embodiments of the invention. For example, as illustrated in block <NUM>, the fracture abundance parameter may be used in fluid flow simulation using the same or a different 3D model of the subsurface formation. Further, as illustrated in block <NUM>, the results of the fluid flow simulation may be used to perform various oilfield operations, e.g., drilling a production and/or injection well, developing a well plan, determining a well trajectory, managing production, mine planning, civil engineering (e.g.: slope stability, tunneling), geotechnical ground control applications, etc. In addition, a directly-calculated fracture abundance parameter such as the directly-calculated P<NUM> value described herein may be used in some embodiments to better constrain and validate a fracture network and/or define correction factors to correct a relative P<NUM> value inferred from borehole data and interpolated in a 3D grid.

It will also be appreciated that the information generated during the various operations described above may also be visualized, e.g., within a graphical tool provided in an E&P platform, including both visualization of a generated fracture network as well as visualization of fracture abundance parameters calculated therefor. Further, it will be appreciated that the various operations may be performed by different tools, and that the operations need not be performed by or within a single tool.

Now turning to <FIG>, as noted above a fracture network is defined in block <NUM> of <FIG> in various manners in different embodiments. <FIG> illustrates the sequence of operations that "grows" or expands three-dimensional fractures from two-dimensional polylines and represents those three-dimensional fractures using one or more geometric primitives. In particular, in block <NUM> fractures are input as a set of two-dimensional polylines, i.e., lines defined by two distinct points and thus having a length along an axis extending between those points, for example. The polylines may be generated, for example, from seismic data, from outcropping data, from borehole images, and, as discussed in greater detail below, from geomechanical data, among other sources. In some embodiments, for example, the 2D polylines may be defined within a common plane such as a two-dimensional map representing a planar slice taken through the subsurface formation, although the invention is not so limited.

Next, in block <NUM>, each 2D polyline is grown or expanded in a predetermined direction and with a predetermined shape and aspect ratio, i.e., a ratio that controls the amount of growth in the predetermined direction relative to the length of a polyline. The predetermined shape may be selected from different potential shapes capable of representing a fracture. In the embodiments, the shape is rectangular or elliptical, although examples not forming part of the present invention may not be so constrained. Rectangular shapes may be favored for performance reasons, while elliptical shapes may be favored for accuracy as many fractures have a profile more closely matching that of an ellipse.

It will be appreciated that in some embodiments, expanding or growing a polyline along a predetermined direction may be considered to include expanding or growing the polyline in two opposite directions, e.g., equidistant from the polyline, or in some instances, different distances from the polyline. Further, expanding or growing a polyline along a predetermined direction generally results in the polyline being expanded within a plane that contains the polyline, referred to herein as a containing plane for the polyline.

In different embodiments, a single direction, shape and aspect ratio may be used to grow all polylines, while in other embodiments the direction, shape and/or aspect ratio may be varied for different polylines. Furthermore, the direction, shape and/or aspect ratio may be manually input by a user in some embodiments, while in other embodiments, one or more of these inputs may be determined programmatically. In some embodiments, for example, an optimizer may apply different inputs to generate different three-dimensional fracture networks that may each be used to determine different fracture abundance parameters, and these different parameters may be used in fluid flow simulations and matched against collected data to determine the combination of inputs that best matches observed data. In addition, it may even be desirable to utilize a randomized approach to generate directions, shapes and/or aspect ratios for different polylines. As such, it will also be appreciated that the respective containing planes of different polylines may in some embodiments extend in a same direction or different directions relative to a common plane within which the polylines are disposed, and further, in some embodiments, the respective containing planes may be substantially orthogonal to such a common plane, e.g., being substantially vertical relative to a substantially horizontal common plane.

Next, in block <NUM>, each grown 2D polyline is transformed into a plurality of geometric primitives, namely triangles in the present invention, to represent the grown shape. Then in block <NUM>, the geometric primitives may optionally be output for visualization or other purposes. For example, in one embodiment, the geometric primitives may be output in a TSURF file format for import and display in the PETREL E&P platform.

With further reference to <FIG> further illustrate the transformation of 2D polylines into 3D sets of geometric primitives, and in the case where the 2D polylines are provided in the form of a horizontal digitized fracture map such as may be generated from seismic or outcrop data, and grown in a vertical direction. <FIG>, in particular, illustrates a fracture map <NUM> including a pair of intersecting 2D polylines <NUM>. <FIG> respectively illustrate the transformation of polylines <NUM> into rectangular, elliptical and circular shapes, respectively.

As illustrated in <FIG>, for example, polylines <NUM> may be grown into rectangular shapes <NUM>, each formed from four points or nodes <NUM> and two triangles <NUM>. Likewise, as illustrated in <FIG>, polylines <NUM> may alternately be grown into elliptical shapes <NUM>, each formed from <NUM> points or nodes <NUM> and <NUM> triangles <NUM>. <FIG> also illustrates circular shapes <NUM>, which are a special case of elliptical shapes where the radius An (semi-axis longest of an ellipse) and Bn (semi-axis shortest of an ellipse) are equal.

As noted above, the amount of growth may be constrained by an aspect ratio, and as such, The starting and ending points of a 2D polyline and the segment length may be used to expand a fracture according to an aspect ratio (Asp) as follows: <MAT> where an is the fracture length divided per two and Asp is the input aspect ratio.

Now returning briefly to blocks <NUM>-<NUM> of <FIG>, once a fracture network is defined, a fracture abundance parameter may be calculated for a fracture network in a number of different manners consistent with the invention. Further, in some embodiments, a 3D observation grid may be generated after defining a fracture network based upon the minimum and maximum longitude/latitude and depth of the fracture network, and using a selected unit cell size to provide the desired resolution for fracture abundance parameter calculations.

<FIG> next illustrates an example sequence of operations <NUM> suitable for implementing blocks <NUM> and <NUM>. In these examples, it assumed that the fracture network is derived from 2D polylines in a horizontal fracture map and grown in a vertical direction, and that the fracture abundance parameter calculated is a P<NUM> fracture density value for each cell within an observation grid.

The calculation of P<NUM> in the 3D grid may in some embodiments be performed by column along the vertical Z-axis, and in some embodiments, may organize or store geometric primitives from the fracture network in an octree or other spatially-organized data structure to optimize calculations. Sequence of operations <NUM> begins in block <NUM> by creating a new 3D grid property for a P<NUM> value. For each cell i,j,k (block <NUM>), an Area_Sum variable is initialized to zero (block <NUM>). Next, for each fracture in the fracture network (block <NUM>) and for each geometric primitive in the fracture (block <NUM>), the sequence calls the ClipInsideCell function (block <NUM>) to create a clipped primitive including only that portion of the primitive that is inside the current cell. Block <NUM> then calls an AreaFromProjections function (discussed in greater detail below) on the clipped primitive to calculate the area of the primitive, and the result of this function is added to the Area_Sum variable (block <NUM>). This process is then repeated for every primitive in every fracture, resulting in Area_Sum storing the combined areas of the fractures within the cell. As such, block <NUM> divides Area_Sum by the Unit Volume of the cell, resulting in the determination of the P<NUM> value for that cell. Each cell in the grid is thereafter processed in a similar manner, and the result is returned in block <NUM> as a 3D matrix of P<NUM> values.

In addition, as noted above, a generated fracture abundance parameter may be used for visualization, among other purposes. <FIG>, for example, illustrates two visualizations <NUM>, <NUM> of a 3D volume or grid used to calculate a P<NUM> fracture density value for the 2D polylines <NUM> illustrated in <FIG> and grown into rectangular shapes as illustrated in <FIG>. In visualization <NUM>, the distribution of the P<NUM> fracture density parameter throughout the full grid is illustrated, while in visualization <NUM>, the distribution is limited to a horizontal layer through the grid. Shadings or colors (mapped in legend <NUM>) denote the varying values of the P<NUM> parameter.

Now turning to <FIG>, in some embodiments of the invention, geometric primitive projection may be used to reduce the computational resources and/or the latency associated with fracture abundance parameter determinations, e.g., to implement the AreaFromProjections function discussed above. As noted above, fracture abundance parameter determinations may be based in part on determining a combined sum of the areas of geometric primitives such as triangular elements within each cell or volume of a three-dimensional observation grid. Consequently, in many embodiments, a fundamental operation that generally consumes a significant portion of the employed computational resources is evaluating the area of a triangle in a cell, with the understanding that the triangle can cut the cell and thus project at least partially outside of the cell. Conventional approaches to area determinations may use simplified fracture geometries, e.g., usually vertical and with a rectangular shape, in order to reduce computational resources for the entire grid (which in some instances may contain millions of cells). As such, the accuracy that may otherwise be achieved using more complex and more realistic shapes to model fractures (e.g., elliptical shapes and/or shapes that project in non-vertical directions) may need to be sacrificed in order to achieve practical runtimes on moderately powerful computer systems.

Some embodiments consistent with the invention, on the other hand, may incorporate primitive projection to accelerate the determination of the area of a geometric primitive such as a triangular element inside a cell, such that the evaluation of fracture abundance may include at least the operations of defining a fracture network within a plurality of cells of a three-dimensional model of a subsurface formation using a plurality of geometric primitives, determining an area of the plurality of geometric primitives within at least a subset of the plurality of cells by summing areas of individual geometric primitives within each of the subset of cells, including determining an area of a first geometric primitive among the plurality of geometric primitives within a first cell in the subset of cells by projecting the first geometric primitive onto each of first, second and third orthogonal planes respectively aligned with faces of the first cell to define respective first, second and third projections and calculating areas of each of the first, second and third projections, and determining a fracture abundance parameter for the fracture network from the determined area of the plurality of geometric primitives. In addition to reducing computation time and/or computational resources, the herein-described technique may also in some embodiments shift the barrier to evaluate any triangular element of any orientation in a 3D Cartesian space, such that rapid calculations may be made of planar triangular elements of any orientation, and generally without involving classical heavy trigonometric algorithms to calculate the area. Further, in some embodiments the area of any subsurface structure capable of being represented by triangular elements, e.g., faults, fractures, horizons, etc., may be determined in a fast and efficient manner using the herein-described techniques, so the herein-described techniques may also be used for evaluating a subsurface formation by in part defining a subsurface structure within a plurality of cells in a three-dimensional model of the subsurface formation using a plurality of geometric primitives, determining an area of the plurality of geometric primitives using projection in the manner described herein, and determining a subsurface structure parameter for the subsurface structure (e.g., a fracture abundance parameter for a fracture network, or another parameter suitable for the particular subsurface structure being modeled) from the determined area of the plurality of geometric primitives.

Primitive projection, in this regard may be considered to refer to an operation that projects a shape of any arbitrary orientation within a three-dimensional cell onto a plane that is aligned with a face of a regular cubic cell. It will be appreciated that a face of a regular cubic cell is generally parallel to a plane formed by two of the three axes of a three-dimensional Cartesian coordinate system, e.g., where points or nodes are identified by (x, y, z) values on mutually-orthogonal X, Y and Z axes, each cell will generally have two faces parallel with each of XY, YZ and ZX planes defined by the X, Y and Z axes. Projection onto a plane aligned with a face of a regular cubic cell may therefore include projection onto a plane that is either coextensive with or parallel to a face of a cell, and thus coextensive with or parallel to the XY, YZ or ZX planes defined for a grid of regular cubic cells.

Projection may be further explained within the context of <FIG>, in particular illustrates the determination of the length L of an arbitrarily oriented one-dimensional line within a two-dimensional plane through projection of the line onto the two orthogonal axes X and Y. The nodes or endpoints of the line have coordinates (x<NUM>, y<NUM>) and (x<NUM>, y<NUM>), and it may be seen that projecting the line onto the X and Y axes generates two projections having lengths of Lx = (x<NUM> - x<NUM>) and Ly = (y<NUM> - y<NUM>). Through the application of the Pythagorean theorem the length L of the line is related to the lengths of the projections by the relationship L<NUM> = Lx<NUM> + Ly<NUM>, so the length L may be determined by taking the square root of the sum of the lengths of the projections.

<FIG> illustrates an extension of this principle into a three-dimensional space, where the area A of an arbitrarily oriented two-dimensional planar shape <NUM> within a three-dimensional volume may be determined through projection of the shape onto the three orthogonal planes XY, ZX and YZ planes to form three projections <NUM>, <NUM> and <NUM> having areas Axy, Azx and Ayz respectively. The area A of shape <NUM>, in particular, is related to the areas of projections <NUM>, <NUM> and <NUM> based upon the relationship A<NUM> = Axy<NUM> + Azx<NUM> + Ayz<NUM>, so the area A may be determined by taking the square root of the sum of the areas of the projections.

A net effect of projecting a shape onto a plane aligned with a face of a regular cubic cell is that all points or nodes of the projection is the reduction of the 3D problem into a simple 2D problem, thus simplifying the determination of the area of a projection into a less computationally-expensive operation. The location of each point or node of a projection of a shape thus may be represented by the other two coordinates. Consequently, assuming that shape <NUM> is defined by three points (x<NUM>, y<NUM>, , z<NUM>), (x<NUM>, y<NUM>, z<NUM>) and (x<NUM>, y<NUM>, z<NUM>), the projections onto the three planes XY, ZX and YZ may be considered to be defined by points:.

Consequently, instead of performing complex computations to determine the area of a triangle inside a three-dimensional cell, the area may be determined by projecting the triangle onto the three principal planes and sum the squares of the resulting projected areas.

It should also be appreciated, however, that clipping may also be performed in connection with projection in order to determine the area of a primitive within a cell, generally prior to projecting the primitive. <FIG>, for example, illustrates the clipped portions of both the primitive <NUM> and each projection <NUM>-<NUM> using darker shading than the portions that fall outside of the cell. Consequently, the coordinates of each point of the primitives may be compared with the coordinates of the boundary of the cell to replace any point disposed outside of the boundary of the cell with one or more points along the boundary of the cell. For example, as illustrated by the XY projection <NUM>, assuming that the projection without clipping would have three points P<NUM>, P<NUM> and P<NUM> and with point P<NUM> lying outside of the cell, this point may be replaced by points P<NUM> and P<NUM> to create a projected shape defined by points P<NUM>, P<NUM>, P<NUM> and P<NUM>. In addition, while various area calculations may be used to determine the areas of clipped projections, in one embodiment the areas of clipped projections may be determined by splitting a clipped projection into a plurality of triangles and summing the areas of the triangles forming the clipped projection.

Now turning to <FIG>, this figure illustrates an example sequence of operations <NUM> for performing projection-based fracture abundance parameter calculations in the manner discussed above, and in particular to determine a P<NUM> fracture density value for each cell in an observation grid. Sequence <NUM>, in the illustrated embodiment, utilizes an octree or other data structure to spatially organize geometric primitives, here triangles, of a defined fracture network, and thereby facilitate lookup of the triangles that are within each cell. Various types of data structures and/or lookup algorithms, including algorithms based on other types of binary and/or spatially-partitioned trees, may be used in other embodiments.

Block <NUM> initially builds an octree of the fractures in the fracture network. Then, for each cell C of the grid (block <NUM>), the octree is accessed to generate a set F of all fractures that at least partially intersect the cell C (block <NUM>). A combined area variable AC is then reset (block <NUM>), and each fracture f in set F (block <NUM>), and each triangle t defining fracture f (block <NUM>) is processed by projecting the area of the triangle t onto planes XY, ZX and YZ (block <NUM>, <NUM> and <NUM>), with the areas of the projections stored in Axy, Azx and Ayz respectively. The areas are then summed and a square root is taken of the sum, with the result added to the combined area variable AC (block <NUM>). Then, for each cell C, the combined area variable Ac is divided by the unit volume of the cell Vc to generate the P<NUM> fracture density value for the cell (block <NUM>). As a result of sequence <NUM>, therefore, a P<NUM> fracture density value is generated for each cell of the grid.

Now turning to <FIG>, in some embodiments fracture abundance evaluation may be based at least in part on geomechanical simulation based on mechanical properties of a subsurface formation. For example, using an engine such as the FAULT MODELER engine available from Schlumberger Ltd. and its affiliates, fracture density and/or fracture height may be estimated from geomechanical simulation using mechanical properties of a subsurface formation, e.g., as collected from well logging. For example, a balance energy operation may be used to estimate a P<NUM> fracture density and fracture height from well log data in some embodiments.

<FIG>, in particular, illustrates an example sequence of operations <NUM> for evaluating fracture abundance based on geomechanical simulation based on mechanical properties. First, in block <NUM>, a P<NUM> fracture density and fracture height for a subsurface formation may be determined from mechanical properties derived from well log data, e.g., one or more of Young's modulus, Poisson's ratio, friction coefficient, cohesion, fault dip, effective vertical stress, fluid pressure, and crack surface energy, among others. Such calculations may be based on a balance energy approach, and may be performed, for example, using geomechanical simulation functionality available in the PETREL, TECHLOG or FAULT MODELER software available from Schlumberger and its affiliates, although it will be appreciated that the invention is not limited to use with such software.

Based upon the results of block <NUM>, an input file is created including 2D polyline representations of fractures associated with fracture height, and taken along a vertical trace corresponding to the well from which the geomechanical properties were obtained that extends along the well (block <NUM>), and then each 2D polyline is grown in a similar manner to that described above in connection with <FIG>, but with the direction of growth being horizontal for a vertical well (block <NUM>). An aspect ratio may also be used to control the growth, and various shapes, including rectangular or elliptical shapes. Again similar to <FIG>, the grown shapes are transformed into geometric primitives (block <NUM>), and then similar to <FIG>, a fracture abundance parameter for the fracture network, e.g., a P<NUM> fracture density value for each cell, is determined (block <NUM>). Fluid flow simulation (block <NUM>) and/or an oilfield operation (block <NUM>) may then be performed, again similar to <FIG>.

In addition, visualization may be also be performed after various operations from <FIG>. For example, <FIG> illustrates a visualization <NUM> of fracture growth around an example well based upon geomechanical properties, while visualization <NUM> illustrates a vertical section of an observation grid containing the directly-calculated "true" P<NUM> value, with legend <NUM> mapping the shading to different ranges of P<NUM> values.

Therefore, some embodiments of the invention support the evaluation of true P<NUM> fracture density based on a 1D geomechanical method. By estimating the true P<NUM>, reservoir simulations in some embodiments may be better constrained, particularly for vertical wells that may run parallel to the fracture network, and may not detect fracture at in some circumstances leading to an optimized production of the reservoir through flow simulation. In addition, in some embodiments, the results of the aforementioned evaluation may be used for other purposes, e.g., for derisking, drilling and fracture connectivity. This is particularly true in certain situations. For example, as illustrated in <FIG>, well log data from a vertical well A may detect one fracture along layer <NUM> while nothing is observed in layers <NUM>, <NUM> and <NUM>. Well B may detect a fracture on layer <NUM>, with nothing observed in in the other layers. Further, well C may not detect a fracture at all. The herein-described techniques may therefore address these situations and estimate 3D fracture abundance with better prediction along wells, which may be used to populate a 3D grid for fracture network generation instead of a direct conventional measure of fracture abundance such as a P10c, value.

Various modifications may be made in other embodiments. For example, it will be appreciated that three-dimensional fracture abundance evaluation in some embodiments may use area determination operations other than the herein-described projection-based area determinations, as well as that projection-based area determinations may have other applications beyond that of fracture abundance evaluation. Further, as noted above, various operations may be used to generate or define a fracture network, so the invention is not limited to the particular geomechanical simulation-based approach disclosed herein, and further, other three-dimensional fracture abundance evaluation approaches beyond those described herein may be used to evaluate fracture abundance using mechanical properties, as long as they fall within the scope of protection that is defined by the appended claims.

Claim 1:
A computer implemented method of evaluating fracture abundance in a subsurface formation, the method comprising:
a) defining (<NUM>) a fracture network within a plurality of cells of a three-dimensional model of the subsurface formation using a plurality of geometric primitives, where the fracture network and the three-dimensional model are generated based on data obtained from a survey of the subsurface formation, wherein defining (<NUM>) the fracture network includes expanding, within a respective containing plane, each of a plurality of two-dimensional polylines (<NUM>) representing fractures in the subsurface formation, wherein:
a two-dimensional polyline comprises one or more line segments that lie within a same plane;
each of the plurality of geometric primitives is a triangular element; and
expanding each of the plurality of two-dimensional polylines includes expanding, in the direction of the respective containing plane, first two-dimensional polyline among the plurality of two-dimensional polylines either into a substantially rectangular shape represented by first and second triangular elements defined by four nodes or into a substantially elliptical shape represented by twelve triangular elements defined by thirteen nodes;
b) determining (<NUM>) an area of the plurality of geometric primitives within at least a subset of the plurality of cells by summing areas of individual geometric primitives within each of the subset of cells, wherein the determining the area of the plurality of geometric primitives within the subset of the plurality of cells further includes:
organising (<NUM>) the plurality of geometric primitives within a spatially-organized data structure;
accessing (<NUM>) the spatially-organized data structure when summing areas of individual geometric primitives within each of the subset of cells to determine which of the plurality of geometric primitives are at least partially within each of the subset of cells; and
clipping individual geometric primitives that are partially within each of the subset of cells; and
c) determining (<NUM>) a fracture abundance parameter for the fracture network from the determined area of the plurality of geometric primitives.