Missakeswal oil field is located 60 km SE of Islamabad, lies on the Eastern extremity of the Potwar sub-basin, Punjab Province, Pakistan (Figure 1). The estimated terrain elevation of area is 500 m above mean sea level [1]. Potwar basin is an important host to hydrocarbon reserves, lies at the northern side of the Upper Indus basin [2]. Geographically study area bordered by Main Boundary thrust in North, Jhelum strike slip fault in East, Kala bagh strike slip fault lies towards West and Salt range in south. Pakistan has high potential of hydrocarbons in its northern (Potwar, Kohat) and southern (Badin, Mari etc.) parts. The topography Potwar plateau is undulating comprised of E-W trending low hills, parallel ridges and dissected valleys. Chorgali and Sakesar Formations predominantly comprised of limestone of Eocene age are the well-known oil producing reservoirs [1]. Due to presence of detachments at different levels, complex structures, anomalous evaporates thickness and technical drilling problems; so far most of the wells were unsuccessful in this area [3-5]. Two-dimensional seismic reflection data reveals the presence of normal faults in basement, and reverse faults in the cover sequence, which can be seen on the interpreted sections [5-7].

Figure 1: Satellite image of Pakistan showing Missakesal area (Khan et al., 2008).

This study is intended to evaluate the hydrocarbon potential by using seismic reflection, well log and velocity data, will be helpful in future for prospects generation and hydrocarbon exploration in Potwar Basin. The average, interval and RMS (Root Mean Square) velocities will helps in the study of subsurface lateral variations, time to depth conversion and their correlation with synthetic seismogram on depth section affirms interpretation. Exploration companies integrate geophysical, geological and well data to carry out reliable subsurface structural interpretation. Precise description of the subsurface geological horizons (Chorgali, Nammal and Sakesar), structures and seismic attributes provide an important insight to hydrocarbon exploration. The basement with normal faults and cover sequence with fault bounded anticlines displaying Jurassic rifting, Pangea splitting and associated diapirism of Himalayan orogeny in the Potwar plateau. The main discoveries of the Misskeswal area can be related with salt cored anticline, which are formed because of crustal shortening in the cover sequences. The calculated average energy and instantons phase attributes helps in study of reflection continuity, faults, lithological dissimilarities and data type. Similar approach used in the proposed study will also aid in delineation and exploration of regional subsurface structures in the Potwar basin and nearby areas.

Geological settings

The Himalayan collision system represents an active collision Orogen between Indian and Eurasian subcontinents. The collision is active since about 55 Ma [8] involves continuous uplifting erosion and deposition of the sediments. The Indus basin, including the Kohat- Potwar (study area) depression, belongs to the category of extracontinental downward basins, which account for 48% of the world petroleum resources [9]. The tectonic depression of Potwar sub-basin formed as result of continent-to-continent collision at NW margin of the Indian plate. The foreland zone; comprised of Salt range, Kohat and Potwar plateau and Hazara Ranges, the area is bounded by Panjal- Khairabad fault in North Salt Range thrust in South, whereas eastern end up with NS trending left lateral Jhelum fault [10]. Potwar sub-basin is predominated by structural traps and mostly seismic data used for the delineation of these structures. Missakeswal area is well known for its hydrocarbon potential, some of the exploration trends and oil field locations of the Potwar plateau is shown in Figure 2.

Figure 2: Geological map of Potwar plateau, oil fields and wells are higlighted (William et al., 1998).

Potwar plateau comprised of 5000-1000 m thick sediments of Precambrian to recent age, these sediments range from fluvial to marine environment of deposition [11,12]. The basement of Precambrian is mainly dominated of slates and phyllites. The Salt Range formation unconfirmably overlying the basement, composed of metamorphic and Volcanic rocks of Indian Shield [13], and overlaind unconfirmably with Khewra Sandstone of Early Cambrian. Contrary to the well data obtained from Adhi and Missakeswal Kussak (sandstone and shale), Jutana (dolomite) of early Cambrian and middle Cambrian Baghanwala (shale, sandstone, siltstone) lies conformably over Khewra sandstone [5]. The Nilawan group comprising of Dandot (sandstone) and early Permian Tobra (conglomerates), whereas Warcha and Sardhai formations are eroded in study area.

