## Abstract

The integrated study of seismic attributes and inversion analysis can provide a better understanding for predicting the hydrocarbon-bearing zones even in extreme heterogeneous reservoirs. This study aims to delineate and characterize the gas saturated zone within the reservoir (Cretaceous C-sand) interval of Sawan gas field, Middle Indus Basin, Pakistan. The hydrocarbon bearing zone is well identified through the seismic attribute analysis along a sand channel. The sparse-spike inversion analysis has efficiently captured the variations in reservoir parameter (P-impedance) for gas prospect. Inversion results indicated that the relatively lower P-impedance values are encountered along the predicted sand channel. To further characterize the reservoir, geostatistical techniques comprising multiattribute regression and probabilistic neural network (PNN) analysis are applied to predict the effective porosity of reservoir. Comparatively, the PNN analysis predicted the targeted property more efficiently and applied its estimations on entire seismic volume. Furthermore, the geostatistical estimations of PNN analysis significantly predicted the gas-bearing zones and confirmed the sand channel as a major contributor of gas accumulation in the area. These estimates are in appropriate agreement with each other, and the workflow adopted here can be applied to various South Asian regions and in other parts of the world for improved characterization of gas reservoirs.

## 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).

## 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].

## 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.

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.

### 3.1 3D seismic attributes

The attributes selected for reservoir structural and stratigraphic evaluation are trace envelope, root-mean-square (RMS) amplitude and sweetness. These are extracted at constant time slice (45 ms below the C-sand top horizon) within the reservoir interval. The resultant attributes are also displayed on the 3D map along with the corresponding vertical section of the arbitrary seismic line.

The trace envelope of the original seismic trace is an analytical signal and can be presented as a complex-valued function *x*(*t*) given by the modulus at each time *t*.

where Re*x*(*t*) is the real part of the analytical signal from seismic trace and Im*x*(*t*) is the imaginary part computed by taking Hilbert transform of the real part [37]. The RMS amplitude can be computed for a set of *n* samples {*s*
_{1}, *s*
_{2}
*…s*
_{
n
}} of stack data as the sum of square root of all trace values *x* squared [38] and are presented in equation 2.

Sweetness *s*(*t*) can be computed by dividing the trace envelope *e*(*t*) with the square root of instantaneous frequency *f*
_{a} of the stacked seismic data [39].

### 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:

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).

#### 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*
_{1}, *x*
_{2}, …, *x*
_{
n
}), model *m* = (*r*
_{1}, *r*
_{2}, …, *r*
_{
n
}) and noise “*n*” is presented as follows:

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):

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:

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:

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:

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:

The solution for *l*
_{1}-norm or *l*
_{2}-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*
_{1}-norm for *P*-impedance *ξ*(*t*) can be expressed as follows:

From the *l*
_{1}-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.

### 3.3 Geostatistical techniques

Geostatistical methods provide linkage between petrophysical rock properties and the seismic data. The methods develop a statistical relationship between seismic-derived and petrophysical properties using linear regression approach [43]. A geostatistical approach is adopted to estimate the effective porosity (porosity) from seismic data as internal attributes and inversion results as external attributes. 3D time slices from seismic attributes and estimated properties are incorporated to better predict the lateral extend of gas-bearing zones at the targeted depth (constant time slice) within the reservoir interval. The geostatistical methods applied in this study include multiattribute regression and probabilistic neural network (PNN) analysis for predicting the targeted rock property (porosity) of the reservoir.

In the multiattribute analysis, numerous seismic attributes are integrated simultaneously to estimate their relationship with the targeted parameter by following the linear regression approach [13]. The attributes that give the optimal correlation with the porosity are selected and cross-plotted to estimate the porosity distribution for the particular time sample. Mathematically, modeling of porosity *ϕ*(*t*) parameter via *n* number of extracted attributes

where
*ϕ*(*t*) as follows:

The mean-squared error *E*
^{2} is calculated by the solution of derived weight of attributes using the least squares method [13]:

The application of the neural network algorithm gives the robust estimation of input attributes and targeted property by following the nonlinear approach [13,44]. Therefore, the PNN analysis as the most effective neural network technique is applied to predict the porosity. For a given series of training samples {*A*
_{1i
}, *A*
_{2i
}, *A*
_{3i
}, *ϕ*
_{
i
}}, PNN linearly combines with the log values of the training set and expresses in new output log value having vector *x* as follows:

where *D* (*x*, *x*
_{
i
}) represents the distance between the training points (*x*
_{
i
}) and point (*x*) to be estimated and can be defined as follows:

