Ms-p-nioc

Analysis of Seismic Data Acquisition Techniques, Parameters and Results in High Velocity Layers of Iranian South West Hydrocarbon Prospects By: SeyedMohsen SeyedAli A THESIS SUBMITED IN FULFILMENT OF THE REQUIREMENT FOR THE DEGREE OF MASTER OF SCIENCE

Supervisors: Dr. Mohammad Ali Riahi Dr. Alireza Alizadeh Attar Co-Advisor: Roohollah Noemani TEHRAN MARCH 2009

Copyright © 2009 SeyedMohsen SeyedAli All Rights Reserved

This thesis is financially supported by Research and Development of National Iranian Oil Company, to which I am greatly indebted.

To My Parents

ABSTRACT Acquiring good images is a debatable challenge in seismic acquisition design. The problem is more robust in areas where high degree of anisotropy and complex structures are present. In excess of these problems, high attenuating formations will reduce the energy content of received signals at the surface. Hence, having reliable images is a challenge in such areas. Conventional methods of data acquisition will result in all above mentioned problems. So it is necessary to utilize more recent approaches in data acquisition strategies to overcome such challenges. One of the most problematic lithologies in seismic acquisition and processing operations are salt layers and salt bodies. They cause several problems due to their intrinsic properties such as high attenuation and high velocity. In addition their steep flanks may cause to energy travels distant from the source location which it causes lost of high amount of energy as a result of geometrical spreading and in some situations there is the probability that the energy never be recorded due to inappropriate estimation of Xmax or migration aperture in acquisition design. Multi-component seismic recording captures the seismic wavefield more completely than conventional single-element techniques. In the last several years, multicomponent surveying has developed rapidly, allowing creation of converted-wave or P-S images. These make use of downgoing P-waves that convert on reflection at their deepest point of penetration to upcoming S-waves. Survey design for acquiring P-S data is similar to that for P-waves, but must take into account subsurface VP/VS values and the asymmetric P-S ray path. P-S surveys use conventional sources, but require several times more recording channels per receiving location. This study analyzes the effectiveness of multi-component data acquisition in areas where high velocity layers such as salt bodies and carbonates exist.

i

ACKNOWLEDGMENT

First of all I would like to express my special appreciations to my thesis supervisor Dr. Riahi for his great encouragements and guidance during my study. I thank to my supervisor Dr. Attar whose ideas and valuable suggestions was the basis for this study. I thank to my co-advisor Mr. Noemani-Rad who provides me with all required facilities for this project. I also would like to express my best regards to Professor J.L. Mari who showed me the way of thinking about mathematical problems in geophysics. I acknowledge to NORSAR Innovation AS for providing me with trial license of NORSAR2D software.

ii

Table of Contents ABSTRACT ...................................................................................................................i ACKNOWLEDGMENT ...............................................................................................ii Table of Contents..........................................................................................................iii List of Figures................................................................................................................ v CHAPTER 1 .................................................................................................................. 1 INTRODUCTION ......................................................................................................... 1 1.1. Framework and objectives this thesis .................................................................2 1.2. Methodology .......................................................................................................2 1.3. Review of some affecting parameters in seismic data acquisition: ....................3 1.3.1. Attenuation:..................................................................................................3 1.3.2. Arrays...........................................................................................................5 1.3.3. Statistics .....................................................................................................10 1.4. Forward Modeling ............................................................................................10 CHAPTER 2 ................................................................................................................ 12 REVIEW OF SOME PREVIOUS SEISMIC DATA ACQUISITION JOBS IN DEZFUL EMBAYMENT AND GEOLOGICAL STUDIES..................................12 2.1- Geological study of Dezful Embayment ..........................................................12 2.1.1. Stratigraphy of the Gachsaran Formation ..................................................14 2.1.2- Structural Geology of Gachsaran Formation.............................................17 2.1.2. Review of seismic challenges of Gachsaran formation .............................18 2.2. Review of some previous seismic data acquisition jobs in Dezful Embayment ..................................................................................................................................21 2.2.1. Seismic data acquisition on Gachsaran Formation (1971-1972) ...............22 2.2.2. Seismic Data Improvements Experiments (1972 and 1973)......................24 2.2.3. Conclusions................................................................................................25 CHAPTER 3 ................................................................................................................ 27 FORWARD MODELING AND RAY TRACING APPROACH .......................27 3.1. Introduction.......................................................................................................27 3.2. Ray Tracing Approach......................................................................................28 3.2.1. Fundamentals of Ray Tracing....................................................................28 3.2.2. Ray Tracing Techniques ............................................................................32 3.2.3. Limitations of ray tracing...........................................................................34 3.3. Modeling ...........................................................................................................36 3.3.1. Introduction................................................................................................36 3.3.2. Modeling Procedure...................................................................................37 3.3.3. Simple Model.............................................................................................38 3.3.4. Complex Model (Based on a real field data in Dezful Embayment).........44 CHAPTER 4 ................................................................................................................ 54 COMMON CONVERSION POINT MAPPING FOR 2D-3C SEISMIC DATA ACQUISITION SURVEY DESIGN IN ANISOTROPIC MEDIA ............................ 54 4.1. Introduction.......................................................................................................54

iii

4.2 Review of seismic anisotropy ............................................................................55 4.2.1 Anisotropy versus heterogeneity.................................................................56 4.2.2 VTI (Polar Anisotropy)...............................................................................57 4.3. Definition of Seismic Anisotropy Parameters ..................................................58 4.3.1. Elastic Tensor and Anisotropy Parameters ................................................58 4.3.2. Explanation of Thomsen Anisotropy Parameters ......................................61 4.4. Conversion Point Mapping ...............................................................................61 4.4.1. Introduction................................................................................................61 4.4.2. Theory ........................................................................................................63 4.5. Methodology .....................................................................................................64 4.6. Results and Discussion .....................................................................................66 CHAPTER 5 ................................................................................................................ 71 CONCLUSIONS AND RECCOMENDATIONS....................................................... 71 5.1. Conclusions.......................................................................................................71 5.1.1 Affecting parameters according to previous jobs .......................................71 5.1.2. Forward modeling approach ......................................................................72 5.1.3. Common conversion point mapping..........................................................73 5.2. Recommendations.............................................................................................73 References.................................................................................................................... 74

iv

List of Figures Figure 1- 1: Schematic of a 2-D array ...........................................................................7 Figure 1- 2: Response curve of an array with twelve receivers and 0.5m receiver interval ...........................................................................................................................8 Figure 1- 3: Response curve of an array with twenty four receivers and 0.5m receiver interval ...........................................................................................................................8 Figure 1-4: Midpoint coverage for: a. a single source into a linear geophone array (after Vermeer, 1998a), b. a linear source array into a linear geophone array, and c. a single source into an areal geophone array. .................................................................10 Figure 2- 1: Tectonic map of Arabian and Iran plates and main structural subdivisions of Zagros orogenic system (Abdolahiefard, 2006). .....................................................13 Figure 2- 2: Schematic Profile of subsurface layers in Dezful Embayment, Izeh and High Zagros zones. Gachsaran Formation is shown by yellow color (Sherkati et al. 2004). ...........................................................................................................................15 Figure 2- 3: Stratigraphy column of Cenozoic Era in Lurestan, Dezful Embayment and Fars area. ...............................................................................................................16 Figure 2- 4: Location of some Iranian oil field in Dezful Embayment and Abadan plane.............................................................................................................................17 Figure 3- 1: Seismic Model..........................................................................................28 Figure 3- 2: Elements of kinematic ray tracing ...........................................................31 Figure 3- 3: Depiction of image and normal incident rays ..........................................33 Figure 3- 4: Un-smoothed model section (Left section) and smoothed model section (Right section)..............................................................................................................35 Figure 3-5: Comparison of PP and PS results. Left figure shows the result of a PP data acquisition, as it is obvious in this figure, gas cloud has obscured the crest of the reservoir. Right figure shows the result of PS data acquisition that has resolved the top of the reservoir and a largr fault in the flank[3]. ...........................................................37 Figure 3- 6: Three layer flat model ..............................................................................38 Figure 3- 7: Representation of traced rays...................................................................39 Figure 3- 8: PP (left) and PS (Right) sections..............................................................39 Figure 3- 9: Incoming reflection angles of PP rays (blue) and PS rays (red) ..............41 Figure 3- 10: Up-going reflection angles of PP rays (blue) and PS rays (red) ............41 Figure 3- 11: Variation of An-elasticity factor versus offset (PP rays (blue) and PS rays (red)).....................................................................................................................42 Figure 3- 12: Geometrical Spreading versus offset (PP rays (blue) and PS rays (red)) ......................................................................................................................................43 Figure 3- 13: AVO analysis (PP rays (blue) and PS rays (red)) ..................................43 Figure 3- 14: Model Blocks .........................................................................................45 Figure 3- 15: Main model ............................................................................................46 Figure 3- 16: Results of normal inident ray tracing on the model illurtrated in figure (3-15) with Station spacing = 250 m ...........................................................................47 v

Figure 3- 17: Reflection from steep flanks ..................................................................47 Figure 3- 18: Results of normal inident ray tracing on the model illurtrated in figure (3-15) with station spacing = 50 m ..............................................................................48 Figure 3- 19: Results of image ray tracing on the model illurtrated in figure (3-15) with - migration aperture = 1km ..................................................................................49 Figure 3- 20) Results of image ray tracing on the model illurtrated in figure (3-15) with - migration aperture = 4km ..................................................................................50 Figure 3- 21: Results of image ray tracing on the model illurtrated in figure (3-15) with - migration aperture = 8km ..................................................................................50 Figure 3- 22: shot gather at X=10 km..........................................................................51 Figure 3- 23: Receiver locations versus mid-point, red crosses show the selected .....52 Figure 3- 24: Non- hyperbolic behavior at far offsets for CMP location of X=8 km..52 Figure 3- 25: Off-end survey. Snapshot from the source-receiver pair where source located at X=5 km and receiver at X= 15 km. .............................................................53 Figure 4- 1: Anisotropy versus Heterogeneity.............................................................56 Figure 4- 2: Velocity distribution in VTI media can be represented by a cylinder. The velocities on all directions within the horizontal plane are identical in amplitude......57 Figure 4- 3: Distinction between Ray (group) and Phase angles.................................60 Figure 4- 4: The conversion point traces a trajectory in the multi-layered model for a certain offset, instead of locating on a vertical axis (Stewart et al., 1999). .................62 Figure 4- 5: Depiction of common conversion point and different angle in an anisotropic media .........................................................................................................65 Figure 4- 6: Variation of displacement of common conversion point in aniotropic media (VTI) versus offset of isotropic conversion point for three values of delta* and epsilon equals to 0.2.....................................................................................................66 Figure 4- 7: Variation of conversion point in aniotropic media (VTI) versus ray angle for three values of delta* and epsilon equals to 0.2......................................................67 Figure 4- 8: Variation of displacement of conversion point with variation of δ * and ratio of isotropic conversion point offset to depth .......................................................68 Figure 4- 9: Variation of displacment of conversion point with variation of epsilon ( δ * = 0.2) ......................................................................................................................68 Figure 4- 10: Variation of displacement of conversion point in weakly Anisotropic medium (blue dots: exact equation results, red dots: linearly approximated results)..69 Figure 4- 11: Variation of displacement of conversion point in strongly anisotropic medium (blue dots: exact equation results, red dots: linearly approximated results)..70

vi

vii

CHAPTER 1 INTRODUCTION A good survey design achieves the geophysical objective while minimizing the cost and time of acquisition and processing. Basic components of the seismic method of hydrocarbon exploration are acquisition, processing, and interpretation. These three rather compartmentalized areas are inextricably related. No amount of clever processing can overcome some deficiencies in acquisition. No magical insight by the interpreter can remedy some acquisition and processing mistakes. Survey design is most effective when processing and interpretation are considered in the design. Emphasis in the pioneer stage of seismic exploration was on acquisition. Later, emphasis shifted to data processing. More recently, interpretation workstations have been a focus of development (Stone, 1994). Survey design is state of the art. Selecting the most efficient acquisition planning plays an important role in optimizing the cost and the results. By sure, the best design is the one that brings about the high resolution images of underground whilst providing us with required data for reservoir characterization. Because the final target of surface seismic surveys is exploring boundaries and its characterization. Delineation of reservoir boundaries will be more difficult when complex structures like salt bodies and tilted blocks exist. In such situations it is necessary to have a good comprehension of wavefield state in media. The problem is stronger where sever lateral velocity variation, high overburden and anisotropy are present. One of the best tools of having such understanding is forward modeling through which we can

