Predicting the gas resource potential in reservoir C-sand interval of Lower Goru Formation, Middle Indus Basin, Pakistan

1 Introduction

From the last few decades, advancements have been seen in hydrocarbon exploration through the invention of prolific geophysical interpretation techniques for hydrocarbon identification [1,2,3]. Among the various seismic interpretation methods, seismic attribute analysis has proven to be the most beneficial technique to detect the reservoir structural and stratigraphic variations away from the well [4]. Seismic attributes are basically any measured information derived by either direct measurements or logical reasoning from the seismic data [5]. The 3D seismic attribute analysis is used to map features from basin to reservoir scale that enhances the spatial prediction of reservoir properties from the seismic data [6,7]. These attributes generate strong reflection coefficient along the lithological interfaces, indicating the bright spots that are usually associated with hydrocarbon (gas) accumulation in sandy reservoir [8,9].

However, besides this, the analyst does not merely rely on the seismic attribute analysis and attempt to enhance the interpretation resolution by inversion of given seismic data to any other reservoir property. Poststack seismic inversion is one of the advanced geophysical techniques that can enhance the interpretation accuracy and has the ability to estimate subsurface rock and fluid models of varying physical characteristics [10]. In the current study, the poststack inversion technique based on the linear programming sparse-spike algorithm is used to obtain full band acoustic impedance (P-impedance) cube. The algorithm of sparse-spike inversion is very vigorous and commonly used to generate the Earth’s P-impedance from the 3D seismic data cube with minimum number of spikes to obtain the absolute subsurface information [11,12]. The postinversion results can be processed further through geostatistical methods to predict reservoir parameters such as lithology, porosity, and fluid saturations in the inter-well region [12,13,14,15]. In addition, the statistical methods aim to search the linear or nonlinear operator that can predict well-logs from the neighboring seismic data. These methods enhance the prediction power by analyzing the nonlinear seismic attributes away from the wells along with the low-resolution seismic reflection data [13].

In the recent studies, Indus Basins of Pakistan have been extensively evaluated for hydrocarbon potential and petroleum prospects [16,17,18]. The authors investigated that the widely distributed reservoir unit in lower and middle Indus Basins is the Lower Goru Formation of Cretaceous age. Despite the reservoir unit is distributed widespread in the area, it has heterogeneous intrinsic properties and variable thickness due to diverse depositional environment. Moreover, the sand of reservoir unit appears as the composite reflection on conventional seismic amplitude data due to intercalations of shale beds, resulting in originating obscurity in the seismic data interpretation [19,20]. Therefore, it is necessary to discriminate the interbedded shale layers from sand bodies, and the objective could be achieved by improving the resolution of the seismic data. Multiple studies attempted previously to characterize the reservoir unit based on petrophysical, rock physical and seismic interpretation [20,21,22]. A recent study by Ali et al. [23] evaluated the reservoir properties of Lower Goru Formation using different poststack inversion algorithms, and another recent study by Ashraf et al. [24] carried out the electrofacies classification using the 3D seismic attribute analysis, but none included essential integrated study of 3D seismic attributes and poststack inversion analysis with efficiency of geostatistical methods for predicting the promising reservoir feature (gas sand channel) in the area.

This study aims to delineate and characterize the gas resource potential within the reservoir interval of the study area using the 3D seismic attribute analysis and poststack inversion technique along with the geostatistical methods. Moreover, to predict the reservoir feature like sand channel that majorly takes part in hydrocarbon accumulation and improves the resource potential. The data used for the interpretative analysis comprised a 3D seismic cube and wireline logs of four wells (Sawan-01, Sawan-07, Sawan-08 and Sawan-14) from the study area of Sawan gas field, Middle Indus Basin, Pakistan (Figure 1).

Figure 1

The study area of Sawan gas field is highlighted on the regional structural map of Pakistan (modified after Anwer et al. [20]) along with the base map.

2 Structural and stratigraphic evolution of the region

The study region of Sawan gas field is a part of Middle Indus Basin of Pakistan, covering the hydrocarbon producing Formation of Cretaceous age (Lower Goru). The study area was discovered in 1997 and located 300 km north-east of Karachi in Thar Desert within Khairpur District of Sindh province [25]. The prolific reservoir unit of Lower Goru Formation is the C-sand interval with approximately 103 m thickness at the average depth interval of 3,000–3,500 m in the area [26].