where *σ*
_{
j
} is the smoothing parameter that responds uniquely for every attribute. The network analysis aims to minimize the validation error for obtaining the smoothing parameter. Thus, for the *K*th target sample, the validation result can be expressed as follows:

In this study, the attributes acquired from multiattribute regression are used for the PNN training set. The PNN analysis selected the optimal attributes using stepwise regression, and after the validation procedure, the PNN analysis derived the statistical relationship among seismic attributes and targeted (porosity) log at well locations. The derived network is then applied to the available seismic data cube.

## 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].

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).

#### 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].

#### 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].

### 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.

#### 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.

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).

### 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).

# | 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 |

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.

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.

## 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
*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.

### 5.1 Significance of the study to gas reservoirs of South Asian Regions

The gas-bearing reservoirs have been investigated in various Basins of Asia by using seismic attributes and poststack inversion techniques. The study performed by Anees et al. [46] is based on paleo sand channels identification in tight gas sand reservoir of Lower Shihezi Formation of northern Ordos Basin, China. Their study successfully recognized the horizontal distribution and geometry of the multiple gas-bearing sand channels through the application of seismic attributes in reservoir. The study by Ali et al. [47] also presented the substantial results in predicting the hydrocarbon (gas) saturated zones by utilizing the seismic attributes and poststack inversion analysis in the secondary reservoir of Eocene carbonates of Lower Indus Basin, Pakistan. Another study by Ali et al. [23] effectively evaluated the reservoir properties of gas sand in Lower Goru Formation using different poststack inversion algorithms along with the geostatistical method. The study by Sinha and Mohanty [48] also presented the reliable estimation of reservoir parameters through poststack inversion algorithms for delineating the gas saturated zones in Krishna Godavari Basin of India. The Cretaceous sand/shale intervals of this Basin have similar characteristics to the lithological intervals of Indus Basin in Sothern Pakistan [49]. Similar results are achieved in the study by Karbalaali et al. [50] from the inversion analysis for Ghar Formation of Persian Gulf Oilfield, Iran. Their study demonstrated the cross-plots of inverted P-impedance from poststack inversion results for predicting the productive hydrocarbon (gas) zones.

Porosity is a most important parameter for prediction of fluid (hydrocarbon) in pore space of reservoir rock volume. However, at the same time, its estimation is difficult due to its variable types and lithological heterogeneity of the reservoir. The porosity along with other parameters (i.e., water saturation, shale volume) is usually measured at particular well locations and represents only vertical variations in a limited area. The aim of seismic inversion is not only to determine the P-impedance model but also to provide the estimation of reservoir parameters laterally away from the wells by in cooperating the geostatistical approach [51,52]. In this study, the geostatistical analysis provides the stepwise estimation of gas reservoir and spatial distribution of its parameters over the whole reservoir surface. This integrated study approach may be considered as advantageous work to delineate the gas saturated zones in Basins of South Asian Regions and in other parts of the world with similar reservoir characteristics.

## 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.

## Acknowledgments

The authors would like to give special recognition to the Directorate General of Petroleum Concessions (DGPC) Pakistan, for the release of 3D seismic and well log data to accomplish the research work. We also deeply acknowledge the Department of Earth Sciences, Quaid-i-Azam University for providing basic research requirements for research students and access to use the relevant software.

## References

[1] Hong Y, Wang H, Tian W, Mo A. Identification and Prediction of the High Heterogeneous Channel Sand in Southern Turgay Basin. Pittsburgh, Pennsylvania: AAPG Annual Convention and Exhibition; 2013. p. 10515.Search in Google Scholar

[2] Li Q, Wu W, Yu S, Kang H, Tong L, Cao X, et al. The application of three-dimensional seismic spectral decomposition and semblance attribute to characterizing the deepwater channel depositional elements in the Taranaki Basin of New Zealand. Acta Oceanol Sin. 2017;36:79–86. 10.1007/s13131-017-1113-0.Search in Google Scholar

[3] Rizwan M, Akhter G, Mustafa A, Bin Nisar U, Ashfaq K. Amplitude versus offset (AVO) modelling and analysis for quantitative interpretation of porosity and saturation: a case study for Sawan gas field, middle Indus basin, Pakistan. Geofis Int. 2018;57(2):151–60.Search in Google Scholar