1

analyze the ray path (ray bending, transmission and reflection) and all other related attributes. In addition it enables us to make choices of data processing steps.

1.1. Framework and objectives this thesis This project is subdivided into four main parts: 1- Review of seismic acquisition techniques in complex structures and analysis of important parameters on data quality (attenuation, static variations, array designing and etc.). 2- Review of some selected previous jobs on study area (Dezful Embayment). 3- Introducing multi-component seismic data acquisition technique in high velocity complex structure areas and representation of forward modeling approach as a clue of success in such regions. 4- Analysis of CCP locating and its importance as the base step of multicomponent seismic data acquisition design in anisotropic media.

1.2. Methodology After some studies and analyses of seismic data acquisition challenges and their results in Iranian SW hydrocarbon oilfields (Dezful Embayment), some typical models will be generated respecting to general structures and petrophysical properties of this area and different techniques of ray tracing will be used to extract several attributes from output data using NORSAR2D software. These attributes will be analyzed to obtain the most efficient techniques and parameters in designing surveys in this area. As this study concerns with acquisition of converted waves in anisotropic media, and most of the commercial softwares neglect the effect of anisotropy in investigation of reflection points, I have generated a MATLAB code to inspect the errors erecting due to neglecting the anisotropy factors in VTI media.

2

1.3. Review of some affecting parameters in seismic data acquisition: 1.3.1. Attenuation: Reduction in amplitude or energy caused by the transmitting media or system is called attenuation. It usually includes geometric divergence effects as waves spread out from a source as well as conversion of energy into heat (absorption) and other factors affecting amplitude, such as transmissivity losses and mode conversion. Spherical Divergence is the decrease in wave strength (energy per unit area of wavefront) with distance as a result of geometric spreading. A spherical wave traveling through the body of a medium continually spreads out so that the energy density decreases. For a point source the energy density decreases inversely as the square of the distance the wave has traveled; this means that the amplitude deceases linearly with the distance traveled. For energy which travels along a surface, the analogous term is cylindrical divergence, where the energy varies inversely as the distance and the amplitude as the square root of the distance (Sheriff, 2002). Another type of energy loss happens at interfaces, while the wave-front hits the interface. At the incident point location several phenomena such as mode conversion, transmission and reflection cause to initial energy breaks down to some other types which will result in attenuation of energy content in each type. Energy loss due to conversion of energy to heat (energy absorption) is another type of attenuation that is represented by attenuation factor (α) which is the exponential decay constant of amplitude of plane wave is traveling in a homogenous medium. The quality factor Q and its inverse Q-1, sometimes called the internal friction or dissipation factor and ζ is the logarithmic decrement. These quantities are related as follows:

1 αv ξ = = Q πf π

(1-1)

Where v is the velocity and f is the frequency. Attenuation coefficient and logarithmic decrement can be expressed as below:

3

ξ = ln(

A1 αv )= A0 f

(1-2)

Where A0 is initial amplitude and A1 is the recorded amplitude. Quality factor is the most common measure of attenuation and is a dimensionless parameter. As an intrinsic property of rock, Q is the ratio of stored energy to dissipated energy. Media with Q values higher than 10 assumed to have low energy absorption coefficient. In some cases that the absorption coefficient takes higher values Hamilton (1992) expressed this equation: 1 = Q

αv α 2v 2 πf − 4πf

(1-3)

However under low-loss assumption (Q>10) the term

α 2v 2 is negligible and may be 4πf

dropped. Some other experimental equations relating quality factor to media velocities have expressed by others such as

1 1000 2 =( ) by Udias, 1991 and Qp Vp

Qs 4 Vs 2 = [ ] by Waters, 1981. 3 V Qp p Among these attenuation parameters geometrical spreading plays the most important role while AVO variation analysis and comparison of P and S-waves. Although the logarithmic decrement( ξ ) for S-waves in some cases are higher than corresponding P-waves, but as the effect of this attenuation is applied exponentially as expressed in equation (1-1) and according to Waters (1981) that shows the ratio of quality factors for P and S-waves can be expressed respect to

Vs ratio, and respecting to the fact that Vp

this ratio increases with depth and tends to unity in highly overburdened layers their effects does not make too much differences in numerical modeling. Another measure of attenuation is Loss-Angle (φ) which can be expressed as the cotangent inverse of quality factor. Mc Donald et al. (1958) examined the dependency of geometrical spreading and energy absorption on wave frequency. Due to their experiments on a compact shale formation in Colorado State, geometrical spreading is more effective in low

4

frequencies and short distances respect to absorption attenuation. By increasing the distance and frequency, absorption increases until becomes more important. Above explanations and terminologies are summarized in table (1-1). Table 1-1: Attenuation terminologies and equations

Term

Symbol

Equation

Geometrical Spreading

G

I2 r = ( 1 )2 I1 r2

Quality Factor

Q

2πE ∆E

Absorption Coefficient

Α

Logarithmic Decrement

ξ

ln( A1 A0 )

Loss Angle

Φ

cot-1(Q)

(f

v) ln( A1 A0 )

1.3.2. Arrays Since early in the history of reflection seismic programs, it have tended to use multiple geophone sensors to record each trace of information. Typically traces are recorded sparsely (between 10 and 60 meter group intervals are typical). This has been due largely to recording channel limitations and the need to record a broad range of offsets. In the youth of seismic, it has been recognized that an organized spatial distribution of analogue recording elements would form a wavenumber domain filter when summed together. This lead to a popular practice of designing arrays to attenuate coherent noise patterns (those with wavelengths generally shorter than the array length). As broader bandwidth became a requirement for sharper imaging of subtle stratigraphic plays, arrays fell into disfavor. Many geophysicists perceived that conventional arrays were so long as to attenuate some of the high frequencies of the far-offset, shallow reflections. This led to a period of insanity where many geophysicists believed it wise to group or “pot” the geophones to eliminate that nasty array effect (Cooper. N, 2002)

5

1.3.2.1. Basic Theory of Array filtering: The filtering of a seismic function S, a succession of seismic impulses arriving spread out over time, by an array consisting of '2n+1' detectors is a spatial filtering over distance 'x'. At a given instance 't', the seismic function S(x,t) is written as S(x). If we call the filter operator w(x), then at every time sample t, the filtered trace is written as: G(x) = S(x)*w(x)

(1-4)

The symbol * represents the convolution operator. An array physically implemented by a series of detectors wired together at separation of 'e'. The array filter is an analog filter. In seismic acquisition survey design array filtering is sometimes called filtering by multiple geophones. Mathematically, the signal G(x) which is summation of signals S(x) at '2n+1' receiver locations separated by 'e' is written as: G(x) = S(x – ne) + S(x – (n – 1)e) + … + S(x) + … + S(x + ne)

(1-5)

In which G(x) = S(x)*w(x)

(1-6)

While w(x) = [D(x - ne) + … + D(x) + D(x + ne)]

(1-7)

And 'D' represents Dirac function. The filter w is an analog filter in wavenumber space (k). For an in-line array consisting of N=2n+1 detectors, Fourier Transform of the filter is written as: +∞

W (k) = FT (w(x)) =

∫ w( x)e

−∞

Where x = (N-1)e,

6

−2πik (( N −1) e )

(1-8)

Solving equation (1-12) numerically in discrete domain will result in: W (k) = AeB

(1-9)

In which A represents the module and B represents the phase. A and B are:

A=

sin Nπke , B = i(N-1)(πke) sin πke

(1-10)

If we plot the absolute value of A respect to k, for constant values of k and N, the response of the filter will be appeared. In figure 1-1 a schematic sketch of a 2-dimensional array has been represented. This array contains nine receivers which are distributed by distance e. All the recorded signals in these receivers will be stacked together and the result will be placed in the center location of the array. Figure 1-2 shows the response curve of an array with twelve receivers and 0.5m receiver interval, as we can see array filtering favors the low wavenumbers associated with reflected arrivals to the detriment of the higher wavenumbers associated with groundrolls. Increasing the number of receivers per array boosts the efficiency of array filtering. Figure 1-3 shows another array containing twenty four receivers with the same receiver interval.

e

n=4

Figure 1- 1: Schematic of a 2-D array

7

Figure 1- 2: Response curve of an array with twelve receivers and 0.5m receiver interval

Figure 1- 3: Response curve of an array with twenty four receivers and 0.5m receiver interval

8

Noise tests are essential in this step of acquisition design to acquire reliable information about the origin and characteristics of noises such as groundrolls to design an optimum array. Array design in 3-D surveys requires careful consideration. If a receiver array is used, it should be as simple as possible. Circles, box patterns, or in-line arrays are commonly used today because of the ease and efficiency of their layout. A designer may place a major emphasis on a complicated receiver array, but this message may be hard to convey to the field personnel who may not understand the necessity of carrying out instructions accurately. In such cases, the designer must work with the field personnel directly to achieve meticulous placement of the arrays. In 3-D surveys, source energy arrives from many directions. The array response varies dramatically depending on the azimuth between point sources (e.g., dynamite) and receiver arrays. If geophone arrays are laid in-line, the response from a point source fired from an inline position will be attenuated, while a broadside source into such an array will not be affected by the array. Because of this variation, many companies prefer no arrays at all when recording 3-D data with single-hole dynamite source points, and use omnidirectional arrays such as circles or clustered arrays. Field tests for linear versus podded geophone arrays have indicated that there may be situations where arrays can attenuate reflections more than groundroll does (Wittick, 1998). The tests included inline as well as cross-line recording, i.e., the wave-fields hit the geophone groups linearly as well as broadside. Arrays in the source line direction should complement arrays in the receiver line direction as much as possible (Vermeer, 1998a). Figure 1-4 shows the different array combinations that may be considered. Figure 1-4a is a single source point being recorded by a linear receiver array, the series of midpoints indicates possible smearing in the receiver line direction. Figure 1-4b shows the areal midpoint sampling that might happen with a combination of linear source and receiver arrays. This areal midpoint coverage is summed as one trace. Figure 1-4c indicates the midpoint coverage that can be achieved with a single source and an areal receiver array. The greater the number of elements in a receiver array, the larger is the areal midpoint coverage.

9

Figure 1-4: Midpoint coverage for: a. a single source into a linear geophone array (after Vermeer, 1998a), b. a linear source array into a linear geophone array, and c. a single source into an areal geophone array.

As we can see in figure 1-4, using distributed shots (shot arrays) can be a solution to reduce back-scattered noise.

1.3.3. Statistics A major goal of reflection processing is to provide reflectivity images of the correct reflector geometry. This goal can be prevented by the statics problem. The statics problem is defined to be static time shifts introduced into the traces by, e.g. nearsurface velocity anomalies and/or topography. These time shifts distort the true geometry of deep reflectors. Uneven topography can produce moveout delays in the shot gathers that can also be interpreted as undulations in the reflector, even if the reflector is flat. After an elevation statics correction (i.e., time shifts applied to traces) the data appear to have been collected on a flat datum plane. Static shifts introduced by topographic variations fall under the class of field statics, and those due to nearsurface lithological variations that occur within a cable length fall under the class of residual statics. Correcting for static shifts in the traces can make a significant difference in the quality of a migrated or stacked image. It is easy to determine elevation static corrections, but not so easy to find the residual static corrections. One means is to determine the near-surface velocity distribution by refraction tomography.

1.4. Forward Modeling By properly modeling the subsurface, geophysicists can:

•

Determine optimum acquisition parameters 10

•

Select proper processing techniques

•

Produce superior seismic results

Relatively simple geometries are effective on flat layers and gently dipping surfaces. In the presence of complex structures wavefield behavior becomes difficult to analyze. Modeling would seem to be a way to improve the survey design. Model based design using ray tracing technique can be a solution for analyzing the problems of conventional designs and finding the optimum parameters of the survey. The specific sequence for using a model in survey design would be: 2- Build a model of significant horizons using geological and previous seismic data and filling in the model by appropriate petrophysical and seismological properties. In this step we can examine several models. 3- Ray trace the model at an arbitrary source spacing and determine the optimum locations of receivers. 4- Shooting the obtained survey and extract useful attributes 5- Analyzing the attributes to find the best solution for proper imaging. The goal of this project is to utilize the ray tracing technique to find a method for having more energy at receiving stations. Several attributes are candidates to be analyzed for this purpose.

