A New Approach for the Efficient Construction of 3D Geological Models for Reservoir Applications

Thomas J. Lasseter, S2S Systems
15325 189th Ave NE, Woodinville, Washington 98072, U.S.A.

1. Introduction

Currently the principal problems in the construction of 3D models for reservoir applications are:

This paper describes the details of a new gridding method that is an extension of the POSC-Rescue multiblock approach. The details of the method and its advantages will be described and discussed.

2. The Model Description

Figure 1 shows a model built with a new method that we will refer to as the XY-orthogonal truncated grid method.

The model was built from two auto-picked horizons. Fault edges for approximately 30 faults were interactively picked on each horizon. The interpretation required one to two hours. Construction of a model that has approximately three million cells required two minutes of computational time on a computer with a Pentium III 1 GHz processor and 512 MB of RAM. Larger models require more memory, but we hope to substantially reduce memory requirements with future improvements.

The XY-orthogonal truncated grid model has two major groups of elements: 1) the boundary and fault sections that form the framework and 2) the surface and volume elements that define the blocks inside these framework elements.

The model consists of a set of sections that divide the model laterally into blocks. Horizons across the model divide the model vertically. The volume between a pair of horizons is a unit. It is assumed that horizons and units extend all the way across the model. Horizons are divided into blockHorizons in each block. A pair of blockHorizons forms a blockUnit as shown in Figure 2. Each block contains a blockUnit for each unit, though some blockUnits may be zero-thickness or null.

Units may be divided into subunits defined by interval thicknesses and layering style (proportional, top-truncated, base-truncated, top-and-base-truncated). BlockUnits are similarly divided into blockSubUnits.

Each blockUnit has a blockUnitGrid composed of IJK cells where IJ corresponds to the original XY-orthogonal grid. The IJ values range over a bounding box around the blockUnitGrid. K is the vertical index. Each cell in the blockUnitGrid is full, empty or truncated. If truncated by one or more bounding sections, the cell geometry is defined by a set of polygons that define the faces of the truncated cell.

Figure 3 shows sections that are defined by sectionLines and sectionEdges. SectionLines are picked along existing sections or are computed from sectionEdge end-points. SectionEdges are typically picked along horizons or fault gaps (if fault sections). The sectionGrid consists of a grid of rowLines and columnLines. XYZ or XYT coordinates are specified at each point in the sectionGrid. To greatly simplify and speed-up key parts of the construction and display processes, rowLines are at constant Z (depth) or T (time) values and sectionLines and columnLines are monotonically increasing in Z or T. In some cases, section intersections cannot be defined by sectionLines and must be trimmed by sectionTrimEdges where they intersect other sections as shown in Figure 3.

A block is defined by a set of blockSides that are the contiguous portions of sections as shown in Figure 4. A block may be open at the top and/or the bottom. A blockSide is defined by a set of blockSideTrimEdges on a section. A blockSideTrimEdge is all or part of a SectionLine, sectionEdge, or boundary edge of a Section.

For each block, we intersect each horizon against the blockSides to create a blockHorizon as is shown in Figure 2.

The complications in using this gridding method come in defining the polyhedral cells on the sections bounding each blockUnit.

In Figure 5, we show a view looking vertically down on two section grid cells with a unit on one side of the fault crossing a unit on the other side. Each section grid cell is divided by I and J vertical planes into a set of IJ-subpolygons. Each such IJ-subpolygon is divided by K horizonEdges or subhorizonEdges on each side of the fault into a set of IJK-subpolygons for each side.

Each of these final IJK-subpolygons is the face of a truncated IJK-polyhedral cell belonging to a block contiguous to the fault section.

In Figure 5, we see cases where we have two IJK-subpolygons with the same IJK belonging to different section grid cells. These are separate faces of the same polyhedral conforming to the exact geometry of the section.

Because each section grid cell contains sets of IJK-subpolygons on each side, the problem of computing the connected area is simply a matter of computing the overlap of polygons between the two sides. It is often the case as in Figure 5 that there is more than one polygon overlap between the same two blockUnitGridCells. For reservoir simulation these should be combined into a single transmissibility value between the two connected cells.

Many of the elements here are identical to those defined in the POSC-Rescue specification for a multiblock model (see www.posc.org/rescue). The key elements missing from the Rescue model are explicit definitions of the elements defining the polyhedral grid cells. It is left to the downstream users of the modeling information to generate these polyhedral cells if they wish to do so. The absence of the polyhedral cells in the geological model representation reduces our ability to directly link and iterate between the geological model and an upscaled simulation grid built from it.

3. Advantages of the XY-orthogonal Truncated Grid

In Figure 6, we compare the approach we are advocating in this paper, the XY-orthogonal truncated grid method with the “standard” gridding approach: the corner-point grid. Two different 2D models are displayed: 1) a cross-sectional grid with an antithetic fault on the left and 2) a horizontal grid with a set of interconnected faults on the right. The figure lists the more obvious advantages and disadvantages of each method in our estimation.

3.1 Ease of Definition. We need only define the fault framework and the geometry of the XY-orthogonal grid (XYorigin, rotation-angle, NX, NY, DX, and DY). The resulting truncated grid is fully defined by this information and is unique. In the corner-point grid case, there are many grids that are possible and significant user interaction is required to define the one desired. In Figure 6a, for example, it is not at all obvious how the corner-point grid on the right should be constructed given the fault framework. If the model is 3D and the faults are non-vertical, this model could not be accurately defined with a corner-point grid.

3.2 Congruence with seismic and property grids. Most 3D seismic data and property grids are uniformly spaced orthogonal grids in the XY plane. The XY-orthogonal truncated grid remains congruent with these grids so that seismic attributes and computed properties can be displayed on the model grid without resampling or rescaling. The corner-point grid generally requires interpolation of seismic and property grids.

