A new method for determining lower density layer in prospection of hydrocarbon

Öz It is known that densities in formations are usually assumed to be constant for gravity model calculations. This also implies that formations are homogeneous and isotropic. However, the formations are usually heterogeneous and densities vary depending on heterogeneity. For this reason, densities should be taken into account as variables. Some scientists consider densities as variables in each formation in model calculations. In other words, density is defined as a function of the required parameters. In fact, functional change is regular. However, density is an irregular variable that depends on the change boundaries of seismic velocity. In this study, it is aimed to take density into account as a variable by using detected seismic velocity boundaries at which seismic velocity changes for each formation. In addition to main formations in model geometry in 3D inversion calculations, another formation was defined. This additional formation has been described by using a combination of all of the change boundaries of seismic velocity present in each formation in a specific order. The density calculated for the additional formation estimated the variation of density between the change boundaries of seismic velocity. This variation is added to the mass densities that are calculated for the description number of each zone. So, lower-density layer comprising oil can be determined by this method. The reliability of the results of the method depends on the reliability of seismic velocity boundaries. Moreover, the increasing number of seismic velocity boundaries leads to the increasing resolution of density variations. Bazı bilim adamları, 3 boyutlu gravite model hesaplamalarında, yoğunlukları her formasyon içinde değişken olarak ele alırlar. Yani yoğunluğu parametrelere bağlı bir fonksiyon olarak tanımlarlar. Bir yeraltı tabakası içindeki yoğunluk değişimi derinlikle orantılı olarak bulunur. Bu çalışmada, her formasyon içinde tespit edilen sismik hız sınırları kullanılarak, yoğunluğun değişken olarak göz önüne alınması amaçlanmıştır. Sismik hız sınırlarının izlediği yol, yoğunluk değişiminin bir göstergesidir. 3B ters çözüm hesaplarında model geometri içindeki ana formasyonlara ek olarak bir formasyon daha tanımlanmıştır. Bu ek formasyon tanımı, her formasyon içinde mevcut olan sismik hız sınırlarının tümü kesintisiz kullanılarak yapılmıştır. İşte bu ek formasyon için hesaplanan yoğunluk, sismik hız sınırları arasındaki yoğunluk değişim miktarı olarak kabul edilmiştir. Bu değişim, ana formasyonlar için hesaplanan yoğunluklara bir düzen içinde ilave edilerek, yoğunluğun derinlikle değişimi ayrıntılı olarak saptanmıştır. Bu çalışma, Adıyaman, Diyarbakır ve Gaziantep bölgesine ait sismik ve açılan kuyulara ait verilerin bir kısmının TPAO’dan alınmasıyla düşük hızlı yer altı modeli oluşturularak yapılmıştır. Bu çalışma sonunda sismik hız sınırlarının ekstra bir kütle olarak alınmasıyla yoğunluğun derinlikle nasıl değiştiği saptanmıştır. Böylece hidrokarbon içeren düşük yoğunluklu tabaka tespit edilmeye çalışılmıştır. Hidrokarbon aramalarında bu yöntem kullanılarak; daha az kuyu açılarak sonuca gidilebilir. Bu çalışmada, başlangıçta yoğunluklar sabit olarak dikkate alınmıştır. Fakat her tabaka içindeki yoğunluklar değişken olarak hesaplanmıştır.


