Introduction

Carbonate rocks are considered as a major host rock for hydrocarbon reservoirs worldwide. It is approximated that 60% of the world’s oil and 40% of gas reserves are contained by carbonate reservoirs (Carrera et al. 2018). These reservoirs behave drastically different from siliciclastic reservoirs due to their depositional environment and complex diagenetic history (Anselmetti and Eberli 1993; Lucia 1995, 1999). Diagenetic (cementation, dissolution and dolomitization), differential compaction and faulting and solution collapse activities alter mineralogy and texture of rocks framework which ultimately cause reservoir rock to exhibit different petrophysical and elastic properties (Zhao et al. 2013; Eberli et al. 2003).

Lucrative exploration and effective exploitation of hydrocarbons from carbonate reservoirs demand development of new technologies and integration concepts for accurate prediction of reservoir litho-facies and petrophysical properties (e.g., clay volume and porosity) distribution within intricate corridors. Conventional qualitative interpretation techniques using well observation and seismic attributes revering geological (depositional environment and complex structures) models do not provide adequate control in effective appraisal and exploitation of these distinct reservoirs (Al Moqbel and Wang 2011). Often, over-simplistic reservoir descriptive models from qualitative interpretation lead to wrong well placement and inefficient appraisal and development plans (Vasudevan 2020). Hence, main challenge for the development of these reservoirs is associated with robust quantitative seismic-based reservoir characterization workflow to elucidate reservoir complexities. Quantitative interpretation can provide a myriad of geological and geophysical data to identify zones with better reservoir quality for right drilling location and provide cost-effective appraisal and development plans.

This study is the first ever attempt made to adopt an improvised quantitatively seismic reservoir characterization workflow to evaluate mature tight carbonate gas reservoir (of Eocene age) field (KDT) in Middle Indus Basin of Pakistan. Several vertical and horizontal wells have been drilled based on the conventional structural-based 3D seismic interpretation of subsurface dome feature. Most of the wells are single string producers (from S2 formation), whereas some are from dual producers (from S2 and S1) formations. The estimation of key reservoir properties, litho-facies and petrophysical properties (clay volume and porosity) are much more challenging because of rock complexity and heterogeneity.

Main objective of this study is to accurately characterize reservoir litho-facies (tight limestone and shale) and petrophysical (clay volume and porosity) properties. Expectations are the reliable estimation of these properties which in turn will contribute immensely toward designing a subsequent reservoir management plan by understanding spatial distribution of the productive part of reservoir and propose future field development plan.

Formation (S1) under investigation is early Eocene with overall 60-m-thick package of net gas pay, exhibiting high porosity (up to 25%) and low permeability. S1 comprised of alternating clean porous to tight, argillaceous stylonodular to variably dolomitic claystone, and dolomitized limestone stringers within claystone intervals. Several vertical and horizontal wells have been drilled for hydrocarbons, yet further development of S1 formation depends on the quantitative studies presented here.

Quantitative seismic reservoir characterization approach employed integrates litho-facies classification, AVO attribute analysis and pre-stack simultaneous inversion technique, supported with the customized well-log and seismic data (gathers) conditioning. Litho-facies classification is estimated at well-based estimates of petrophysical properties (such as gamma ray, volume of clay, water saturation and mineralogical volumes). Standard AVO attributes (intercept and gradient) are calculated as a convenient reconnaissance inspection of amplitude behavior (reflection coefficients) of the possible AVO (class I) anomaly in the area. Unfortunately, AVO analysis has led many geoscientists to conclude that it has of little practical importance due to seismic data processing issues and often been ignored to invest the effort required to understand the technique. Here, we used AVO analysis as a helping reconnaissance tool to investigate the prevailing AVO classes, fluids and lithology characteristics. Finally, geophysical data are transformed into rock elastic properties through an inverse problem. Various studies have been performed to prove the effectiveness of seismic inversion methods to characterize carbonate reservoirs and pointed out new potential hydrocarbon areas of good quality (Liu and Wang 2017; Karimpouli et al. 2013; Benabentos et al. 2007). Hence, pre-stack simultaneous inversion (Hampson et al. 2005) is applied to estimate elastic attributes (p-impedance, s-impedance, VpVs ratio and density). The interdependency of p-impedance and density attributes with petrophysical properties (porosity, water saturation and clay volume) in carbonate reservoirs is very well studied (Garia et al. 2019; Cataldo and Leite 2018; Huuse and Feary 2005). Therefore, we used p-impedance and density attributes to map the defined litho-facies (tight limestone and shale) distribution. Probable litho-facies maps are generated using Bayes theorem approach. A good match of litho-facies is achieved after cross-validation at well locations. Petrophysical properties (clay volume and porosity) are derived from Probabilistic Neural Network (PNN) approach using well logs, seismic attributes (such as instantaneous, windowed frequency, filters and derivatives) and attribute volumes from seismic inversion (p-impedance and density). We preferred neural networking approach over conventional methods (such as linear or multi-linear regression) to handle reservoir nonlinearity and achieved better correlation of clay and porosity prediction at the well locations. Results proved the validity and robustness of quantitative seismic reservoir characterization framework employed to litho-facies (tight limestone and shale) and predict petrophysical properties (volume of clay and porosity), and helped us to effectively position and plan of high-technology horizontal wells as a part of redevelopment program of the field. This study is a part of ongoing quantitative seismic reservoir characterization work of producing carbonate reservoir field for additional infill wells to improve hydrocarbon production.