The unconfirmable boundary exposed the area from upper Permian carbonates to lower Paleocene [14]. Hangu (sandstone), Lockhart (limestone) and Patala (shales) are well developed Paleocene sequence and Nammal, Sakesar and Chorgali Formations of Eocene age conformably overlies Paleocene strata. The Mesozoic sediments are eroded and Oligocene sediments have not been found in the region [2]. Stratigraphy of the Potwar plateau is well developed along the Salt range; a generalized stratigraphic chart is given in Table 1.

Table 1: Generalized stratigraphic Chart of the study area.

Potwar marine facies has great potential of hydrocarbon generation. Recent discoveries in Potwar result in delineation of deep subsurface crest [4]. The gray shales of the Datta (Jurassaic age) and Patala Formations (Paleocene age) are potential source rocks in Salt Range Potwar- Foreland Basin (SRPFB) [1]. The Cambrian, Permian, Jurassic, Paleocene and Eocene reservoirs are producing oil in Salt Range Potwar-Foreland Basin (SRPFB). Thick layers of evaporite and shale have good sealing potential for Infra-Cambrian reservoir. Most of the fields discovered in Kohat-Potwar geological province are overturned faulted anticlines, popup structures or fault-block traps. The latest trapforming thrust event began at approximately 5 to 2 Ma [15].

Presently two-fold division is envisaged for the Potwar sub-basin. Based on the deformation style; (1) the Northern Potwar Deformation Zone (NPDZ) and (2) the southern Potwar platform zone [1,4,8].

The Northern part of Potwar is more intensely deformed, known as the North Potwar Deformed Zone (NPDZ), characterized by East West, tight and complex folds, overturned to the south and sheared by steep-angle faults [10]. The central Potwar comprised of fault bounded structures mainly thrust and back thrusts, while asymmetric anticlines with single faults have been observed in different locations. Main structure in the Potwar area based on seismic interpretation is Salt cored anticlines bounded with faults on limbs and Triangle zones [16].

Work Flow and Methodology

Seismic, well log, Velocity and Navigation data of fifteen seismic lines were acquired from DGPC (Directorate General of Petroleum concession), Government of Pakistan. Missa-keswal area is projected in Universal Transverse Mercator (UTM, Zone 43) projection system, Base map generated using K-tron Precision Matrix is shown in Figure 3.

Figure 3: Base map of the study area displaying wells and seismic lines (Khan, 2000)

Different modules of K-tron software used to accomplish various processes; workflow to carry out this research work is shown in Figure 4. Interpretation has been carried out by using X-works software, provides an interface for marking horizons and faults, velocity and depth data for contouring and for further analysis such as crustal shortening and 2D seismic modelling. VAS (Velocity Analysis System) was used to calibrate velocity data which is used in time to depth conversion [17]. Crustal shortening and 2D modeling performed to confirm compressional regime and structural interpretation. The seismic velocities used in different applications like time to depth conversion, contouring and Modeling. Using depth section and velocities, we can also depict structures, lithological boundaries and the lateral variation in the subsurface area. Seismic attributes analysis affirms interpretation along with that also provide information about bedding continuity, lithology and helps us to predict the probable reservoir of the area.

Figure 4: Integrated Work Flow used in study.

Seismic interpretation

Initially different reflector was marked based on prominent reflections; well data used as secondary information from subsurface depths utilized to pick prominent reflections. Each reflector is marked with different color; so, that they can easily distinguished, main reflectors of on the sections are Chorgali, Lockhart and Sakesar Formations (Figure 5). Processing of Seismic velocities (conversion, interpolation and smoothing) carried out by using X-works software, which then used in different applications like time to depth conversion, contouring and Modelling [18]. Depth sections give real picture of the subsurface and the horizontal variation in the horizon, which occur due to the lateral variation of the velocity, because only time section is insufficient to delineate the results (Figure 6). Average velocities are multiplied with time of each reflector and the divided by two to acquire depth section (Figure 7). Horizon interpolated average velocity contour map of Chorgali Formation generated to understand the spatial velocity variations provide true depth imaging of the formation (Figure 8).