[4] Ashraf U, Zhu P, Yasin Q, Anees A, Imraz M, Mangi HN, et al. Classification of reservoir facies using well log and 3D seismic attributes for prospect evaluation and field development: a case study of Sawan gas field, Pakistan. J Pet Sci Eng. 2019;175:338–51. 10.1016/j.petrol.2018.12.060.Search in Google Scholar

[5] Coren F, Volpi V, Tinivella U. Gas hydrate physical properties imaging by multiattribute analysis-black ridge bsr case history. Mar Geol J. 2001;178:197–210. 10.1016/S0025-3227(01)00156-6.Search in Google Scholar

[6] Taner MT, Berge T, Walls JD, Smith M, Taylor G, Dumas D, et al. Well log calibration of Kohonen-classified seismic attributes using Bayesian logic. J Petro Geo. 2011;24(4):405–16. 10.1111/j.1747-5457.2001.tb00683.x.Search in Google Scholar

[7] Maleki S, Ramazi HR, Gholami R, Sadeghzadeh F. Application of seismic attributes in structural study and fracture analysis of DQ oil field. Iran Egypt J Pet. 2015;24:119–30. 10.1016/j.ejpe.2015.05.001.Search in Google Scholar

[8] Chopra S, Marfurt KJ. Emerging and future trends in seismic attributes. Lead Edge. 2008;27(3):298–318.10.1190/1.2896620Search in Google Scholar

[9] Sarhan MA. The efficiency of seismic attributes to differentiate between massive and non-massive carbonate successions for hydrocarbon exploration activity. NRIAG J Astro Geophy. 2017;6(2):311–25. 10.1016/j.nrjag.2017.06.003.Search in Google Scholar

[10] Krebs JR, Anderson JE, Hinkley D, Neelamani R, Lee Sunwoong, Baumstein A, et al. Fast full-wave field seismic inversion using encoded sources. Geophys WCC177–WCC188. 2009;74(6):177–88. 10.1190/1.3230502.Search in Google Scholar

[11] Li Q. LP sparse spike impedance inversion. CGG Veritas: Hampson-Russell Software Services Ltd, CSEG; 2001.Search in Google Scholar

[12] Maurya SP, Singh NP. Characterising sand channel from seismic data using linear programming (L1-norm) sparse spike inversion technique: a case study from offshore. Can Explo Geophy. 2019;50(4):449–60. 10.1080/08123985.2019.1606206.Search in Google Scholar

[13] Hampson D, Schuelke J, Quirein J. Use of multi-attribute transforms to predict log properties from seismic data. Geophy. 2001;66:220–36. 10.1190/1.1444899.Search in Google Scholar

[14] Kumar R, Das B, Chatterjee R, Sain K. A methodology of porosity estimation from inversion of post-stack seismic data. J Nat Gas Sci Eng. 2016;28:356–64. 10.1016/j.jngse.2015.12.028.Search in Google Scholar

[15] Zahmatkesh I, Kadkhodaie A, Soleimani B, Golalzadeh A, Azarpour M. Estimating Vsand and reservoir properties from seismic attributes and acoustic impedance inversion: A case study from the Mansuri oilfield, SW Iran. J Pet Sci Eng. 2018;161:259–74. 10.1016/j.petrol.2017.11.060.Search in Google Scholar

[16] Zaigham N, Mallick K. Prospect of hydrocarbon associated with fossil-rift structures of the southern Indus Basin. Pak AAPG Bull. 2000;84(11):1833–48. 10.1306/8626C3A7-173B-11D7-8645000102C1865.Search in Google Scholar

[17] Droz L, Bellaiche G. Seismic Facies and Geologic Evolution of the Central Portion of the Indus Fan. Seismic Facies and Sedimentary Processes of Submarine Fans and Turbidite Systems. New York: Springer; 1991. p. 383–402.10.1007/978-1-4684-8276-8_21Search in Google Scholar

[18] Quadri VN, Shuaib M. Hydrocarbon prospects of Southern Indus Basin. AAPG Bull. 1986;70:730–47.Search in Google Scholar

[19] Ali A, Hussain M, Rehman K, Toqeer M. Effect of shale distribution on hydrocarbon sands integrated with anisotropic rock physics for AVA modelling: a case study. Act Geophys. 2016;64(4):1139–63.10.1515/acgeo-2016-0041Search in Google Scholar

[20] Anwer HM, Alves TM, Ali A, Zubair. Effects of sand-shale anisotropy on amplitude variation with angle (AVA) modelling: the Sawan Gas Field (Pakistan) as a key case-study for South Asia’s sedimentary Basins. J Asian Earth Sci. 2017;147:516–31. 10.1016/j.jseaes.2017.07.047.Search in Google Scholar