Geological and stratigraphic setting

Exploration history of KDT field can be traced back to 1959 and is considered one of the matured hydrocarbons (gas) field in Lower Indus Basin of Pakistan. KDT field is located to east of Kirthar Range and south of the Sulaiman Range and extends east to Pakistan–India border and southward to Lower Indus Basin (Fig. 1). Both Kirthar and Sulaiman Ranges are fold belts, which were deformed and uplifted during the Himalayan orogeny (Shah and Malik 1995). This orogeny was result of collision of Indo-Pakistan plate with the Eurasian plate in mid-Eocene. Major deformation of Kirthar and Sulaiman fold belts began in early Oligocene with the resulting foreland basin depocenter located adjacent to fold belt. This foreland basin effectively overlies the pre-collisional Middle Indus Basin. Two principle lows are the axial depocenters of Kirthar and Sulaiman Foredeeps. Kirthar Foredeep is oriented north to south and is parallel and adjacent to mountain front of the Kirthar Range. Pre-collisional Middle Indus Basin comprises a thick sequence of Mesozoic and Lower Tertiary sediments, the youngest being of Eocene age. These sediments thickened from east to west, measuring 8–10 km beneath the frontal folds of the Sulaiman Range. Field is covered with alluvium, and there is no evidence of surface structure. However, based on seismic and gravity data, three distinct subsurface features (domes) have been identified. These separate subsurface domes lie on the east, west and southern directions of block boundaries. All the three domes have commercial gas discoveries from different stratigraphic compartments. This study focuses only on the east dome. HL limestone formation (83-m-thick) is first reservoir zone (not under investigation here) of the field containing locally abundant of argillaceous lamina along with mix minerals of argillaceous/chert (Fig. 2a). S1 reservoir unit of Eocene age is the second gas producing zone after shallower HL formation (Siddiqui 2000a). S1 unit is overlain by GS (shale) formation of early Eocene age and is capped by LA (shale) formation of Eocene age. LA formation is predominantly shale with abundant of limestone stringers and very local dolomites in the upper and lower sections. GS formation is widely pervasive consisting of roughly 50 m of thick clean shale section, separating S1 and S2 formations. S2 (limestone) formation is the third gas producing reservoir zone (which is also not under investigation here) roughly 295 meters thick. Upper 10 m section of S2 contains frequent shale stringers followed by the gas producing reservoirs. S2 is locally very argillaceous with abundant of disseminated pyrite and traces of dolomite and calcareous shale. Regional extensive GS formation provides a good quality top seal over S2 formation. SR and LG (shale) formations of early cretaceous are main source rocks in the area. Discoveries in and around the field support the presence of mature source and migration of hydrocarbon into the available traps. Depositional environment of S1 formation is marine found in the deeper part of lagoon, offshore and on immediate flank (Siddiqui 2004). It forms a layer cake reservoir with variable thickness ranging up to 65 m. Its lithology is predominantly limestone with minor thin shale stingers and rare dolomite occurring in basal part of the section. Simplistically, S1 sequences can be considered to comprise of five main lithotypes: (1) clean porous limestone, (2) clean tight limestone, (3) argillaceous stylonodular limestone, (4) variably dolomitic claystone and (5) dolomitized limestone stringers within claystone intervals (Fig. 2b). Formation thin section exhibits dirty white to milky white and light gray in parts benthonic foraminifera–echinoderm lime wackestone forming the more common lithology type with common algal lime packstone with some argillaceous content (Fig. 2c).

Fig. 1
figure 1

Location map of the study area lying in Middle Indus Basin of Pakistan

Full size image
Fig. 2
figure 2

a Stratigraphic column of the study area, b sequence stratigraphy of S1 limestone and c thin section from S1 reservoir showing pelagic wackestone—probably deposited in a very low energy deep-water base-of-slope/basin-floor environment with poor fracture and local vuggy, matrix micro- and intraparticle porosities and mixed Ovulites-foram wack–packstone—probably being deposited in the flanks of small Ovulites shoals

Full size image

Porosity development in S1 reservoir lays mainly over crest of structure where most wells of the field have been drilled. Much of the structural closure displays relatively high seismic amplitude corresponding with good porosity in the reservoirs, while deeper on the flanks seismic amplitude decreases markedly due to porosity reduction associated with diagenesis (secondary cementation).

Geophysical approach

