Second Quarterly Report
June 7, 1996 to August 31, 1996
A Project of:
Fundamental Geoscience
Research and Development
BDM-Oklahoma
U.S. Department of Energy
National Oil and Related Programs
Contract Number #G4S51728
Prepared by:
William S. Dershowitz
Paul R. LaPointe
Herbert H. Einstein
Violetta S. Ivanova
Golder Associates Inc.
Redmond, Washington
October 24, 1996
963-1357.521
| Table 2-1 | Fracture Sets Classification in Terms of Fracture Orientation |
| Table 2-2 | Stratigraphic Record of the Yates Field |
| Table 2-3 | Tract 17 and 49 Area Windows |
| Table 2-4 | Yates Field Unit Window Area Coordinates |
| Table 2-5 | Field Data on Fracture Properties in the Yates Field |
| Table 2-6 | Block Size Algorithm Test Case Theoretical Values |
| Table 2-7 | MDS and CH Block Size Algorithms |
| Table 2-8 | Drainage Volume Algorithm, Test Case 1 |
| Table 2-9 | Expected Errors for Drainage Volume Algorithm, Test Case 1 |
| Table 2-10 | Drainage Volume Algorithm, Test Case 2 |
| Table 2-11 | Fracture Property Classes |
| Table 2-12 | ISIS/N User Interface |
| Table 2-13 | ISIS Data Format (.ISI) |
| Table 2-14 | Borehole Data Format (.SAM) |
| Table 2-15 | ISIS/N Verification Case |
| Table 2-16 | Calculation of Volume and Horizontal Extent for Hypothetical Model |
| Table 3-1 | Schedule Revisions |
| Table 3-2 | Milestones and Deliverables |
This quarterly progress report describes activities during the period June 7, 1996 to August 31, 1996. Work was carried out on the following tasks:
The major activity in the quarter was intensive development of the neural network code for fracture set identification, and field data collection in support of development of the hierarchical fracture model. Major efforts were also made for development of technologies for reservoir compartmentalization analyses, and for updating of the WWW site.
The following project deliverables were scheduled or submitted during this quarter.
Deliverable |
Scheduled Date |
Date Submitted |
Data Available on WWW |
96.06.30 |
96.06.15 |
First Quarter Progress Report |
96.06.15 |
96.06.26 |
The following issues were addressed during the quarter:
The following tasks were active during the quarter:
This section describes progress during the quarter for each of the active tasks.
During the quarter, Marathon Oil Company collected fracture and production data from project study site and provided this data to Golder Associates to provide the basis for the initial data warehouse. The following additional data was donated to the project:
This data was assembled and processed into an information warehouse within the World-Wide-Web (WWW) server application (Task 5.1.2). Figures 2-1 through 2-5 present examples of the data provided by Marathon during the quarter.
2.2.2 Task 1.2.1: 3D Hierarchical Fracture Model
During the quarter, MIT continued its research on the mathematical description of rock fracture systems based on geologic fracture genesis processes. The general 3D hierarchical model was further developed, so that it will soon be suitable for verification with field data from the Yates Field. A visit to the Midland office of Marathon Oil and the Yates Field at Iraan, TX, made it possible to collect data on the geology and fracturing of the region. This geologic information will be used in the modeling of the fracture systems in the Yates reservoir, based on the 3D hierarchical model. The data gathered included both surface outcrop observations and measurements of fracturing, and subsurface fracture data from core and well logs.
2.2.2.1 Classification of Fracture Systems Formed in Different Geologic Processes
A fracture set is defined as an assembly of fractures that have been created by common geological mechanisms, and share common geological characteristics, such as orientation, mineral infillings and surface morphology. Fracture systems (or networks) consist of one or more fracture sets. Structural geologic mechanisms, such as folding, faulting and crustal extension, cause development of a multitude of diverse fracture sets in complex relationships to the geologic processes and to one another. (Figure 2-6) illustrates some fracture systems, typical for several geologic mechanisms.
2.2.2.1.1 Fracture Sets Associated with Folding
Folds are geologic structures that form due to bending, buckling or plastic flow of rocks. Buckling can be caused by a compressive force, acting parallel to the surface of the rock. Bending may result from solution of underlying beds or from differential compaction of sediments. The type of folding and the associated fracture sets, characteristic for a given reservoir, depend on the rheological and mechanical properties of individual strata, the bed thickness, and the depositional order and variety of strata composing the reservoir. For example, several massive dolomite layers might deform as a single unit if there are no appreciably thick shale units separating them. If relatively incompetent layers, such as a shale, separate more competent beds, then each competent layer may behave more independently, and shear stresses that would otherwise produce fracturing are dissipated in slip in the shales.
2.2.2.1.2 Fracture Sets Due to Folding and Faulting
Several models have been proposed for how fractures form when rock is folded. Within the petroleum industry, models such as those proposed by Stearns (1971), Stearns and Friedman (1972), Suppe (1983) and Mitra (1993) are widely used. Physical model and seismic studies by Dula (1991) and Withjack et al. (1995) have been used to develop models to predict fault and joint orientation due to crustal-scale tectonic processess such as passive margin rifting. While there are diverse causes of fracturing and jointing, the patterns that result share an important mathematical aspect; they are hierarchical.
Figure 2-6 shows fracture patterns developed in many reservoirs. Fractures produced during concentric folding include bedding-plane faults, tension joints and conjugate pairs of shear fractures on the fold crest and within the flanks (Figure 2-6a). Such fractures are commonly found in foreland areas of thrust systems.
Figure 2-6b illustrates the fracturing style found in many strike-slip systems. Strike-slip faulting is caused by lateral compression. Sometimes strike-slip re-activates pre-existing fractures (joints or older faults). Simple or compound zones of overlap between en echelon strike-slip faults usually contain a multitude of secondary fractures: splay cracks, steps, bends, and other oblique fractures (Figure 2-6b).
During horizontal extension, steeply dipping or vertical sets of numerous subparallel joints form in competent rocks of all lithologies and tectonic settings. When some joints propagate, they inhibit the growth of other joints in their vicinity. Joints in sedimentary and other layered rocks usually occur in sets of parallel fractures, oriented perpendicular to the bedding. The spacing between joints of a given set is relatively constant within a single layer and is a function of the rock type, bedding thickness and mechanical contrast with the rock above and below the bed. Joints are usually more closely spaced in thinner layers. In the case of interbedded rocks of different lithologies, bed-perpendicular jointing may occur in only one type of rock. For example, Figure 2-6c shows joints in competent siltstone, propagating across some thin lenses of incompetent shale, but terminating at the thick shale beds.
A slight rotation of the orientation of the remote tensional stress may cause the propagation of en echelon cracks from the ends of the parent joints. The echelon cracks propagate along twisting surfaces into the new principal stress plane. Echelon cracks may grow to lengths that are relatively great compared to their width, which generally is only a small portion of the length of the parent joint.
A fourth class of fracturing mechanisms involve vertical geological processes, such as diapirism and force-folding. Typical structures, which are formed by a vertical active force and have a circular shape in plan view, are diapiric and collapse structures.
Salt domes are common diapiric structures, which form in result of accumulation and rising of salt underlying sedimentary rocks. As the salt dome grows, stretching of the surrounding and overlying sedimentary strata causes the formation of a complex network of radial and concentric normal and thrust faults, dipping downwards toward the salt stack (Figure 2-6d).
Collapse structures are common in areas containing soluble rock. In this case collapse occurs in empty spaces created by solution of limestone, gypsum, and salt or by the differential compaction of shale masses. Typical fractures are normal faults, dipping towards the center of the collapse.
2.2.2.2 Summary of Fracture System Geometry
Despite the great variation of geologic processes and conditions, natural fracture systems are composed of fracture sets that form a hierarchical pattern. This underlying similarity makes it possible to develop a method for creating three-dimensional fracture models based upon hierarchical processes. Incorporation of these hierarchical processes is an important advance in creating more useful fractured reservoir models.
2.2.2.2.1 Fracture Sets Classified According to their Hierarchical Domains
A fracture set generally develops within a portion of a geologic formation, roughly enclosed between several planar or non-planar surfaces. For every fracture set it is important to identify where exactly it forms in relation to other geologic structures. Two major types of fracture sets can be distinguished: (1) primary sets which are created either in intact rock, or without relation to other fracture sets in previously fractured rock (for example, subparallel joint sets in a rock formation subjected to remote tension); (2) secondary sets which develop in close relation to primary sets (for example, oblique fractures within simple and compound fault zones). The region of existence of primary sets is confined by the boundaries of the geologic structure, within which a particular stress field exists. The region of existence of a secondary fracture set is defined by one or several primary fracture sets and the local stress field to which the primary set gives rise.
2.2.2.2.2 Types of Fracture Sets According to the Orientation
Fracture set orientations are summarized in Table 2-1. Most fracture sets consist of subparallel fractures. Some fracture sets are made up of non-parallel fractures, usually with orientations symmetric around a plane or an axis.
Table 2-1 Fracture Sets Classification in Terms of Fracture Orientation
| I. FRACTURE SET ORIENTATION RELATED TO A PLANE OR SURFACE | ||
| Class A: | Description: | The fractures within the set are subparallel to a plane or surface. |
| Example: | Bedding-plane fracturesin concentric folds, parallel to the sedimentary beds. | |
| Class B: | Description: | The subparallel fractures of the set are oblique (or perpendicular) to a plane or surface. |
| Example: | Antithetic faults, oblique to normal faults. | |
| Class C: | Description: | The orientations of subparallel fractures within a set are symmetric around a plane. |
| Example: | Longitudinal tensile joints on an anticlinal crest, symmetric around the fold axial plane. | |
| Class D: | Description: | Sets of conjugate pairs, symmetric to a plane or surface. |
| Example: | Conjugate shear fractures on an anticlinal crest, symmetric around the fold axial plane. | |
| Class E: | Description: | Special case: sets of non-parallel connected fractures, perpendicular to a plane or surface. |
| Example: | Syneresis or dessication, perpendicular to the topographic and contact surfaces. | |
| II. FRACTURE SET ORIENTATION RELATED TO AN AXIS | ||
| Class F: | Description: | Fracture sets, striking radially form an axis. |
| Example: | Normal faults and some thrust faults next to and above a salt dome, striking radially from the salt dome axis. | |
| Class G: | Description: | Fracture sets, concentric around an axis. |
| Example: | Cone sheets, concentric around the magma chamber axis. | |
2.2.2.2.3 Types of Fracture Sets According to Fracture Intensity
The intensity of a fracture set is a measure of the ratio of fractured and intact rock. It is convenient to distinguish three types of fracture intensity: (a) low; (b) high, due to a small number of large fractures; and (c) high, due to a large number of small fractures. The fracture intensity of a set is low only if the set consists of a small number of small fractures. An example is a set of small cracks in the rock between outcrop-scale fractures in a tensile field. Rocks may be highly fractured either through a system of large fractures, or though a very large number of small fractures. An example of the former is a set of bedding-plane shear fractures in a concentric fold.
2.2.2.3 Geologic Stochastic Model of 3D Fracture Networks
During the quarter, the 3D Hierarchical Model (Ivanova 1995 and Ivanova et al. 1995), model was enhanced with some new procedures, the most important being the Fractal Line Tessellation. This is expected to be particularly useful for modeling the San Andres formation and similar fractured reservoirs because it reproduces the clustering of fractures into zones.
2.2.2.3.1 Basic Concepts of the 3D Hierarchical Model
The 3D Hierarchical model is characterized by the following basic features. The model is:
In the current version of the model, the fractures are convex polygonal planar objects of discontinuous rock, randomly located in three-dimensional space. The pole orientation and the coordinates of the center and the vertices of an individual fracture are determined indirectly when this fracture is generated as a member of a fracture set.
2.2.2.3.2 Generation of a Single Fracture Set with the 3D Hierarchical Model
Geometrically, a fracture set is a collection of fractures with related orientations, sizes, and locations. A fracture set is characterized by the following parameters: (1) mean orientation; (2) probability density function (PDF) describing the distribution of fracture orientations around the mean orientation; (3) size variation of the fractures which belong to that particular set; (4) density of fracturing as it varies in space.
A fracture set is generated by applying a sequence of three stochastic processes in a modeling volume, representing the rock mass in which the fracture network developed. The modeling volumes of independent sets are defined by the boundaries of geologic structures in which the particular primary fractures are created. Dependent sets are usually generated within modeling volumes, determined by the randomly located independent fracture sets. Figure 2-7 illustrates the modeling volumes for secondary and tertiary fractures, enclosed between the primary and secondary fractures, respectively, in strike-slip fault zones.
The primary process consists of a homogeneous, isotropic or anisotropic, Poisson plane network, representing the principal tensile and shear planes defined by the existing stress field. The variation of plane orientations is described with PDF's. Figure 2-8 shows some examples of how the stress field determines the fracture orientations, and the choice of appropriate PDF in these cases.
The secondary process, subdivision of each plane into a fractured region and its complementary region of intact rock, is accomplished by a homogeneous line network and a non-homogeneous polygon marking procedure. The secondary process produces sets of fractures which have certain size and shape variation, and are arranged in clusters. The secondary stochastic process represents the intensity of fracturing, which can be defined as the ratio of fractured to intact rock.
The subdivision of every plane is obtained with the Fractal Line Tessellation (Figure 2-9). First, the original plane is divided into large polygons by a small number of Poisson lines and some of the polygons are marked as fracture zones. The probability of marking a polygon as fractured, P, is calculated for every polygon as:
P[ a polygon is fractured] = P[location, size, shape, termination, ...]
Next, a Poisson line tessellation of small density and a marking process is performed only in the fracture zones. This procedure is repeated in a number of iterations. The parameters of the Fractal Line Tessellation (number of iterations, density of Poisson lines in every iteration, and the probabilityof marking polygons as fracture zones) can be chosen to represent low or high intensity of fracturing, as well as different fracture persistence. The Fractal Line Tessellation leads to a relatively fast subdivision of a plane into the desired fractured and intact areas.
The tertiary stochastic process involves random shifting and rotation of (a portion of) the fractured polygons in the vicinity of their original positions. The tertiary stochastic process represents possible variations of the stress field, for example, rotation of the principal stress axes, which leads to some specific fracture patterns such as fractures arranged en echelon.
2.2.2.3.3 Generation of a Fracture System with the 3D Hierarchical Model
A natural rock fracture system is reproduced in three dimensions through a superposition of independent and dependent fracture sets. Dependence on previously generated sets can be obtained, for example, (1) by specifying the marking probability for a polygon as a function of the distance to the primary fractures, and (2) by defining a function for termination of fractured polygons from a secondary set at the intersections with fractures from a primary set.
2.2.2.3.4 Applicability of the 3D Hierarchical Model for Representation of Geologic Fracture Networks
The 3D Hierarchical Model has been chosen for simulation of the fracture system in the San Andres formation because of its capability to reproduce important fracture system characteristics such as hierarchical fracture genesis and clustering. The model combines the advantages of the purely geometric and the purely mechanical models.
The conceptual model duplicates natural fracture systems on the basis of the available geologic information. The implemented stochastic processes are inherently related to the underlying mechanics and geology, which makes the model promising as a tool for simulation of geologic fracture systems.
2.2.2.4 Geology of the Yates Field
Significant efforts were made by the project team during the quarter to develop a quantitative understanding of the geology of the Yates field study site as it relates to the implementation of the 3-D Hierarchical Fracture Model. This section summarizes the geological setting of the Yates Field, and some preliminary observations on the fracture geology.
2.2.2.4.1 Regional Geology of the Permian Basin
A generalized column of the stratigraphy and lithology of the Permian Basin of west Texas and southwest New Mexico region is shown in Figure 2-10. In earlier studies of the regional geologic setting of the Yates field, performed by Marathon Oil geologists (Craig 1958), the geology of the Permian Basin is discussed in terms of five major phases.
2.2.2.4.1.1 First Phase: Upper Cambrian and Early Ordovician
More than 3,000 feet of sediments were deposited during the Upper Cambrian and Lower Ordovician. The Upper Cambrian sea advanced from the south and southeast and deposited conglomerates, coarse sandstones, and sandy carbonates on a broad, low relief shelf of Precambrian rocks. The Lower Ordovician sediments (Ellenburger Group) are mostly evenly bedded lithographic to coarsely crystalline limestones and dolomites, slightly cherty and sandy at some locations. The upper surface of the Ellenburger group has undergone several episodes of erosion.
2.2.2.4.1.2 Second Phase: Middle Ordovician through Mississippian
Relatively uniform deposition of marine sediments occurred within the limits of the broad Tobosa Basin (Figure 2-11). The limits of the basin are defined on the east and west by two broad shelf areas (Concho arch and Diablo Platform), and on the north by a landmass of Precambrian rocks (Pedernal Massif). The Middle Ordovician (Simpson) Group consists mostly of fine sandstones and green shales with minor limestone, totaling more than 2,000 feet near the center of the Tobosa basin. The rocks of the Montoya Formation (Upper Ordovician) are mostly cherty carbonates with some green Sylvan Shale. Silurian and Devonian rocks consist of more than 1,500 feet of limestone and dolomite, with some green shale and chert. A black, carbonaceous shale unit (Woodford), with maximum thickness of 600 feet near the center of the basin and gradual increase in silt and sand content near the rim, comprises the boundary between Devonian and Mississippian rocks. Mississippian rocks, composed mostly of limestone, interbedded and in some places overlain by shale, are up to 1,200 feet thick. The sediments from this phase were greatly reduced in thickness and geographic extent by erosion at three regional unconformities.
2.2.2.4.1.3 Third Phase: Pennsylvanian
The first period of regional folding and faulting to affect the region took place in the Pennsylvanian. Figure 2-12 depicts the framework of basins and shelves which developed during the Pennsylvanian. Due to substantial upwarping of the ancestral Central Basin Platform, the regional sedimentary basin was subdivided into the Midland and the Delaware basins (the three elements together constitute the Permian Basin). This framework was bounded by contiguous shelves on the north, east and west, and a trough on the south. Early Pennsylvanian (Bend) sediments consist mostly of coarse sands and shales. Along the east flank of the Central Basin Platform and on other sea floor highs, some dark, shaly limestones were also deposited. The basin was filled with 10,000 feet of clay, sand, and fine sand detrital material from the Llanoria landmass. These thick Bend sediments were folded and thrust northward during the Late Pennsylvanian. Post-Bend Pennsylvanian (Strawn, Canyon, and Cisco) rocks are mostly shales, siltstones, and limestones, and some sandstones, with a maximum total thickness of 2,000 feet. For the first time in the Paleozoic, limestone ridges were formed during the Post-Bend Pennsylvanian at the shelf margins along the eastern flank of the Midland basin, and along the Central Basin and Diablo platforms. These limestone deposits were formed as carbonate buildups in the shallow water down shelfward from the zone of surf and do not contain evaporitic material. During the Post-Bend, shales and thin, dark, shaly limestones were deposited in the basins. Due to their rapid subsidence throughout the Pennsylvanian, the basins remained only partially filled with sediments. A regional unconformity affects the entire area, except the deepest basins, and separates the Pennsylvanian rocks from the Permian.
2.2.2.4.1.4 Fourth Phase: Permian
Deposition of sediments into the basins developed during the Pennsylvanian continued in the Permian (Figure 2-12). Sedimentation was influenced by the gradual drying up and constriction of the seas in the Late Permian. Permian rocks total more than 20,000 feet in thickness and can be grouped into four series: Wolfcamp, Leonard, Guadalupe, and Ochoa.
Wolfcamp
Typical Wolfcamp sediments in the region are sands and clays, shed from the Llanoria Landmass. These sediments spilled into the basins, thinning northward and changing facies to dark, basin limestones and fine clastics. An uplift of the Eastern shelf caused the deposition of some carbonate buildups in its northern half. By the end of the Wolfcamp, the Permian Basin has become structurally stable. The region as a whole gradually sank, thus accommodating the Permian seas and their sedimentation. An erosional unconformity, separating the Wolfcamp from the Leonard, is present in the Central Basin Platform and the shelves, and absent in the two basins.
Leonard
Progressive constriction of the sea in the Midland Basin area continued during the Leonard. Facies typical of barrier reef environments are present in Leonard rocks. Reefs around the basin rims consist of irregularly bedded, fossiliferous limestone and dolomite, grading basinward into silliceous shale, clay shale, sandstone and thin-bedded, dark, shaly and cherty limestone, reaching a total thickness of more than 4,000 feet in the middle of the Delaware Basin. The rocks in the back-reef are composed of green, red and gray shales, thin-bedded limestones, and dolomites, with some anhydrite. An unconformity on the Central Basin Platform and the shelf areas separates the Leonard from the overlying Guadalupe rocks.
Guadalupe
The Guadalupe Series of the Permian Basin can be divided into two units. The lateral relationships between the lithofacies of the Lower and the Upper Guadalupe are shown in Figures 2-13 and 2-14, respectively. The vertical relationships among the major facies are shown in Figure 2-15. An erosional unconformity locally separates the Lower and Upper Guadalupe.
The Lower Guadalupe consists of three major facies (Figure 2-13): (1) the Cherry Canyon Sandstone of the Delaware Basin: a 2,000 to 2,500 feet thick sequence of dark, fine-grained basin sandstones and siltstones; (2) the Goat Seep Reef of the Delaware Basin margins: up to 1,200 feet thick limestone and dolomite reef and carbonate buildup; and (3) the San Andres of the Central Basin Platform: approximately 1,000 feet thick lagoonal or back-reef dolomites, grading landward into progressively more anhydritic, shaly, and sandy carbonates. Shrinkage of the Midland Basin continued during the Lower Guadalupe, but its rate was exceeded by the rate of sediment influx. The sediments in the Midland Basin, analogous to the basin, reef and back-reef types of the Delaware Basin, exhibit less sharply defined facies boundaries.
The Upper Guadalupe also exhibits the three major facies of the Permian: reef, fore-reef, and back-reef (Figure 2-14). Uplift of the southern shelf of the Midland basin shallowed the Shefield Channel and initiated the growth of the Capitan reef, consisting of up to 1,200 feet thick, irregularly bedded, fossiliferous limestone and dolomite. Fore-reef and basin rocks (the Bell Canyon Formation of the Delaware Basin) are dark siltstones and fine-grained sandstones with a few thin, dark basin limestones. The carbonate reef growth was confined to the periphery of the Delaware Basin, and the Midland Basin became part of the back-reef environment. Its subsidence controlled the thickness of the evaporitic, clastic and carbonate sediments which accumulated in the lagoonal and back-reef area into a sequence of five formations (the Artesia Group) with distinctive lithologies: (1) the oldest, Grayburg: up to 300 feet thick, sandy or silty and locally oolitic dolomite, with some bentonite beds and frosted quartz grains, and small amounts of anhydrite; (2) Queen, which filled structural and erosional lows on the Grayburg: a 300 feet thick sequence of interbedded red and gray sandstones and sandy dolomites, with abundant round, frosted quartz grains and some anhydrite and salt; (3) the thickest Seven Rivers: about 1,400 feet thick in the Midland Basin, thinning to 400 feet on the Central Basin Platform, mostly calcium sulfate minerals (gypsum and anhydrite), locally present salt and polyhalite, with common thin interbeds of siltstone and some basal silty dolomite; (4) Yates: up to 350 feet thick gray and red, fine quartz sandstone and siltstone,gray and red shale, a few thin beds of dolomite, and locally present anhydrite and salt; (5) the youngest, Tansill: 60 to 320 feet thick dolomite, changing to anhydrite and salt shoreward from the reef.
Ochoa
The Ochoa Series is underlain by a disconformity and a bed of maroon shale, up to 40 feet thick, which may represent an old soil. At the end of the Guadalupe, the barrier Capitan Reef stopped growing and all of the basin area became back-reef to reef or carbonate buildup, extending across the Hovey Channel and partially obstructing the flow of sea water into the Delaware Basin. The Ochoa Series consists of four formations: (1) the oldest Castille: evaporitic in nature, up to 2,100 feet thick in the Delaware basin, mostly “banded anhydrite” (alternating laminae of anhydrite and calcite), with minor thin-bedded salt; (2) Salado: evaporitic, up to 2,200 feet thick in the Delaware Basin, composed primarily of salt with small amounts of anhydrite, and locally important dolomite, potash materials and magnesite; (3) Rustler: 375 feet thick sequence of evaporitic and clastic rocks, deposited in an environment of fluctuating sea level, composed of basal thin-bedded sandstone, conglomerate, and variegated shale, overlain by upper and lower dolomite separated by a section of anhydrite and salt; (4) the final Permian deposit, the Dewey Lake: a marine sequence of orange-red siltstones and very finely-grained sandstones, partly cemented by anhydrite and gypsum, 281 feet thick in the Midland Basin. Between the Salado and the Rustler, there was a period of emergence and mild folding, which led to the deposition of Rustler sediments on upturned and eroded surface of Salado evaporites.
2.2.2.4.1.5 Fifth Phase: Post-Permian
After the withdrawal of Permian seas, a broad, flat plain near the sea level was left in the region. The plain remained stable and little changed, except for some subarial erosion during the Early and Middle Triassic, until the Late Triassic deposition of stream and lake sediments. More than 2,000 feet of variegated clay and silt, with lenticular bodies of sand and gravel, were deposited on a broad alluvial plain in warm and fairly moist climate during the Late Triassic. Most of the rocks in the Late Triassic sequence, conglomerates, sandstones and red shales, were derived from eroded Permian and Paleozoic rocks in the highlands. Some Triassic rocks were removed during the prolonged period of erosion during the Jurassic.
The stratigraphy of the Cretaceous in the Permian basin is very similar to that of the Upper Cambrian, due to the slow advance of the seas from the southeast over the flat plain. Initial deposition of coarse basal clastics (sands and gravels, Trinity group) and some finer clastics was followed by deposition of limestones (Washita-Fredericksburg group), which was predominant during the Cretaceous. The Cretaceous rocks, varying in total thickness from 800 to 1000 feet, were deposited on a shelf dipping gently to the southeast. At the end of the Cretaceous and during the Early Tertiary, the area west of the western edge of the Delaware Basin was strongly uplifted, folded, and faulted to produce a system of small mountains and detritus-filled valleys. Pliocene shales, sandstones, and conglomerates were laid by streams on a broad alluvial plain eastward from the mountains. Pleistocene deposits are mostly sand and gravel, formed as alluvial fans. Pliocene and Pleistocene deposits together total up to 5,000 feet in thickness. Increased rainfall, combined with the structural fracturing of the beds overlying the Salado salt, caused solution of the salt along the flanks of the Central Basin Platform.
Table 2-2 Stratigraphic Record of the Yates Field
2.2.2.4.2 Geology of the Yates field
The stratigraphic record in the Yates Field area, with an indication of rock types and thicknesses, is summarized in Table 2-2.
The total thickness of the Upper Cambrian and Lower Ordovician conglomerates, sandstones and carbonates within the limits of the Yates field has been estimated to 1,100 ft, based on data from deep wells. Some of the Lower Ordovician carbonate beds were probably eroded in the field area before deposition of Middle Ordovician rocks.
Some of the rocks deposited during the second phase (Middle Ordovician through Mississippian) have been completely eroded from the field area which lies well within the Tobosa Basin. The present rocks, variegated shales and sandstones of the Middle Ordovician, and dolomites of the Devonian and Silurian, are separated by regional unconformities. The total thickness of all pre-Pennsylvanian (first and second phase) rocks under the Yates field is estimated to be about 2,700 feet.
The total thickness of Pennsylvanian rocks under the Yates field is estimated to be from 900 to 1,250 feet, based on data from deep wells. The Yates field is located at the southernmost tip of the Central Basin Platform. In the field area, and over the southern portion of the Central Basin Platform in general, the upper layers consists of up to 320 feet thick sandstone and red and green shale interbedded with limestone. Beneath the sandstone-shale sequence, which may not persist over the entire field, the Pennsylvanian rocks consist of uniform fossiliferous, cherty limestone with some minor sandstone and gray shale, approximately 700 to 950 thick in the field area.
In the Yates Field, Permian-aged San Andres, Grayburg, and Queen Formations, totaling approximately 5,100 feet in thickness, constitute the main oil reservoir. The reservoir seal is provided by the Seven Rivers anhydrite. Wolfcamp and Leonard rocks, predominantly limestone with some shale, deposited along the Central Basin Platform, total between 2,500 and 3,500 feet in thickness. The Lower Guadalupe consists of back-reef deposits (San Andres), mostly dolomite. The San Andres is the top unit of 4,000-foot-thick platform carbonates, deposited on the eroded surface of pre-Pennsylvanian rocks. The top surface of the San Andres (plotted in Figure 2-16) is marked by a pronounced regional unconformity which caused the high solution porosity in the upper layers of the formation. In the vicinity of the field, the Grayburg dolomite is from 0 (absent by non-deposition at some locations) to 300 feet thick, averaging about 50 feet within the productive limits of the field. The Queen sandstones and the basal layers of the Seven Rivers anhydrite thickened into topographic depressions and thinned over topographic highs on the Grayburg surface. The Seven Rivers anhydrite, forming the seal of the main Reservoir, ranges from 310 to 950 feet in thickness, averaging about 400 feet within the productive limits of the field. The middle and upper Seven Rivers, the Yates sandstone, and the Tansill dolomite, have relatively constant thickness and lithologies over the field area. The Salado salt is absent over most of the Yates Field structure, except for a stratum of anhydrite 30 to 40 feet thick, surrounding the field anticline. The salt originally covered the field, but was later removed due to solution of the subsurface, preceding the advance of the Cretaceous seas, and due to physical erosion by stream action. The gentle folding that occurred in the interval between Salado and Rustler deposition produced the present Yates structure. The Salado layers are upturned and truncated beneath the Rustler formation, which is absent over most of the field area.
Triassic clastics were deposited over the entire area of the Yates Field, but were later eroded during the Jurassic and are missing over the apex of the field structure. Rocks, belonging to the three major groups of the Early Cretaceous (Trinity clastics, and Washita-Fredericksburg limestones) are also present over the field.
2.2.2.4.3 San Andres Formation in Tract 17 and 49 Areas
Since the Unit horizon extends down from the Seven Rivers anyhdrite to an elevation of 800 feet above sea level, there are no field data available from deeper strata. The elevations of the top surface of the San Andres (which is the formation of interest for fracture modeling) and the formation contacts above it have been determined in field surveys. Figure 2-17 depicts a structural map of the San Andres (top surface) in Tracts 17 and 49. The structural apex is in Tract 49, whereas a local peak is located in the southeast portion of Tract 17.
The stratigraphy of the tracts has been interpreted, based on measurements of porosity, gamma ray, and other log data. In Tract 17 the San Andres formation consists of massive dolomite, interbedded with some shale layers. Figure 2-18 shows the shale distribution, based on interpretation of gamma ray response in Tract 17 area. The lighter colors (higher radioactivity) most likely indicate shale beds. In Tract 49, the San Andres formation is composed almost entirely of massive dolomite. Figure 2-19 illustrates the variation of rock matrix porosity in the Tract 49 area, as it has been directly measured in core rock. The distribution of calcite content in core rock (mostly crystalized in fractures) in the San Andres formation in Tract 49 is shown in Figure 2-20. Most of horizon depth data, gamma ray response, matrix porosity, and other field measurements are available at the project web site. Figures 18 to 20 have been generated with a numerical model, called Stratamodel, a commercially available visualization system which is installed at the Marathon Oil Technology Center in Littleton, CO.
2.2.2.5 Fracture Data and Modeling in Permian Rocks of the Yates Field
The main objective of the pending fracture modeling is to realistically reproduce the fracture system, significant for oil accumulation and transport, in the top portion of the San Andres formation. The rock in the San Andres, down to a depth of 800 feet above sea level, contains the portion of the Main Reservoir, from which the Yates Field Unit currently obtains most of its oil production. According to Marathon Oil geologists, the San Andres was probably fractured due to differential subsidence of sediments in deeper strata, prior to deposition of the upper layers. This accounts for the less intense fracturing of the Grayburg and Queen formations, which overlie the San Andres and are also part of the Main Reservoir.
2.2.2.5.1 Selected Area for Fracture Modeling and Verification of the 3D Hierarchical Model
The area of Tracts 17 and 49 has been selected for verification of the 3D Hierarchical Model via simulation of the fracture network in the San Andres. The coordinates of Tract 17 and 49 within the local system of the Yates Field Unit are defined in Table 2-3. The structure of the San Andres has its apex and a local peak within Tracts 49 and 17, respectively. The modeling volume for fracture generation for either tract is defined by a horizontal bottom plane at elevation 800 feet above sea level, by four vertical planes coinciding with the boundaries of the tract area, and by a convex upward surface defined by the top of the San Andres.
Table 2-3 Tract 17 and 49 Area Windows
| Tract 17 Area Window | Tract 49 Area Window |
|
| Northern limit (Ymax) | 114,696 |
109,966 |
| Southern limit (Ymin) | 112,967 |
108,123 |
| Eastern limit (Xmax) | 110,141 |
114,762 |
| Western limit (Xmin) | 106,271 |
112,876 |
All values for coordinates are from the local survey coordinate system used in the Yates field area. The unit of measurement is “vara” (33.333 in). Conversion to Texas Central State Plan can be accomplished using the transformations provided in Table 2-4.
Table 2-4 Yates Field Unit Window Area Coordinates
| Scale factor: | 2.7771151 feet/vara | |
| Clockwise rotation: | 1.162667 degrees | |
| Point of rotation: | Yates (vara) | Texas Central |
| X=112,475.983 -> | X=1,498,103.600 | |
| Y=107,276.091 -> | Y=441,536.200 |
2.2.2.5.2 Fracture Size Distribution in the Yates Field Rocks
There are three important scales of fracture size in the fracture system of the Yates field. These include microfractures which help to drain the matrix, but do not play a major role in regional permeability; joints and on the order of several to several hundreds of feet; and major faults or fracture zones. The contribution of the small-scale fractures has been incorporated by Marathon geologists into a 1% increase of the rock matrix porosity.
The joints and smaller faults interconnect and form the 3D network into which oil can drain and be produced. The hierarchical model is used to reproduce the network of intermediate-scale fractures as a stochastic component of the reservoir.
In the current model of the Yates field, the large scale zones are modeled deterministically. This third significant fracture system is composed of fractures with sizes comparable to the entire Yates Field area. Their horizontal extent has been interpreted in a previous study Marathon Oil (Curran, 1996, personal comm.). Figure 2-21 shows a map of the large-scale fractures in the Tract 17 and 49 area. Little is known at present about the vertical extent of these field-scale fractures, and whether they are shear faults or extensional joints.
2.2.2.5.3 Field Reconnaissance on Fracturing of Permian Rocks in the Study Area
Figure 2-22 illustrates the locations in the Tract 17 and 49 area, in which field information on fracture properties has been collected. The available information is summarized in Table 2-5.
2.2.2.5.4 Relation to Fracture Outcrops in Cretaceous Rocks in the Yates Field
During a visit to the Yates field at Iraan, TX, several surface exposures of fracturing in Cretaceous rocks in outcrops and road cuts were observed by researchers from MIT and Golder Associates. Fracture and bed orientations (dip and strike) and spacing were measured at five locations near Tracts 17 and 49. The predominant outcropping sets are (Figure 2-25a and b): (1) bedding plane fractures, parallel to the beds which are approximately horizontal or slightly dipping; and (2) two or three sets of vertical or steeply dipping fractures, cutting across one or several beds. Some fractures are solution enhanced (Figure 2-25b), a few of them to openings as large as small caves. Solution-enhanced fractures tend to be sub-vertical and strike north-south. The moldic porosity in the top Cretaceous layers (Figure 2-25c) may be similar to that in the top portion of the San Andres. Most of the Cretaceous fractures are joints (i.e. opened in tension); at one locationin Tract 49 some faults associated with local folding were observed (Figure 2-25d). Some fractures are filled with precipitated calcite.
The fracturing in the Cretaceous rocks is interesting because these rocks were deposited in very similar environment to those in the San Andres, and were subjectd to many of the same geological processes. Because of the similar lithology of the rocks in the Cretaceous and in the San Andres, the fracture systems in the former may indicate what type of fracturing has occurred in the latter.
Table 2-5 Field Data on Fracture Properties in the Yates Field
| Locations | Description | Comments |
| YU1711 and YU1755 (within Tract 17); | Detailed information on orientation (dip and azimuth), depth, truncation, and mineralization of fractures, intersected by the well logs. | Excel files (available on the web page for well YU1711). Interpretation of various field measurements by Brendan Curran of Marathon Oil. |
| YU2511 (near Tract 17); | ||
| YU5127, YU4903, and | ||
| YU4007 (within Tract 49); | ||
| YU3728 and YU4959 (next to Tract 49); | ||
| YU2440 and YU2414 (next to Tract 17); | Interpretation of FMI/FMS logs in the form of stereoplots and other plots of fracture properties. | Example in Figure 2-23. The data need to be digitized prior to use in numerical modeling. |
| YU1711 (in Tract 17); | High resolution water injection profiles. Information on locatio of conductive fractures within the wells. | Example in Figure 2-24. Excel files for YU1711, Yu2511, Yu2414, and YU5127. Other data have to be digitized. |
| YU2414, YU2511, and YU15c1 (near Tract 17); | ||
| YU5127 and YU4903 (within Tract 49). |
Since rock fracturing is controlled by the geologic and material characteristics of the host rock, the project team members are studying the geological models of the Yates field and of tracts 17 and 49. For this purpose, a visit to the Marathon research center in Colorado is scheduled for the next quarter. The Stratamodel for the Yates Field is currently updated with all available field information. Analysis will focus on (1) the distribution of shale lenses in the San Andres dolomite, since shale may inhibit the propagation of fractures across dolomite beds; (2) the variation of the thickness of dolomite beds, since it may control to a great extent the spacing of both bed-perpendicular and bedding-plane fractures; (3) the variation of rock porosity and hence permeability, which will facilitate the determination of the percentage of conductive fractures; (4) the level of calcite precipitation, since calcite may have sealed some fracture sets at depth; and (5) the effect of the anhydrite overlying the main reservoir.
The most important objective for the next quarter is to establish relations between the geology of the Yates Field and the geometric characteristics of the fracture system. Since many fractures are steeply dipping, most of them are not intersected by the vertical wells. It will be necessary to relate those subvertical fracture sets to the curvature of the San Andres flexure. Bedding-parallel fractures can be modeled as subparallel to the San Andres flexure, if it can be established that the deformed dolomite beds preserved their thickness. Besides fractures related to the flexure, there might be some secondary fractures related to the field-scale discontinuities. The limited number of water injection profiles will be best applied for comparison to flow simulations through the synthetic networks, generated with the 3D Model.
2.2.3 Task 1.3.1 Implementation of Block Size Algorithms
During the quarter, a range of new algorithms were implemented and evaluated for estimation of rock block volume. The assessment of rock block volume is essential for efficient development and production of fractured reservoirs.
2.2.3.1 Alternative Block Size Algorithms
Reservoir simulation in dual porosity systems requires calculation of several parameters having to do with matrix block shape and size. The fracture surface area of matrix blocks within a simulation grid cell influences the rate and quantity of fluids that can move between the matrix and the fracture system. The Z-dimension of matrix blocks influences gravity drainage mechanisms. The shape of the matrix blocks influences the choice of sugar cube, matchstick or slab idealization. A realistic description of block size and shape in a way that can be implemented in existing dual porosity simulators will benefit not only the thermal simulation for TAGS processes, but also non-thermal simulation of injection or production in fractured reservoirs. Two algorithms have been developed and tested to compute matrix block shape and size that are :
The first algorithm is a fast computational method to compute blocks based upon fracture spacing distributions in several directions. Its main advantage is that the calculation is very fast. Its disadvantage is that it assumes that block dimensions are uncorrelated. This algorithm is referred to as the multidirectional spacing distribution algorithm.
The second algorithm is based upon graph theory and is based upon computing the convex hull of points lying on fractures bounding or partially bounding a matrix block. This algorithm is referred to as the convex hull algorithm.
2.2.3.1.1 Multi-Directional Spacing Distribution Algorithm
Figure 2-26 illustrates the multi-directional spacing (MDS) algorithm. For each realization of the discrete fracture model, a series of randomly-located lines in user-specified directions are generated. The location of fractures intersected by each line is recorded. This leads to a spacing frequency distribution in several directions. Typically, the directions include the vertical direction, in order to calculate the vertical dimension of blocks for gravity drainage considerations, and in two or three orthogonal directions that relate to simulator grid layering geometry and the fracture system.
The spacing probability distributions are multiplied together using Monte Carlo sampling techniques to produce a frequency distribution of block volumes and surface areas. This is carried out by selecting X, Y and Z spacing values at random with selection probability proportional to their frequency, and multiplying them together to create a prismatic block.
The volume of the block for the i-th Monte Carlo computation is the product of the spacing values:
The surface area of the i-th block is given by:
The block shape is given by the ratios of the mean X, Y and Z spacings:
The vertical dimension of the blocks is taken as the mean value of the Z spacing distribution:
2.2.3.1.2 Convex Hull Algorithm
The convex hull (CH) algorithm is more computationally intensive, but measures the actual dimensions of the blocks, rather than reconstructing blocks stochastically from spacing frequencies. Thus, any correlations among block dimensions or non-prismatic block shapes do not present problems. The algorithm's accuracy is governed by two factors: whether in fact the matrix blocks are convex; and how many points are required to accurately characterize the convex block. The algorithm as implemented allows the user to specify the number of points for characterizing the convex block, but computation time increases as the number of points increases.
Figure 2-27 illustrates the CH algorithm. From a random point within a discrete fracture model, lines are projected out in user-specified directions until they intersect a fracture. The intersection point coordinates for all of the line intersections constitutes a sample of the matrix block. As is shown in the 2D illustration (Figure 2-27a), some of these intersection points may not all be on the block exterior. Next, a convex hull (Figure 2-27b) is calculated for the set of points. A convex hull has the property that it includes all of the points, and that it is convex. In two dimensions, the hull is a polygon. This polygon has the additional property that any line joining two arbitrary points within or on the boundaries of the polygon lies entirely on the boundary or in the interior of the polygon. In practice, this means that the polygon has no "holes" or embayments, and it is a single polygon, not several isolated ones. In three dimensions, the polygon becomes a polyhedron.
It is a simple matter to compute the volume and the surface area of the convex hull. The vertical dimension of the matrix block is taken to be the difference between the minimum and maximum Z coordinate of the convex hull. The shape of the matrix blocks approximated as convex hulls can be determined by visually displaying the hulls, or by computing the aspect ratios for the X, Y and Z directions from the hull coordinates.
The convex blocks that are calculated are estimated by projecting the lines out from randomly selected points within the model. Since a larger block has a proportionately greater probability of having a randomly located point within it, the algorithm will overestimate the number of large blocks and underestimate the number of small blocks unless a check is made to see if an estimate has already been made for the block. The current algorithm checks to see if each new block has previously been estimated, and if so, it retains the block with the larger volume. This is because a convex hull estimate of a convex block is always equal to or less than the volume of the actual block.
2.2.3.2 Block Scale Algorithm Test Cases
Three test cases were devised to assess how well each of these two algorithms correctly estimated the matrix block size distribution. Cases 1 and 2, which are not geologically realistic, serve as verification and evaluation cases for the algorithms. It is not possible to calculate the true block volume distribution for Case 3, but this case illustrates how the algorithms might be applied in more realistic DFN models. Case 1 (Figure 2-28a) consists of a 200 m by 200 m by 200m cube dissected into 8000 10 m cubes. Each cube is a matrix block, and is completely bounded by fractures. Case 2 (Figure 2-28b) also consists of cubes, but of varying sizes and with dimensions that are partially correlated. The distribution of block volumes is bi-modal. These two test cases are geologically unrealistic, but clearly illustrate how well the algorithms perform. Case 3 (Figure 2-28c) is more geologically realistic. In this model, fractures do not completely isolate blocks. There are three sets in approximately three orthogonal orientations, much like the fracture pattern that develops in carbonates or other layered rocks.
Table 2-6 summarizes the values for Z dimension, specific surface area, P32 and volume for the three test cases. Both the specific surface area and P32 quantify how much fracture surface area per volume of rock there is in the model. The specific surface area is calculated by summing all of the surface area for all of the identified blocks, and dividing by the total volume. In test cases 1 and 2, this means that almost every fracture will be part of two blocks. P32 on the other hand, is the total amount of fracture area divided by the volume of rock, and is independent as to whether these fractures form blocks. For test cases 1 and 2, the specific fracture surface area should be approximately twice the value of P32. For simulation purposes, it is necessary to use the surface area that relates to the matrix block idealization, rather than the absolute P32 value.
Table 2-6 Block Size Algorithm Test Case Theoretical Values
| Test Case | Z Dimension (m) |
Surface Area/Rock Volume (m2/m3) | P32 (m2/m3) |
Volume (m3) |
| Case 1 | 10.0 |
0.6 |
0.312 |
1,000 |
| Case 2 | mean = 22.22 |
0.27 |
0.146 |
mean = 10,974 |
| Case 3 | not calculated |
not calculated |
0.112 |
not calculated |
Table 2-7 presents the results for MDS and CH analysis of Test Cases 1, 2, and 3. The MDS algorithm matches the correct block dimensions exactly for Case 1 when the sampling directions were chosen to be perpendicular to the block faces. In many situations involving fracturing of layered rocks, principal fracture orientations should be parallel and perpendicular to bedding, so it may be possible to choose sampling orientations that give very good matrix block volume estimates. The MDS algorithm produces reasonable estimates for Case 2, but has a tendency to somewhat overestimate the number of very small matrix blocks. Nevertheless, the MDS algorithm reproduces the shape and bimodality of the correct distribution reasonably well (Figure 2-29).
Table 2-7 MDS and CH Block Size Algorithms
| TEST CASE | Z Dimension (m) |
Surface Area/Rock Volume (m2/m3) |
P32 (m2/m3) |
Volume |
| Case 1 - MDS | 10.0 (constant) |
0.60 |
0.312 |
1,000 (constant) |
| Case 1 - CH | 14.1 (mean) |
0.63 |
0.312 |
823 (mean) |
| Case 2 - MDS | 23.4 (mean) |
0.25 |
0.146 |
6,653 (mean) |
| Case 2 - CH | 29.3 (mean) |
0.32 |
0.146 |
9,154 (mean) |
| Case 3 - MDS | 11.1 (mean) |
0.36 |
0.112 |
5,484 (mean) |
| Case 3 - CH | 43.0 (mean) |
0.36 |
0.112 |
10,642 (mean) |
The CH algorithm provides better results for Case 2 and worse results for Case 1. For Case 1, most convex hulls will contain a smaller volume than the cubes in which they are calculated. This is due to the fact that the convex hull is based upon a finite number of points lying on the edge surfaces of the cubes. This means that corners of the cube may be effectively "trimmed off" by calculating the hull, which in turns leads to an underestimate of block volumes. For Case 2, the median and mean block volumes are more accurately estimated by the CH algorithm than by the MDS algorithm. The inherent ability of the CH algorithm to account for the correlation in block dimensions appears to offset the tendency to underestimate the volume of cubical matrix blocks.
Case 3 is a more realistic DFN model. As horizontal and vertical cross-sections through one DFN model realization illustrate (Figure 2-30), blocks are not completely isolated on all faces by fractures. There are blocks of many different sizes, but in horizontal cross-section, the blocks are on the order of 40 m, and in vertical cross-section blocks are on the order of 10 m.
In summary, the three test cases suggest that both the MDS and CH algorithms provide reliable and consistent estimates of fracture surface area, at least for simple fracture geometries. The CH algorithm appears to provide better estimates of the mean volume of matrix blocks when block dimensions are partially correlated. Since jointing in many sedimentary rocks is characterized by pseudo-periodic spacings (e.g. La Pointe and Hudson, 1985), it may be preferable to use the CH algorithm to estimate block volumes. On the other hand, the geometric construction of a convex hull from a sparse data set creates hulls with slightly greater average Z-dimensions than the MDS algorithm. In both test cases 1 and 2 the MDS algorithm provided more accurate estimates of the Z-dimension. Thus, both algorithms have proven useful and necessary to provide estimates of matrix block parameters, and neither alone is completely satisfactory.
2.2.4 Task 1.3.2 Drainage Volume Analysis
Research on tributary drainage volume initiated during the first quarter of the project was continued during this quarter. Tributary drainage volume is related to both block size and compartmentalization. Tributary drainage volume is the estimated volume of matrix that a fracture system intersected by a well or perforated zone can access. In those reservoirs in which the matrix provides the storage, and the fracture system the reservoir permeability, the amount of hydrocarbon that can be produced is a function of the fracture network connectivity and geometry. In the TAGS process, there is no need to use elaborate calculations to determine the amount of matrix volume that has been effectively heated by steam injection. Studies performed by Marathon Oil (Kenyon, 1996, personal comm.) show that the volume of heated matrix can be directly calculated from the volume of injected steam. For the TAGS process, it is more important to estimate the shape of the condensation front through time, which is the reason for development of the steam compartmentalization algorithm described in the previous section.
Nonetheless, the mobilization and drainage of oil due to pressure depletion occurs in many fractured reservoirs. Improving forecasts of the potential volume of reservoir that a well can drain through the connected fracture system is useful for planning field development in many fractured reservoirs.
The algorithm developed to compute the tributary drainage volume is divided into two steps:
Step 1 utilizes the algorithm developed for compartmentalization analysis, as described in Section 2.2.6 below. As in the compartmentalization algorithms, it is possible to process the fractures in a DFN model to identify those fractures which are connected to a well or a portion of a well.
There are different ways that might be used to accomplish Step 2. If the fracture network were very dense, then the volume of the matrix accessed by the fracture network will be closely approximated by the volume of the convex hull enclosing the network. In essence, this means that the compartment volume is equal to the tributary volume. However, there is no reason to expect that all fracture networks will be sufficiently dense that the convex hull method produces the best estimate. In less dense fracture networks, some or maybe a significant portion of the volume inside the convex hull will be too far from any of the fractures to be easily produced.
A better estimate which does not include matrix than might not be efficiently produced through pressure depletion drainage is to surround each fracture in the network with a polygon that is calculated from the area of the fracture and the distance away from the fracture over which drainage might be effective (Figure 2-31). This leads to a prism that encloses the fracture. For pressure depletion mechanisms, the fracture forms the midplane of the prism.
An obvious problem with this algorithm is to avoid double-counting the volume where there is overlap between the prisms. Calculations based upon solid geometry to compute the volume while accounting for the overlap are highly time-consuming for the number of fractures that might commonly be encountered in a fracture network. A simpler method has been devised which is much more efficient, though not as numerically exact.
The prism enclosing each fracture is approximated by a series of cells that are referenced to a global origin. All of the cells whose centers lie within the prism are said to belong to the prism. Initially, all cells have the value of "0". Each fracture is processed in turn, and all cells belonging to the fracture prism are incremented by the value "1". If a cell belongs to more than one fracture prism, then its counter is incremented to a number equal to the number of prisms it belongs to. After all of the fractures are processed, the total number of non-zero cells multiplied by the cell volume is taken as the estimate of the tributary volume, or:
where
is the Dirac
function which is equal to 1.0 if the cell has a value greater than 0.0,
and is set to 0.0 if the value in the cell is 0.0, n is the number of cells,
and dx, dy, dz are the cell dimensions.
An example of this calculation for a simple, verifiable test case is shown in Figure 2-32. The model consists of three horizontal fractures with partial horizontal overlap. Expected an actual results for this test case are shown in Table 2-8.
Table 2-8 Drainage Volume Algorithm, Test Case 1
Case 1 |
Expected Volume (m3) |
Estimated Volume (m3) |
Percent Difference |
1.0 |
600 |
600 |
0.0% |
2.0 |
1000 |
1000 |
0.0% |
8.0 |
3400 |
3400 |
0.0% |
These series of verification cases are simplistic, but demonstrate that the algorithm correctly computes the volume and does not double-count when the tributary fracture volumes overlap. The reason for the exact match between expected and predicted volumes is that the fractures are horizontal, making the prisms exactly approximated by a series of stacked cells, and the cell size (1.0 m cubes) were evenly divisible into the prism thickness. For example, if a cell size of 3 m were used, then the algorithm would be expected to underestimate the volume. The predicted underestimates and the actual results are shown in Table 2-9.
Table 2-9 Expected Errors for Drainage Volume Algorithm, Test Case 1
Case 1 |
Expected Underestimated Volume (m3) |
Estimated Volume (m3) |
Percent Difference with True Volume |
1.0 |
567 |
567 |
-5.5% |
2.0 |
972 |
972 |
-2.8% |
8.0 |
3240 |
3240 |
-4.7% |
These cases verify that the algorithm functions as intended, and shows that a poorly chosen cell size can lead to a mis-estimation on the order of 5%.
Test Case 2 is defined for demonstration purposes, because there is no analytical solution for the case of incompletely constrained blocks in discrete fracture networks. Test Case 2 is defined using the fracture geometry for fractures connected to a vertical well drilled in the center of Block Size Test Case 3 (Figure 2-28). For this analysis we assume that the prism surrounding the fracture is of variable thickness, similar to the increasing volume of matrix produced through time as the reservoir is depleted. The measuring cells are cubes with edges of length 1.0 m. Results are shown in Table 2-10.
Table 2-10 Drainage Volume Algorithm, Test Case 2
Case 2 |
Volume from Drainage |
1.0 |
2908 |
2.0 |
5735 |
8.0 |
21840 |
Table 2-10 shows that the tributary volume does not scale linearly with increasing prism thickness, due to the overlap of the drainage regions for each fracture.
2.2.5 Task 2.1.1 Fracture Sets Analysis
During the quarter, we completed development of neural net technology for fracture set identification. The technology is based on the use of probabilistic neural networks to classify fractures. The algorithm shows the following advantages over conventional approaches:
This following section describe the algorithm implemented, the user interface, and a verification test case.
The basic concepts behind the use of neural networks for set definition were provided in our First Quarterly Report (Dershowitz, 1996), and will not be repeated here. The development below assumes familiarity presented in the First Quarterly Report.
The probabilistic neural network (Specht, 1990) is based on a combination of probability theory and Bayesian statistics, and was developed primarily for solving multivariate classification problems (Masters, 1993). The definition of fracture sets can be considered a classical classification problem, with the following special attributes:
The probabilistic neural network algorithm is designed to provide a classification which minimizes mis-classification of fractures to the wrong set. The classification system which has the minimum "cost" of mis-classification is termed "Bayes Optimal" (Parzan, 1962).
The algorithm for the implemented probabilistic neural network is illustrated in Figure 2-33. The algorithm proceeds as follows.
The user specifies weightings wi for each of the classification variables, for use in the utility function for evaluation of set classification.
where V(c) is the utility for classification c, Wi is the weighting for variable I, and Di(j|c) is the euclidian distance from the data point j for its classification c,
The classes of fracture properties which can be used in this algorithm are provided in Table 2-11.
Table 2-11 Fracture Property Classes
| Property Class | Description | Example |
| Real | Real valued number | Trace length, aperture |
| Integer | Integer valued number | JRC, RQD |
| Orientation | Trend on [0,360] and dip on [0,90] for the dip vector (D) or pole vector (P). For calculation of spherical angles, the minimum angle of either the upper or lower hemisphere orientation vector is used. The default is lower hemisphere | Fracture orientation, striation orientation, foliation orientation. |
| Vector | Similar to orientation, but uses only lower hemisphere values | |
| Class | Membership in a group, as a logical (0,1) value | Rock type, fracture termination mode |
| Ordinal | Positive, integer value | Fracture Set Number |
The user interface for ISIS/N is illustrated in Figure 2-34. The user interface is designed based on an "object oriented", tree structure model. This model is illustrated in Figure 2-35. Each fracture can have any number of properties, defined according to the class-types of Table 2-11. The procedure for set definition is described below. The menu options are listed in Table 2-12.
Table 2-12 ISIS/N User Interface
| Primary Menu Item/ Command | Secondary Menu Item/ Command | Action |
| File | New | Create a new analysis (the previous analysis will be kept open) |
| Open | Open an .ISI data file, and the .SAM file describing the boreholes and traceplanes the data was collected from | |
| Close | Close the current analysis | |
| Save | Save the fracture set definitions as .ISI, and .ORS files, and save the statistical reports as .STS files. Save using default file names | |
| Save As | Save the fracture set definitions as .ISI, and .ORS files, and save the statistical reports as .STS files. Save using user provided file names | |
| Print the contents of the current window, which can be either (a) data, (b) the analysis object tree, or (c) stereoplot visualizations | ||
| Print Preview | Display on-screen a preview of the items to be printed | |
| Print Setup | Setup the printer | |
| Exit | Leave ISIS/N | |
| Edit | Undo | Undo the previous text entry |
| Cut | Cut the selected text field | |
| Copy | Copy the selected text field | |
| Paste | Paste the copied text field at the selection | |
| Terzaghi | Carry out a Terzaghi correction on the selected data | |
| Define New Set | Define the properties, default values, and weights to be used for the next set in the current analysis | |
| Neural Net | Select Data Set | Select a subset of the currently open data file for analysis |
| Generate Training Set | Generate training sets based on the specified set statistics | |
| Load Training Set | Load a previously defined training set | |
| Train | Run the neural network on the training set to create a neural network | |
| Classify | Classify the fractures to sets using the neural network developed from the training set | |
| View | Toolbar | Display the toolbar on screen |
| Status Bar | Display the status bar on screen | |
| Stereoplots | Display a window containing a stereoplot of the current data | |
| Window | New Window | Create a new analysis (as with the New menu item) |
| Cascade | Cascade the windows | |
| Tile | Tile the windows | |
| Arrange Icons | Arrange icons neatly for the minimized windows | |
| <File Name> | Switch to the main analysis window for the analysis of the file name displayed | |
| Help | Help Topics | Provide context sensitive help |
| About ISIS | Provide information about ISIS. |
Table 2-13 ISIS Data Format (.ISI)
| # Any line beginning with a # symbol is a
comment. BEGIN FIELDS Count = 3 # Count is the number of fields. BEGIN FIELD Name = "Orientation" Type = ORIENTATION # data type orientation is expressed as two real values # theta (trend) and phi (plunge) > Pole = TRUE # Pole = TRUE if data is pole trend, plunge, FALSE if data is dip-dir,dip Corrected = TRUE # Corrected = TRUE if a Terzaghi correction has been made END BEGIN FIELD Name = "Survey ID" # Files can be linked to borehole and traceplanes through # an integer Survey ID in the .SAM file to facilitate Terzaghi correction. Type = INTEGER END BEGIN FIELD Name = "Tracelength" Type = REAL END END #now the data #TREND PLUNGE SID TRACELENGTH BEGIN DATA Count = 200 211 29 1 3.96 224 34 1 7.05 201 13 1 2.22 <<additional data records>> 341 36 1 4.73 314 31 1 7.73 END |
Table 2-14 Borehole Data Format (.SAM)
| # Any record beginning with a # is a
comment #Borehole BEGIN borehole name = "Borehole NE-1X" survey_id = 1 origin = 80 -80 0 scan_trend = 0 scan_plunge = 0 scan_length = 160 radius = 0.12 END #Traceplane BEGIN traceplane name = "Tracemap XJX-43-2K" survey_id = 2 origin = 80 80 0 scan_trend = 0 scan_plunge = 0 scan_length = 160 tran_trend = 0 tran_plunge = 90 tran_width = 160 END |
The verification case was defined by generating two overlapping Fisher distributed fracture sets using the statistics given in Table 2-15. The stereoplot before fracture separation by ISIS/N is provided in Figure 2-39. The statistics for the fracture sets following neural network analysis, and the stereoplots for the fractures assigned to the sets are provided in Figure 2-40.
Table 2-15 ISIS/N Verification Case
Set |
Orientation Distribution | Mean Pole Trend |
Mean Pole Plunge |
Dispersion |
1 |
Fisher | 36.67 |
89.45 |
10 |
2 |
Fisher | -0.28 |
59.65 |
20 |
2.2.6 Task 2.1.4 Compartmentalization Analysis
A major effort was made during the quarter to develop new technologies for compartmentalization analysis. Compartmentalization in reservoirs where fractures dominate permeability often leads to wells that produce at different rates and volumes than expected. Within compartments, pressure communication can be nearly instantaneous, while nearby wells in different compartments may have very little communication. This situation can lead to several undesired consequences, among them, unanticipated interference among wells, reduced recovery efficiency, and increased production uncertainty. Estimating the degree of compartmentalization is important at all stages of field development in order to properly engineer the field and to provide realistic recovery estimates and rates for financial decisions.
Compartmentalization may be due to several factors: fault offset of the producing horizon; the existence of high permeability sub-vertical faults that form barriers to lateral migration; the reduction in permeability of fractures due to mineralization; and the natural geometrical clustering of joint and smaller fault networks when matrix permeability is low. In terms of thermal recovery projects, there is an additional type of compartmentalization that is dynamic. Injected steam loses heat to the matrix and condenses. This condensation "front", which may be less of a definable surface than it is a series of fingers reaching out from the well along fracture paths, forms the steam "compartment". Knowledge of the shape and extent of this compartment is crucial to optimizing production. It is a combination of the natural geometrical compartmentalization of the fracture networks, and the ability of the fractures in the networks to exchange heat with the matrix. The geometry of the fracture networks places constraints on the steam compartments; the steam compartments can never be bigger than the networks themselves.
It is usually possible to identify large-scale fault-bounded reservoir compartments from seismic or production histories. It is far more difficult to assess the compartmentalization due to joint network geometry and connectivity, for which seismic is of little use. Joint network compartmentalization is often suspected when static and dynamic calculations of recoverable oil or gas do not agree, and their is no evidence for fault-offset or others types of fault-related compartmentalization in a fractured reservoir.
For this reason, two approaches have been developed and are being evaluated to address the more general problem of joint network compartmentalization for non-thermal applications, and also to extend the calculation for steam injection applications.
Figure 2-41 shows a hypothetical discrete fracture network. In this model, there are clusters of fractures, despite the apparent high density of fracturing. Each well shown in the model intersects a different cluster. In other words, these wells could only communicate with each other through the matrix. If this were a field, it would be important to determine the volumes of these compartments, which ultimately relates to the amount of hydrocarbon that a well can produce, and the horizontal dimensions of these compartments, which relates to optimal well spacings and the efficiency of primary or secondary recovery for a specified well pattern.
This task is focused on development of the technology for analysis of compartmentalization in fractured rock masses, where the compartmentalization is controlled by such fracture location heterogeneity. This task is related to Tasks 1.3.1 and 1.3.2, which address similar issues of fracture geometry, and Task 2.1.2, which deals with fracture spatial location distributions.
2.2.6.1 Algorithm to Compute the Volume and Horizontal Extent of Joint Network Compartments
The computation of the volume and horizontal extent of joint network compartments is a three-step process:
The identification of the individual fracture networks is a straightforward process for which robust algorithms have been developed. After the DFN model has been created, the intersections among the fractures is calculated. This leads to a symmetrical matrix of "1"'s and "0"'s, where "1" indicates an intersection. This connectivity matrix is then searched beginning with fracture number 1 in a standard tree search. This search process is repeated for all fractures that do not belong to any previously identified cluster until all fractures in the model have been processed. The results of this first step are saved as a list of fractures belonging to each cluster.
The second step is the calculation of the bounding surface for each cluster. Each network may be irregular. While it would be possible to compute the "bounding box" for a network, and use this box volume and horizontal cross-section as surrogates for compartment volume and horizontal extent, this would lead to an overestimation in most cases of both volume and cross-section. This in turn would produce overestimates of the ultimate recovery from wells, and suggest greater well spacings and recovery efficiencies than would actually be the case (Figure 2-42).
To reduce the potential for overestimates, it is necessary to calculate a bounding surface that better approximates the outer limits of the network. A convex hull, as described in Section 2.2.3.1.2 above meets these requirements.
Figure 2-43 shows the three-dimensional convex hull calculated using the QuickHull algorithm and the Qhull software package (Barber et al., 1995). As this figure shows, the shapes of the discrete fracture network (DFN) of Figure 2-41 would not be well-approximated by a bounding box for many of the networks.
The final step is to compute the volume and cross-sectional area of each hull. This is easily done. First, Qhull contains options to compute the volume of the hull, and to output the coordinates of vertices belonging to the hull. From these vertices, it is possible to compute another convex hull which represents the horizontal extent (Figure 2-44).
For vertical wells, the probability of whether the well will intersect the cluster (given that the well will be drilled deeply enough to intersect the cluster if the cluster lies under the well) does not depend upon the Z-coordinates of the cluster. The intersection probability is only a function of the X- and Y-coordinates. The shape and area of this horizontal extent is easily determined by calculating the two-dimensional convex hull of the (X,Y) data points belonging to the three-dimensional convex hull that bounds the network cluster. This algorithm can easily be extended to non-vertical wells by projecting the coordinates of the three-dimensional hull onto a plane that is perpendicular to the well trajectory.
Results for the volume and horizontal extent for the hypothetical model are shown in Table 2-16.
Table 2-16 Calculation of Volume and Horizontal Extent for Hypothetical Model
| Volume (m3) | 156,789 |
| Horizontal Area (m2) | 4369.0 |
| Bounding Box (m) | 81.0 x 82.4 |
| Hull/Bounding Box | 65% |
2.2.6.2 Algorithm to Compute Steam Compartmentalization
The compartmentalization defined by the condensation front of injected steam is a dynamically-changing process. As the matrix heats up with time due to steam injection, the condensation front moves downstream into new fractures. Not all of the matrix encompassed by the fracture compartment defined by fracture geometry is immediately heated by injected steam. Factors such as the conductivity of the fractures, their surface area, and the adjacent rock temperature gradient influence this process.
While it is not possible to carry out a detailed, dual-phase, dual-porosity simulation of this process for large fracture networks, it is possible to capture first-order effects of this process for the purposes of estimating steam compartmentalization.
The algorithm builds upon the geometrical clusters identified in the previous steps. To evaluate the rate and shape of compartment growth, a well is inserted into one the of geometric network clusters. Heat "particles" are injected into the well in the same way that tracer particles are injected in particle-tracking schemes in groundwater flow modeling codes. Miller et al. (1995) describe the MAFIC code which carries out this type of particle-tracking simulation for tracers. The equations for tracer transport and heat transport (as hot water) are identical, except for the need to have heat capacity for the fracture and a heat capacity for the effective bulk matrix. To simplify the process, diffusive transport within the fracture flows is neglected. We assume that the dominant mechanism for heating the matrix is conduction of heat from the fracture to the matrix, and that other heat transfer mechanisms are less important.
In this initial effort, overburden heat losses are ignored. Eventually, it may be necessary to modify the algorithm to represent conduction of heat into semi-infinite media to account for this heat loss mechanism.
This algorithm represents a first stage in modeling steam compartment formation through time. Clearly, flow of condensing steam is not the same as flow of an incompressible, single phase fluid, but the single-phase approximation used in this algorithm is a reasonable starting point to learn about heat transfer in a complicated discrete fracture network. In reality, as steam condenses in the fractures, the velocity of flow changes due to steam collapse and density contrasts near the periphery of the steam zone. It may be possible to use an effective viscosity approach to model the moving steam in a way that accounts for the change in steam quality with distance. The effective viscosity would depend upon steam quality, which in turn, is a function of the heat balance along the fracture.
2.2.7 Task 4.1.1 Fracture Image Data Acquisition
During the quarter, Marathon collected and processed fracture image data, and provided the data for posting on the WWW server.
2.2.8 Task 4.1.2 Well Testing Data Acquisition
During the quarter, Marathon collected and processed well test and hydraulic response data, and provided the data for posting on the WWW server. The well test and hydraulic responses provided by Marathon during the quarter are summarized in Figures 2-45, 2-46, and 2-47.
2.2.9 Task 5.1.2 WWW Site Updates
During the quarter, Golder Associates updated the World Wide Web (WWW) Server for the project. The update included the following postings::
2.2.10 Task 5.2.1 Progress Reports
The quarterly progress report for the period March 7, 1996 to June 6, 1996 was prepared during the quarter, and was posted to the WWW.
The following significant management activities were carried out during the quarter:
No project deliverables were prepared or delivered during the quarter.
This section presents revisions to the schedule and deliverable list provided with the first quarterly report. These revisions reflect changes expected due to logistical constraints for coordination between Golder, Marathon, and MIT. No change has been made to the project work scope.
The following revisions have been made to the project deliverable schedule:
The revised schedule for Task 2 is provided in Table 3-1. The schedule for other project activities remains as in our previous report.
Table 3-1 Schedule Revisions
| 2. Technology Development | |||||
| 2.1 | Fracture Data Analysis Technology | 96.03.07 | 97.01.31 | Technology Available on WWW
Research Report |
97.01.30 97.03.15 |
| 2.1.1 | Fracture Set Analysis | 96.03.07 | 96.07.31 | ||
| 2.1.2 | Spatial Location Analysis | 96.07.01 | 97.01.30 | ||
| 2.1.3 | Hydraulic Parameter Analysis | 96.09.01 | 96.11.30 | ||
| 2.1.4 | Compartmentalization Analysis | 96.09.01 | 96.11.30 | ||
3.2 Milestones and Deliverables
The following project milestones during the quarter are described in Table 3-2. Milestones and deliverables during the quarter were met.
Table 3-2 Milestones and Deliverables
| Milestone or Deliverable | Scheduled Date | Delivery Date |
| Data Available on WWW | 96.06.30 | 96.06.15 |
| First Quarter Progress Report | 96.06.15 | 96.06.26 |