The study area is bounded by several structural highs and lows within the tectonically controlled basement of the Precambrian age [27]. The Jacobabad–Khairpur and Mari–Kandhkot highs cover the southern and Sargodha High covers the northern side of the study area (Figure 1). The Indian Shield and Kirthar range bound the area from eastern and western sides, respectively [18]. The horst structure of Jacobabad–Khairpur high divided the Indus Basin into Central and Southern Basins and Sawan area is situated on the southeastern flank of Jacobabad–Khairpur high. The uplifting of Jacobabad–Khairpur high significantly developed the structural traps of reservoir quality in Sawan area and in other adjacent exploration fields during the Cretaceous time [18,28].

The regional tectonic and stratigraphic information revealed that the Paleogene uplifts of Ranikot Formation were initially encountered as pinch-out with Jacobabad–Khairpur high. The existence of carbonates and associated marls indicated the development of complex basinal and subbasinal paleo-topography during the early Eocene epoch and latterly filled with the Gazih Shale [27,28].

The tectonic collision of Indian and Eurasian plates had contributed in the formation of deep basement-rooted and shallower wrench faults in the area passing through the whole Cretaceous sequence [29,30]. The crustal thrusting developed the forebulge of paleo-highs during the Mesozoic and Tertiary periods in the northwest parts of the study area, and this was the last episode of hydrocarbon migration and trap formation in the Lower Indus Basin [28]. The stratigraphic sequence (Jurassic to Quaternary age) ranging from Chiltan limestones at the base to alluvial rocks at the top covers in-between the Sembar and Lower Goru Formations of Cretaceous age. The Sembar Formation is a proven source rock of Indus Basin and mainly comprised black shale with minor filling of sandstone, argillaceous limestone and dark siltstone [31]. The Sembar Formation is thermally mature due to its deep burial from western edge and less mature and comparatively shallower from eastern edge in the Indus basin with variable thickness ranging from 0 to 260 m [32]. The tectonic uplifting of Jacobabad–Khairpur high developed the favorable depositional scenario for the main reservoir unit; the proximal reservoir sand was accumulated in deep structural areas, whereas nonreservoir shale layers were positioned up-dip in the distal side to develop the major trap for hydrocarbon accumulation [33,34]. The Lower Goru member of Goru Formation covered the Sembar Formation and was considered as the major reservoir unit of the study area [35]

The sequence of Sembar–Goru Formations has various progrades passing from shallow distal to proximal facies [34]. The composition of Lower Goru Member from upper side is mainly shales, whereas it comprises medium to coarse-grained sandstones with thin layers of limestone and shale from its lower part. The subordinate portion of the reservoir interval can be classified into distinct sand intervals: A, B, C and D intervals (Figure 2). Petrographically, A and B intervals have mostly quartz-rich arenites, whereas the major portion of D-interval consists of black shales. The C-sand interval mostly comprises pore lining iron chlorite cement and partially altered volcanic rock fragments. The sands of the C-sand interval can be classified as sublitharenites, lithic arkoses and litharenites, whereas other diagenetic components include carbonates, glauconite quartz cement and chlorite cement [35]. The C-sand interval covers the main gas reservoirs in the research area, accommodating high average porosity of around 17% within the depth range of 3,000–3,500 m. The Upper Goru Member consists of reservoir transgressive shales, siltstones and marl, which acts as a regional seal and extend from shoreface to foreshore facies of the prograding shoreline [34].

Figure 2

Stratigraphic chart showing the petroleum play of the study area. The principal gas-producing reservoir zone (C-sand interval) is highlighted with red color (modified after Krois et al. [53]).

3 Dataset calibration and methodology

The proposed methodology for predicting and characterizing the gas saturated zones within the C-sand interval of Lower Goru Formation of the study area using the 3D seismic and wireline log data is presented in Figure 3.

Figure 3

Generalized workflow of integrated study applied to predict the gas resource potential in the area.