Every quantitative seismic reservoir characterization project presents a unique problem requiring a tailored solution following the integrated approach. Quantitative interpretation (QI) can reduce exploration and development risk by extracting additional high value information by integrating the data from various disciplines. To observe seismic amplitude variation in terms of rock properties, we transformed geophysical seismic attributes into rock litho-facies and reservoir properties by posing as an inverse problem. Seismic inversion uses seismic data, well logs, re-classified litho-facies, seismic velocities and interpreted horizons to generate the rock properties volumes (P-impedances, VpVs ratio and density). Data quality (both seismic and well) is seen as an important requisition for reliable well to seismic calibration. Litho-facies (tight limestone, gas bearing limestone and shale) are defined at well location based on core sample and area geological (depositional and sedimentological) information using well-log-based estimates of petrophysical properties (such as gamma ray, volume of clay and water saturation). As a reservoir reconnaissance, hydrocarbon anomalies are assessed by amplitude variation with offset (AVO) attributes (such as AVO synthetic gathers, angle versus reflection amplitude gradient curves and intercept versus gradient plot). AVO attributes calibrate hydrocarbon anomalies effect modeled from well logs to seismic data. This calibration is also the main objective of seismic inversion to identify similar reservoir seismic signatures. Finally, an impedance model (low-frequency model) is made using calibrated wells, wavelets, seismic velocities and horizons used to simultaneously invert seismic data (angle gathers) to a rock property volume. Inverted elastic properties (impedance and density) attributes are used to estimate probability density functions (pdfs) of the reservoir facies and supervised Bayesian classification scheme generated probability cubes for the distribution of different litho-facies. Additionally, we used Probabilistic Neural Network (PNN) to estimate reservoir petrophysical properties (clay volume and porosity) using well logs, seismic attributes (instantaneous, windowed frequency, filters, derivatives, integrated, time) and attribute volumes from seismic inversion (p-impedance and density). Quantitative estimation of litho-facies and petrophysical properties helped us to facilitate the interpretation in locating sweet spots. Figure 3 shows an integrated quantitative geophysical approach process flow diagram employed in this study.

Fig. 3
figure 3

Integrated quantitative geophysical process flow diagram

Full size image

Data quality control

Quality and consistency of 3D seismic and well-log data play an important role in reservoir characterization. To have suitable inputs, both conditioning and calibration of seismic and well data are performed as required to meet the objectives of study.

Seismic data conditioning

Reservoir characterization necessitates conditioning of seismic data to be noise free, minimizes the amplitude damaging effects and time aligned across all the angles (Durrani et al. 2017; Singleton 2009). Poorly processed seismic data only focused to improve data image quality for structural and stratigraphic interpretation can severely influence the results of rock properties inversion algorithm and ultimately impact final interpretation. Various studies have shown that reservoir compliant (amplitude friendly) processed data pose a significant effect on the behavior of data being interpreted in AVO (intercept and gradient) attributes domain (Verma 2015; Coleou et al. 2013). Therefore, customized AVO compliant seismic data conditioning workflow is applied to have suitable input seismic data for better results (Fig. 4). Each stage of the workflow has been tested with several methods and parameters to have seismic data which are:

Fig. 4
figure 4

Customized AVO compliant seismic data conditioning workflow

Full size image

  • Noise and multiple free—removed random noise

  • Having true amplitude—compensate for the amplitude loss

  • Wideband and stabilized low frequencies through Fresnel frequency filtering

  • Correctly imaged and time alignment to mitigate for residual move-out

Processed seismic data can never fully satisfy these criteria; however, we tried to be as much closer to the above-mentioned objectives as possible. Because inversion algorithm tries to model seismic data by a convolution of a reflectivity series with a wavelet, so if seismic data contain noise—whether random or multiple—then inversion will attempt to model the noise. Seismic amplitudes also represent subsurface geology—if they do not, the convolution assumption is violated through processes such as automatic gain control (AGC) which is avoided in this case. Similarly, if data are not imaged correctly—for example by a less than optimum migration, then the resultant impedance volume can also be incorrectly imaged. Bandwidth of data also affects the results in two ways: High frequencies in seismic data allow inversion to resolve thinner layers, and low frequencies in seismic data reduce uncertainty associated with generating an earth model to infill these frequencies. Figure 5 shows comparison of before and after conditioning of seismic (gathers) data.

Fig. 5
figure 5

a Unconditioned gathers and b conditioned gathers. Reservoir target zone is drawn with red marker (horizon). Reservoir zone shows an area of significant improvement after data conditioning

Full size image

Figure 5a shows unconditioned gathers are spurious and contain noise (multiples) and reflection events are not time aligned, whereas after conditioning, reflection events are flat, sharp and continuous with negligible effect on amplitudes (Fig. 5b).

Well data conditioning