11

CHAPTER 2 REVIEW OF SOME PREVIOUS SEISMIC DATA ACQUISITION JOBS IN DEZFUL EMBAYMENT AND GEOLOGICAL STUDIES This chapter is subdivided into two main parts: 1- Geological review of south west regions of Iran with special attention to problematic layers and formations from geophysical point of view. As it is mentioned earlier high attenuation lithologies such as salt and anhydrite due to some intrinsic properties cause to energy is not recorded at the surface with acceptable energy content. Hence, as one of the most developed formations which plays the role of cap rock for plenty of reservoirs in Dezful Embayment is Gachsaran formation, this will be discussed in more details. 2- Review of some selected seismic data acquisition jobs which have done in this region, and discussion about the results and parameters.

2.1- Geological study of Dezful Embayment The Zagros basin is defined by a 7–14 km thick succession of cover sediments deposited over an extraordinary wide and long region along the north–northeast edge of the Arabian plate (Fig. 2.1; Bahroudi, 2004; Sepehr, 2003). Since the end of Precambrian times, this region has evolved through a number of different tectonic settings. Like other fold–thrust belts, it shows intense shortening close to the suture which becomes less intense toward the foreland. It is dominated by NW–SE trending

12

folds and thrusts, which swing into a NE orientation as they approach the Strait of Hormuz. The thick sedimentary sequences of the Zagros basin contain rocks ranging in age from Cambrian to Recent. The Zagros basin was part of the stable super continent of Gowndwana in Paleozoic times. The geological evidence suggests that the region was a passive continental margin in the Cambrian-Carboniferous and subsequently underwent rifting in the Permo-Trias and collision in the Late Tertiary (Berberian & King, 1981; Beydoun, Hughes Clarke, & Stoneley, 1992; Sepehr, 2003; Stocklin, 1974).

Figure 2- 1: Tectonic map of Arabian and Iran plates and main structural subdivisions of Zagros orogenic system (Abdolahiefard, 2006).

13

Sherkati et al. (2004) suggest a Middle to Post Miocene shift of sedimentary depocenter to the southwest allowed rapid subsidence and thick accumulation of the Fars group in the Dezful Embayment. These consist of the Lower Fars Gachsaran, the Middle Fars Mishan, and the Upper Fars Aghajari Formations. Each of these was restricted to a new foreland basin in front of the evolving Zagros Fold-Thrust Belt (ZFTB) (Bahroudi, 2004). Near the end of early Miocene, the evaporatic environment

in the Zagros basin led to sedimentation of Mid Miocene Gachsaran Formation. In Farsi, “Gach” covers anhydrite, gypsum and salt and, although the Gachsaran Formation contains mainly of evaporites, it also contains marls, limestones, and shales. Salt of the Gachsaran Formation was considered as the thickest basin-centre facies (Falcon, 1958; Gill & Ala, 1972; Kashfi, 1980). Gachsaran Formation is one of the most important seals for hydrocarbon reservoirs (especially in the Asmari Limestone) in one of the richest habitats for hydrocarbon in the world (Beydoun, 1991; Beydoun, Hughes, & Stoneley, 1992; Colman-Sadd, 1978; Dunnington, 1968; Edgell, 1996; James & Wynd, 1965; Motiei, 1993; Murris, 1980; O’Brien, 1957; Stöcklin, 1968). Gachsaran formation is particularly deposited in the area of Izeh, Lurestan, and Dezful Embayment. The Early Cambrian Hormuz salt located at the base sedimentary column and the Gachsaran Formation higher in the Zagros stratigraphic column as thickest and most wide-spread evaporatic units have acted as important detachment horizons in the southwestern region of the ZFTB (Sherkati et al., 2005). O'Brien (1950) introduced diapirism in the Gachsaran Formation to explain decoupling between competent and incompetent units. In addition, the divergence within the Gachsaran Formation used to be interpreted by thrusting of the upper members of the Gachsaran Formation over the incompetent members.

2.1.1. Stratigraphy of the Gachsaran Formation Detailed stratigraphic studies of the Gachsaran Formation report rhythmic bedding (Dunnington, 1968; James & Wynd, 1965; Stöcklin, 1968) with 22 distinct cycles in the Pusht-e-Kuh (Lorestan) province (Gill & Ala, 1972) and 14 cycles in the Kirkuk embayment (see Dunnington, 1968). Each rhythmic unit consists of bluishgreen marl, overlain in turn by limestone, dolomite and anhydrite, and with/without bedded salt (Dunnington, 1968; Gill & Ala, 1970; Motiei, 1993; Slinger & Crichton, 1959). 14

Figure 2- 2: Schematic Profile of subsurface layers in Dezful Embayment, Izeh and High Zagros zones. Gachsaran Formation is shown by yellow color (Sherkati et al. 2004).

15

Figure 2- 3: Stratigraphy column of Cenozoic Era in Lurestan, Dezful Embayment and Fars area.

The geologists of Iranian Oil Operating Companies (I.O.O.C.) for the first time reduced the order of the complexities of the Gachsaran Formation in 1946 (Setudehnia, 1977). James and Wynd (1965) modified the subdivisions of former workers in seven members; however, a distinct location as type section was not introduced for the Gachsaran Formation. The stratigraphic description of the Gachsaran Formation is based on results of the Gachsaran Oilfield wells in the Dezful Embayment. Member 1, also called Cap Rock, consists of interbedded anhydrite and limestone associated with bituminous shale. This member lies conformably on the Asmari Formation (Setudehnia, 1977) and acts as an important seal for the Asmari reservoir. Member 2 is mainly composed of salt with some anhydrate and limestone intercalations. Member 3 consists of thick anhydrate with subordinate salt. Members 4-6 are composed of mainly anhydrate with intercalation of marl, salt, and limestone. Thick salt beds are the main constituents of the Member 4, while, exposure of the salt beds at surface is rare (Motiei, 1993). Member 7 is overlain conformably by the Mid Miocene Mishan Formation. This member consists of alternating anhydrate, marl and limestone (Setudehnia, 1977). Previously, Mbr7 was considered as Stage 3, Mbr 6 as Stage 2 and Mbrs 1-5 as Stage 1 (Motiei, 2003).

16

Figure 2- 4: Location of some Iranian oil field in Dezful Embayment and Abadan plane

In geophysical point of view, Gachsaran usually considered as geological unit, mainly composed of 4 main lithologies. (Table 2-1) Table 2-1: List of main lithologies of Gachsaran Formation in geophysical ponit of view (Tamimi.N, 2007). Lithology

Abbr.

Constituent

Anhydrite

An

Anhydrite + Gypsum

Marl

Ml

Gray Marl + Red Marl

Claystone

Clst

Red Claystone + Gray Claystone + Shale

Salt

Salt

Salt

Limestone

Lst

Limestone

2.1.2- Structural Geology of Gachsaran Formation Because of tectonic mobility of members 2 to 5, they act as an incompetent unit with disharmonic folding between members 1 and 6.

17

The maximum expected thickness of the Gachsaran Formation is 1800 meters based on drilling in the Gachsaran Field. Whereas, drilling of some other fields showed that the thickness of Gachsaran Formation exceeds 3000 meters. For example, in the well Zeloi 5 the Gachsaran Formation outcropped at surface and its measured thickness is 3447 meters. This abnormal increase in thickness is supposed to be related to tectonism as the thickness of the Gachsaran Formation usually varies over short distances (Bahroudi and Koyi, 2004). Bahroudi and Koyi (2004) consider salt of the Gachsaran Formation as the thickest basin-centre facies, whereas, salt deposited throughout the local basins (with minimum thickness variations) and with other incompetent sediments are displaced by tectonic activities in later stages. Layer parallel thrusts were forced to propagate as ramps above the anticlines. Of significance is that; the thickness of the Gachsaran Formation steadily reduces towards the SW and it changes laterally to a more competent unit in the Abadan Plain. Therefore, the Gachsaran Formation is not acing as a well-developed detachment horizon in the Abadan Plain. However, few observations that confirm thrusting in this otherwise regional detachment level may also be due to lower strain in the frontal ZFTB compared to more extensive shortening in the Dezful Embayment.

2.1.2. Review of seismic challenges of Gachsaran formation Tamimi. N, 2007, reviewed the velocity challenges and their causes in Gachsaran formation. The preceding explanations are based on his researches and analysis. The natural complexities of petroleum reservoirs systems continue to provide a challenge to geoscientists. Reliance on classic methods and current technologies without adequate understanding about problems is becoming less satisfactory. Understanding and characterization of challenges is an inevitable solution to decrease the risk of exploration in new prospects. Seismic imaging as the most important geophysical method is commonly used to evaluate subsurface. Different seismic methods aim at evaluation and subsurface from two points of view: Geological Structure and Geological Stratigraphy. In order to achieve above-mentioned goals, quality of seismic image plays an important role.

18

Obviously, quality of acquisition and processing directly affect the quality of seismic data and consequently the quality of interpretation. Geophysical studies of overlying layers are as much important as reservoir intervals. Plastic and incompetent formations (e.g. evaporites) are the major group of overlying cap rocks. Plastic formations quiet often make serious challenges in exploration of new prospects due to high structural deformation, abrupt lateral velocity variation, non-ideal elastic properties and high wave energy attenuation. Gachsaran Formation, as an overlying formation of most important oilfields in Dezful and Kirkuk embayments, is selected for this study. Gachsaran Fm. shows all the above-mentioned limitations and challenges related to plastic and incompetent formations. The tectonic incompetency of the Gachsaran Formation is a remarkable indication for a seismic interpreter in the study area. In some cases, which relevant well data such as geological markers and well velocities are note available, the Cap Rock (Gachsaran Member 1) is detectable based on divergency of reflectors (Abdolahie Fard, 2006). The top Cap Rock is usually conformable with the lower reflectors, whereas, the incompetent members show chaotic reflector pattern. The chaotic pattern is seen in the southern flank of the main anticlines, where incompetent units are moved from the crest of underling anticlines and accumulated there. The incompetent members of the Gachsaran Formation are mainly combination of anhydrate, shale, and salt. The dominant rock governs the velocity of the whole swelling part. Thickness of salt within the Gachsaran Formation increases in the SW of the Mountain Front Fault, whereas shale content of the Gachsaran Formation is more than salt in SW of the Dezful Embayment. In area where quantity of anhydrate or salt is high, then, high interval velocity is expected for the Gachsaran Formation. On the contrary, dominancy of shale leads to lower interval velocity. Therefore, high velocity anomalies within the Gachsaran Formation cause velocity pull-up effects on the underling reflectors in the seismic sections processed in time domain. In contrast, lower interval velocities show opposite results. Both cases affect the seismic quality of surrounding area especially the southern flank of the main anticlinal structure. Usually conventional time-migrated seismic sections are distorted and obscure due to accumulation of incompetent units of the Gachsaran Formation. Lateral velocity variations lead to complicated seismic ray-path and such seismic imaging is proved to be more challenging as precise velocity analysis and advanced seismic processing 19

