1. Introduction
The detailed study of the geometry of structural traps is one of the most common methods used when attempting to determine the position and intensity of fracture development within a reservoir. Fracture prediction is commonly based on geometric and kinematic models such as analysis of fold curvature (Lisle, Reference Lisle1994, Reference Lisle, Cosgrove and Ameen2000; Stewart & Podolski, Reference Stewart, Podolski, Coward, Daltaban and Johnson1998; Fischer & Wilkerson, Reference Fischer and Wilkerson2000; Hennings, Olsen & Thompson, Reference Hennings, Olson and Thompson2000; Bergbauer & Pollard, Reference Bergbauer and Pollard2004; Bergbauer, Reference Bergbauer, Jolley, Barr, Walsh and Knipe2007) or seismic-based techniques (Mueller, Reference Mueller1992; Gray, Roberts & Head, Reference Gray, Roberts and Head2002; Hall, Kendall & Barkved, Reference Hall, Kendall and Barkved2002; Hall & Kendall, Reference Hall and Kendall2003; Masaferro et al. Reference Masaferro, Bulnes, Poblet and Casson2003), but far less commonly with geomechanical modelling (Bourne et al. Reference Bourne, Rijkels, Stephenson and Willemse2001; Dee et al. Reference Dee, Yielding, Freeman, Healy, Kusznir, Grant, Ellis, Lonergan, Jolly, Rawnsley and Sanderson2007; Olson, Reference Olson, Lewis and Couples2007; Wilkins, Reference Wilkins, Jolley, Barr, Walsh and Knipe2007; Smart, Reference Smart2009).
‘Curvature analysis’, in which it is assumed that regions of maximum curvature on a surface coincide with regions of maximum fracturing, is such a method. Curvature can be measured in several ways, e.g. by an inscribed circle or as the second derivative (Rijks & Jauffred, Reference Rijks and Jauffred1991; Antonellini & Aydin, Reference Antonellini and Aydin1995).
For instance, a reservoir characterization mainly based on rock typing and its extension to uncored wells, a fracture characterization and modelling study, and on the construction of a full field geostatistical fine grid model has been done by Beicip-Franlab (Galard et al. Reference Galard, Zoormand, Ghanizadeh, Daltaban and Camus2005).
On the other hand, ‘strain analysis’, in which it is assumed that regions of maximum strain coincide with regions of maximum fracturing, is another, relatively new method of fracture prediction in a reservoir which, with the advances in computer processing technology, can now be extended into three dimensions.
The Gachsaran anticline is one of the most important giant oil fields in the Zagros fold–thrust belt, and consequently it has excellent seismic coverage and over 350 wells. These data together with the production index data provide an ideal opportunity to compare the two fracture prediction techniques. This comparison is the focus of this paper.
In this study, a curvature analysis, based on the variation of the second derivative of the structural surface that defines the reservoir, was carried out using new techniques and software for determining the isodip contours, and the locations of maximum and minimum curvature.
2. Geological and tectonic setting
The Gachsaran anticline is situated in the Iranian section of the Zagros fold–thrust belt, a belt which stretches from the Anatolian fault in eastern Turkey to the Minab fault near the Makran region in the SE of Iran. This fold–thrust belt is part of an active foreland basin complex which was initiated in Late Tertiary time (Stocklin, Reference Stocklin, Burk and Drake1974; Berberian & King, Reference Berberian and King1981) by the collision between the Iranian and Arabian plates (Fig. 1). Geologically, the belt forms the northeastern border of the Arabian plate; geographically it is located in SW Iran.
Figure 1. The Iranian section of the Zagros fold–thrust belt showing some of the main belt parallel faults (the Main Zagros Fault and the High Zagros Fault) and some major faults that traverse the belt (e.g. the Kazerun Fault). The belt formed as a result of the collision between the Arabian and Iranian plates. The location of the Gachsaran anticline in the upper part of Dezful South is shown in the rectangle.
The Main Zagros Fault, which marks the northeastern margin of the belt, represents a suture zone, (Fig. 1). This zone is now seismically inactive and the seismicity has migrated southwestwards towards the frontal part of the belt (Berberian, Reference Berberian1983). The belt is divided into two NW–SE-trending domains by the High Zagros Fault. These are the Imbricate Zone or High Zagros located between the Main Zagros Fault and the High Zagros Fault, and the Folded Zagros Zone, which runs southwest of and parallel to the High Zagros Fault (Fig. 1).
In addition to this parallel structural division of the belt that is based on the dominant mechanism of shortening (thrust dominated compared to fold dominated domains), the belt has also been divided along its length into the Lurestan, Dezful Embayment and Fars regions (Fig. 1). These regions are separated by major N–S- and E–W-trending faults across which major sedimentary facies changes occur (Motiei, Reference Motiei1994). Indeed, these faults have controlled the location and orientation of some of the sedimentary basins and have had a marked impact on sedimentation in these regions (Setudehnia, Reference Setudehnia1978). Based on these divisions, the Gachsaran anticline is located in the Dezful Embayment region and more accurately in Dezful South near the Zagros mountain front (Fig. 1).
This field is an elongated doubly-plunging anticlinal structure with dimension of 65 km length and 4–8 km width, and it is considered one of the most important productive oil fields in the Oligocene–Lower Miocene carbonate horizons (Asmari Formation) and the Middle Cretaceous carbonate horizons (Sarvak Formation) (Fig. 2; Motiei, Reference Motiei1994). Because of the large amount of data available on the Asmari reservoir horizon, this study focuses exclusively on these Tertiary carbonates.
Figure 2. The Mesozoic–Cenozoic stratigraphic column of the Dezful Embayment area showing the position of the two major carbonate reservoirs in the Gachsaran anticline field, namely the Asmari and Sarvak formations.
3. Curvature analysis
Attempting to map the position of fractures and fracture intensity in all parts of the field on the basis of well data alone, even when many wells have been drilled as in the reservoir under consideration, is extremely difficult. Many wells have no image logs. Cores, when available, are often not oriented, and often wells have not been drilled in key regions of the reservoir. Consequently, other techniques for determining the state of fracturing in a reservoir have been developed. Two of the most used are curvature analysis and strain analysis (see Section 4).
Although it is known that fracture intensity on a field scale is dependent on both structural position within the reservoir and lithology, curvature analysis takes no account of lithology and determines fracture position, orientation and intensity solely on the basis of the variation in curvature of the surface defining the reservoir bed. The mapping and analysis of bedding plane dip variations has been used to predict fracture density in different parts of an anticlinal reservoir (Lisle, Reference Lisle1994). The basis of this method is that fold-related fractures are created as a result of the flexing (buckling) of the bedding, which is assumed to behave as an elastic plate. Areas of maximum curvature are assumed to be areas of maximum fracturing. Such a relationship has been noted by Rijks & Jauffred (Reference Rijks and Jauffred1991) in the North Sea and by Antonellini & Aydin (Reference Antonellini and Aydin1995) in North America.
The maximum and minimum curvature at any point occurs parallel to the directions of maximum and minimum dip variation, respectively. Fracture orientation is taken to be parallel to the direction of minimum curvature, and maximum curvature is used as a proxy for fracture density. The absolute amounts of the minimum and maximum curvature and of the ratio of the maximum to minimum curvature are well known for the various folds and fold-related structures (cylindrical folds, domes, saddles, etc.) (see e.g. Masaferro et al. Reference Masaferro, Bulnes, Poblet and Casson2003).
For the curvature analysis reported in this study, the latest sub-surface structural contour map has been used (Fig. 3). This has been constrained by being tied to well tops.
Figure 3. Well-tied sub-surface contour map of the Oligocene–Lower Miocene carbonate horizon (the Asmari Formation) that represents the upper reservoir unit in the Gachsaran field. Lishtar as a closure in the NW of the main closure is commonly considered a part of the Gachsaran oil field.
Regions of artificially high curvature may appear near faults as a result of ‘smoothing’ across the fault edge. In such a situation the curvature will reflect the amount of displacement on the fault and will be unrelated to the folding. In any curvature analysis it is clearly important to eliminate such fault effects (Fig. 4).
Figure 4. The curvature distribution map for the top of the Asmari Formation within the Gachsaran field. Note that the effects of faulting have been eliminated.
Gaussian curvature is the product of the minimum and maximum curvature at any point. Unlike simple curvature it has the dimensions of 1/area, so special consideration has to be given to the interval over which the curvature is to be measured. Conventionally the curvature is taken to be positive for convex regions, negative for concave regions and zero for planes. Surfaces which have a Gaussian curvature are non-cylindrical. Cylindrical folds have a Gaussian curvature of zero and the value of this curvature for a conical fold is dependent on the wavelength of measurement (Stewart & Podolski, Reference Stewart, Podolski, Coward, Daltaban and Johnson1998).
One of the objections to the use of curvature analysis for fracture predictions is that it considers only the final shape of the structure and does not consider the tectonic events the reservoir has passed through or the detailed structural evolution of the fold. For instance in a box-shaped fold the curvature in the limbs is low, almost zero (i.e. there are no dip changes over significant distances), and so it is to be expected that the fracture density in these regions will be low. However, when the evolution of such folds (which includes limb rotation and hinge migration) is considered, it is clear that most of the deformation occurs in the limbs and that these will therefore be areas of high fracture density.
4. 3D strain analysis
3D strain analysis of a reservoir is based on knowledge of the tectonic evolution of the structure and of the evolution of the basin containing the structure. This method depends upon being able to determine the structural evolution of the reservoir. If this can be achieved then the reservoir can be restored back to its original geometry (i.e. pre folding and thrusting). This type of strain analysis shows the total deformation of a particle in 3D regardless of whether or not the deformation was in the main transport direction (Kane et al. Reference Kane, Williams, Buddin, Egan and Hodgetts1997). Strain analyses can be divided into two types, namely ‘dilation analysis’ and ‘finite strain analysis’.
4.a. Dilation analysis
This method analyses the dilation, i.e. the ratio of the change in length, area or volume of the strained object, relative to the unstrained values. It can be applied to lines, areas and volumes and can be positive or negative. Consequently the volumetric dilation is calculated by: volumetric dilation = (1 + e1) (1 + e2) (1 + e3) − 1 (where e is elongation).
4.b. Finite strain analysis
The finite strain method can only be applied to volumes. Analysis is carried out on the change in position of the vertices of each tetrahedra. By default the X,Y,Z position of each vertex is expressed relative to the central ● position of the tetrahedra. The change in position of the vertices of the strained, relative to the unstrained tetrahedra (not the centre of the tetrahedra) allows the calculation of the volumetric dilation and principal strain values/axes. These axes are then positioned in the centre of the tetrahedra for strain visualization purposes (Fig. 5).
Figure 5. Strain tetrahedra with strain ellipsoid. The XYZ axes are equivalent to eigen values, which are equivalent to the principal strain axes. Principal strain is a measure of the new magnitude of the xyz axes of the strained object relative to its unstrained position.
Strain orientations were also calculated from the coordinate position of the strained object relative to the coordinate position of the unstrained object. Therefore, with calculation of the magnitude and the ratios of the principal strain axes it is possible to determine the plane strain ratio and the nature of the strain (plane, constructional or flattening).
The same structural map used in the curvature analysis discussed earlier was used for the strain analysis method in order to facilitate the comparison between the two techniques of fracture analysis. However, because of the complex and time consuming processes involved in the strain analysis, a coarser grid was used. As noted earlier, corrections need to be made for the ‘smoothing’ that occurs at faults because of the need to generate a continuous surface to describe the reservoir surface.
Although strain analysis can be carried out in one, two or three dimensions, it is better to apply strain analysis in three dimensions, i.e. to consider the vertical variations within the different reservoir zones. In addition, it is also important to take into account other reservoir data such as that coming from the ‘mud loss model’ and the ‘facies model’. A volumetric structural model that does this for the Gachsaran field is provided as a final deformed object.
4.c. Flexural slip restoration
One of the most important stages in any strain analysis is the decision relating to the selection of the restoration mode and the algorithms to be used. These depend to a great extent on the knowledge of the regional geological setting and the structural evolution.
In parallel folds developed in bedded sediments the flexural slip unfolding algorithm can be applied. The algorithm works by restoring a folded template and its template-parallel slip system to a horizontal or regional datum. It is therefore ideally suited to use on folds such as the Gachsaran anticline that characterize this part of the Zagros belt. This restoration removes the bedding parallel slip associated with the folding (e.g. Griffith et al. Reference Griffiths, Jones, Salter, Schaefer, Osfield and Reiser2002).
There are some principles for this algorithm: (a) line length of the template bed is preserved in the unfolding direction; (b) line length of all beds parallel to the template bed is maintained in the unfolding direction; (c) surface area of cylindrical folds with one fold generation is preserved; (d) volume (area in 2D) of the fold is maintained; (e) relative bed thickness is constant (discrete slip between beds will change thickness measured at specific points along the template bed).
As sedimentary beds were deposited horizontally in sedimentary basins, they are restored to this orientation. It is acknowledged that the model would display some ‘noise’, especially near the boundaries. Consequently, it is important that errors linked to this are removed.
It has been established that in the region of the Gachsaran fold no phase of tectonic extension has occurred since the Oligocene–Lower Miocene carbonate sedimentation. Thus, the final strain state in the model is the cumulative dilation model for the Asmari Formation in the Gachsaran anticline.
Finally, a 3D strain model based on flexural slip restoration is created (Fig. 6). This figure illustrates some dip and strike sections of the 3D model, and many spatial variations are observable. For example, computed strain values are generally larger in the tight folded central part of structure as well as the highly deformed steep flanks.
Figure 6. A 3D strain model of the Asmari reservoir within the Gachsaran field that is visualized in dip and strike sections to be more perceptible. It is created based on flexural slip restoration, and great variations in modelled strain are noticeable both vertically and laterally.
In addition, a maximum principal strain axes model has been determined (Fig. 7), which is another essential piece of data for determining the orientation and distribution of fractures. It can be seen that the maximum principal strain axes increase sharply in the fold limbs (Fig. 7).
Figure 7. Maximum principal strain axes model of the Asmari Formation in the Gachsaran anticline. Note that noise has been removed already.
5. A comparison of the curvature and 3D strain analysis methods
Curvature analysis examines the curvature variations on a surface and is essentially 2D (Fig. 4), and in this regard is different from the model analysed in the strain analysis method in which a 3D stain analysis yields a 3D fracture pattern. Figure 8 shows the result of a curvature analysis of the Asmari Formation in the Gachsaran anticline that is projected onto the base of the formation together with a few vertical sections that give an indication of the 3D strain distribution within the structure.
Figure 8. A curvature map for the top of the Asmari Formation accompanied by some sections of absolute volumetric strain; for better visualization the curvature map is projected on to the base of Asmari Formation. Vertical variation in the strain model is clear, which leads to some considerable dissimilarity with curvature trends.
A comparison of the average strain intensity and curvature maps for the top of the Asmari Formation is illustrated in Figures 9 and 10. Overall, they are generally quite similar but there are some big differences in the detail. For example, closer inspection of the models shows that in a section normal to the fold axis, the strain model shows a peak in intensity near the main flexure of the anticline and that this decreases towards the limbs reaching a minimum at around 45 degrees after which it begins to increase again.
Figure 9. 3D strain average map calculated for the Asmari reservoir interval overprinted with the production index (PI) bubbles for all the wells with PI data (the bubble radius is an index of PI). A good correlation, especially for highly productive wells, is observable where larger circles are mostly concentrated in high strain areas, mainly in the southern flank. There is a poor correlation on the NW flank of the anticline (rectangle) that is thought to be affected by local fault activity.
Figure 10. Curvature map calculated for the top of the Asmari Formation overprinted with the production index (PI) bubbles for all the wells with PI data (the bubble radius is an index of PI). In general, a moderate correlation between PI and curvature is observable.
In contrast, in the curvature model, the intensity decreases continuously and becomes zero when it passes through a region of no dip change (Figs 9, 10). In the strain analysis, the 3D strain values are seen to change vertically and in some areas their values can be reversed, i.e. extension at the reservoir top and compression at its bottom.
Another major difference between the two methods is their sensitivity to faults. It was noted earlier that the smoothing effect used in curvature analysis can result in faults generating significant local curvature patterns unrelated to folding (Fig. 10).
6. A comparison of the curvature and 3D strain method predictions with the production index
The average strain distribution and the curvature distribution maps for the top of the Asmari Formation are shown in Figures 9 and 10, respectively, overprinted with the production index (PI) bubbles for all the wells with PI data in the Asmari reservoir of the Gachsaran field (the bubble radius is an index of PI). Accessible production data as a well-known formation productivity test and a suitable indicator for fracture existence were limited but fairly well distributed. Thickness corrections have been carried out on these data via dividing the PI by the thickness of the producing zone. It is considered a modification of PI data to highlight the effect of fracturing.
A better correlation is found between strain intensity and fracturing (high PI) than between curvature and fracturing. In almost all the areas with a high PI value the strain intensity is also high. In one area located on the northwest flank of the anticline (shown by a rectangle in Fig. 9), a poor correlation is observed. It is suggested that this is because of local fault-related fracturing producing fractures not predicted by the strain model.
7. Conclusion
The main conclusion reached in this comparison of the fracture predictions of a curvature analysis and a strain intensity analysis are that in two important traits the strain analysis performs better. The first is that it is a 3D analysis capable of determining the strain distribution throughout the reservoir, i.e. within the various units of the reservoir compared to the curvature analysis which considers the geometric properties of a single surface (2D) and takes no account of the geometry in the third direction. In addition, there is a better correlation between 3D strain and the production data, in particular the production index, which is directly related to fracture intensity.
Acknowledgements
The authors are grateful to M. Pourkermani, P. Omodi and J. Cosgrove for sharing their experience and advice. J. Cosgrove is gratefully acknowledged for his constructive review, which considerably improved the final version. We acknowledge our colleagues at NIOC who participated in regional studies and field campaigns: M. G. Hasan Goodarzi, M. Valinejad, M. H. Norouzi, F. Beik, K. Shah Hoseini. We thank also the NIOC Research and Technology Directorate for its financial support of our research project.