Well-log data are routinely being used for petrophysical evaluation, rock physics analysis, wavelet estimation, low-frequency model building, seismic velocity calibration and time-to-depth conversion. It is often the case that well-logs data do not satisfy certain criteria for accurate rock properties estimation for various reasons: poor logging conditions, different tool measurements and invasion of borehole fluids into formations and alteration of formation properties. Therefore, well-logs data need to be conditioned, corrected and synthesized to provide complete and reliable input for petrophysical interpretation, litho-facies analysis and calibration of elastic logs (density, p-sonic and often s-sonic logs) required for seismic inversion. Figure 6 shows the workflow adopted to condition well-log data which include the following measures.

Fig. 6
figure 6

Well-logs data conditioning workflow

Full size image

  • Depth shifting for inter- and intra-well consistency

  • Removing of sonic cycle skipping and spikes through normalization

  • Reconstruction of density log

  • Covered the entire vertical interval of interest

  • Covered reservoir and non-reservoir rocks

  • Represented the true, undisturbed rock properties

Three wells have been used for data conditioning to achieve high-quality and consistent dataset for the study. Log curves are analyzed for unreliable intervals and edited using a range of manual, empirical and petrophysical log editing techniques. Measured well-logs quality especially the density and sonic logs are conditioned and checked using multi-well cross-plots and histograms (Fig. 7). Plots show target log data of reservoir zone (S1) after conditioning which have low scatter distribution covering the whole range of sonic and density values.

Fig. 7
figure 7

Multi-wells (five wells data—only three wells names are disclosed due to data confidentiality) cross-plot distribution of density, sonic and neutron porosity logs after conditioning of the target reservoir zone (S1)

Full size image

For shear wave (Vs) prediction in reservoir (S1) and non-reservoir zones (GS and S2), we utilized techniques of multi-linear regression and cross-validation from a nearby well with measured shear wave log information to develop relationships using suite of other well-log curves. Those relationships are then used to predict pseudo-shear wave log curves in wells. Multi-linear Eqs. (1 to 3) are as follows:

$$ \begin{aligned} & {\text{S - wave}}\, \left( {\text{Vs}} \right)\, {\text{for}}\, {\text{S1 }}\,{\text{unit}} \\ & = 99.55 + 0.26 * {\text{Vp }}\left( {{\text{us}}/{\text{ft}}} \right) - 7.19 * {\text{Density }}\left( {{\text{g}}/{\text{cc}}} \right) \\ & \quad + 272.83 * {\text{Nphi}} - 2.85 *10^{ - 04} * {\text{LLD}} + 0.253 * {\text{MSFL}} + 0.317 * {\text{GR}} \\ \end{aligned} $$
(1)
$$ \begin{aligned} {\text{S - wave}}\, \left( {\text{Vs}} \right)\, {\text{for}}\, {\text{GS}}\, {\text{unit}} \\ & = 186.61 + 0.52 * {\text{Vp }}\left( {{\text{us}}/{\text{ft}}} \right) - 39.546 * {\text{Density }}\left( {{\text{g}}/{\text{cc}}} \right) \\ & \quad + 91.63 * {\text{Nphi}} + 0.11 * {\text{LLD}} - 0.72 * {\text{MSFL}} + 0.004 * {\text{GR}} \\ \end{aligned} $$
(2)
$$ \begin{aligned} {\text{S - wave}}\, \left( {\text{Vs}} \right)\, {\text{for}}\, {\text{S2 }}\,{\text{unit}} \\ & = 215.60 + 0.418 * {\text{Vp }}\left( {{\text{us}}/{\text{ft}}} \right) - 51.10 * {\text{Density }}\left( {{\text{g}}/{\text{cc}}} \right) + 111.399 * {\text{Nphi}} \\ & \quad + 0.17 * {\text{LLD}} - 0.17 * {\text{MSFL}} + 0.005 * {\text{GR}} \\ \end{aligned} $$
(3)

Log quality control and assessment of the workflow helped to make sure logs are run and calibrated under control (Fig. 8).

Fig. 8
figure 8

Composite display of multi-wells (KDT-04, KDT-07 and KDT-08) well logs (gamma ray (GR), compressional wave (DT) and shear wave (DTS)) after conditioning covering the reservoir (S1) and non-reservoir lithology (G1 and S2) highlighted with formation

Full size image

Facies re-classification

Facies classification is very important in reservoir characterization. In static reservoir model, distribution of reservoir properties depends on the facies re-classification, which is generally achieved from geophysical measurements such as well logs and seismic data (Lindberg and Grana 2015). In this study, facies are defined based on core sample analysis constrained with area depositional and sedimentological information. Log facies are re-arranged at well location using well-log-based estimates of petrophysical properties computed in formation evaluation analysis (such as gamma ray, volume of clay, water saturation and mineralogical volumes). Deterministic cutoff method is used to classify log facies based on a limited number of well logs defined as follows:

  • Shale—characterized with high gamma ray response (> 50 API) and contains volume of clay > 0.25

  • Gas bearing limestone—defined with volume of clay < 0.25 and water saturation < 0.70

  • Tight limestone—defined with volume of clay < 0.25 and water saturation > 0.70