[21] Akhter G, Ahmed Z, Ishaq A, Ali A. Integrated interpretation with Gassmann fluid substitution for optimum field development of Sanghar area, Pakistan: a case study. Arab J Geosci. 2015;8(9):7467–79. 10.1007/s12517-014-1664-8.Search in Google Scholar

[22] Azeem T, Yanchun W, Khalid P, Xueqing L, Yuan F, Lifang C. An application of seismic attributes analysis for mapping of gas bearing sand zones in the Sawan gas field Pakistan. Acta Geod Geophys. 2016;51:723–44. 10.1007/s40328-015-0155-z.Search in Google Scholar

[23] Ali A, Alves TM, Saad FA, Ullah M, Toqeer M, Hussain M. Resource potential of gas reservoirs in South Pakistan and adjacent Indian subcontinent revealed by post-stack inversion techniques. J Nat Gas Sci Eng. 2018;49:41–55. 10.1016/j.jngse.2017.10.010.Search in Google Scholar

[24] Ashraf U, Zhu P, Yasin Q, Anees A, Imraz M, Mangi HN, et al. Classification of reservoir facies using well log and 3D seismic attributes for prospect evaluation and field development: A case study of Sawan gas field, Pakistan. J Pet Sci Eng. 2019;2019(175):338–51.10.1016/j.petrol.2018.12.060Search in Google Scholar

[25] Seeber L, Quittmeyer RC, Armbruster JG. Seismotectonics of Pakistan: a review of results from network data and implications for the Central Himalayas. Geol Bull Univ Peshawar. 1980;13:151–68.Search in Google Scholar

[26] Afzal J, Kuffner T, Rahman A, Ibrahim M, Seismic and well-log based sequence stratigraphy of the early Cretaceous, Lower Goru “C” sand of the Sawan gas field, middle Indus Platform, Pakistan. In Proceedings, Society of Petroleum Engineers (SPE) Annual Technical Conference, Islamabad, Pakistan. 2009.Search in Google Scholar

[27] Ahmed AR, Ahmad M, Rehman A. Comparison of core/log and well test permeabilities – a closer look sawan tight sands. Annual technical conference. Soc Pet Eng. 2010;142836:1–6. 10.2118/142836-MS.Search in Google Scholar

[28] Kazmi AH, Jan MQ. Geology and Tectonics of Pakistan. Karachi, Pakistan: Graphic Publishers; 1997.Search in Google Scholar

[29] Kazmi AH, Rana RA. Tectonic map of Pakistan. Geol Sur Pak, Quetta Scale. 1982;1(2).Search in Google Scholar

[30] Kazmi AH. Active fault system in Pakistan. In: Farah A, DeJong KA, editors. Geodynamics of Pakistan. Quetta: Geol. Sur. Pak; 1979. p. 285–94Search in Google Scholar

[31] Khalid P, Qayyum F, Yasin Q. Data-driven sequence stratigraphy of the cretaceous depositional system, Punjab Platform, Pakistan. Surv Geophys. 2014;35:1065–88. 10.1007/s10712-014-9289-8.Search in Google Scholar

[32] Wandrey CJ, Law BE, Shah HA. Sembar Goru/Ghazij composite total petroleum system, Indus and Sulaiman-Kirthar geologic provinces, Pakistan and India. US Geol Sur Bull 2208-C. 2004. http://pubs.usgs.gov/bul/b2208-c.Search in Google Scholar

[33] Azeem T, Yanchun W, Khalid P, Xueqing L, Yuan F, Lifang C. An application of seismic attributes analysis for mapping of gas bearing sand zones in the Sawan gas field Pakistan. Acta Geod Geophys. 2016;51:723–44. 10.1007/s40328-015-0155-z.Search in Google Scholar

[34] Berger A, Gier S, Krois P. Porosity-preserving chlorite cements in shallow-marine volcaniclastic sandstones: evidence from cretaceous sandstones of the Sawan gas field. Pak Am Assoc Petro Geol Bull. 2009;93:595–615. 10.1306/01300908096.Search in Google Scholar

[35] Munir K, Iqbal MA, Farid A, Shabih SM. Mapping the productive sands of Lower Goru Formation by using seismic stratigraphy and rock physical studies in Sawan area, southern Pakistan: A case study. J Pet Explor Prod Technol. 2011;1:33–42. 10.1007/s13202-011-0003-9.Search in Google Scholar