3.3 Simplified operations. Surface-to-surface operations are much simpler because all surfaces use the same grid. Vertical slices along grid-lines are can be computed and displayed on the fly. Interpolation of top and bottom unit surfaces to generate intermediate surfaces is simple and much faster than the case with a corner-point grid.

3.4 Process Iteration. If we upscale by combining rows and columns, the simulation grid maintains a simple relationship to the original XY-orthogonal grid. Simulation results can thus be easily mapped back onto the fine-grid and compared directly with seismic attributes and fine-scale properties We can make the process of iterating from seismic to simulation and back again finally practical.

Figure 6. Comparison of Gridding Methods

4. Model Construction Details and Issues

In the following sections will be discussed some of the key issues in the model construction procedure.

4.1. Gapping Horizon Grid Back from Faults. An important detail is how to efficiently intersect the horizon grids against the bounding sections. In Figure 7a we show points on a blockHorizonGrid for one block. We have chosen to gap back the grids a distance of one grid cell distance for the associated faults. In general, we define a gapping distance on each side of each fault. Any point that is within the gapping distance of any fault is eliminated. In Figure 7a, the solid circle points are the remaining points and the open circle points are gapped points. Using the remaining points, we compute appropriately interpolated/extrapolated values for the gapped points. There are many choices for the algorithm to be used; the primary criteria are efficiency and smoothness. Because we want to intersect the blockSides against the blockHorizonGrid, we must extrapolate points across the anticipated intersection with the blockSide. These extrapolated points are shown as open square points.

Figure 7. Model Grid Constuction Details

4.2. Handling Dying Faults. An important problem we have to address is how to handle dying faults that do not necessarily form a closed block. This is accomplished by adding auxiliary sections that extend or intersect these dying faults and form a closed block. Figure 7a shows an example where an auxiliary section (shown as a dashed line) has been added. Because the block is “open” across this auxiliary section, additional points shown as solid square points are used in the interpolation. The weight associated with these points should be reduced if they are partially screened from the point being interpolated by a dying fault.

4.3. Defining Grid Cells Around Laterally and Vertically Dying Faults. At faults we have shown that cells are truncated. How do we handle grid cells at the ends of dying faults, across auxiliary sections, and vertically when faults do not extend all the way through the model?
A dying fault is first trimmed back to the last row or column of the blockHorizonGrid that it crosses. For each cell split by an auxiliary section, we reassemble the two pieces and assign it to the blockHorizon owning the larger piece.

For a fault that does not extend all the way through the model, the active area of the fault is defined by a trim-loop. Outside this trim-loop, there is no displacement and the fault does not effectively exist. Where the fault splits a cell and part of or the entire cell is outside the active area trim-loop, we reassemble the two pieces with no displacement and assign it to the blockUnit owning the larger piece.

5. Interfacing to Property Modeling

Because the 3D blockUnitGrids need to be populated with properties, we need to be able to efficiently interface with property models and property modeling packages. Because most property models operate on an XY-orthogonal grid, the simplest approach is to have the property modeling package construct a flattened unfaulted grid with the same XY-grid as is being used here. Because there is then a one-to-one cell mapping between the faulted model and the property model, data defined on the faulted model such as log-curve averages, seismic attributes, etc. can be output easily to the property-modeling package. The resulting 3D property model can then be mapped back to the faulted model. Because faults may divide a single cell into multiple polyhedral pieces, each will have the same property value. If this is unacceptable, than a more elaborate grid-to-grid transform would be required, the details of which go beyond the space limits of this paper.

6. Output For Upscaling and Reservoir Simulation

The reservoir model is organized into a set of blocks. Each block contains a set of blockUnits. Each blockUnit has a 3D blockUnitGrid. These blockUnitGrids share the same global XY grid. Each cell in each blockUnitGrid is either: complete (a hexahedral), null, or a polyhedral.

In the terminology of reservoir grid description, we now have a multiblock grid description: each cell can be indexed by block number and IJK indices within the block. Where the simulation/upscaling package does not understand multiblock grid descriptions, these indices must be mapped to a single domain grid index I’J’K to reindex cells from different blocks that have the same IJK. Note that the K index remains the same; only IJ needs to be reindexed. Visually we can think of this reindexing process as pulling the blocks apart in the I and J directions enough to eliminate any overlap in IJ. While some null cells will be generated in this process, the number should be minimal. This approach in fact has fewer nulls than the multiblock approach where many nulls may exist in the bounding box around a given block.

Because simulation/upscaling packages do not generally allow for the input of polyhedral cells, the volumes of these cells must be computed and provided as “cell volume overrides”. Non-neighbor connections across faults or vertically across truncated cells must also be provided. A significant remaining problem is how to upscale non-neighbor connections. The advantage of this method is that the cell is orthogonal in XY which means that X and Y permeability values refer to absolute coordinates as opposed to local coordinates which is the case in corner-point gridding. Transmissibility calculations can thus be computed more accurately in this method. It is preferable however to provide the polyhedral geometries to the upscaling/simulation package so that the transmissibilities can be recomputed there in case the user wishes to change permeability values. This will require tighter integration between the modeling and upscaling packages than is currently available.

7. Future Developments

Specific areas we are working on or plan to work on include:

8. Conclusions

This paper describes and discusses the advantages of the XY-orthogonal truncated grid method. Performance and functionality in building large complex models has been clearly demonstrated. Because the method extends the POSC-Rescue format, its general usefulness has already been well established in commercial applications. The extensions though allow it to be more broadly used in iterating the complete seismic-to-simulation process. While more work needs to be done, the method offers significant benefits in terms of simplicity, functionality, and performance for 3D reservoir modeling.