Figure 9 shows facies classification are linked with petrophysical (porosity and water saturation) and elastic properties (p-impedance). We can see from cross-plots in 3D space that only tight limestone (blue) and shale (green) could be discriminated clearly from well logs and gas limestone (red) is quite challenging due to similar petrophysical (porosity and mineralogy) and acoustic (p-wave and s-wave velocity) properties.

Fig. 9
figure 9

Facies classification of three wells. 3D cross-plots of p-impedance versus water saturation (Sw) and effective porosity (Phie) of the target reservoir (S1) color coded with lithology log

Full size image

Seismic–well tie and wavelet estimation

Accurate well to seismic tie and reliable wavelet estimation are essential steps before AVO analysis and seismic inversion and provide confidence on quantitative interpretation. However, tying a well with seismic has never been a straightforward practice and can be affected by poor well and seismic data quality and geological complexities (Ziolkowski et al. 1999). Well to seismic tie involves forward modeling procedure by generating synthetic trace using sonic and density logs, then matching with the observed seismic reflection data. This procedure generates a relationship between logs (measured in depth) and seismic data (measured in time). For good well to seismic tie, additional factors need to be considered are check shot correction, wavelet estimation using suitable window and identifying horizon to be picked for the initial model. Seismic wavelet provides an important link between earth reflectivity (we are trying to determine) and seismic data (which we have). It requires reliable sonic and density logs for reflectivity generation. In this study, wavelets are estimated by two methods: (1) statistical and deterministic method described by Oldenburg et al. (1987) and (2) multi-wells deterministic wavelet. First, commonly practiced statistical method is opted to estimate wavelet from 3D seismic data which uses auto-correlation process and assumes a user-defined constant phase. Missing well-log time series are replaced by maximum entropy deconvolution of seismic records about the stationarity of seismic records. Secondly, multi-wells deterministic wavelet is estimated in the time window (800 ms of length) which relied on known well-logs reflectivity and derived wavelet by minimizing the differences between synthetic seismogram and seismic trace at the same well location. Quality and stability of the deterministic wavelet are vetted with quality control (QC) checks: seismic-to-synthetic match, wavelet shape, residuals and spectral analysis. Optimum wavelet is determined to be the one which, when convolved with well-log reflectivities, results in a synthetic that best matches with the input seismic data in a least square way (Fig. 10).

Fig. 10
figure 10

a Blue trace beside the p-wave seismic (red trace) is the synthetic zero offset p-wave response calculated from logs (Vp and density) information extracted along the well bore, b estimated deterministic wavelet and maximum cross-correlation coefficient (80%) and c well to seismic tie shown over K-04 well

Full size image

Agreement between the datasets (well and seismic picked horizons) ensures that the extracted wavelet and reflectivity data from logs produce a representation of the subsurface compatible with seismic data. As field of study has good well control, prior information is used to overcome any issue causing a mismatch between synthetic and real seismic at well control tied to the appropriate events. All three wells were tied with good cross-correlation (> 80%) observed at well location with seismic data (Fig. 11).

Fig. 11
figure 11

Multi-wells well to seismic tie analysis—arbitrary seismic line passing through three wells

Full size image

AVO attribute analysis

Assessing amplitude variation with offset (AVO) behavior in pre-stack seismic (gathers) data is an archetypal direct hydrocarbon indicator (DHI). Under different geological conditions, DHI principle relies on change in amplitude of seismic (reflection) data with angle/offset is due to fluid and/or change of reservoir lithology. Different simplified AVO approximations (Bortfeld 1961; Aki and Richards 1980; Shuey 1985; Fatti et al. 1994) to the Zoeppritz (1919) equation have been published to estimate common attributes intercept (A) and gradient (B) based on the fact that seismic amplitudes at the boundaries/interfaces (above and below) are affected by variations of physical properties. Based on Aki and Richards’s approximation, Ostrander (1984) established that AVO analysis can be used to distinguish hydrocarbon (gas)-related amplitude anomalies from non-hydrocarbon bearing types. Later, Rutherford and Williams (1989) identified three major AVO classes (I, II and III) of hydrocarbon anomalies based on acoustic impedance contrast with the surrounding rock (e.g., shale).