The wireline log data of density (RHOB), sonic (DT), gamma ray (GR), caliper (CALI), resistivity, spontaneous potential (SP) and neutron porosity (NPHI) logs are used for this study that lies within the seismic domain of study area. An integrated method is adopted using the 3D seismic and well-log data to establish the tie for identification of key stratigraphic surfaces (i.e., Lower Goru, Top and Base of the C-sand interval) on the basis of strong reflection continuity (Figure 4a). The 3D view of the arbitrary seismic line passing through the four wells location is displayed along with time amplitude map of C-sand top horizon (Figure 4b). The petrophysical analysis is carried out to transfer the wireline log data into predictable parameters such as volume of shale, porosity, water saturation and pay zone [36]. The porosity of the shale is eliminated using estimated shale volume to correctly calculate the effective porosity. The general workflow adopted for the petrophysical analysis to quantify the reservoir parameters is presented in Figure 5.

Figure 4

(a) Interpreted seismic section showcasing the horizons of interest along arbitrary seismic line. (b) 3D view of arbitrary line passing through the four wells over the time amplitude map of C-sand top horizon.

Figure 5

The workflow adopted for petrophysical analysis in this study.

3.2 Poststack inversion procedure

The poststack inversion procedure involves the synthetic seismogram generation (wavelet estimation and well to seismic data calibration), initial model/low-frequency model, inversion analysis and linear programming sparse-spike inversion.

The necessary transformation of seismic amplitudes to impedance values in inversion can be done by the process of convolution. The process generates the synthetic seismic trace (S _t) by convolving the estimated wavelet with Earth’s reflectivity series [40,41], which can be expressed as follows:

(4) S t = W t × R t + N t ,

where W _t is the extracted statistical wavelet at a constant phase, R _t is a reflection coefficient series and N _t is the random noise. The statistical wavelet used for the correlation between traced and inverted reflectivity from seismic data at Sawan-08 well location is extracted within the window length of 1,700–2,500 ms in the seismic time domain with 200 ms wavelength, and synthetic seismogram is then generated by the convolution of the extracted statistical wavelet with reflectivity series. The seismic data lack the low-frequency component due to its band-limited nature, which is essential for estimating the reservoir properties. The P-impedance logs generated from the wells data must be added to provide the missing frequency components (low and high) to seismic data for realistic model building [41]. The low-frequency model of P-impedance is estimated within the time window of 1,700–2,500 ms for the application of linear programing sparse-spike inversion (Figure 6).

Figure 6

Low-frequency model of P-impedance estimated for the application of sparse-spike inversion.

3.2.1 Linear programing sparse-spike inversion algorithm

The low-frequency model provides the missing frequency component for analysis and also constrains the inversion results to control the nonuniqueness of the solution. After setting the input parameters for the inversion process, the linear programming algorithm is applied on the input data in two steps. In the first step of quality control (QC) test, inversion analysis is performed between P-impedance log and extracted composite traces that are inverted for P-impedance around well location. In the second step, the technique transformed the entire seismic section into P-impedance volume.

The inversion technique based on the linear programming theory, which was developed by Sacchi and Ulrych [42] and later modified by Li [11], is adopted for this study. The purpose of inversion methodology is to generate a series of reflection coefficients as a bulk subsurface model. The objective could be achieved by using seismic data (S _t) and source wavelet (W _t) as input parameters and by estimating P-impedance from that model. The relationship among the data d = (x ₁, x ₂, …, x _n ), model m = (r ₁, r ₂, …, r _n ) and noise “n” is presented as follows:

(5) L m + n = d ,

where L is the operator to measure the misfit between d and m. Generally, the d is known, and the subsurface m vector can be defined by a probability P(m|d) and can be expressed in Bayes’ formula in equation (6):

(6) P ( m | d ) = P ( d | m ) × P ( m ) P ( d ) ,

where P(d) is a normalization factor from prior knowledge of data vector, P(m) is a prior knowledge of the model and P(d|m) is a likelihood function between model and the observations. Thus, this additional information of the prior model is obtained from the wireline log data and is used to estimate the difference between synthetic and seismic traces [11]. For the model based on the posterior probability P(m|d), we could use the maximum-a-posteriori (MAP) solution that maximizes the P(m|d). The objective function Y can be obtained as follows:

(7) Y = − log ( P ( m | d ) ) = − log ( p ( d | m ) ) − log ( p ( m ) ) .