techniques are needed. For example, pre-stack seismic migration is required to image and correctly position the steep flanks of the underling anticlines. This technique is considered to be the appropriate method for imaging targets in the presence of overburden, especially in the presence of a salt body (Oezsen, 2004). However, the quality of seismic data depends to the acquisition parameters. Due to structural complexity, it is necessary to use larger offsets (distance between seismic source and farthest receiver) to reach a correct illumination of the whole salt mass (Mougenot and Al-Shakhis, 1999). In addition, 2D seismic technique is unable to image properly anomalies within the Gachsaran due to scattering ray paths and interferences with outof plane reflections. In such complicated case, 3D seismic migration is needed to avoid interferences of side effects and focus seismic energy. High fluctuations of layer thickness have been recognized in geological sections and seismic profiles especially in northern and eastern area of Dezful Embayment. In some cases, the variation of Gachsaran thickness reaches to thousands of meters just for 1km away. The high variations in thickness effectively cause complexity of study of wave propagation through the medium. Moreover, steep reflectors do not reflect wave significantly and it makes some problems in acquisition and processing steps. It resembles to the case of salt flanks and maybe has the same solution. Also, the existences of thrust faults in Gachsaran exert abnormal lithostatic pressure on shales and because of low permeability; it causes high pore pressure in marly and shaly layers in Gachsaran. High pore pressure decreases p-wave velocity significantly and add more geophysical complexity to Gachsaran. Such complexities disable us to introduce a confident velocity for Gachsaran Formation. As mentioned previously, Gachsaran was precipitated where affected by high tectonic activity. Folding and incompetent layer causes high lateral lithology variation. Ideally, high lateral lithology variation leads to high lateral velocity variation. In complex situations, high lateral velocity variations cause non-hyperbolic components to be added to the CMP gathers (Alaei, 2005). In order to get the correct laterally varying velocity field, analysis in the pre-stack domain can resolve the non-hyperbolic moveout component. A variety of techniques have been proposed for using pre-stack reflection multi-coverage surface seismic data to perform velocity analysis in complex areas with lateral velocity variation (Bishop et al., 1985; Faye and Jeannot, 1986; Etgen, 1990; Van Trier, 1990; Stork and Clayton, 1991; Stork, 1992; Kosloff et al.,

20

1996). High velocity variation not only reduces quality of depth-migrated image but also the reliability of depth-converted seismic image. Beside high velocity variation, the most noticeable fact about seismic of Gachsaran which directly affect quality of acquisition of seismic data is high attenuation in incompetent layers. Such phenomena highly attenuate reflected waves and amplitude of reflection in Gachsaran is very low. Weak reflection can not be used as good criteria in making velocity model from seismic reflection. Also, interpretation of reflectors in and below Gachsaran is very difficult due to this high attenuation. Sometimes in thick Gachsaran, the problem becomes critical. It has been guessed that main lithology which causes high attenuation is Salt. Thick layers of salt in Eastern and Northern areas of Dezful Embayment made too many problems for interpreters, processors and field men. In addition, salt exhibits plastic deformation and creeping. Plastic deformation characteristic make nonhyperbolic reflection where thick layers of salt exit. Interpretation and making velocity model from velocity spectrum in non-hyperbolic condition is too hard and erroneous. Marl is another lithology in Gachsaran and shows low seismic velocity. Thin to medium beds of marl in anhydrite and salt layers with high velocity add more difficulty to wave propagation and velocity model. In geological term, marl is made of limestone and shale which are different in velocity (limestone characterized by high velocity and shale exhibits low velocity medium). Thus seismic characteristic of marl depends on its constituent and we can not consider marl as unique lithology.

2.2. Review of some previous seismic data acquisition jobs in Dezful Embayment Several seismic data acquisition surveys have been performed in various fields of Dezful Embayment. Some of them which are more problematic are selected and their acquisition parameters and conclusions are analyzed.

21

2.2.1. Seismic data acquisition on Gachsaran Formation (1971-1972)

2.2.1.1. Introduction Reflection profiles shot in this Area over surface outcrops of the Gachsaran Formation had not often produced usable seismic sections, whatever field techniques were used. Geologically the Gachsaran Formation of the Miocene Fars Group in the “oilfields area” is described by James and Wynd (Stratigraphic Nomenclature of Iranian Oil Consortium Agreement Area, 1965) as having “a mobile and saliferous nature” and containing Anhydrites, Salt, Marl and thinly bedded limestones. Geophysically these rock types can be expected to transmit energy at high velocity and being incompetent may give great dispersion and refraction off strata which dips steeply or is overturned, with faulting adding further complications. Consequently where Gachsaran is thick, energy absorption and dispersion can be expected, preventing deeper reflections from being recorded at the surface. The results and analysis of previous jobs showed that Surface noises to be in the range of 40 meters to 90 meters and wave-length with sixty meters predominating. Noise creates bad effects particularly near steeply dipping surface beds. Large shot patterns and possibly smaller geophone station intervals gave better results, but as single cover and six-fold stacking usually gave poor or no results, it was difficult to make comparisons. Sometimes six fold cover shows some improvement over single cover and the one twelve-fold stacked line was a little better. Surface patterns proved best when used with an offset of about 1200 meters using a star pattern of large area with 150 lbs. charge divided into separate 2.5 lb charges; there seemed to be no advantage when reducing the sizes of individual charges below 2.5 lbs. Geophone station intervals of 50 m. or 60 m. were normally used, though 25 m. gave better definition on some lines, 24 to 48 geophones per trace, 8 meters or more apart with 30 meters between strings looked better than patterns with smaller area, Of the various recording filters used the high cut filter did not seem critical so long as it was above 50 Hz., but on the low cut side, since some reflections of interest have a 10 to 16 Hz. content, a setting of 8 Hz. seemed suitable.

22

2.2.1.2. The 1971 Experiment This series of experiments was designed to analyses the coherent noise, compare different field parameters by altering only one at a time and using these results, to design suitable field configurations for shooting the line a number of times. The noise test analyses showed coherent noise to be almost entirely in-line and not as strong as might be expected. Response curves for the shot and geophone patterns showed good theoretical attenuation and it was concluded that coherent noise was not responsible for the poor results usually obtained when shooting on Gachsaran outcrops. Comparisons of charge size and patterns were made using the field monitors and the in-line patterns were seen to be better than the “star patterns. Large charges in-line seemed to be best and offset comparisons were made using these parameters and 1225 meter for offset was selected as giving the best overall data. The processed section suggested, however that the larger in-line charge sizes were very much better and that a longer offset might also have given improvement. A number of field filter tests suggested an 8 to124 Hz range to be reasonable and the processed data confirmed this to be so. From these experiments it was decided to shoot the line using 50 meter station intervals, then using 25 meter station intervals each with appropriate geophone patterns. Finally the better of these two techniques would be selected and the necessary additional shots made to convert the original sixfold cover into twelve-fold cover. From the field monitors the 25 meter interval was thought to be best and this technique was used to produce the twelve-fold stacked section. A general comparison of the six fold versions obtained so far showed that none of them had recorded any useable data under the thick Gachsaran. The use of deconvolution did not greatly improve the data, but did enable better application of residual statics by “sharpening” the data.

2.2.1.3. The 1972 Experiments The experiments showed that although coherent noise could be troublesome from some shot points it was not the main cause of the poor data usually obtained. Random noise combined with absorption of input energy by the highly folded and faulted

23

rocks together with rough topography would appear to be responsible for lack of results when shooting on Gachsaran. Useful reflections were obtained on some parts of one of the lines where nothing had been seen previously and poor reflections had been greatly improved on this line. A number of factors combined to give this improvement some of which were in the field recording technique and some in the data processing. The field parameters which helped to give great improvement were the increase in area and number of elements of the multi-element shot and geophone patterns, the use of longer offsets combined with the “right” spread length and the use of end-on stacking, all these being combined and used to produced greater multiplicity of cover. In the processing of the data, residual statics particularly the method using a dip selective approach produced the greatest improvement, but the ability to apply time variant processes such as deconvolution, filtering and equalization were also useful.

2.2.2. Seismic Data Improvements Experiments (1972 and 1973)

2.2.2.1 Introduction During the seismic season September 1972 to May 1973 a number of field and processing experiments were carried out with the aim of improving the data quality of seismic reflection work. The experiments were aimed at following up the leads given by the 1971-1972 ’Gachsaran Experiments’ and areas with poor or bad data quality were chosen for the field work. In order to investigate charge parameters a number of comparisons were made in Deh-Dasht area including the use of tapered patterns. Further experiments were made by altering shot patterns during the shooting of twelve fold lines in order to obtain two sets of six fold sections for comparison. No stacked comparisons were made between different geophone patterns, since operationally it is easier to compare shot patterns and the results should be valid for the geophone patterns as well. Residual statistics were used on most lines where the signal to noise ratio was good enough to the operator be effective.

24

2.2.2.2. Improvements in poor data areas Although the Gachsaran Experiment had been aimed primarily at obtaining data under outcrops of Gachsaran, the parameters which contributed to the improvement also improved poor data where other formations outcrop. Usually lines shot on Bakhtyari Formation and alluvium outcrops obtain fair to good results, but this area is an exception. Increasing the fold of coverage and charge size while other parameters kept constant resulted in noticeable improvements in data quality.

2.2.2.3. Deh-Dasht Experiment The experiments in this area included an in-line noise spread together with two charge patterns shot into part of the noise spread. The noise test was analyzed in the field for frequency and wavelength and results were used to select charge patterns for the comparison. In this area the main noise had wavelength of between 140 meter and 220 meter and appear to be associated with the near surface refractions. The charge comparisons showed that near traces (1-24) were not contributing much, if any useful signal, but data quality was so poor that this type of experiment has limited value in such an area. It was concluded that where long offsets were needed to get better results, but at the same time fairly shallow data was required. Charge tests show that in this area using a specific amount of dynamite spread over a large area might produces better results. Mixing was also tried over poor data, but this did not help very much and could not be used when steep dips occurred.

2.2.3. Conclusions

Considerable improvement was achieved during the 1972-1973 season by increasing the fold of multi-cover together with the use of off-end shooting with longer offsets and larger shot and geophone patterns. In addition due to mountainous nature of this

25

area (ZFTB) residual static correction seems to be very crucial for obtaining better images.

26

CHAPTER 3

FORWARD MODELING AND RAY TRACING APPROACH

3.1. Introduction Forward seismic modeling is the process, through which a subsurface geologic (acoustic impedance) model, in one- two- or three-dimensions, is transformed into a synthetic seismic record of one, two or three dimensions. Vertical depths within the geologic model are converted to two-way transit time. Acoustic impedance contrasts within the geologic model are converted to reflection amplitudes. Often, the relationships between geologic models and corresponding synthetic seismic records can be readily deduced through visual examination. Synthetic seismic records are typically generated both before and after the acquisition of seismic field data. Synthetic seismic records generated before field acquisition are used to determine if an intended geologic target will generate an interpretable signature on output processed reflection seismic data. Synthetics can also be a valuable tool with respect to the design of an acquisition program. Synthetic records generated after acquisition and processing of seismic field data facilitate the interpretation of the processed field data. Synthetics make the confident correlation of the observed reflections and geologic interfaces possible, and verify that the seismic responses of interpreted

27

conceptual geologic models are consistent with the actual seismic data. They are essential interpretation tools.

3.2. Ray Tracing Approach In this study, NORSAR2D is used for forward modeling purposes via 2D ray tracing. In preceding part, the ray tracing parameters, requirements and theories (Advanced mathematical theories have been excluded) have been introduced for better understanding of this concept.

3.2.1. Fundamentals of Ray Tracing

3.2.1.1. Seismic Model Seismic model consists of three basic types of elements which are represented in figure 3-1.

Figure 3- 1: Seismic Model

1- Interfaces, representing the discontinuity surfaces between the rock properties.

28

2- Blocks (layers) defined by the volumes (areas) between the interfaces. 3- Material properties, assigned to each block. In this context, the material properties are seismic P- and S- velocities and densities, and optionally the P and S Q-factors characterizing the an-elastic damping of the medium. The interfaces, as well as the material properties, shall be represented by continuous and smooth functions. The only discontinuities in material properties occur when crossing an interface.

3.2.1.2. Ray Code Ray Codes are the instructions telling the ray tracing process what to do after a ray has hit an interface in the model. The ray may be reflected back to the incident block or transmitted into the adjacent block. The mode of the ray may also change from P to S (i.e., from longitudinal to transversal particle displacement, so-called mode conversion), or vice versa. A ray may initially start off as P or S.

3.2.1.3. Ray tracing with given initial point, ray direction, and ray code 3.2.1.3.1. Introduction Ray tracing is a process by which one can calculate important quantities tied to seismic wave propagation through a layered medium (i.e., a seismic model as described above). Ray tracing may be classified as an approximate solution of the seismic wave equation, valid for high frequencies. For a comprehensive description of the ray methods, is referred to Cerveny (1985). The ray theoretical aspects discussed here pertain to so-called geometrical rays, i.e. rays following the geometrical law of reflection or transmission (Snell’s law) at all interfaces. Another ray family are diffracted rays, which follow Keller’s law of edge diffraction at a pre-defined diffraction point, i.e. a point where two of the model interfaces intersect. Diffracted rays consist of three geometrical portions; a geometrical ray from the ray start point (shot) to the diffraction point, a geometrical ray from the diffraction point to the end point (receiver), and finally, a geometrical ray path describing a so29