Preliminary analysis of high-quality and conditioned seismic data (gathers) is performed by AVO analysis to have a convenient inspection of amplitude behavior (reflection coefficients) in seismic (gather) data at different angles of incidence. We calculated standard AVO attributes (intercept and gradient) using exact solution of Zoeppritz (1919) equation to explain the increase in amplitude from shale limestone interface with offset. Zoeppritz equation uses layers of rock properties (p-wave velocity, s-wave velocity and bulk density) to explore AVO behavior of single interface between the two geological formations (shale and limestone). Using the conditioned seismic gathers and the velocity model was calibrated using well logs, a set of AVO attributes are generated to analyze the possible anomalies present in the S1 and S2 reservoirs. Figure 12a exhibits AVO response from synthetic seismic gathers at three wells (K-04, K-07 and K-08) location, AVO gradient curves at S1 and S2 levels, and intercept versus gradient cross-plot. Figure 12b shows AVO class I anomaly for S1 reservoir with positive intercept indicates higher relative impedance than overlying unit with decreasing gradient with offset/angle falling in the IV quadrant. Figure 12c exhibits cross-plot between intercept and gradient at S1 samples from three wells showing class I response. Figure 13a shows practical plot of intercept versus gradient for all the time samples at all trace location within specified window. Deviation effect from the background trends shows the effects of changes in reservoir properties on AVO classes’ response. The degree of shift in trends from the background might be due to the rock stiffness, fluids content and porosity. Different trends are highlighted with different colored ellipses indicating the top and base of AVO classes in cross-plot. Figure 13b exhibits data from the selected ellipses projected onto seismic cross section. The red oval represents data points for AVO class I response of target S1 reservoir (in red), gold oval represents AVO class IV or base of reservoir S1 and blue oval shows base of AVO class III. Figure 13c shows time slice of intercept and gradient product (AxB) AVO attribute at S1 level.

Fig. 12
figure 12

a Comparison of AVO response from synthetic seismic angle gathers at three wells location, b angle versus reflection amplitude curves showing class I response (a positive normal reflection coefficient and a negative gradient) related to S1 level with shale as caprock and c intercept versus gradient cross-plot for identifying the magnitude and class of seismic reflection. Polarity convention in this figure display denotes that an increase in acoustic impedance signifies a peak

Full size image
Fig. 13
figure 13

AVO attributes analysis: a intercept versus gradient cross-plot from seismic data highlighting effects of changes in reservoir properties on AVO response, b seismic cross section along the AVO captured (ovals) anomalies: AVO class I (in red), AVO class IV (in gold) and base of AVO class III (in blue), and c intercept and gradient product (AxB) time slice at target reservoir (S1) level

Full size image

Simultaneous pre-stack inversion

Seismic inversion converts seismic reflection data into subsurface quantitative rock property model consistent with seismic data. The estimated quantitative rock properties increase interpreter’s ability to discriminate between different lithology and fluids and, hence, result in the detailed reservoir characterization. Various seismic inversion methods are frequently being used to invert seismic attributes (e.g., p- and s-impedance, density, Vp/Vs and lame parameters) using seismic and well-log data (Inichinbia et al. 2014; Wagner et al. 2012; Goodway et al. 1997; Russell and Hampson 1997). In this study, we used Hampson et al. (2005) approach which is an extension work of Simmons and Backus (1996) and Buland and Omre (2003). The approach utilizes angle-dependent pre-stack seismic data with preserved amplitude variation information unlike the post-stacked data where amplitude variation information is with angle (or offset) are lost. Pre-stack inversion involves the following steps; (1) Establish initial p- and s-wave impedance through the extrapolation of low-frequency model (LFM) from well-log data, (2) calculate reflection coefficients of each layer with Fatti’s AVO approximation (Fatti et al. 1994), (3) obtains synthetic angle gather by convolving reflection coefficients with the stable wavelet, (4) compare response with the observed seismic data to calculate errors between them and (5) then minimize the error between synthetic gathers and observed seismic data using conjugate methods. Finally, pre-stack data are used to simultaneously invert for p-impedance, s-impedance and density from a single inversion process. Linear relationship between logarithm of p-impedance correlated with both s-impedance and density makes simultaneous inversion accurate and stable.

Inversion parameter quality control (QC)

Setting of inversion parameters is an iterative and continues process until the best set of quality control (QC) has been obtained. Standard inversion QC includes seismic signal-to-noise level, well-log normalized correlation and standard deviation of the inversion result and elastic parameter sparseness. In each QC plot, optimum parameters are selected with low residual observation. Residual is difference between input seismic and synthetic seismic generated from reflectivity and seismic wavelet. Hence, residual data are inspected for any coherent information that should be included in the inversion results. The best point to optimize inversion parameters is at well control, where well logs define the ground truth for p-impedance (Fig. 14).

Fig. 14
figure 14

a Pre-stack inverted elastic logs p-impedance (Zp), s-impedance (Zs) and density comparisons at three well locations: actual (in blue) versus inverted (in red) along with low-frequency log (in black), b cross-plot of actual versus inverted p-impedance and c cross-plot of actual versus inverted density

Full size image

After testing and QC at well location using optimal parameters, simultaneous pre-stack inversion is applied over 3D seismic survey. Figure 15 shows inverted p-impedance and density section representing physical layer properties of the subsurface and changes in elastic properties (impedance and density) can often be related to changes in specific reservoir properties such as porosity, lithology and saturation. Figure 16 shows inverted p-impedance and density time slices corresponding to delineate S1 reservoir properties. Both p-impedance and density maps indicate the presence of new sweet spot (blue ellipse).

Fig. 15
figure 15

Simultaneous pre-stack inversion results: a p-impedance cross section passing through three wells corresponding to delineate reservoir (S1) and b inverted density along the same cross section. Base map (upper left corner) shows profile of the cross section

Full size image
Fig. 16
figure 16