Introduction
Some scientists take the densities to be variables in each formation for 3D modeling using mathematical functions.In other words, density is defined as a function of the required parameters.The gravity anomaly from the whole mass body is an algebraic sum of the contributions of all vertical rectangular prisms at appropriate depths and distances from the observation point.This procedure is widely used in gravityanomaly forward modeling and inversion 2,3,6, 14-19,34,38,3942.For a rectangular prism, a closedform equation for the gravity anomaly is derived by 1,34 when the density contrast is a constant, by 16,38 when the density contrast is a quadratic function of depth and by 18 when the density contrast varies with depth following a cubic polynomial law.
Historically, increase of density and decrease of porosity with depth is of primary interest because of the mechanical compaction arising from the overburden and diagenesis resulting in reduced porosity and vertically layered structure 5-8,16,17,20,24-26,31,33, 36-38,41,46-[49.However, because of complicated geological and geochemical processes in the diagenesis of rocks, metamorphism, intrusives, extrusive volcanics, and facies changes, the density contrast of earth material can also depend arbitrarily on horizontal positions 32,46,49.For instance, as a sediment ages, organic matter combines with mineral constituents largely by physical forces.Changes in the density distribution of sediments can be caused by changes in oxidation or reduction by the surface charges that bind the components into a composite aggregate 4.
Specifically, mechanisms that cause variability in density contrast include dipping layered intrusions 40, folded sedimentary formations, exhumation, overpressure, salt that can result in off-normal compaction curves in sediments, fan development 9, no uniform stratification, physical and chemical cementation 44 and gradual horizontal change in density between two rock types caused by metamorphism 22,35,40.The density contrast of earth material also can depend arbitrarily on horizontal position 32,46,48,49.
Using line integrals is an efficient method to calculate the gravity anomaly for a given density-contrast model 43, 47-49.It is obtained a line integral for irregular 2D masses of constant density contrast for calculating the gravity anomaly 27.It is studied line integrals systematically for irregular 2D masses by defining 2D vector gravity potential and obtains line integrals when the density contrast is depth-dependent or varies vertically and horizontally 47-49.
In this study, the variable densities are defined by using a new method.In this study, process is made by using an algorithm 10,12,13.
As it is mentioned in the summary, an additional formation was also defined in addition to the main formation in the model geometry in 3D inversion process.The amount of density calculated for this additional formation shows the change of the density between the boundaries of seismic velocity with depth.
The regions between seismic velocity boundaries are once again re-defined as the number of rotations.This density calculated for this additional formation is the change amount of common density among seismic velocity boundaries.Density changes in formations are found in detail by adding this amount to the densities calculated for each main formation as the number of rotations.In other words, density is defined as a variable according to the geometry of seismic velocity boundaries.
In the reduction of seismic velocity with depth, density also decreases.In this method, low velocity layers are added to the inversion calculations by changing the direction of the definition of the mass geometry.

3D gravity algorithm
3D model geometry has been triangulated to describe the upper surfaces of the masses.Three points were used to define planes, and the method offered the most convenience for 3D modeling.Even for very complicated mass shapes, a good description can be obtained by increasing the number of triangles.The gravity effect of the bodies was calculated first for the tetrahedron expanded by an "observation" point P to each triangle and then adding them up in a certain sequence.28 Figure 1, Eqs. (1), ( 2) and (3); see for details 11,13,29.The surface of a uniform 3D body can be well approximated as a polyhedron composed of plane triangles to any degree of detail.Moreover, this parameterization is flexible and efficient.
Figure 1 shows the basic uniform tetrahedron, expanded from the "observation" point  at  to an arbitrarily oriented planar triangle, or A-B-C in an earth-oriented Cartesian coordinate system (, , ) with origin  and  pointing vertically downward.This does not reduce generality.First, the gravitational potential  of the tetrahedron is derived for its apex , and then the gravity effect, , is obtained as the vertical component of the potential gradient, by vertical differentiation of .The effects  and  of the polyhedron are the sums of the basic tetrahedral effects.With a consistently defined sequence of computational steps, the partial effects were automatically calculated with the correct sign, i.e., positive for "inside" triangles and negative for "outside" triangles; more specifically, "in" and "out" signify the geometrical relation of the observation point  and the polyhedron.The far-side basic tetrahedral effects were added, while the near-side tetrahedral effects were subtracted such that only the effects of the intervening polyhedron remained.The calculations were also correct if P was enclosed in a polyhedron 13.
Integration of the tetrahedral potential effect, , in arbitrary orientation is awkward, but the orientation is irrelevant for the potential.
Therefore a suitable coordinate transformation is carried out (Figure 2): the , ƞ,  system is defined such that the triangle is in the  − ƞ plane and one edge ( − ) is parallel to .The coordinate transformation is carried out by vector operations (see 13).A FORTRAN code that performs the triangulation and the integration has been developed in the previous studies by ÇAVŞAK 10,12,13.Evaluating the complicated term is unnecessary, and the simpler expression eases the evaluation of  and increases the numerical accuracy.In contrast, when completely written, the expression is fairly complex.
Gravity potential 10 is gravity constant and ℎ is height of tetrahedron.
The gravity effect of the polyhedron is given by the vertical derivative of the potential effect.
Figure 2 shows the parameters of a tetrahedron in Eq. (1).

Geology of the study area
The investigation area is located at the boundary of Southeast Anatolian Thrust Belt and Tauride Orogenic Belt 30.The Southeast Anatolian Orogenic Belt has developed as a result of geological events that had occurred during the closure of the southern branch of Neotethys which was bordered by Taurus at north and Arabian Platforms at south between Late Cretaceous-Miocene time period.The evolution of this belt especially consists of the movement of nappes relatively to the south, Arabian Plate between Late Cretaceous -Miocene 45.
Southeast Anatolian Orogenic Belt consists of three different tectonic units which are E-W trending and separated from each other by north dipping main thrust planes 45.These tectonical units from north to south are the Nappe zone, Accretional Prism and the Arabian Platform (Figure 2).

3D model calculations
3D deep seismic reflection studies that have been done in southeast of Turkey (Figure 3) by TPC (Turkish Petroleum Corporation) have been benefited.The velocity information of the 3D layer and the information of the opened wells which are on the profiles were obtained from TPC.On the initial model; 3 wells on the  = −450  profile, 4 wells on the  = 0  profile and 3 wells on the  = 450  profile were opened during TPC operation in the region (Figure 5).The velocity and depth information of these 10 wells was used in the calculations.Densities of the layers were calculated from the velocity values by using Eq. ( 4) 21.Seismic sections (Figure 4) and depth information of the wells were evaluated together in creating model.Because of the gravity anomalies cannot be obtained, Bouguer gravity data was generated from the initial model with forward modeling.Firstly, gravity values were obtained by utilizing forward modeling (Figure 5).In this solution, an underground model is used.According to the algorithm used in this study, model definition was performed by defining the bottom surface first, then the top surface of every mass.The letters K, L, M and N represent the layer boundaries from the surface to the bottom (Figure 6).The numbers I, II, III and IV written in red are the layer numbers from the surface to the bottom (Figure 6).According to the method used in this study, previously surface L and then surface K of Şelmo formation were defined.Later, surface M and then surface L were defined for the Midyat formation.Then, surface N and then surface M were defined for the Germav formation.Thus, the definition process is completed (Figure 6).Also, Kastel formation is bedrock.Arbitrary noise (white noise) is added by chance to the calculated gravity values for inversion process.Then, density is calculated by inversion calculations for the additional mass constituted by using seismic velocity boundaries.
The new differences in density found from the inversion solution results for masses in Figure 6 are given in Table 3.New densities belonging to the layers were found by using Eq. ( 5).The densities were calculated as a result of the inverse solution Table 6.Thus, an attempt was made to determine real densities by adding the aforementioned change in amount of density, calculated with two different identifications for additional mass, to the main formation densities by the number of definitions for each region between seismic velocity boundaries.
Locations of the opened wells are shown on the each profile (Figure 6 and 7).
For layers at which the velocity increases downward; first, the sub-surface, then the top surface are defined in definition.For layers at which the velocity decreases downward; first, the top surface, then the sub-surface is defined.The letters A, B, C, D, E, F, G and H represent the surfaces of seismic velocity boundary from surface to the bottom (Figure 7).The numbers 1, 2, 3, 4, 5, 6 and 7 written in red are the layer numbers from surface to the bottom (Figure 7).
In the first identification of the additional mass, it was accepted that the velocity increases downward in Germav formation.The definition of the masses between seismic velocity surfaces is made by using B-A-C-A-D-A-E-A-F-A-G-A-H-A surfaces continuously in ascending order.As it can be seen here, the masses are described once more starting in the amount of rotations in such a way that one is within the other.In other words, while the number 1 mass is defined 7 times, the number 2 mass is defined 6 times and the last number 7 mass is defined only once (Figure 7).The inverse solution technique is applied by giving the defined values of seismic velocity boundaries to the algorithm 10 as additional mass in addition to the defined values of the model.The differences in density of the masses in the model are calculated Table 4.
Table 4: Density differences found from the inverse solution made after the first definition of additional mass.

Mass
Nu.
Calculated Density Difference g/cm 3 I -0.3043II -0.0436 III -0.0860Additional Mass -0.0026 Masses have been turned a few times from the defined surface in the defining process of seismic velocity boundaries.Therefore, differences in real density for the first definition are calculated by using Eq. ( 6). Here; ∆  : Difference in real density calculated for the places between each seismic velocity boundary calculated Table 6.
In the second identification of the additional mass, it was accepted that the velocity decreases down in Germav formation.The definition of the masses between seismic velocity surfaces is made by using B-A-C-A-D-A-E-A-G-H-F-H-E-H surfaces continuously in ascending order.The number 1 mass is defined 4 times, the number 2 mass is defined 3 times and the last number 7 mass is defined 3 times (Figure 7).
The differences in density obtained from the inverse solution made after this second identification are shown below Table 5.Also in this last inverse calculations made by changing the data calculated by the straight analysis, in addition to the digitized values of the underground model, also the digitized values of seismic velocity boundaries are given to 3D gravity algorithm 10,12,13, as the additional mass, inversion technique is applied and the density differences are calculated between the masses of the subsurface model Table 6.It has been restored digitized circus several times in the digitization process of the seismic velocity boundaries.Therefore, the true density differences are calculated using the Eq. ( 6).
Actual density differences for identification in the third inversion; If the differences in density for the second definition are calculated;   6.
As shown in Table 6, using seismic velocity boundaries; the first mass is divided into two layers, the second mass is divided into two layers and the third mass is divided into three layers (Figure 7).The density of each layer is found from the inversion calculations.It has been neglected that the seismic velocity decreases with depth in the mass III in the first identification.It has been included to account that the seismic velocity decreases with depth in the mass III in the second identification.
Then again using (Eq.4), real densities from density differences are calculated Table 6.
As it can be seen in Table 6, layers densities are obtained as a result of two different definitions that are made for the additional mass.It has been found that the layer 7 has the lowest density from the result of the inversion calculation in the second identification.

Conclusions
In this study, different densities for both main formations and additional mass are calculated as a result of inversion calculations made by using two different methods employed in the definition of additional formation.The layer 7 (Figure 7) in the third layer, since the density is low, the layer is though as probably oil bearing reservoir rock.The cost can be minimized by using this method and opening less exploratory wells in the hydrocarbon exploration areas.
Irregularity in the geometry of seismic velocity boundaries refers to the irregularity shown in the density of the variable.This situation varies because of the practice of taking the density as a variable in a mathematical expression.Mathematical definition of the density as a variable is actually an expression of the regular density change in every formation.This requires that seismic velocity boundaries have to be parallel to each other, i.e. parallel to the ground.Mathematical definition of irregular changes of seismic velocity boundaries does not have a practical side as in inversion process.This method which takes the density as a variable depending on seismic velocity boundaries represents the density variations in a mass very well.

Figure 2 :
Figure 2: Explanation of the nomenclature used in describing the tetrahedron, after coordinate transformation from , ,  to the system , ƞ,  (see text) with the triangle, lying in the bottom  − ƞ plane.
The intensive tectonical activity that occurred at the end of Cretaceous and Miocene and caused the settlement of allochthonous units on the region has also given rise to the development of marine deposition at the same period and to the closure of the basins.Allochthonous units have taken their recent positions at the end of Upper Miocene as gravity sliding and overthrust sheets.Late Cretaceous Koçali complex has tectonically been settled on Karadut complex (lower allochthonous series) at bottom and on the Kastel formation forming the uppermost horizons of the Arabian Plate.These are uncomformably overlain by Upper Maestrichtian-Paleocene transgressive deposits (Terbüzek formation, Besni formation and Germav formation) of the Arabian Platform belonging to the Autochthonous Series.The region has been transgressed at the beginning of Eocene and Lower Eocene-Lower Miocene transgressive deposits (Gercüş formation, Midyat formation, Gaziantep formation and Fırat formation) have uncomformably overlain units at bottom.

Figure 7 :
Figure 7: Seismic velocity boundaries in the model geometry for the additional mass definition.

Table 1 :
Thicknesses, velocities and densities of layers in model.The differences in density were found by subtracting the reference density from the mass densities.The densities belonging to the model are shown in Table2.These are reliable density parameters.
Figure 4: Seismic sections of model.

Table 2 :
The densities and density differences of masses.

Table 3 :
The density differences of masses calculated with inversion solution.

Table 5 :
Density differences found from the inverse solution made after the second definition of additional mass.

Table 6 :
Densities calculated from the three inverse solutions; (a): The layers separated with seismic velocity boundaries, (b): Density calculated by without using seismic velocities, (c): Calculated density without taking into account the velocity decrease downwards in the Germav formation (III) and (d): Calculated density with taking into account the velocity decrease downwards in the Germav formation (III).