[36] Ali A, Kashif M, Hussain M, Siddique J, Aslam I, Ahmed Z. An integrated analysis of petrophysics, cross-plots and Gassmann fluid substitution for characterization of Fimkassar area, Pakistan: A case study. Arab J Sci Eng. 2015;40:181–93. 10.1007/s13369-014-1500-1.Search in Google Scholar

[37] Taner MT, Schuelke JS, Doherty R, Baysal E. Seismic attributes revisited. SEG, Expanded Abstr. 1994;1104–6. 10.1190/1.1822709.Search in Google Scholar

[38] Koson S, Chenrai P, Choowong M. Seismic attributes and their applications in seismic geomorphology. Bull Earth Sci Thail. 2014;6:1–9.Search in Google Scholar

[39] Radovich BJ, Oliveros RB. 3D sequence interpretation of seismic instantaneous attributes from the Gorgon Field. Lead Edge. 1998;17:1286–93. 10.1190/1.1438125.Search in Google Scholar

[40] Cooke D, Cant J. Model-based seismic inversion: comparing deterministic and probabilistic approaches. CSEG Rec. 2010;35:28–39.Search in Google Scholar

[41] Oliveira SAM, Lupinacci WM. L1 norm inversion method for deconvolution in attenuating media. Geophy Prospec. 2013;61(4):771–7. 10.1111/1365-2478.12002.Search in Google Scholar

[42] Sacchi M, Ulrych T. High-resolution velocity gathers and offset space reconstruction. Geophy. 1995;60(4):1169–77.10.1190/1.1443845Search in Google Scholar

[43] Doyen PM. Porosity from seismic data: A geostatistical approach. Geophy. 1988;53(10):1263–75. 10.1190/1.1442404.Search in Google Scholar

[44] Calderon JE, Castagna JP. Porosity and lithology estimation using rock physics and multi- attribute transforms in Balcon field, Colombia. Lead Edge. 2007;26:142–50. 10.1190/1.2542439.Search in Google Scholar

[45] Hart BS. Channel detection in 3-D seismic data using sweetness. AAPG Bull. 2008;92(6):733–42. 10.1306/02050807127.Search in Google Scholar

[46] Anees A, Shi W, Ashraf U, Xu Q. Channel identification using 3D seismic attributes and well logging in lower Shihezi Formation of Hangjinqi area, northern Ordos Basin, China. J Appl Geophys. 2019;163:139–50.10.1016/j.jappgeo.2019.02.015Search in Google Scholar

[47] Ali A, Younas M, Ullah M, Hussain M, Toqeer M, Bhatti AS, et al. Characterization of secondary reservoir potential via seismic inversion and attribute analysis: A case study. J Pet Sci Eng. 2019;178:272–93. 10.1016/j.petrol.2019.03.039.Search in Google Scholar

[48] Sinha B, Mohanty PR. Post stack inversion for reservoir characterization of KG Basin associated with gas hydrate prospects. J Ind Geophys Union. 2015;19(2):200–4.Search in Google Scholar

[49] Gupta SK. Basin architecture and petroleum system of Krishna Godavari Basin, east coast of India. Lead Edge. 2006;25(7):830–7. 10.1190/1.2221360.Search in Google Scholar

[50] Karbalaali H, Shadizadeh SR, Riahi MA. Delineating hydrocarbon bearing zones using elastic impedance inversion: A Persian Gulf example. Iran J Oil Gas Sci Tech. 2013;2(2):8–19. 10.22050/IJOGST.2013.3534.Search in Google Scholar

[51] Leite EP, Vidal AC. 3D porosity prediction from seismic inversion and neural networks. Comp Geosci. 2011;37(8):1174–80.10.1016/j.cageo.2010.08.001Search in Google Scholar

[52] Rijks EJK, Jauffred JCEM. Attribute extraction: an important application in any detailed 3D interpretation study. Lead Edge. 1991;10:11–9. 10.1190/1.1436837.Search in Google Scholar

[53] Krois P, Mahmood T, Milan G. Miano Field, Pakistan, a case history of model driven exploration Proc. Pakistan Petroleum Convention. Islamabad: Pakistan Assoc. Petro. Geo; 1998. p. 111–31.Search in Google Scholar

**Received:**2020-01-23

**Revised:**2020-05-16

**Accepted:**2020-06-09

**Published Online:**2021-01-27

© 2021 Muhammad Rizwan Mughal and Gulraiz Akhter, published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.