called shadow boundary. The ray theory described below applies separately to each of the geometrical ray portions.

3.2.1.3.2. Geometrical ray tracing Here we will divide ray tracing in two basic parts, that is: • Kinematic ray tracing, • Dynamic ray tracing. Kinematic ray tracing calculates the location of the ray paths and the travel-time along these ray paths. For a number of applications in seismics only the travel-times are of interest. The kinematic ray tracer only needs the wave velocity in the model blocks, i.e., P-velocity and/or S-velocity depending on which wave modes are of interest. Dynamic ray tracing also calculates the dynamic properties of the seismic wave field, such as the geometrical spreading factors, wavefront curvature, and amplitude coefficients along the ray paths. In addition to the seismic P- and S- velocities, the dynamic ray tracer makes use of the density functions. By using a full dynamic ray tracing, one has enough information about the wave field to calculate synthetic seismograms in receivers located close to the rays.

3.2.1.3.3. Kinematic ray tracing When tracing a single ray into a model, the following elements must be given (Figure 3-2): • The initial point of the ray • The initial direction of the ray (ray tangent) • A ray code (i.e., rules determining which ray to follow next when hitting an interface (reflected or transmitted P or S)).

30

Figure 3- 2: Elements of kinematic ray tracing

The ray path and travel-times along a ray within a continuous block of the model are calculated by solving a system of differential equations (the kinematic ray tracing system). The velocity function in the block determines the ray behavior, for example, a constant velocity (homogenous layer) gives straight ray paths, whereas a linear velocity gives circular rays. When a ray is about to leave a block, the intersection point with the interface is calculated. At the interface, the ray bends (changes direction) according to the ray’s angle of incidence at the interface, and the velocity contrast across the interface. This change in ray direction is calculated by the wellknown Snell’s law, and the ray path can be calculated through the next block. Thus, the ray path is calculated successively for block after block by alternating use of a block tracer, an intersection algorithm, and Snell’s law. At each ray or interface intersection point the ray code is checked to determine what ray type to follow through the next block (reflected or transmitted, P or S). The ray code also contains a stop criterion for finishing the ray tracing, for example, when a certain travel-time is exceeded.

3.2.1.3.4. Dynamic ray tracing When the ray path and travel-times have been calculated by the kinematic ray tracer, one may optionally use a dynamic ray tracer to calculate dynamic ray quantities (ray attributes). This calculation is performed along the already found ray path, using a similar system of equations as in the kinematic ray tracer. In this approach of ray

31

tracing several factors such as geometrical spreading factors and wavefront conversion are included.

3.2.1.3.5. Ray tracing from shot to given receivers (two-point ray tracing) The ray tracing method described in the previous section is based on initial conditions making possible a “direct” calculation of the ray path and its associated attributes; when the start point and start direction of the ray are given. However, in many (most) cases we have a number of specified receiver positions in which we want the ray paths to end, that is, we want to find ray paths between given shots and receiver positions. Unfortunately, for a general structural model with curved interfaces and/or varying layer velocities, there is no a priory way of knowing what initial ray direction will make the ray arrive in a pre-specified receiver point. As a consequence of this, a search procedure must be used in order to calculate the proper ray path. This “trial and error” procedure can be designed in various ways - we use a technique referred to as the “shooting method”. The procedure consists of “shooting” out a fan of rays having a large coverage of ray directions, and on basis of those rays passing close to some specified receiver point, repeating the shooting for nearby directions until the ray hits within a certain capture radius of the receiver.

3.2.2. Ray Tracing Techniques

3.2.2.1. Normal Incident Ray Tracing In this technique, rays trace with Normal Incidence Paths with respect to selected interfaces. The travel-times may simulate interpretations of stacked un-migrated data. The survey is zero offset, i.e., for each shot there is one and only one associated receiver situated at the same position. The shot-receiver points, also termed SR points, are regularly distributed along a horizontal line. The SR points may be projected vertically to a reference interface. Rays are traced from the selected reflectors to all shot-receiver points.

32

3.2.2.2. Image Ray Tracing The travel-times from ray tracing in this module may simulate interpretations of migrated data. In the survey, rays are traced from start points regularly distributed along a horizontal line. The rays start off in a vertical direction. Two-way travel-times are calculated every time an interface is intersected. Rays are traced until they reach the model box or encounter supercritical reflections. The important criterion in this technique is Migration Aperture. As this technique simulates the post stack migrated section selection of this parameter plays very important role on quality of migrated resulting section. Migration aperture in ray modeling can be expressed as the search distance respect to initiation of image ray in which traced normal incident rays from each horizon are taken into account for migration. In fact, while an image ray traces the incident points of the ray-path and interfaces records in a dataset. Then after at each incident location Normal Incident ray will be traced and the resulting travel time doubles and records as the zero-offset migrated time of that location in Image ray source location. As it is depicted in figure 3-3 choosing insufficient migration aperture value may leads in lack of sampling in some locations, especially where steep flanks are present which cause to ray travels to long distance.

2.75 km

Normal Incident Rays

Image

Figure 3- 3: Depiction of image and normal incident rays

33

3.2.2.3. Common Shot Ray Tracing This technique simulates ordinary horizontal surveys. Shots are regularly distributed along a horizontal line. The receivers are specified as fixed or relative to the shots. This general definition is valid for all shots, hence each shot has a related sequence of receivers. The receivers all lie along a horizontal line, which may differ from the shot line. Both shots and receivers may be projected vertically to a reference interface. Rays are traced from the shots to the related receivers.

3.2.3. Limitations of ray tracing 3.2.3.1. Frequency Limitations As mentioned in the previous section the ray tracing technique is based on an approximation to the general wave equation - strictly valid for high frequency signals only. Mathematically, this means that the theory is valid for indefinitely high frequencies, or put in another way, for indefinitely short wavelengths. What does this really mean in practice? For real applications this makes restrictions with respect to the seismic model in which we want to do ray tracing. This means that the seismic wavelength L must be considerably shorter than the “length of the smallest details in the model”. In practical terms, L must be considerably smaller than quantities such as: • The radius of curvature of the interfaces • The length of the interfaces • The block size • Measures of the inhomogeneity of the material property in a block (i.e., expressions of the type

v ). | grad (v) |

It means that for example having a media with compressional velocity of 2500m/s, after ray tracing performed and pulses created, if we choose a wavelet of 20 Hz frequency, the seismic wavelength will be 125 meter. It means that in this analysis we will not be able to image features smaller than this value. In practice, the model dimensions must be at least four times larger than wavelength. A common way of

34

expressing the smoothness criterion in ray tracing is to say that you should prepare a macro model, i.e., a model containing only the large scale features (implicitly referring to the seismic wavelength dimension).

3.2.3.2. Lack of smoothness A very typical pitfall when building a seismic model for ray tracing is that the model interfaces are too detailed, that is, not smooth enough. Especially, when importing interfaces from seismic interpretation systems or digitizing tables, they tend to be much too oversampled, and to lack the necessary smoothness to directly fulfill the ray tracing criteria. In most cases the small insignificant “ripples” on the interfaces can be avoided. This problem causes to algorithm fall into mistake while estimating angle of reflection or transmission on an interface. This error affects on reflection and transmission coefficient values and overally will result in unexpected amplitudes in final section. Figure 3-4 shows two sections obtained by shooting on the same model with same parameters in one smoothing ignored and in another one model has been smoothed. Sudden variations of amplitude in un-smoothed model section are obviously seen in the left section.

Figure 3- 4: Un-smoothed model section (Left section) and smoothed model section (Right section)

35

3.3. Modeling

3.3.1. Introduction In this part of project I have tried to find some solutions for seismic data acquisition under complex structures while high velocity layers with high degree of attenuation present. First of all an approximate model of study area has been built. For this target several models and time interpretations have been reviewed and the common features extracted. In addition petrophysical properties have been modeled by statistical techniques. Relevant data from different wells which are drilled in this area have been utilized and in some cases core data and laboratory measurements have been applied. The main characteristics of this area (Dezful Embayment) are the presence of thrust systems, existence of plastic layers which are highly attenuative such as Gachsaran formation as the cap rock in many cases, existence of steep flanks and lateral property variations. One of the probable solutions for seismic data acquisition in such area is multicomponent data acquisition. There are several reasons to recommend this technique: 1- As it is mentioned in chapter 2, one of the most important problems in areas in which thick sediments of Gachsaran formation are present is the record of acceptable energy. Hence, as S-waves has lower velocity than P-waves, according to Snell's Law they will be reflected with lower angle, so they can be recorded in shorter offsets, this fact not only may cause to records of higher amplitudes, but also can provide us with some records of energy from reflection points whose P-reflections never could be recorded in the surface, because they travel to very long offsets which maybe far beyond the length of our spread. 2- S-waves are not affected by gas or other kind of fluids. Based on reservoir engineering, after production of oil for a long time and declining the pressure of the reservoir, gas expels from the oil and forms a gas cap at the crest of the reservoir. This gas affects the resolution of P-waves and may obscure the structural features and variations in the crest of anticlines in images resulted

36

from this kind of data (PP). But as S-waves are not affected with the gases, they can bring out better images from structural features (figure 3-5). In addition V p Vs ratio can be used as a reliable indicator of fluid content. 3- In the case of thrusting systems, in some situations P-reflections from underlying horizons hit the thrust with critical or over critical angles. In such circumstances the incident rays will be reflected down again and no energy can be recorded from underlying layers of the thrust, while S-converted waves are capable of passing through the thrust and bringing energy to receivers.

Figure 3-5: Comparison of PP and PS results. Left figure shows the result of a PP data acquisition, as it is obvious in this figure, gas cloud has obscured the crest of the reservoir. Right figure shows the result of PS data acquisition that has resolved the top of the reservoir and a largr fault in the flank[3].

3.3.2. Modeling Procedure The modeling has been performed in two steps: 1- Construction of a simple horizontal layers model to compare some attributes of PP and PS reflections. 2- Construction of a complex model based on common structural and petrophysical features of Dezful Embayment for analyses of acquisition parameters in such areas. In each step after the construction of the models, ray tracing is performed and from the results, some attributes have been extracted and analyzed.

37

3.3.3. Simple Model At the first step of this study a simple model has been utilized. The geometrical and petrophysical properties of this model is based on simplification and averaging of geometrical and properties of several well logs, core data and interpreted sections in Dezful Embayment. Heterogeneity and anisotropy have been neglected in this model, but plasticity which is one of the main properties of attenuating formations of this area has been considered. The model is depicted in figure 3-6.

Figure 3- 6: Three layer flat model

Properties used for this modeling are represented in table 3-1. Table 3-1- Block properties Thickness

V p (km/s)

Vs (km/s)

Density(gr/cm3)

Qp

Qs

(km) Block 1

2.500

5.3

2.944

2.19

22

18

Block 2

.491

5.85

3.900

2.68

35

30

Block 3

1.009

6.2

4.133

2.75

30

28

For the purpose of our analysis Shot Gather Ray Tracing has been used. One shot has been laid at x = 20 km and the receivers have been spread with split spread geometry with maximum offset of 10 km on each side. Figure 3-7 shows the overlaid PP and PS reflection of traced rays.

38

Figure 3- 7: Representation of traced rays

As it is illustrated in figure 3-7 the difference between the last CDP (common depth point) and last CCP (common conversion point) on the second reflector which are recorded by the same receivers is about 3km. This is another advantage of multicomponent data acquisition by which it is possible to record data from deeper and/or more distant targets using the same length of spread (Xmax). Recorded events from ray tracing which play the role of reflectivity series, can be easily converted to seismic amplitudes by convolving them into a wavelet. The 50 Hz Zero-Phase Ricker Wavelet has been utilized in this part. PP and PS resulted seismic sections after AGC correction have been separated and represented in figure 3-8. It is expectable to record converted waves with a time delay respect to P-waves. Moreover as it is seen in the figure there is no energy at zero offset for converted waves.