Figure 5: Interpreted time section of seismic line GNA-19.

Figure 6: Seismic velocities of GNA-19 combinedly displayed at different CDP’s.

Figure 7: Depth Section of seismic line GNA-19.

Figure 8: Chorgali Formation Horizon Velocity contour map.

Crustal shortening

When two continental crusts collide and neither sub ducts, the material pushed towards the Earth’s surface resulting in the formation of mountains like Mount Everest. This causes the crust to reduce in size [19]. The crustal shortening of all the three formations of seismic line GNA-14, is determined by loading the depth section in X-Works that gives the length of each segment of horizon (Figure 9). The horizon divided into different segments because of reverse faulting. The crustal shortening calculated by using the formula given below;

Figure 9: Geological cross-Section for crustal shortening with length of each faulted part.

CS = L – Lo (1)

L is the total length of the horizon obtained by adding all the segments of the horizon. Lo the length of the seismic section calculated by multiplying the CDP (Common Depth Point) number with the CDP interval. Crustal shortening performed to confirm compressional regime and structural interpretation. The horizon length and their corresponding crustal extension are given in Table 2.

Table 2: The horizon lengths and their corresponding crustal shortening.

Volume of shale

Missa-01 well illustrates reservoir and source zones, the prominent shale signature is corresponding to the location of the Patala shales at depth 1750 and 2250 M as given by the information from the well tops, while the zone which can be seen above shale is reservoir at depth of about 1800-2020 m, which corresponds approximately to the depth of Chorgali and Sakesar Formation. The response of Gamma ray log along with clear description of shale and clean line is shown in Figure 10.

Figure 10: Gamma Ray log (red) with calculated Volume of Shale curve (blue). [Cut-off line (sky blue), shale line (green) and sand line (yellow)].

Seismic modeling

1D- forward modeling: An artificial seismic reflection record, prepared from velocity log data by convolving the product of acoustic and density logs (reflectivity function) with seismic data derived wavelet is termed as Synthetic Seismogram [20]. Synthetic seismogram of Missa-02 well generated by using Wavelets software [21].

In this procedure, the petrophysical logs; Sonic (DT) and Bulk Density (RHOB), which provide the inverse velocity and density information of subsurface layers respectively, this velocity and density data used to compute a series of reflection coefficients called reflectivity series. Source Ricker wavelet with a dominant frequency of 35 Hz is generated. The reflectivity series then convolved with the source wavelet shown in Figure 11. Synthetic seismogram matched with the interpreted seismic depth section at the well point to correlate the succession of reflectors as shown in Figure 12. It also used to calibrate our seismic velocities. 1D modeling is important to study the area at different points and it does not give information of lateral variation that is computed by using 2D modeling.

Figure 11: A Ricker Wavelet Generated with frequency 35 Hz used as a source for Seismic Modeling.

Figure 12: Geological Depth Section of GNA-19 Correlated with Synthetic Seismogram generated from sonic log.

Depth section based seismic model (2D Seismic Modeling)

The reverse process of seismic interpretation is now as seismic modeling, interpretation is carried out by marking the horizons based on continuity of reflector on the Seismic section whereas, in 2D seismic modelling the synthetic seismic section, is generated by the convolution of Ricker wavelet of 35 Hz with the Geological section of seismic line, here we are using GNA-19 at 25 trace intervals as shown in Figure 11.

The time sections are converted to depth sections using velocity data and reflection coefficient were assigned to each horizon as an attribute. The geological section is convolved Ricker wavelet at every 10th CDP interval to generate a 2D synthetic seismic section as in Figure 13. This modeling is based on structural data, so the generated synthetic section corresponds to the subsurface structures. 2D Modeling is important in the areas or points where they arise lateral variation that may give rise to focusing problem, shadow zone, diffractions, and beneficial for adjustment of acquisition parameters foe neo surveys and seismic interpretation affirmation etc. [22].

Figure 13: Geological cross-section with generated Synthetic Seismic Model.