As the constant term of log(P(d)) can be omitted for simple operation, the solution can be obtained from the maximization of objective function Y in least square sense [11]. The prior knowledge about the model is generally known as a global constraint S(m). Afterward, the principle of maximum entropy can be used to compute the prior probability [42].

Consider the continuous model parameters (m) with the probability density function p(m). The entropy factor h is presented as follows:

(8) h = ∫ P ( m ) log [ p ( m ) ] d m .

Equation (8) gives the uncertainty associated with the probability density function p(m). If we have available information of global constraint S(m) for parameter (m), then the corresponding probability distribution for maximum entropy will be expressed as follows:

(9) P ( m ) = A e − S ( m ) ,

where A is the constant for normalization. The p(d|m) as a function of discrepancy between the model and the observation in generalized version can be written as follows:

(10) p ( d | m ) = p 1 − 1 / p 2 σ p Γ ( 1 / p ) exp − 1 p | d − L m | p ( σ p ) p .

The solution for l ₁-norm or l ₂-norm can be estimated if one selects p = 1 or 2 in equation (10), respectively. In this way, error between the observation and the model parameters can be reduced that maximizes the objective function Y [11]. The solution of l ₁-norm for P-impedance ξ(t) can be expressed as follows:

(11) ξ ( t ) = ξ ( 0 ) exp 2 ∫ 0 t R ( t ) d t .

From the l ₁-norm solution in equation (11), it can be observed that the sparseness in R(t) may give a blocky structure in ξ(t) [11,42]. The workflow followed for linear programming sparse-spike inversion for estimation of P-impedance within subsurface is displayed in Figure 7.

Figure 7

Sparse-spike inversion procedure adopted to estimate the P-impedance.

4 Results

4.1 Petrophysical analysis

Petrophysical interpretation of the C-sand interval for Sawan-8 well is carried out to delineate the potential gas bearing zones in reservoir (Figure 8). The results of the calculated reservoir parameters are presented in Table 1. The Sawan-08 well encountered with the reservoir interval (C-sand) at the depth level of 3,267 m, having 95 m gross thickness. A 44 m thick net-pay sand zone is delineated within the C-sand interval, which shows prolific results for hydrocarbon (gas) potential. The pay zone exhibits the rock layer of high effective porosity (15%) and saturation with hydrocarbon (58%) of economical level. High porosity values are good indicator of substantial hydrocarbon saturation within productive zones of reservoir [36].

Figure 8

Petrophysical interpretation for the estimation of reservoir parameters of the C-sand interval in Sawan-08 well.

Table 1

Petrophysical analysis of reservoir C-sand interval for Sawan-08 well used in this study.

Reservoir properties	Sawan-08
Depth range (m)	3,267–3,362
Gross thickness (m)	95
Reservoir thickness (m)	61
Pay thickness (m)	44
Average porosity (%)	16.1
Average effective porosity (%)	15
Average water saturation (%)	42
Average hydrocarbon saturation (%)	58
Average volume of shale (%)	15

4.2 3D seismic attribute analysis

4.2.1 Trace envelope

The relatively higher envelope values in time slice indicated toward a channel feature of reservoir sand. The sand channel is laterally extended from the south-west (SW) to the north-east (NE) direction covering the productive wells (Sawan-01, Sawan-07 and Sawan-08) of the Sawan gas field. The section display of the envelope clearly shows the amplitude anomalies at reservoir interface that highlighted the favorable lithological changes for hydrocarbon accumulation in the C-sand interval (Figure 9).

Figure 9

(a) Constant time slice map of trace envelope attribute indicating the predicted sand channel. (b) The corresponding seismic section of the trace envelope attribute along arbitrary line. (c) 3D map of trace envelope attribute of the targeted area.

4.2.2 RMS amplitude

The application of RMS amplitude slice shows the high values along the sand channel. The prevalence of the high amplitude values within the C-sand interval indicates the existence of relatively high porous sand unit (Figure 10). The subsurface lithologies with high porosity significantly contributed toward the accumulation of hydrocarbons and usually related to high values of RMS amplitudes [9].

Figure 10

(a) Constant time slice map of RMS amplitude attribute indicating the predicted sand channel. (b) The corresponding seismic section of the RMS amplitude attribute along arbitrary line. (c) 3D map of RMS amplitude attribute of the targeted area.