Figure 3- 8: PP (left) and PS (Right) sections

Having these reflection amplitudes, we are ready to extract some attributes for further analysis. The attribute which are used in this study are: 1- Incoming ray angle at first point of reflection. 2- Outgoing ray angle at first point of reflection. 39

3- An-elasticity Factor which is the quantity telling the loss of wave energy to heat. This factor is frequency dependent (see chapter one). If, t = *

s( X r )

ds , in which "s" represents the ray path from source to receiver v( s )Q( s ) s( X s )

∫

and Q is quality factor, then the frequency dependent amplitude factor is: *

A( f ) = e −π × f ×t , hence, the more the an-elasticity factor, the less the amplitude

factor. 4- Total Geometrical Spreading which is the measure of total effect of geometrical spreading on rays. The more this value, the lower the attenuation amount of the energy. 5- Amplitude Coefficients variations for P-P and PS reflected rays versus offset. In below figures, red dots represents the events belong to PS reflections and blue dots related to PP reflections. Figures 3-9 and 3-10 illustrate the incoming and upcoming angles in reflection points. These graphs are very useful when we are going to analyze reflections in angle domain. By the way, in this study we are just going to represent the reasons which cause to receiving converted waves have more energy respect to corresponding Pwaves. Lower reflection angles which are subsequent of lower converted wave velocities, confirm existence of smaller ray-path for these rays which will result in lower geometrical spreading and absorption attenuation in a plastic media. An-elasticity factor is a property which is not in favor of PS reflection in comparison with PP ones. In most of cases quality factor of media for P-wave is higher than corresponding S-wave, by the way according to Waters (1981) and other experimental practices,

Qs V V ratio is proportional to [ s ]2 ratio, and as this ratio ( s ) tends to unity Vp Qp Vp

versus depth,

Qs also tends to unity at deep targets. In this model, as illustrated in Qp

table 3-1, because of plastic nature of Gachsaran formation quality factors are higher for P-waves.

40

Figure 3- 9: Incoming reflection angles of PP rays (blue) and PS rays (red)

Figure 3- 10: Up-going reflection angles of PP rays (blue) and PS rays (red)

Figure 3-11 shows the variation of An-elasticity factor versus offset.

41

Figure 3- 11: Variation of An-elasticity factor versus offset (PP rays (blue) and PS rays (red))

The next parameter which plays very important role in amplitude variations is Geometrical Spreading. As mentioned earlier, S-waves undergo lower geometrical spreading attenuation effect as a result of their lower velocity. For having more precise comparison of PP and PS result it is necessary to compare the recorded amplitudes from the same reflection point, not from the same offset. As it is depicted in figure 3-7, the last receiver records data corresponding to P and S reflections from two different points that are about 3km far from each other. Obviously the traveled path will be much more for that S-ray. Figure 3-12 shows the variation of geometrical factor versus offset from the first horizon. Although the above mentioned criterion has not been included, however Swaves show better result. According to the figure, as the offset increases, the effect of geometrical factor decreases (because of inverse proportionality of energy and distance). It is clear that all over the 10 km offset, the geometrical factor is higher for S-wave. The last comparison is AVO comparison. Figure 3-13 shows the variation of amplitude versus offset for both PP and PS reflections.

42

Figure 3- 12: Geometrical Spreading versus offset (PP rays (blue) and PS rays (red))

In the most of offsets the recorded amplitude of S-reflections are higher than its corresponding P-reflections. The novelty of this figure is that several factors including absorption, an-elasticity, reflection coefficients and geometrical effects have been considered for computations.

Figure 3- 13: AVO analysis (PP rays (blue) and PS rays (red))

According to the first step of modeling in this study it can be concluded that multicomponent data acquisition technique, not only can result in better images of subsurface features in areas with high attenuating layers, but also this technique will

43

provide us with some auxiliary data (such as

Vp Vs

ratio and S-wave polarizations)

which can be used for pore fluid identification, fracture detection and anisotropy modeling.

3.3.4. Complex Model (Based on a real field data in Dezful Embayment) 3.3.4.1. Introduction Structural interpretations of reflection seismic data in foreland fold belts are dependent on structural styles (Bally, 1983; Mitra and Fisher, 1992). The quality of Reflection seismic profile from the Zagros foreland fold thrust belt of southwest Iran, like other mountain belts (e.g. Canadian Rockies) are typically not good enough to draw the complete configuration of thrusted folds especially in deeper levels. The reliance on structural styles is necessary in view of the fact that the quality of the seismic images of fold and thrust structures is poor. Parts of the stratal geometry maybe clearly shown while other parts show either a lack or a confusing overabundance of reflection signals (Lingrey, 1991). Structural problems in such areas generate some seismic signatures that are not taken into account in processing. Ray tracing techniques enable us to find better understanding of wave propagation and ray paths in underground horizons. This matter can help us in setting appropriate acquisition and processing parameters. Moreover, some seismic signatures such as diffraction patterns can guide the interpreter to find the more exact locations of thrust faults and other ambiguous seismic features. In the second step of modeling, different method of ray tracing techniques are applied to investigate the problems of imaging under complex structures and thrusting systems. Three different techniques of ray tracing have been used for specific purposes. The techniques are: 1- Normal Incident Ray Tracing Technique 2- Image Ray Tracing 3- Common Shot Ray Tracing 44

3.3.4.2. Model Building The first model of this part has been built based on some depth migrated sections and the properties obtained from well logs and distributed base on geostatistical methods. This model shows three anticlines which are related together by two thrust faults (figure 3-14).

Gypsum and Anhydrite

Figure 3- 14: Model Blocks

As it is illustrated in the figure overlying layer mainly contain gypsum and anhydrite. Due to the plastic and incompetent nature of such lithologies they can flow during folding and their structure does not obey from underlying layers. In addition lateral velocity changes cause to ray path bending. To respect these properties the overlying layer which plays the role of cap rock has been subdivided into some narrower members with different petrophysical properties. Figure 3-15 shows the resulting model.

45

Figure 3- 15: Main model

3.3.4.3. Ray Tracing As it stated before, three different techniques of ray tracing are utilized for analysis in this part. The first approach is Normal Incident Ray Tracing. In figure 3-16, 251 sourcereceiver pairs are spread on the surface from the location of 2.5 km in model reference coordinates until 15 km. Normal rays traced from each horizon and after some basic corrections the following section has been obtained. The red rectangle in figure 3-16 may be interpreted as diffraction pattern, but the fact is that it is nothing but a signature that is as a consequence of existence of steep flanks. Figure 3-17 shows how the normal incident rays go travel far from its reflecting point. Red arrow shows a possible simple path for traveling of the ray from location of 10 km and recording at some locations around 6 km after about 3 seconds.

46

Figure 3- 16: Results of normal inident ray tracing on the model illurtrated in figure (3-15) with Station spacing = 250 m

Figure 3-18 shows the result of denser sampling (station interval equals to 50 m). It is obvious that no improvement has been met.

Figure 3- 17: Reflection from steep flanks

In such conditions increasing the number of channels does not make any difference, because this energy must be focused just by migration operator. A probable solution in such cases is the use of multi-component data acquisition, because reflected Swaves can be recorded at much closer distances.

47

Figure 3- 18: Results of normal inident ray tracing on the model illurtrated in figure (3-15) with station spacing = 50 m

To correct the above problems especially for better positioning of thrust faults, Image Ray Tracing approach has been utilized. In the first attempt a dense sampling spread over 12.5 km (from 2.5 to 15 km) has been chosen. The migration aperture for migration of normal incident rays is equal to 1 km. As it is shown in figure 3-19 neither the left flank of the structures nor the thrust planes are been recorded. In such cases may seem that PSDM in the only solution of imaging of such complex structures, but increasing the migration aperture which implies the increase of spread length can be a solution to this problem.

48

A

B

Figure 3- 19: Results of image ray tracing on the model illurtrated in figure (3-15) with migration aperture = 1km

Figure 3-20 shows the same survey except for migration aperture value which has been increased to 4 km. There are some improvements in left structure (A) and moreover the left thrust fault is completely imaged. For the right structure (B), there are some improvements in illumination of the top of the reservoir, but the relating thrust can not be seen at all. This is due to high degree of steepness of this flank. In such cases reflecting P-waves achieve critical angle while hitting the flank of such over-turned structures and can not refract into the structure and arrive to surface. But S-waves which travel with lower velocities will hit the structure with lower angle, and easily bring the energy to the surface. In the next survey (Figure 3-21), migration aperture has been increased to 8 km, and it has resulted in better illumination of overturned flak layers which end on thrust plane. We can conclude that in thrusting systems long offset survey result in higher resolution images, but when attenuating layers such as salt and anhydrite present, P-S data acquisition is preferred (see 3.3.3)

49

A

B

Figure 3- 20) Results of image ray tracing on the model illurtrated in figure (3-15) with migration aperture = 4km

B

Figure 3- 21: Results of image ray tracing on the model illurtrated in figure (3-15) with migration aperture = 8km

The last ray-tracing mode applied, is the shot gather offset modelling. This resembles the real acquisition. Shot gathers, or common shot gathers are traces recorded at a 50

series of receivers propagating from a single source. This is the way data are collected and its results could be useful in evaluation and design of seismic surveys (Fagin, 1991). At this stage, the focus is more on the events and not on the image from subsurface. One of the important applications of ray-tracing methods is the modelling and identification of specific events on seismic records (Carcione et al., 2002). At first a survey with following parameters has been designed. Table-3.2 Spread Type First shot coordinate Shot interval Number of shots Max Offset Receiver Interval Number of channels per shot Maximum Fold

Split Spread 5 km 0.25 km 25 6 km 0.25 km 48 24

A section from shot at location of 10 km is illustrated at figure 3-22.

Figure 3- 22: shot gather at X=10 km

The red arrow in figure 3-22 shows the effect of fault bend folding. In such cases CMP sorted data do not show the hyperbolic behavior. To show this event, the

51

reflections from top of the first reservoir (the first layer underlying the gypsum and anhydrite in figure 3-14) has been selected and illustrated in figure 3-24.

Figure 3- 23: Receiver locations versus mid-point, red crosses show the selected

Figure 3- 24: Non- hyperbolic behavior at far offsets for CMP location of X=8 km

The non-hyperbolic-behavior in this section can be regarded to the complexity of structures. In such situations higher order velocity analysis is recommended, because simple second order equation can not take this behavior into account. In addition results of off-end spread survey confirms the effect of offset in recording the energy from steep flanks, so as single-ended spread creates longer offsets, in such

52

areas this type of spread is recommended. Figure 3-25 shows that most of the energy from the crest of the reservoir is recorded at offset of 10 km.

Figure 3- 25: Off-end survey. Snapshot from the source-receiver pair where source located at X=5 km and receiver at X= 15 km.

Another recommendation which takes more time and expense is using pre-stack depth migration. Comparing all analysis in this chapter, it can be concluded that P-S data acquisition is a volunteer for imaging of complex structures especially in reservoirs having gas caps like Iranian south west hydrocarbon prospects. As P-S data acquisition is recommended for this region, the basic part of such surveying which is positioning of common conversion points have been discussed in the next chapter as the start point of 3 component data acquisition.

53

CHAPTER 4 COMMON CONVERSION POINT MAPPING FOR 2D-3C SEISMIC DATA ACQUISITION SURVEY DESIGN IN ANISOTROPIC MEDIA

4.1. Introduction Converted waves are being more and more widely used in seismic exploration, for they can provide high-quality images where conventional images are poor (Stewart et al., 1999). It has been proven that the location of the conversion point is critical to P-S seismic survey design (Lawton and Cary, 2000). Also the P-S conversion point location is of great importance to P-S data binning and Common Conversion Point (CCP) stacking and NMO corrections. In isotropic media, the location of the P-S conversion point can be calculated by using the Snell’s law. It then can be simplified by using asymptotic binning (Tessmer and Behle, 1988). However, it is well known that the earth’s crust is inhomogeneous and anisotropic. Almost all conventional seismic surveys and processing procedures have not taken anisotropy into consideration. This isotropy assumption could lead to large errors in NMO correction, stacking and migration, which have been shown by many authors (e.g. Alkhalifah and Larner, 1994). Anisotropy has less influence on P-waves than on S- (or converted) waves. In an anisotropic case, due to the difference between the group velocities and phase