a p-impedance and b density time slice at reservoir (S1) level

Full size image

Litho-facies estimation

Spatial distribution, connectivity and heterogeneity of rock properties in the static reservoir model depend on the litho-facies classification, which is quantitatively achieved from geophysical measurements such as well-logs and seismic attributes data. Litho-facies help to provide well-constrained stratigraphic framework to facilitate the interpretation in locating potential hydrocarbon areas (Zhao et al. 2014: Saussus and Sams 2012). Generating realistic litho-facies maps requires integration of various amount of information from different relevant domains such as subsurface geological understanding, well data (litho-facies), petrophysics and seismic attributes (e.g., p-impedance, density or VpVs). Such sound integration is optimally achieved through a workflow that is driven in Bayesian context and allows for the estimation of litho-facies using seismic inversion results. Figure 17 shows the workflow followed for litho-facies classification using well data and seismic attributes (p-impedance and density) from seismic pre-stack inversion. First, we created facies classification using the well-logs curves corresponding to the input inversion results (p-impedance and VpVs for example). In addition, we need numbered coded lithology log which identified zones within the well-log curves themselves. We used rock physics cross-plot (p-impedance density) and defined/highlighted different lithology zones or litho-facies (hydrocarbon sands, shales and shalier sands) color coded with lithology (based on petrophysical cutoff) log in generalized satisfactory way. Instead of defining sharp cutoff for the zone using lithology log, we statistically identified the probability distribution function (pdf) for each zone—pdf of each zone has maximum value within the cluster of points and falls off gradually as we had gone farther away from the clustered points (defined for the zone). We cross-validated and checked with zone with seismic, and this probability distribution is mapped to the inversion using “Bayesian classification.” In this way, we not identified best litho-faceis zone but also the probability of occurrence of each litho-facies. From the well-log data, we extracted the training data which are basically a cross-plot of elastic attribute color coded with lithology log. From that extracted data, we calculated probability distribution functions (pdfs). Good separation of those pdfs indicates high ability to separate the defined classes where they overlap shows the non-uniqueness or difficulties (hydrocarbon sand facie in our case which we did not incorporate) which is helped by the mathematical process. From the probability distributions, we calculate at every single point the probability of each class (shale and tight limestone). The mathematical tool which is used to calculate the probability distribution is Bayes theorem. The probability of the certain class given a particular set of seismic attributes is calculated as follows:

$$ p \left( {c_{i} |X} \right) = \frac{{\left[ {p\left( {X |c_{i} } \right)* p\left( {c_{i} } \right)} \right]}}{p\left( X \right)} $$
(4)

where \( c \) is a class (e.g., limestone or shale), \( X \) represents a seismic attributes vector (e.g., Zp, density), \( p\left( c \right) \) is a prior probability for class (e.g., probability of getting limestone in general), \( p(X|c) \) is the probability of attributes \( X \) knowing class \( c \) (e.g., distribution of Zp, density in limestone), and \( p\left( X \right) \) is the probability of attributes \( X \). Instead of defining arbitrary sharp cutoffs for lithological zones, multi-dimensional probability distribution functions (pdfs) are computed from p-impedance versus density attributes cross-plot using kernel density estimation method which finds continuous density function using the following equation:

$$ F\left( {Xx} \right) = \frac{1}{{nh \mathop \sum \nolimits_{i = 1}^{n} K \left[ {\frac{{\left( {x - x_{i} } \right)}}{h}} \right]}} $$
(5)

where \( K \) is the kernel function, \( x_{i} \) represent the seismic attributes (e.g., Zp) and \( h \) is the parameter which controls the smoothing of the distribution.

Fig. 17
figure 17

Schematic workflow for the litho-facies prediction from seismic inversion attributes

Full size image

Next, pdfs are fitted to each cluster of data points on the cross-plot which represents the mathematical measure of the overlap and degree of separation between the classes (Fig. 18a). Cross-validation is exercised by calculating probability of occurrence through confidence matrix which provides the information about the classifier success (Fig. 18b). The confidence matrix exhibits the probability of getting the target facies (tight limestone and shale) and shows the actual tight limestone (82.50%) and shale (94%) presence on upscaled log. Lithology discrimination using volume of clay (Vcl) as cutoff is consistent to all reservoir. Hence, lithology separation between tight limestone and shale in p-impedance and density cross-plots is quite clear except some of the middle part which has overlapping mixed acoustic (p- and s-wave velocity and density) and mineralogy (e.g., shaly limestone/argillaceous). Figure 18c shows good prediction of litho-log classification calibrated between well logs with no abnormal mismatch, error or confusion between the two defined classes. Finally, the pdfs are applied sample by sample to the input 3D inverted elastic attributes and the likelihood of each litho-facies given the inverted attributes is calculated at each sample location. This has yielded a 3D volume for each defined litho-facies type (Fig. 19). Figure 19 shows litho-facies section provides us a good idea about the reservoir (S1) gross thickness (increasing toward northeastern side) consistent with the wells data calibration.