4.2.3 Sweetness

The time slice of sweetness attribute also highlighted the fluid bearing zone as the sand channel body. The presence of the hydrocarbon increases the sweetness amplitude strength and decreases the frequency content. Hence, the distinct lithological contrast at the C-sand interval is the indicator of sweet spots that shows the gas-bearing zones (Figure 11). The high sweetness values are the evidence for the presence of reservoir quality lithofacies [45].

Figure 11

(a) Constant time slice map of sweetness attribute indicating the predicted sand channel. (b) The corresponding seismic section of the sweetness attribute along arbitrary line. (c) 3D map of sweetness attribute of the targeted area.

4.3 3D Poststack inversion analysis

The synthetic seismogram is generated by the convolution of the extracted statistical wavelet with reflectivity series derived at Sawan-08 well location. The correlation of synthetic seismogram with seismic traces at Sawan-08 well is displayed in Figure 12. The synthetic traces using the statistical wavelet have shown the excellent correlation results with extracted seismic traces where the correlation coefficient value is 0.854. The inversion analysis is then applied to the seismic data using linear programming sparse-spike inversion algorithms.

Figure 12

Right panel shows correlation of estimated seismic traces (blue color) using statistical wavelet with extracted seismic traces (red color) with the coefficient value of 0.854 for synthetic seismogram generation using data of Sawan-08 well.

4.3.1 Sparse-spike inversion

The inversion analysis is performed between the inverted P-impedance and the actual P-impedance log at the Sawan-08 well location (Figure 13). In Figure 13a, log-derived (blue line) and seismic-derived (red line) P-impedance is in good agreement particularly within the C-sand interval. The correlation between the synthetic (acquired after convolution process) traces (red) and seismic traces (black) at well location gives the excellent results (correlation coefficient = 0.996) for P-impedance (Figure 13b). The estimated RMS error between original and inverted P-impedance is 0.08 and 653 (m/s) × (g/cc), respectively (Figure 13c). After performing the inversion analysis, the spatial distribution of P-impedance away from the wells is well captured by the inversion using the linear programming algorithm (Figure 14). Low P-impedance values within the C-sand interval are observed in the range of 8,000–9,700 (m/s) × (g/cc) at two-way time window of 2,240–2,280 ms, which is a good indicative of potential gas-bearing sand zone at this reservoir level. The overlying layer near the C-sand top horizon within purple-blue color range represents the high impedance anomaly having value range of 10,700–11,950 (m/s) × (g/cc) and considered as the seal rock unit of shale lithology.

Figure 13

QC analysis between sparse-spike inversion and low-frequency model of P-impedance at Sawan-08 well. (a) Log-derived (blue line), seismic-derived (red line) and initial models (black line) of P-impedance. (b) Synthetic traces (red) after convolution process and extracted seismic traces (black). (c) Traces of RMS error between original P-impedance and inverted P-impedance.

Figure 14

Sparse-spike inversion-based P-impedance along the arbitrary seismic line. The lithological variations are well captured by the linear programming algorithm. The low P-impedance anomaly is clearly visible in the C-sand interval.

A constant time slice map is extracted from the inversion result to analyze the spatial distribution of P-impedance at the desired reservoir depth level (Figure 15). A relatively lower impedance anomaly is observed along the sand channel, which is assumed to be responsible for prolific hydrocarbon (gas) accumulation in Sawan-01, Sawan-07 and Sawan-08 wells of the study area. P-impedance inversion results also indicated that Sawan-14 well is not situated at the suitable place compared to other three wells in the area. The less prolific inversion results for Sawan-14 well can be attributed toward its minimum production level of hydrocarbon (gas).

Figure 15

Constant time slice map of estimated P-impedance. The spatial distribution of its lowest values indicating toward the hydrocarbon (gas) saturation along the sand channel.

4.4 Geostatistical estimations

4.4.1 Multiattribute regression