54

velocities and the difference between the ray angles and phase angles, there could be a difference in the P-S conversion point location between anisotropic and isotropic cases. The most commonly considered anisotropic medium is the Vertical Transverse Isotropy (VTI) case. When the formation is shown to be VTI, the location of the conversion point will depart horizontally from that in an isotropic case, and it cannot be approximated by the isotropic case. Therefore, the horizontal position of the conversion point on the reflector has to be calculated specifically for VTI media. However, this requires knowledge of VP, VS, and anisotropy parameters. Otherwise, it will lead to problems in NMO correction and stacking which could cause traces that do not contain reflector energy to be summed, and those that do contain reflector energy not to be summed (Tessmer et al., 1990).

4.2 Review of seismic anisotropy The subject of anisotropy has a long history but now has come to be a central feature of geophysics to explore hydrocarbons (Thomsen 2002). "It can be said that any process which involves the concept of scalar velocity field is subject to error if, in fact, the actual velocity is a vector whose magnitude depends upon its direction."(Tsvankin, 1999). By definition, anisotropy means the variation of a physical property depending on the direction in which it is measured. In the field of our interest, seismic anisotropy is defined to be the dependence of seismic velocity on either the direction of travel or the direction of polarization (Sheriff, 2002). Seismic velocity here includes not only the P-waves, but also the other waves (e.g. S-waves). Accounting for the effects of anisotropy of the earth begins with understanding the geologic foundations of anisotropy. Anisotropy in sedimentary rocks may be caused by various factors (Thomsen 1986 and 2002): (a) Thin isotropic or anisotropic layers (b) Preferred orientation of anisotropic mineral grains (c) Preferred orientation of cracks and flat pores (sandstones) These attributes can exist in the rocks independently or in combination. Carbonates can display anisotropy due to the oriented distribution of cracks.

55

4.2.1 Anisotropy versus heterogeneity Anisotropy is easily to be confused with heterogeneity by many people since both terms are used to characterize velocity variation within rock formations. The distinction between these two terms has to be made clear. Heterogeneity is defined as the variation of a certain physical property from point to point. However, the anisotropy means the directional variation of a certain physical property at one point. The figure 4-1 shows the distinction between heterogeneity from anisotropy. However, heterogeneity may be related to anisotropy by using different scales. For example, a rock sample, sandstone or shale, is heterogeneous at the grain scale. The grains can have an order, a fabric or a texture, which makes the rock block appears to be anisotropic. This order may be caused by sedimentation. It is this ordered heterogeneity on the small scale, which causes seismic anisotropy on the large scale. Anisotropy and heterogeneity may coexist. There are four possible conditions: isotropic

and

homogeneous,

isotropic

and

heterogeneous,

anisotropic

and

homogeneous, and anisotropic heterogeneous. The real earth is known as being both heterogeneous and anisotropic.

Heterogeneous and Isotropic

Homogenous and Isotropic

Homogenous and Anisotropic

Heterogeneous and Anisotropic

Figure 4- 1: Anisotropy versus Heterogeneity

56

The standard of large or small scale means whether the seismic wavelength is larger or smaller than the heterogeneity scale. When the seismic wavelength is large compared to the size of the heterogeneous object, the object can be considered as a homogeneous, anisotropic one (Backus, 1962). In this case, the waves obey the laws of anisotropic wave propagation (Thomsen, 2002).

4.2.2 VTI (Polar Anisotropy) In seismic exploration, the most commonly considered type of anisotropy is polar anisotropy, which is also called Transverse Isotropy (TI). It has a symmetry axis and typically the axis is perpendicular to bedding. The velocities on the plane normal to this axis are identical. When this symmetric axis is vertical, the medium is called Vertical Transversely Isotropic (VTI) and isotropy is limited within the horizontal plane. For example, a set of horizontal thin-bed shale belongs to VTI case. The velocity distribution in VTI can be shown schematically by a cylinder (Figure 4-2). This type of anisotropy is important in sedimentary basins, and is also the media that this thesis going to discuss.

Figure 4- 2: Velocity distribution in VTI media can be represented by a cylinder. The velocities on all directions within the horizontal plane are identical in amplitude.

When the symmetry axis is not vertical, the case is termed as Tilted Transversely Isotropic (TTI). When the symmetry axis is horizontal, it is called Horizontal Transversely Isotropic (HTI), or horizontal polar anisotropy. In this case, the isotropy is constrained within the vertical plane. HTI media corresponds to the vertical tilted beds or to the case of vertical fractures in isotropic formation.

57

4.3. Definition of Seismic Anisotropy Parameters 4.3.1. Elastic Tensor and Anisotropy Parameters A linear elastic material is defined as one in which each component of stress σ ij is linearly dependent upon every component of strain ε kl (Nye, 1957). Since each directional index may assume values of 1, 2, 3, (representing directions x, y, z), there are nine such relations, each involving one component of stress and nine components of strain. These nine equations may be written compactly as: 3

3

σ ij = ∑ ∑ C ijkl ε kl , i, j = 1,2,3,

(4-1)

k =1 l =1

where the 3 × 3 × 3 × 3 elastic modulus tensor Cijkl completely characterizes the elasticity of the medium. Because of the symmetry of stress (σij = σji), only six of these equations are independent. Because of symmetry of strain ( ε kl = ε lk ), only six of the terms on the right side of each set of equations (4-1) are independent (Thomsen, 1986). Hence without loss of generality, the elasticity may be represented more compactly with a change of indices, following the Voigt recipe:

(4-2) So that the 3 × 3 × 3 × 3 tensor Cijkl may be represented by the 6 × 6 matrix Cαβ. Each symmetry class has its own pattern of nonzero, independent components Cαβ. In practical cases it is not possible to measure all 36 components, so assuming some symmetry in media cause to simplification of matrix. The simplest type of symmetry which is called Isotropic case assumes the symmetry in all directions and can be expressed by three components which two of them are independent. Equation (4-3) shows the isotropic matrix:

58

Cαβ

⎡ C33 ⎢C − 2C 44 ⎢ 33 ⎢C − 2C 44 = ⎢ 33 ⎢ ⎢ ⎢ ⎣⎢

C33 − 2C 44

C33 − 2C 44

C33 C33 − 2C 44

C33 − 2C 44 C33 C 44 C 44

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ C 44 ⎦⎥

(4-3)

These components can be expressed by Lame Parameters and Bulk modulus: 4 µ and C44 = µ (4-4) 3 The simplest anisotropic case of broad geophysical applicability has one distinct C33 = λ+ 2µ = K +

direction (usually, but not always, vertical), while the other two directions are equivalent to each other. This case- called transverse isotropy. The elastic modulus matrix has five independent components among twelve nonzero components, giving the elastic modulus matrix the form of: ⎡ C11 ⎢C − 2C 66 ⎢ 11 ⎢ C13 =⎢ ⎢ ⎢ ⎢ ⎢⎣

Cαβ

C11 − 2C 66 C11 C13

C13 C13 C 33 C 44 C 44

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ C 66 ⎥⎦

(4-4)

where the three-direction (z-axis) is taken as the unique axis. It is significant that the generalization from isotropy to anisotropy introduces three new elastic moduli, rather than just one or two (Thomsen, 1986). Daley and Hron (1977) give a clear derivation of the directional dependence of the three phase velocities:

[

]

(4-5)

[

]

(4-6)

ρv p 2 (θ ) =

1 C 33 + C 44 + (C11 − C 33 ) Sin 2 (θ ) + D(θ ) 2

ρv sv 2 (θ ) =

1 C 33 + C 44 + (C11 − C 33 ) Sin 2 (θ ) − D(θ ) 2

ρv sh 2 (θ ) = C 66 Sin 2 (θ ) + C 44 Cos 2 (θ )

(4-7)

With, 59

[

[

]

[

]

D(θ ) = (C33 − C44 ) 2 + 2 2(C13 + C44 ) 2 − (C33 − C44 )(C11 + C33 − 2C44 ) Sin2 (θ ) + (C11 + C33 − 2C44 ) 2 − 4(C13 + C44 ) 2 Sin4 (θ )

(4-8) Where ρ is density and θ represents the phase angle which is the angle between the wavefront normal and the unique (vertical) axis (Figure 4-3).

Figure 4- 3: Distinction between Ray (group) and Phase angles

And vertical sound speeds for P - and S-waves are, respectively: v p (0) = α 0 =

C 33

(4-9)

ρ

C 44

v sv (0) = v sh (0) = β 0 =

ρ

(4-10)

And horizontal sound speeds for P - and S-waves are, respectively:

v p (90) = α 0 =

C11

v sv (90) = β 90 =

C 44

v sh (90) =

ρ ρ

C 66

ρ

60

(4-11)

(4-12)

(4-13)

]

0.5

Thomsen (1986) rearranged the elastic components and combined some common expressions in phase velocity equations to derive the anisotropy parameters. ε is used to describe the fractional difference between vertical and horizontal P velocities and γ corresponds to the conventional meaning of “SH anisotropy”:

δ* =

1 2C 33

2

ε=

C11 − C 33 v p (90) − α 0 = α0 2C 33

(4-14)

γ=

C 66 − C 44 v sh (90) − β 0 = 2C 44 β0

(4-15)

[2(C

13

+ C 44 ) 2 − (C 33 − C 44 )(C11 + C 33 − 2C 44 )

]

(4-16)

4.3.2. Explanation of Thomsen Anisotropy Parameters The Thomsen parameter epsilon and delta are dimensionless parameters that describe the deviation of a media from isotropy. When the medium is isotropic, both parameters are zero. If the absolute values are around 0.10 or smaller, the medium is weakly anisotropic. If the absolute values are larger than 0.20, the medium is strongly anisotropic. In most rocks epsilon is from -0.05 to 0.30 and delta is between -0.20 and 0.30.Epsilon is the fractional difference between the P-velocity along the axis of symmetry and the velocity orthogonal to the axis. Delta controls the anisotropy in directions close to the axis of symmetry, and in particular the P-wave NMO velocity.

4.4. Conversion Point Mapping

4.4.1. Introduction Tessmer and Behle (1988) showed that in multi-layered strata, the actual P-S conversion point is not a constant offset from the source, but traces a trajectory that moves towards the receiver as the depth decreases (Figure 4-4). The conversion points

61

moves towards the receiver when the offset/depth ratio increases. The asymptote of this trajectory, at large offset/depth is defined by equation (4.17): Xc =

X v 1+ ( s ) vp

(4.17)

Where X is the source-receiver offset, and X c is the conversion point offset from the source. As S-wave velocity is smaller than P-wave velocity, in Common Conversion Point mapping, conversion point lays further respect to source position.

Figure 4- 4: The conversion point traces a trajectory in the multi-layered model for a certain offset, instead of locating on a vertical axis (Stewart et al., 1999).

For larger offsets (or shallow depths), Tessmer and Behle (1988) expressed the trajectory equation in single layer explicitly as: 2

⎡ 2 ⎡⎛ X c ( X − X c ) ⎞⎤ (2τ ) 2 X ⎤ + X − X ( X c − )⎥ = 0 ⎜ ⎟ c ⎢ ⎢ ⎥ 2 z 2 ⎦ τ −1 ⎠⎦ ⎣⎝ ⎣

(4.18)

Where X c is the conversion point offset and X is the source receiver offset.

τ represents the P-S velocity ratio. As these equations are just based on geometrical calculation and algebraic approximations, they will not be reliable for anisotropic media. In the next section the theory of anisotropic CCP binning is explained.

62

4.4.2. Theory As discussed in previous sections Thomsen (1986) introduced anisotropy parameters by which it is possible to model the seismic velocities in each direction in anisotropic media. According to equations (4-5) - (4-8) and (4-14) - (4-16), we can write: v 2 p (θ ) = α 0 [1 + εSin 2θ + D * (θ )] 2

(4-19)

α02 α02 * 2 v sv (θ ) = β 0 [1 + 2 εSin θ − 2 D (θ )] β0 β0 2

2

v 2 sh (θ ) = β 0 [1 + 2γSin 2 (θ )] 2

(4-20)

(4-21)

with

β 1 ⎛⎜ 1− 02 ⎜ 2 ⎝ α0 2

D * (θ ) =