Seismic attribute analysis

Seismic attributes are widely used for lithological and petro physical prediction of reservoirs. Various methodologies have developed for their application to broader hydrocarbon exploration and development decision-making [23]. There are no specific rules for computing attributes, a quantity that is directly calculated from seismic data is called seismic attribute [24]. Seismic attributes affirm interpretation along with that also provide information about bedding continuity, porosity, phase, and lithology and helps us to predict the probable reservoir location [25].

Envelope of trace (reflection strength) and average energy/ apparent polarity

Envelop as a measure of local energy level of seismic traces has related to the interval thickness of the corresponding geological unit.

Trace envelop helps in horizon and fault interpretation particularly smaller faults, to recognize the major changes in depositional environment and provide aid in Spatial correlation to porosity and other lithological variations [26].

Envelop trace is computed for seismic line GNA-19 to see the major changes in lithologies. The maximum reflection strength indicated by the thick yellow packages, and it also showing spatial patterns representing changes in the limestone thickness and breakage due to the faults (Figure 14). Weak reflections are due to the interbedded shales [27,28]. Even negative reflection coefficients such as limestone formation overlaid on clayey formation would generate a positive response in this attribute.

Figure 14: An Envelop attribute map of GNA-19 seismic line [Chorgali (Red), Sakesar (Green), Lockhart (Blue)].

Average energy is a post-stack wavelet attribute, which is computed at the peak of envelop attribute within a specified window, in a defined zone the attribute of wavelets are represented by trace envelope minima. Thus, wavelet attributes are computed over each cycle of the seismic trace [29]. Apparent polarity shows the blocky nature as it represents the average value over each cycle in the seismic trace, which can be observed in Figure 15. Each peak in the reflection indicates the cycle. The Polarity base attribute can have used for better understanding in interpretation particularly for fault interpretation [30].

Figure 15: Average energy attribute map of GNA-19 [Chorgali (Red), Sakesar (Green), Lockhart (Blue)].

Instantaneous phase

The information of phase is related with the propagation of seismic wave front, and independent of trace amplitude. It computed from real and imaginary traces as given below;

Ө (t) = tan-1[h(t)/f(t)] (2)

Instantaneous phase attribute used to highlight the deeper horizons that are not visible in the normal amplitude section as shown in the Figure 16. The interpreted horizons lie over the minimum phase regions indicated by orange color, instantaneous phase attributes which changes from -1800 to +1800 [31]. This attribute further confirms the interpretation as the input data is minimum phase. The negative value in color bar is indicating negative phase. It can be observed in comparison to amplitude-based sections that the instantaneous phase shows much deeper horizons. Phase attribute is sensitive to lithological boundaries and good indicator of lateral continuity, so this can be used for marking boundaries in interpretation.

Figure 16: Average energy attribute map of GNA-19 [Chorgali (Red), Sakesar (Green), Lockhart (Blue)].

• Interpreted seismic reflection data of Missakeswal area indicates the presence of normal faults in the basement and salt cored anticlines bounded by reverse faults on both limbs in the cover sequence.

• The study area had experienced crustal shortening, as result of continent-continent collision at NW margin of Indian plate and the tectonic forces acting on North Potwar Deformed Zone (NPDZ) affirmed the compressional tectonic regime in the study area.

• ISO (International Organization for Standardization) Velocity contour maps indicate the velocities vary laterally and horizontally due to compaction, lithological contrast, faulting and different burial depth.

• Synthetic Seismogram matches with the marked reflector, which confirms stratigraphic interpretation, and interpreted 2D model of geological cross-section confirms the type and nature of sub-surface structures.

• The structural and lithological interpretation is further affirmed through seismic attributes maps i.e. envelope, instantaneous phase and average energy. Instantaneous phase attribute provided good results and the precise positions of faults. The overall results indicate the economic viability of Chorgali limestone as a reservoir.

I am thankful to the Director General, Petroleum Concessions (DGPC), Pakistan, for the provision of data for this research work and to Sir Khalid Amin Khan (OIST) for his thoughtful guidance and technical manner, I found it quite helpful for me.