In the application of multiattribute regression analyses, the prediction (training) error is minimized to obtain optimal operator values between the predicted and the actual target log (i.e., porosity). At the same time, most suitable attributes are selected for predictive variable in an unsupervised manner under the stepwise regression process. However, a validation test is required to constrain the number of investigated attributes. Hence, run a cross-validation process that computed the error for all the wells and repeated until the optimal number of attributes is achieved. Afterward, a plot between validation error and number of attributes is generated to identify the suitable attributes for predicting the targeted reservoir parameter. The entire dataset of internal and external attributes is randomly split into 10 folds. Table 2 presents the list of suitable attributes along with their numerical errors for predicting the porosity. The optimal number of attributes is determined from the plot of average error versus number of attributes (Figure 16). Results showed that initial eight attributes are statistically significant and utilized in the prediction process because the validation error starts to rise beyond this number. In Figure 17, the predicted porosity is cross-plotted against the actual porosity for QC check and the predicted data points are closely distributed to the actual data points and least scattered from the best fitted line of cross correlation. The correlation coefficient between the actual and predicted porosity is 0.866 with the average training error of 0.032. However, the validation result decreases to 0.773 with 0.041 error (Table 3).

Table 2

List of multiattributes from stepwise regression for predicting porosity in the reservoir.

#	Target log	Final attribute	Training error	Validation error
1	Porosity	Average frequency (inversion result)	0.04789	0.059123
2	Porosity	X-coordinate (inversion result)	0.040144	0.047124
3	Porosity	Derivative	0.034673	0.046932
4	Porosity	Integrated absolute amplitude	0.033773	0.046933
5	Porosity	Instantaneous phase	0.03304	0.045306
6	Porosity	Integrated absolute amplitude (inversion result)	0.032505	0.043155
7	Porosity	Second derivative	0.031903	0.043799
8	Porosity	Quadrature trace (inversion result)	0.031692	0.042945
9	Porosity	Filter 5/10–15/20	0.031372	0.043216
10	Porosity	Dominant frequency	0.030681	0.044052

Figure 16

Validation plot between average error and number of attributes for all wells, showing that the optimal number of attributes is eight for predicting the porosity. The error using all wells and when a well is removed are indicated by a black curve and a red curve, respectively.

Figure 17

Cross-plot of actual and predicted porosity using multiattribute regression with correlation coefficient of 0.866. Data points from the C-sand interval of each well are shown in one color.

Table 3

Summary of porosity prediction by multiattribute regression and PNN analysis, and the correlation reasonably improved in PNN analysis.

Target log	Technique	Validation result		Training result (QC)
		Cross correlation	Average error	Cross correlation	Average error
Porosity	Multiattribute regression	0.773	0.041	0.866	0.032
	Neural network (PNN)	0.810	0.038	0.896	0.029

4.4.2 PNN

For the neural network analysis, the suitable attributes are selected through the stepwise regression process and further validated to get the final attribute set for neural network training. Afterward, the constrained attributes and targeted log (porosity) at well locations are correlated to derive the statistical relationship. Finally, the derived neural network is applied to the seismic data volume. The attributes obtained in the multiattribute regression analysis are used in the neural network training process. Figure 18 displays the comparison between actual and modeled (predicted) porosity with the average training error at each well location within the C-sand interval under the PNN analysis, and the training results showed the correlation coefficient value of 0.896 with the average training error of 0.029 (Table 3). It is evident from the cross-plot of actual and predicted porosity that the PNN estimation (Figure 19) showed less scatter relative to the best fitted regression line when compared with cross-plot estimation of multiattribute regression for actual and predicted porosity (Figure 17). Higher cross correlations with less average errors indicate that the prediction efficiency improved under the PNN analysis compared to the multiattribute regression analysis.

Figure 18

Comparison between actual porosity (black curve) from well and modeled (predicted) porosity (red curve) for composite trace with average training error at each well location.

Figure 19

Cross-plot of actual and predicted porosity using PNN showing high correlation coefficient of 0.896. Data points from the C-sand interval of each well are shown in one color.

The PNN is approved as a better geostatistical approach based on the training and validation test results and applied on the seismic data to predict the porosity away from the wells within the C-sand interval. The results showed that the high porosity values are prominent in gas saturated zone of the C-sand interval (Figure 20). High porosity anomaly varies between its fractional volume range of 0.14–0.20 v/v in the seismic time window of 22,40–2,280 ms, which coincides with the low P-impedance values predicted in inversion results. The extracted time slice map from porosity volume also showed the distribution of its relatively high values along the identified channel feature (Figure 21). The favorable distribution of porosity in the direction of low P-impedance values has successfully mapped the gas saturated zone in reservoir.