0.5 2 ⎧⎡ ⎫ ⎤ β 0 ⎪⎢ ⎪ 4 ( 1 ) − + ε ε ⎥ 2 * ⎞⎪⎪⎢ ⎪ α 4 δ 0 ⎟⎨ 1 + Sin 2 (θ )Cos 2 (θ ) + Sin 4 (θ )⎥ − 1⎬ 2 2 ⎟ ⎢ ⎥ β ⎪ ⎠⎪⎢ 1 − β 0 ⎥ 1− 02 2 ⎪ ⎪⎢ ⎥⎦ α0 α0 ⎪⎩⎣ ⎭

(4-22) Where α 0 , β 0 are P- ad S- wave velocity along vertical axis in VTI media respectively. The ray angle ( φ (θ ) ) and associated ray velocity ( v φ (θ ) ) can be expressed as:

(tan(θ ) + ( 1 )(dv )) v dθ tan(φ (θ )) = θ tan( ) dv 1− ( )( ) v dθ

(4-23)

And v 2 (φ (θ )) = v 2 (θ ) + (

dv 2 ) dθ

(4-24)

Using Taylor expansion for ε , γ , δ * Thomsen simplified D * (θ ) as: D * (θ ) ≈

δ* Sin 2 (θ )Cos 2 (θ ) + εSin 4 (θ ) 2 2 1 − β0 α0 63

(4-25)

Substitution of equation (4-25) into above equations cause to derivates calculated analytically and equations simplify to: v p (θ ) = α 0 [1 + δSin 2θCos 2θ + εSin 4 (θ )] v sv (θ ) = β 0 [1 +

(4-26)

α02 (ε − δ ) Sin 2θCos 2 (θ )] 2 β0

v sh (θ ) = β 0 [1 + γSin 2 (θ )]

(4-27)

(4-28)

Defining a new parameter δ as below: 1 2

δ = (ε +

δ* ) 2 2 1 − β0 α 0

(4-29)

The associated equations for ray angle of P- and converted-waves can be written as: tan(φ p ) = tan(θ p )[1 + 2δ + 4(ε − δ ) Sin 2 (θ )]

α02 tan(φ sv ) = tan(θ sv )[1 + 2 2 (ε − δ )(1 − 2 Sin 2 (θ sv ))] β0 tan(φ sh ) = tan θ sh (1 + 2γ )

(4-30) (4-31)

(4-32)

According to above equations, it is possible to calculate to ray behavior in a VTI media. In preceding section the methodology of these calculations are described.

4.5. Methodology According to the discussed equations in previous section, I have used an algorithm to find the variation of common conversion point respect to variation of anisotropy parameters, depth of reflection and some other variables. This algorithm has been implemented by MATLAB and results are presented in the next section. Figure 4-5 is the basis of calculations.

64

Figure 4- 5: Depiction of common conversion point and different angle in an anisotropic media

According to figure 4-5 the following procedure has been applied to determine the ccp location and its difference with corresponding isotropic media. 1- A fan of phase angles has been created at an specific depth (Z) 2- Corresponding phase velocities and in consequence ray parameters are obtained 3- Solving the following equation: p(ray parameter) =

Sinθ s for calculation of vs

appropriate phase angles for converted wave corresponding to each P-wave phase angle. 4- Calculating corresponding converted wave velocities 5- Substituting in equations (4-23)

and (4-24) to find group (ray) angles and

velocities 6- Using simple trigonometric relationship for finding the conversion point. 7- Repeating above procedure for isotropic case and obtaining corresponding conversion points 8- Subtracting the results of steps of 6 and 7 to find the difference In the next section this algorithm has been applied and results have been presented.

65

4.6. Results and Discussion The model which is used for this numerical modeling is the one represented in figure 3-6, with the properties of table 3-1. Figure 4-6 show the variation of conversion point from isotropic media with variation of δ * where epsilon is constant and equals to 0.2. This figure is the result of the solution of exact Thomsen equations. As it is shown, by increasing the value of δ * , the variation of conversion point from its corresponding in isotropic media increases to a certain offset. But this trend reverses after that offset which is dependent on the petrophysical properties of the medium.

Figure 4- 6: Variation of displacement of common conversion point in aniotropic media (VTI) versus offset of isotropic conversion point for three values of delta* and epsilon equals to 0.2.

Figure 4-7 shows the displacement versus ray angle of incident P-wave. The general form of these displacement variations respects to variation of δ * has been demonstrated in figure 4-8. It can be deduced that by increasing of the δ * common conversion point displacement has larger rate of increment up to certain offsets and after that the rate of deviation increases for smaller values of δ * .

66

Figure 4-9 represents the effect of epsilon on conversion point displacement which is resulted from solution of exact Thomsen's equations. It is obvious that there is a direct relationship between these two variables.

Figure 4- 7: Variation of conversion point in aniotropic media (VTI) versus ray angle for three values of delta* and epsilon equals to 0.2

67

Figure 4- 8: Variation of displacement of conversion point with variation of isotropic conversion point offset to depth

δ*

and ratio of

Figure 4- 9: Variation of displacment of conversion point with variation of epsilon ( δ = 0.2) *

68

The next analogy examines the effectiveness of linear approximated and exact equations in different media. Figures 4-10 and 4-11 compare the result of linear approximated and exact equations in two different media. The anisotropy parameters are listed in table 4-1.

Strongly Anisotropic Medium Weakly Anisotropic Medium

Epsilon

Delta*

0.3

0,22

0.05

0.03

Table 4-1: Anisotropy parameters for two media one with weak and another with strong anisotropy

Figure 4- 10: Variation of displacement of conversion point in weakly Anisotropic medium (blue dots: exact equation results, red dots: linearly approximated results)

69

Figure 4- 11: Variation of displacement of conversion point in strongly anisotropic medium (blue dots: exact equation results, red dots: linearly approximated results)

The figures 4-10 and 4-11 obviously show that the divergence of results (exact and linearly approximated) for strongly anisotropic medium is much higher than weakly isotropic case. So it is recommended that in the multi-component seismic data acquisition design in media with strong anisotropy like where carbonates and evaporates present, exact equations utilized for common conversion point binning. It is noticeable that solving the exact equations is a time consuming process and results of commercial softwares like NORSAR2D is very close to the results of linearly approximated equations. So care must be taken in such cases.

70

CHAPTER 5

CONCLUSIONS AND RECCOMENDATIONS This study can be subdivided into three main parts: 1- Review of some previous jobs in the study area and analysis of the geological and seismological features of this region to find the data acquisition parameters which are important in recording of the energy. According to this step, the main problems have been detected. 2- Finding the solutions for above problems by optimizing the involving parameters using forward modeling and investigation of multi-component data acquisition effectiveness. 3- Discussing about common conversion point mapping as the key step of multicomponent data acquisition survey design process and analyzing the effect of anisotropy in displacement of common conversion points.

5.1. Conclusions

5.1.1 Affecting parameters according to previous jobs 1- Attenuation is the one of the most affecting parameters on data quality in this region. Due to the presence of steep flanks energy travels to distant locations that even may exceed from the length of acquisition spread. In such cases geometrical spreading and energy absorption are the most attenuating agents of energy respectively. Because of the presence of plastic formations such as Gachsaran formation which is mainly composed of salt and anhydrite, utilizing a method which causes to shorter ray-path can result in better energy records at the surface.

71

2- Arrays are another affecting factor that can improve data quality if designed properly. In this case, experiments show that using source arrays (distributed shots) can attenuate back scattering noises. In addition while using receiver arrays on complex structures it is very important to decrease the size of the arrays to prevent severe static problems. 3- Applying iterative static corrections in processing unit with more reliable velocity model can improve the quality of the data.

5.1.2. Forward modeling approach 1- Forward modeling can be used for better understanding of ray's behavior at underground levels. This technique enables the survey designer to simulate the pseudo real data acquisition for determination of more suitable acquisition parameters. 2- According to the results of ray tracing (forward modeling) by Image Ray Method in the study area on the complex structure, it can be concluded that increasing length of the migration aperture and/or the spread length cause to having better image from study area. Hence using the off end spread is recommended instead of split spread in 2D acquisition surveys. 3- In addition forward modeling by Normal Incident Method showed that denser sampling does not improve the image quality if spread length does meet the appropriate length. 4- Results of simulated field data acquisition by Common Shot Ray Tracing Method illuminated that there are some non-hyperbolic behavior in CMP sorted data at far offsets. Hence, conventional second order velocity analysis will not lead in correct answer in such cases. 5- Having all above mentioned problems, utilizing multi-component data acquisition technique can result in better results. Because in this method not only the attenuation effect are less due to shorter travel path of the rays between source and receivers, but also as S-waves are not affected by gas effects, this method can image in a more efficient manner the target reflectors.

72

5.1.3. Common conversion point mapping Common conversion point mapping is the key step of multi-component data acquisition survey design. According to the results in chapter four, below conclusions can be made: 1- Neglecting the anisotropy factor in common conversion point mapping process cause to severe errors appears in CCP binning that influences on real fold of coverage of the acquired data. 2- According to the results of the solutions of the exact and linearly approximated anisotropy equations of Thomsen (1982), it can be concluded that using linearly approximated equations although leads in faster calculations but generates a large value of error in mapping of conversion points at far offsets.

5.2. Recommendations As it was the first time that forward modeling is applied for analysis of data acquisition parameters, techniques and results in Dezful Embayment area, some recommendations have been presented for next jobs. 1- Recording of S-velocity in well logs for construction of more reliable geological models which must be used in forward modeling technique. 2- Performing three dimensional ray modeling for optimizing 3D data acquisition design parameters. 3- Utilizing more sophisticated softwares than what have been used in our industry for seismic survey design purposes and ray modeling such as NORSAR2D and NORSAR3D.

73

References

1. Abdollahie Fard, I., 2006. Structural models for the South Khuzestan area based on reflection seismic data, Ph.D., Shahid Beheshti University, Tehran. 2. Alaei, B., 2005. Seismic forward modelling of two fault- related folds from the dezful embayment of the iranian zagros mountains, Journal of Seismic Exploration, 14, 13-30. 3. Barkved, O., 2004, Many Facets of Multicomponent Seismic Data, Schlumberger Oilfield Review, Vol. 16, No. 2. 4. Cordsen, A., Galbraith, M., Peirce, J., 2000, Planning Land 3-D Seismic Surveys, Society of Exploration Geophysicists, Tulsa, Oklahoma. 5. Mari, J. L., 1999, Signal Processing for Geologists and Geophysicists, Editions Technip. 6. NORSAR-2D User’s Guide - Version 5.1 rev. 0, 2007. 7. NORSAR-3D User’s Guide - Version 5.3 rev. 0, 2008. 8. Robert Lewis, J. E., 1972, Gachsaran Experiment 1971-1972, Report No. 1201. 9. Robert Lewis, J. E., 1972, Gachsaran Seismic data improvement experiments September 1972 - May 1973, Report No. 1211. 10. Robert Lewis, J. E., 1972, Gachsaran Seismic data improvement experiments September 1973 - May 1974, Report No. 1244. 11. Stone, D. G., 1995, Designing Seismic Surveys in Two and Three Dimensions, Society of Exploration Geophysics, Tulsa, Oklahoma. 12. Sherkati, S., 2004. Tectonic style and kinematics of folding in the Iranian Zagros (Izeh zone): with special emphasis on petroleum systems, Ph.D., Université de Cergy Pontoise. 13. Sherrif, R., 1989, Encyclopedic dictionary of exploration geophysics: Society of Exploration Geophysics. 14. Taksoz, M., 1981, Seismic Wave Attenuation, Society of Exploration Geophysics. 15. Tamimi, N., 2008, Analysis of Seismic Wave Velocities in Incompetent and Plastic Gachsaran Formation, Petroleum University of Technology.

74

16. Thomsen, L., 1986, Weak elastic anisotropy: Geophysics, Vol.51, No.10, 1954-1966. 17. Vermeer, G., 2002, 3-D Seismic Survey Design, Society of Exploration Geophysics, Tulsa, Oklahoma. 18. Walker, D. J., 1976, Experiment Seismic Survey Sulabedar, Report No. 1263. 19. Yang, J., 2003, Numerical and physical modeling of P-S converted waves in VTI media, MSc, University of Calgary.

75