Fig. 18
figure 18

a Accumulative multivariate pdf calculated from Zp versus density cross-plot opted for the litho-facies classification schemes, b confidence matrix provided the probability of each defined facies and (c) QC plot of litho-facies classification calibrated at three wells locations: (1) well logs (actual), (2) upscaled (actual), (3) well logs (computed) and (4) upscaled (computed)

Full size image
Fig. 19
figure 19

Litho-facies (tight limestone and shale) cross section along the three wells

Full size image

Reservoir petrophysical property estimation

Generally, petrophysical properties (clay volume and porosity) are estimated using well logs and seismic attributes. An optimum linear transform is generated well data to relate the petrophysical property from logs and elastic property (e.g., p-impedance and density) obtained from seismic inversion. Early investigations had utilized linear and multi-linear regression relations to establish links between elastic properties and petrophysical properties as suggested by some acoustic propagation theories (Lorenzen 2018; Draper and Smith 1966; Hampson et al. 2001). However, results are considered with skepticism when it came to suggestion that they are applied on seismic data and determine a nonlinear relationship between seismic properties and well-log information. In spite of the obvious success of multi-linear regression relations in many applications, the acute variability associated with the subsurface reservoir properties variations limits the reliability of multi-linear regression results (Russell 2004: Maurya and Singh 2019). However, nonlinear models based on data-driven approaches (such as artificial neural network) are used to handle complex problem without having prior knowledge of the process (Herrera et al. 2006; Masters 1995). Therefore, we preferred to use the Probabilistic Neural Network (PNN) approach to estimate petrophysical properties (clay volume and total porosity) which can solve nonlinear relationships. Neural network approach determined the nonlinear relationship between seismic attributes and well-log properties. This relationship is used to calculate the well-log property at every location throughout the seismic data by using external pre-stack inverted elastic (p-impedance and density) attributes and sample-based internal seismic (instantaneous, windowed frequency, filters, derivatives, integrated and time) attributes. Petrophysical properties results are cross-validated at well locations and found that PNN approach provided reliable estimation of the petrophysical properties (clay volume (VCL) and total porosity) of interests with the task of minimizing the potential error which is otherwise difficult to obtain. Figure 20 shows VCL prediction through PNN with higher correlation coefficient (≈ 99%) achieved from well log. Figure 21 exhibits total porosity prediction through PNN with good correlation coefficient (94%) and an average error of 1.05% estimated from well logs. Figure 21b shows good porosity (22%) limestone (S1) reservoir is sitting on the flank position of the anticline (marked with blue oval).

Fig. 20
figure 20

Application of Probabilistic Neural Network (PNN) volume of clay (VCL) estimation: a PNN training results QC for VCL prediction at wells showing total correlation coefficient (≈ 99%) achieved between actual well-log (in black) and estimated curve (in red), b VCL cross section along three wells (K-04, K-07 and K-08) and c VCL cross section along the three wells with cutoff (VCL < 0.25) for limestone shown in olive green and G1 shale (VCL > 0.25) shown in aqua color

Full size image
Fig. 21
figure 21

Application of Probabilistic Neural Network (PNN) for total porosity estimation: a PNN training results QC for total porosity prediction at wells showing total correlation (≈ 99%) between actual well-log (in black) and estimated curve (in red) with average error of 1.05%, b predicted total porosity cross section passing through three wells and c total porosity data slice at target reservoir (S1) level

Full size image

Conclusion

To characterize tight carbonate gas reservoirs (of early Eocene age), quantitative geophysical approach is implemented to estimate reservoir litho-facies and petrophysical properties from a mature gas field in the Middle Indus Basin of Pakistan. Quantitative seismic reservoir characterization methodology integrated well-based litho-facies classification, AVO attributes analysis and pre-stack simultaneous inversion attributes (p-impedance and density), supported with customized well-log and seismic data conditioning to achieve reliable results. Well-based deterministic cutoff method is used to define log facies at well location using well-log-based estimates of petrophysical properties. AVO attributes (intercept and gradient) conveniently provided inspection of amplitude behavior (reflection coefficients) and the presence of possible AVO (class I) anomaly in Eocene (S1) reservoir. Probable litho-facies (tight limestone and shale) are estimated based on the precise analysis of well-data-constrained litho-facies classification and inverted seismic attributes (P-impedance and density) using Bayes approach. Litho-facies section provides us a good idea about the reservoir gross thickness (increasing toward northeast) and is consistent with the wells data. The spatial distribution of the petrophysical property (clay volume and porosity) is preferably achieved through Probabilistic Neural Network (PNN) approach using well logs and pre-stack inverted attributes (p-impedance and density) constrained with sample based seismic attributes (instantaneous, windowed frequency, filters, derivatives, integrated and time). Accurate estimation of petrophysical properties at Eocene (S1 formation) level helped to identify sweet spot presence which subsequently can help in reservoir field development plan by understanding the spatial distribution of the productive zones of the reservoir in the area.