Figure 20

Section display of computed porosity using the PNN analysis along the arbitrary seismic line. The high porosity zone is clearly visible in the C-sand interval, where the low P-impedance anomaly encountered in the inversion analysis.

Figure 21

Constant time slice map showcasing the porosity distribution within the reservoir interval. The high porosity values indicating the gas accumulation along a sand channel in the reservoir.

5 Discussion

Initially, the analysis starts with the extraction of attributes from the seismic data at the desired depth level within the reservoir C-sand interval. The high amplitude anomalies of attributes predict the sand channel at the targeted reservoir depth level trending in SW to NE direction. Afterward, the reservoir elastic property (P-impedance) is estimated through poststack seismic inversion to predict the promising zone of hydrocarbon (gas) accumulation. The inversion results showed that the low P-impedance values ranging from 8,000 to 9,700 (m/s) × (g/cc) between the time interval of 2,240 and 2,280 ms are good indicatives of hydrocarbon (gas) zone in the reservoir. The spatial distribution of P-impedance revealed that the relatively lower values are observed along the predicted sand channel.

To further characterize the reservoir, porosity is estimated in the inter-well region using multiattribute regression and the PNN analysis within the C-sand interval. The PNN analysis delivered more significant cross-correlation results between predicted and actual porosity compared to the multiattribute regression analysis and therefore applied on entire seismic data to generate the porosity volume. The transect and time slice display of the estimated porosity indicated that the high porosity anomaly coincided with the low P-impedance zone of the C-sand interval. The uniformity in these estimates (i.e., P-impedance and porosity) has efficiently predicted the gas saturated zones in the reservoir. The estimates also highlighted the stratigraphic feature of sand channel, which is identified as the main contributor of gas accumulation in the area.

A cross-plot is generated between the predicted porosity and inverted P-impedance estimated at constant time slice from geostatistical and inversion analysis, respectively (Figure 22). A linear regression approach is adopted to find the correlation between porosity and inverted P-impedance, as the relationship is linear between these two parameters with a negative slope. The result showed a reasonably good correlation coefficient value (i.e., 0.807). The best fitted line developed the relationship between P-impedance and porosity as follows:

where ϕ is the porosity and A is the inverted P-impedance. The porosity can be estimated directly from P-impedance and vice versa by using the estimated relationship (equation (14)) in the study region. The outcomes of petrophysical, seismic attribute and sparse-spike inversion analyses are consistent with geostatistical results in which the spatial distribution of the estimated reservoir parameter (porosity) confirmed the presence of gas sand channel (Figure 22) as a major gas saturated zone within the reservoir. Moreover, the extent of predicted gas sand channel can be further demarcated at the regional level for the C-sand interval of Lower Goru Formation by following the applied methodology of this study for the area.

Figure 22

Cross-plot of predicted porosity versus inverted P-impedance at constant time slice with the significant value of correlation coefficient. The gas sand channel zone is highlighted with a black circle.

6 Conclusions

The integrated use of advanced geophysical techniques (seismic attributes, sparse-spike inversion and geostatistical analysis) effectively estimated the gas-bearing zones of Cretaceous C-sand interval of Lower Goru Formation in Middle Indus Basin, Southern Pakistan. The high amplitude anomalies of seismic attributes predicted the sand channel at the targeted reservoir depth level trending in SW to NE direction. The sparse-spike inversion analysis successfully captured the detailed lithological variations and revealed that the relatively lower impedance values are observed along the predicted sand channel, which is good indicatives of hydrocarbon (gas) saturated zone within reservoir interval. To further characterize the reservoir, porosity is estimated in the inter-well region using multiattribute regression and the PNN analysis. The PNN analysis delivered more significant cross-correlation results between predicted and actual porosity and therefore applied on entire seismic data to generate the porosity volume. The spatial distribution of estimated porosity is in good agreement with its petrophysical estimation in Sawan-08 well. Furthermore, the current study provides the knowledge toward the estimation of reservoir properties through advanced integrated geophysical approach to robustly predict the hydrocarbon bearing zones in reservoir. The methodologies applied in this study can be executed in other global regions beside Pakistan that exhibit the similar geological conditions for better understanding and characterization of gas reservoirs.