Quantitative Characterization of Shallow Marine Sediments in Tight Gas Fields of Middle Indus Basin: A Rational Approach of Multiple Rock Physics Diagnostic Models

Abstract

For the successful discovery and development of tight sand gas reserves, it is necessary to locate sand with certain features. These features must largely include a significant accumulation of hydrocarbons, rock physics models, and mechanical properties. However, the effective representation of such reservoir properties using applicable parameters is challenging due to the complicated heterogeneous structural characteristics of hydrocarbon sand. Rock physics modeling of sandstone reservoirs from the Lower Goru Basin gas fields represents the link between reservoir parameters and seismic properties. Rock physics diagnostic models have been utilized to describe the reservoir sands of two wells inside this Middle Indus Basin, including contact cement, constant cement, and friable sand. The results showed that sorting the grain and coating cement on the grain’s surface both affected the cementation process. According to the models, the cementation levels in the reservoir sands of the two wells ranged from 2% to more than 6%. The rock physics models established in the study would improve the understanding of characteristics for the relatively high V_p/V_s unconsolidated reservoir sands under study. Integrating rock physics models would improve the prediction of reservoir properties from the elastic properties estimated from seismic data. The velocity–porosity and elastic moduli-porosity patterns for the reservoir zones of the two wells are distinct. To generate a rock physics template (RPT) for the Lower Goru sand from the Early Cretaceous period, an approach based on fluid replacement modeling has been chosen. The ratio of P-wave velocity to S-wave velocity (V_p/V_s) and the P-impedance template can detect cap shale, brine sand, and gas-saturated sand with varying water saturation and porosity from wells in the Rehmat and Miano gas fields, both of which have the same shallow marine depositional characteristics. Conventional neutron-density cross-plot analysis matches up quite well with this RPT’s expected detection of water and gas sands.

1. Introduction

Seismic reflections of rock properties play an imperative role in seismic reservoir characterization and require integrating knowledge and concepts from different disciplines [1,2,3,4,5]. The seismic reflection method is a cost-effective and efficient way to characterize the reservoirs representing an indirect measurement acquired by the propagation of sound waves through various subsurface formations [6,7,8]. Understanding the seismic response of various folded structures is crucial to interpret seismic profiles in terms of tectonics and geology [9,10,11]. This type of interpretation of the seismic profiles is called structural interpretation, which implies the mapping of the geologic features of the subsurface to configurations such as folds, faults, and fractures [12,13]. The difference in rock elastic properties in the subsurface determines seismic reflections [14]. Elastic properties of rocks, that is, P- and S-wave velocities and bulk density, depend in turn on the porosity, lithology, rock texture, and pore fluid [15,16]. Acoustic impedance, the product of bulk density and velocity, is a common result of seismic inversion [17]. By itself, acoustic impedance is meaningless to the geologist and engineer. However, once it is interpreted in terms of reservoir rock properties, namely pore fluid, porosity, and lithology, it transforms into vital parameters for reservoir evaluation and drilling decisions [18].

Due to this link between geological properties and the seismic response, it is possible to transform the information from wells using seismic amplitudes through a method or a function [19]. This integration implies a transformation based on inverse and forward modeling methods. The rock physics model (RPM) is a transformation used to interpret a seismic response into reservoir rock properties and vice versa using inverse and forward modeling [15,20]. In other words, inverse modeling aims to reconstruct a geologic model using seismic data. In contrast, forward modeling is applied when core and well data are available to derive seismic signatures of rock properties for a particular geologic model. This type of interpretation of seismic reflections is called a quantitative seismic interpretation that quantifies the petrophysical parameters of the subsurface for each rock unit within the geometrical interpretation framework.

The only way to establish relationships between various rock attributes and rock physics models is via controlled experiments by simple cross-plotting between velocity and porosity [21,22]. Controlled experiments are conducted at data factories which can be categorized into three types of data: laboratory, well (wireline), and computational rock physics. Rock physics diagnostics (RPD) is a methodology for establishing RPMs that describe the relationship between seismic velocity and porosity [23]. To perform RPD, reservoir rock properties should be examined from the wireline logging dataset or laboratory core data analysis [24]. Certainly, the former approach has its advantages, such as in situ measurements of rock properties, a large number of data points, and high vertical resolution. However, one faces difficulties interpreting lateral changes in reservoir rock properties [25]. RPD is useful for understanding the relationships between seismic velocity and porosity and quantifying the cement content from well log data, which is an important aspect of the mechanical strengths of reservoir rocks and fluid flow [19,26]. Moreover, it describes the texture of the rocks, particularly grain size sorting, the position of diagenetic cement, the effect of clay, etc. All these parameters control the response of the elastic properties of pure minerals and their mixtures [14,22].

A number of rock physics models were developed, and relevant studies were conducted to establish and understand the relationship between various rock attributes using laboratory, wireline, and digital rock physics data [24]. According to published literature, several RPD studies have been conducted to predict rock properties quantitatively using the RPM, which matched the well data under examination [27,28,29,30,31,32]. The authors [33] performed RPD to estimate the amount of contact and noncontact cement in sandstone formations from Saudi Arabia. This study combined stiff-sand, soft-sand, and intermediate stiff-sand models to quantify the cement content present in sedimentary rock. Thin section analysis was used to determine the amount of quartz cement. It was found that the thin section was well correlated with the results from RPD. In another study [34], the author presented a RPD study of an offshore gas field to examine wireline data from three wells. They established an RPM as a constant-cement model. A correction of measured shear-wave velocity was also conducted as it appeared unreasonably low. Moreover, they performed perturbational forward modeling to predict the seismic response of various geological scenarios away from well control. RPD was also used in a study for the prediction of the missing S-wave velocity curve [35]. The authors in [36] performed RPD to evaluate the amount of contact cement. Utilizing RPD, we can quantitatively explain the well data.

Numerous retrospective research has been conducted to address the methods for determining and estimating these geological parameters in the Lower Goru sand reservoir. Therefore, it is important to develop the finest RPM and template to identify a previously unidentified, insufficiently unconsolidated reservoir. We illustrate the deterministic methodology to establish RPMs for the well using well data intervals from a clastic depositional environment: tight gas sandstone (TGS) in the Lower Goru fields. Later, the derived reservoir properties can be extrapolated to the prospects away from well control to investigate specific what-if scenarios, such as plausible lithology or fluid variations, thereby building catalogues of RPMs and templates of rock properties. We propose that our rock physics model will provide further insights for accurate reservoir characterization of the Middle and Lower Goru clastic reservoir, Pakistan, and other sites around the world that have the same geology.

2. Geological Setting and Stratigraphy

The study area lies within the N-E pattern dominating the Middle Indus Basin, also known as the Lower Indus Basin. The Middle Indus Basin is situated in a manner that the Indian Shield is abutted to the east; the Sulaiman Fold, Thrust Belt, and Kirther Ranges in the west; Jacobabad-Khairpur High in the south; with Sargodha High in the north [37]. The research area lies close to Jacobabad-Khairpur High (Figure 1). Khairpur High is located in the northwest direction of the study area, which is distinguished by the highest geothermal gradient considered to be effective up to 40.68 °F/328 ft. Khairpur High is dominant in forming structural traps in Miano, Kadanwari, and neighboring areas [13,38,39]. Since the Upper, Middle, and Lower Indus Basins are quite abundant in hydrocarbons and cover rare complete petroleum systems. Many wells in the Lower Indus basin have been drilled [37,40,41]. The strongest rock source is known for being organically abundant, black-colored shales deposited in a quite oxygen-depleted marine environment [42].

In the field of study, the Sembar Formation Lower Cretaceous shale is recognized for the foundation of oil and gas contained in the region due to its thermal maturity and organic richness [43,44]. The Sembar formation is overlaid by the Goru Formation, divided into two distinct sections, the Lower Goru member and the Upper Goru member [44]. The Goru Formation was deposited in the Cretaceous age [37]. The Goru Formation comprises shale, limestone, sandstone, and siltstone [45]. The layered limestone is thin as well as fine-grained within the reservoir rocks. Goru Formation has been deposited on the shelf of the shallow marine environment [12,46]. Limestone and shale have occupied the upper section. Basal and massive sands dominate the lower section of the Goru Formation, and these are the main reservoir rocks in the Middle Indus basin (Dolan, 1990). The lower section of the Lower Goru formation was additionally divided (from bottom to top) into four separate sand intervals: A, B, C, and D (Figure 2) [47,48]. The B, C, and D sand intervals behave as a potential hydrocarbon reservoir in the study area. Two wells, namely, Miano-07 and Rehmat-02, from two different gas fields were examined for rock physics analysis.

3. Materials and Methodology

The study was carried out using well data from Miano-07 and Rehmat-02. In the research, the conventional geophysical logs such as caliper (CAL), gamma-ray (GR), spontaneous potential (SP), neutron porosity (NPHI), density (RHOB), shallow resistivity (LLS), deep resistivity (LLD), sonic (DT), and shear sonic (DTS) were available. The workflow of the technique first applies the conditioning of the well logging data [49,50,51]. Then petrophysical analysis was conducted across the zone of interest to enhance the log data quality and estimate attributes of the reservoir such as porosity (∅), the volume of shale (V_sh), and water saturation (S_w) [52,53]. In the second stage, the RHOB log was corrected. The elastic parameters are computed in a subsequent step. To fine-tune the Gassmann fluid substitution model, we will need to identify its petrophysical characteristics in Step 3. Log tracking was then used to compare the degrees of agreement between the predicted and measured elastic properties. Many diagnostic models of rock physics are contrasted and assessed to characterize the velocity–porosity behavior of a clastic reservoir.

The following steps were followed to validate the rock physics model.

3.1. Density Log Conditioning and Curve Synthesis

After uploading and inspecting the well log data, it was discovered that most log responses are aligned in-depth and did not require a shift in depth. However, when testing the raw log data, it was found that they have CAL log readings, which can affect the log quality. In particular, pad-type devices (e.g., micro-resistivity and density tools) can suffer from inadequate pad contact and may give a substandard quality result. Figure 3a demonstrates the measured RHOB curve response that needs log conditioning at various depths.

The RHOB log plays the most important role in providing a good source of porosity and clay volume data, and the rock physics study is reliant on the RHOB log. Various methods are available for creating a synthetic RHOB; Gardner’s empirical relationship is one of the most famous methods [54]. The main aim of curve synthesis and reconstruction is to shape the logging data at break-out/wash-out intervals on an absent section or where the data is incorrect. Similarly, the RHOB curve was designed at the bad hole interval by curve synthesis, rebuilding the synthetic curve, and applying data from the nearest interval with fine borehole conditions. Gardner et al. [54] investigated primary sedimentary rocks of different lithologies through a wide range of basins, geological ages, and depths and determined relationships between the P-wave velocity and density in rocks of different lithologies. The relationship between the P-wave velocity and density for shale has long been used in seismic analyses:

$${\mathsf{\rho }}_{\mathrm{syn}}=\mathsf{\alpha }\times {\left({\mathrm{v}}_{\mathrm{p}}\right)}^{\mathsf{\beta }}$$

(1)

where V_p is the compressional wave velocity in m/s, ρ_syn is the synthetic RHOB in g/cc, and α = 0.31 when V_p is in m/s or α = 0.23 when V_p is in ft/s and β = 0.25. Major sedimentary rocks generally define a narrow stripe along this trend line, except for coals and evaporate.

In this study, we used the default Gardner parameters to transform velocity data to density. To get a visual feel for the differences between the default Gardner method and the measured RHOB log, Figure 3b is plotted as a typical comparison for the methods derived from velocity logs.

In Figure 3b, the curve in blue color is the evaluated RHOB data, where a few parts of inadequate quality data have been observed against the breakout/wash-out portion. Similarly, for this interval, the synthetic RHOB curves were generated using Gardner’s equation and V_p data in the interval 3525–3580 m. Subsequently, Gardner’s equation is mentioned below as a mathematical expression (Gardner’s Equation (1)).

Figure 3b shows that the log plot verifies the reliability of the synthetic RHOB curve. This also indicates a good match with the measured RHOB data; thus, where the whole condition is fine, it can replace the RHOB data with the portion of the bad hole where the measured data is low. Thus, the calculated RHOB was replaced by the synthesized RHOB curve at 3525–3580 m.

Finally, the synthetic RHOB curve was generated where the measured RHOB curve was not acceptable. Before starting the petrophysical interpretation on Schlumberger Techlog software, all these changes were made. Determining fluids and reservoir parameters were carried out with emphasis on the reservoir, and the closest non-reservoir zones were plotted. V_p, V_s, and ρ were obtained from these zones, which were used later in the rock physics analysis.

3.2. Well Log Response for Identifying Hydrocarbon Reservoir

An alternative procedure for identifying hydrocarbon-bearing zones is mostly based on conventional well log data obtained by the exploratory well-run sondes [55]. The Miano-07 has penetrated Lower Goru sand at a depth of 3420 m, followed by Upper Goru shale at a depth interval of 3320–3410 m (Figure 4). Another well, Rehmat-02, has intruded the Lower Goru at the depth interval of 3420–3650 m. Figure 5 shows the well log response of Rehmat-02 at the targeted sand interval within the Lower Goru Formation. Wireline log reflexes, for example, GR, LLS, LLD, NPHI, and RHOB logs, were used to estimate petrophysical parameters for selected depth intervals in the gas-bearing zones. Figure 4 and Figure 5 show the log responses of P-wave velocity (V_p), S-wave velocity (V_s), P-impedance, and V_p/V_s from each of these wells in the hydrocarbon-bearing regions. The Rehmat-02 well revealed a clean sand reservoir with low RHOB and GR values and a low V_p/V_s ratio. Miano-07 well showed a shaly sand reservoir with alternate high-low GR, low V_p/V_s ratio, and low RHOB responses. Shale and silt lamina produces a relatively high GR response, implying a shaly sand reservoir. The sand reservoir in Miano-07 and Rehmat-02 is part of the Lower Goru Formation of the Early Cretaceous age, deposited in a shelf in the shallow marine environment [12]. Water-saturated reservoirs ranging from 15 to 23% are cemented reservoirs with an effective porosity (∅_e) of 15 to 25% and V_sh ranging from 2.5% in the clean sand reservoir to 17% in the shallow reservoir.

Table 1. Prediction of petrophysical parameters from chosen depth intervals of two wells in the Middle Indus Basin.

Wells	Gas Zones		Petrophysical Parameters
	Depth Interval (m)	Thickness (m)	Vsh (%)	${\mathsf{\varnothing }}_{{e}}$	${{S}}_{{w}}$
Miano-07	3325	1.5	0.45	10	0.55
Rehmat-02	3540	5	0.4	15	0.45

shows the petrophysical parameters of Rehmat-02 and Miano-07 at the targeted depths.

3.3. Rock Mechanical Properties from Well Logs

The significant goal of petrophysical analyses was to obtain proper physical parameters to differentiate fluid and lithology from borehole (log) data. The three fundamental independent elastic parameters that might result from well logs data are V_p, V_s, and ρ. These three fundamental parameters can be used to extract several helpful elastic parameters and have different abilities to discriminate between lithologies and fluids. The elastic properties of rocks can be classified into compressional-type and shear-type parameters [56]. Shear-type parameters are not sensitive to fluid because shear stress cannot transmit in fluids. Examples of shear-type parameters include V_s, SI, and G. Parameters such as V_p, acoustic impedance (AI), Bulk Modulus (K), and Lame’s first parameter (λ), which is sensitive to the fluid and are usually known as compressional-type parameters. Some of the compressional-type and shear-type parameters and their hybrids will be used for the well analysis to describe fluid saturation zones and indicate fluid and lithological changes.

The elastic parameters are estimated by using equations are:

$$\begin{array}{ccc}\mathrm{Poisson}\mathsf{’}\mathrm{s\; ratio} \left(\mathrm{PR}\right)& & \mathsf{\sigma }=\frac{3\mathrm{K}-2\mathrm{G} }{2\left(\mathrm{K}+\mathrm{G}\right)}\end{array}$$

(2)

$$\begin{array}{ccc}\mathrm{Bulk\; modulus} \left(\mathrm{K}\right)& & \mathrm{K}=\mathsf{\rho } \left({\mathrm{v}}_{\mathrm{p}}^2-\frac{4}{3}{\mathrm{v}}_{\mathrm{s}}^2\right)\end{array}$$

(3)

$$\begin{array}{ccc}\mathrm{Shear\; modulus} \left(\mathrm{G}\right)& & \mathrm{G}= \mathsf{\rho }\times { \mathrm{v}}_{\mathrm{S}}^2\end{array}$$

(4)

$$\begin{array}{ccc}\mathrm{Young’s\; modulus} \left(\mathrm{E}\right)& & \mathrm{E}=\frac{9\mathrm{K}\times \mathrm{G}}{3\mathrm{K}+\mathrm{G}}\end{array}$$

(5)

$$\begin{array}{ccc}\mathrm{Compressional\; modulus} \left(\mathsf{\mu }\right)& & \mathsf{\mu }=\mathsf{\rho }\times { \mathrm{v}}_{\mathrm{P}}^2\end{array}$$

(6)

$$\begin{array}{ccc}\mathrm{Lame}\mathsf{’}\mathrm{s\; parameter} \left(\mathsf{\lambda }\right)& & \mathsf{\lambda }=\mathrm{K}-\frac{2\mathsf{\mu }}{3}\end{array}$$

(7)

In the Lower Goru sand, the values of G, E, K, and UCS (unconfined compressive strength) in the Upper Goru shale are increased, followed by an increase in sediments. The value of the Poisson ratio (PR) does differ from 0.29 to 0.31 in the Upper Goru shale, with a decrease to 0.13 in the Lower Goru sand formation. In the Lower Goru sand formation, the gas-producing zone is determined by the low PR values between depths of 3533 m and 3597 m in Rehmat-02 (Figure 6) and similarly, depths of 3333–3339 m in the Miano-07 well (Figure 7).

3.4. Rock Physics Model

We focused on building reliable RPD modeling and RPT analyses for the Lower Goru sand reservoir in the Middle Indus Basin. The model’s friable sand (FS), contact cement (CC), and constant cement (CC) for specific areas of the Miano-07 and Rehmat-02 wells have been implemented to associate seismic and elastic rock properties with geological properties and to estimate the degree to which these geological directions may impact the fluid effect of the seismic parameters and the lithology. We have also employed the Hertz–Mindlin theory (HMT) to operate the modified Hashin–Shtrikman (HS) limits. An effort was made to drive models for clean sands, shaly sands, and shales since those are the lithologies found under the study zones in the Miano-07 and Rehmat-02 wells.

3.5. Rock Physics Diagnostics (RPD)

A standard procedure in the rock physical model is establishing specific trends in porosity–seismic velocity to be determined from the entire data volume and assigning these individual tendencies to correct depth intervals and deposition series. This process is usually referred to as an RPD. The RPD is normally performed on the core and wireline log data [26]. This procedure explains the rock texture: the position of diagenetic cemented and uncemented, grain size sorting, and the influence of shale or clay, etc. The V_p and V_s are usually not well-thought out to be the correct fluid indicator device because of the coupling effect among V_p and V_s via K and G.

In comparison, rock physics trends show up separately in the modulus–porosity or impedance–porosity plane than in the velocity–porosity plane. K is extremely sensitive to pore filling (S_w) deformation due to seismic wave changes resulting from various pore volumes. Any fluid typically does not impact the G. The RPD functional place is the plane of rock physics, which can be (a) modulus–porosity; (b) impedance–porosity; and (c) velocity–porosity plane. Compared to rock diagnostic velocity–porosity planes, we favored elastic moduli porosity or impedance–porosity planes primarily to change the correlated fluid effects in K. In the current research, V_p, K, and G act as a function of porosity for the Miano-07 and Rehmat-02 wells and are overlaid for the three different combined models CC, CC, and FS. These RPD models, with previous knowledge of the quantification of cement, are beneficial to predict the constant cement model with a percentage of cement that will be utilized for searching coordination numbers in HMT to build the dry modulus of rock under the crucial porosity developed by the rock physics template (RPT). Since we do not have the amount of quartz cement in the study in the form of thin sections, the model reported in [57] was applied to quantify cement for RPT generation in the basin.

The cement quantification model is designed to fix log data into a rock modulus–porosity and velocity–porosity system. To define the FS model, HMT and the lower HS bound were implemented. The upper bound HS was applied for cement sand. As a general rule, lithology regulates the elastic properties of rock (texture and composition), porosity (type and quantity), pore fluid, frequency, anisotropy, and depth that include different parameters such as temperature, differential pressure, age, and lithology. In a complex local geological setting, these control parameters are usually different. HMT is utilized to build the dry modulus of high porous rocks, commonly called critical porosity (∅_e). The vital porosity of the sand is endorsed (~40%). The authors in [58] reported that the dry elastic moduli at very high porosity could be extracted using the following relationships:

$${\mathrm{K}}_{\mathrm{HM}}={\left[\frac{{\mathrm{n}}^2\left(1-{\varnothing }_{\mathrm{c}}\right){\mathrm{G}}^2}{18{\mathsf{\pi }}^2{\left(1-\mathrm{PR}\right)}^2}\mathrm{P}\right]}^{\frac{1}{3}}$$

(8)

$${\mathrm{M}}_{\mathrm{HM}}=\frac{\left(5-4\mathrm{PR}\right)}{5\left(1-2\mathrm{PR}\right)}{\left[\frac{3{\mathrm{n}}^2\left(1-{\varnothing }_{\mathrm{c}}\right){\mathrm{G}}^2}{2{\mathsf{\pi }}^2{\left(1-\mathrm{PR}\right)}^2}\mathrm{P}\right]}^{\frac{1}{3}}$$

(9)

where ${\mathrm{M}}_{\mathrm{HM}}$, ${\mathrm{K}}_{\mathrm{HM}}$ = dry rock shear and bulk moduli, respectively; n = coordination number; P = net effective pressure, at ${\varnothing }_c$; PR = Poisson’s ratio for solid phase; and M = shear modulus (mineral modulus) for solid phase. The coordination number (average number of contacts per grain) increases with reduced porosity, resulting from more well-organized packing under increased efficient pressure [59]. High porosity shales have quite a limited number of coordinates and inverses. Cement coordination numbers are estimated to be 8.6 and 4.6 (at crucial porosity) for the reservoir sand and shale, respectively.

The model for rock physics concerns the elastic properties of rock physics. Rocks consist of two fundamental components: pore fluid and mineral frame. The mineral frame consists of several minerals. The conventional way of dealing with this condition is to systematically construct a solitary or efficient mineral whose elastic properties rely on the elastic properties of the mineral components that might be used to measure the dry mineral frame’s elastic properties. The “effective solid” hypothesis is that the effective mineral is a complex consisting of several elastic components of pure minerals with known volumetric fractions. The elastic composite’s K and G have upper and lower boundaries [59]. Established on the HS bounds between zero porosity and crucial porosity, the elastic dry rock modulus can be calculated. For the measurement of dry elastic rock modules with HS bounds for the Miano-07 and Rehmat-02 wells in the Middle Indus Basin, the RHOB and elastic modulus of the various dry solid matrix, for example, feldspar, quartz, and clay are set out in

Table 2. The expected elastic moduli and density of minerals for RPT [59].

Minerals	Bulk Modulus (GPa)	Shear Modulus (GPa)	Density (g/cc)
Quartz	36.6	45.0	2.65
Feldspar	75.6	25.6	2.60
Clay	20.9	6.9	2.58

.

The HS bounds deliver the thinnest variety of elastic moduli without indicating the geometries of components [59,60]:

$${\mathrm{K}}_{\mathrm{HS}\pm }={\mathrm{K}}_1+\frac{{\mathrm{f}}_2}{{\left({\mathrm{K}}_2-{ \mathrm{K}}_1\right)}^{-1}+{\mathrm{f}}_1{\left({\mathrm{K}}_1+\frac{4}{3}{\mathrm{M}}_1\right)}^{-1}}$$

(10)

$${\mathrm{M}}_{\mathrm{HS}\pm }={\mathrm{M}}_1+\frac{{\mathrm{f}}_2}{{\left({\mathrm{M}}_2-{ \mathrm{M}}_1\right)}^{-1}+\frac{2{\mathrm{f}}_1{\left({\mathrm{K}}_1+2{\mathrm{M}}_1\right)}^{}}{5{\mathrm{M}}_1{\left({\mathrm{K}}_1+\frac{4}{3}{\mathrm{M}}_1\right)}^{}}}$$

(11)

where ${\mathrm{M}}_{\mathrm{HS}}$, ${\mathrm{K}}_{\mathrm{HS}}$ = shear and bulk modules, respectively, are calculated via HS bounds; M, K = shear and bulk mineral moduli, respectively, of various components (index mentions to separate phase 1 or 2); and f = fraction of volume of separate phases.

Computing the upper and lower bounds is usually based on material stiffness and softness. The upper bound means the material is rigid, and parameters in (Equations (10) and (11)) will be subscribed with 1. In the current work, the goal of modeling can be achieved using the modified HS bound in which if the material is soft, the bound might be lower, and the parameters can be subscribed with 1. This attaches to end members in the porosity–moduli plane (crucial porosity and zero porosity), i.e., the mineral point and fundamental point are linked by unconsolidated lines, which are modified by lower HS boundaries [58]. CC and CC models with porosity due to cementation take the temperament of V_p and elastic rock modules (K and G). A thick sand interval (Figure 5) was determined in the Rehmat-02 well by low and different GR values along with high V_p (approximately 3 to 3.5 km/s). This sand body is surrounded by shale, contrasting with the sandstone reservoir’s GR and Vp. For Miano-07, we noted the GR log response varying from shale to clean sand in the reservoir zone (Figure 4). In both wells, Miano-07 and Rehmat-02, the shaly sand and the clean sand reservoir zone reflect the same stratigraphic structures, although they are located at various depths and in diverse gas fields. Lower Goru sand in Miano-07 and Rehmat-02 wells consists of shale, cemented sandstones, and friable sands. The hardness of rock depends not only on mineralogy and porosity but also on the rock’s texture (microstructure), i.e., the arrangements on the pore scale of solid phase components. The V_p is graphed in contradiction of porosity for sands that appreciate the microstructure effects. The cement model presupposes that the porosity reduces the unchanging deposition of cement layers at the grain surface in proportion to the initial critical porosity value. Diagenetic cement significantly increases the rigidity of the sand by improving the grain contacts [21]. For diagnostic purposes, the elastic properties are cross-plotted from similar depth intervals of two wells (3533 m and 3597 m for Rehmat-02 and 3335–3339 m for Miano-07) in moduli–porosity and velocity–porosity planes overlaying hypothetical models of rock physics (Figure 8a, Figure 9a and Figure 10a). The cross-plots show that increasing depth of clay content is related to increasing V_p and elastic values. The sorted well data of Miano-07 through the Lower Goru sand was divided into two sections. One is the sand that contains more quartz cement at the contact of grains, and the second is the sand gained from Lower Goru sand with 2–6% cement. A great portion of the scatter is detected in elastic moduli–porosity and velocity–porosity cross-plots for Rehmat-02. This dispersed impact is due to the presence of clay in the sand of the Lower Goru.

The data from the Lower Goru sand of Rehmat-02 was sorted and separated into three parts: friable sand (uncemented sand), sands with large volumes of quartz cement at grain contacts, and sands with 2–6% cement. The upper curve refers to pure quartz cases, whereas the lower curve applies to uncemented rock (Figure 8b, Figure 9b and Figure 10b). We observed that the maximum data value falls on the upper line of the velocity–porosity plane.

The plane reflected all cement deposited exactly on grain contacts with a flat slope for Miano-07 and steep slope for Rehmat-2. The data value drops along constant cement lines, which resembles cement’s uniform deposition on the grain’s surface. For both the Miano-07 and Rehmat-02 wells, the data value dropped at grain contacts and had a steep slope on the upper link of the bulk modulus–porosity plane. The elastic moduli–porosity plane values decreased among the constant cement line and showed a gentle slope towards the steep slope. Microcrack variability is a possible key explanation for the disparity or shift in slope between the porosity of the velocity plane and the porosity of the modulus. Recall that the highest data value is on a constant cement line (6%) that reflects uniformly deposited cement on the grain surface. As observed in Figure 8, Figure 9 and Figure 10, a CC model shows the reservoir zones recognized in the depth intervals: 3533 m and 3597 m for Rehmat-02 and 3335–3339 m for Miano-07. In these RPD cross-plots, the value of the distributions indicates complex geological trends showing lower to higher porosity values. Gas-saturated pay zone grew distinctive and specific elastic moduli–porosity and velocity–porosity patterns for the Miano-07 and Rehmat-02 wells. In the Lower Goru sand reservoir in the Miano-07 well, strong cementation in excess of 6% was noted. In the Rehmat-02, well, the uncemented sand and cemented sand with diverse cementation from 2 to over 6% were noted. Together with the three models of RPD, the CC-CC models were well-thought out to examine the reservoir rock texture and dry rock modulus of these two wells.

3.6. Gassmann Fluid Substitution

Gassmann fluid substitution is usually used for elastic estimation modulus of fluid-saturated rocks such as hydrocarbon or brine. In a RPD and elastic modulus, Gassmann’s equation was used to measure the effect of fluid substitution on the desired saturation. Utilizing the fluid substitution model [61], the elastic modulus of saturated porous rocks was computed, which provided the association between Ks of saturated porous rock, dry rock modulus, pore fluid, and solid matrix.

$${\mathrm{K}}_{\mathrm{Sat}}={\mathrm{K}}_{\mathrm{Dry}}+{\left[\frac{\left(1-\frac{{\mathrm{K}}_{\mathrm{Dry}}}{{\mathrm{K}}_{\mathrm{mineral}}}\right)}{\frac{\varnothing }{{\mathrm{K}}_{\mathrm{fluid}}}+\frac{1-\varnothing }{{\mathrm{K}}_{\mathrm{mineral}}}+\frac{{\mathrm{K}}_{\mathrm{Dry}}}{{\mathrm{K}}_{\mathrm{mineral}}}}\right]}^{}$$

(12)

$${\mathsf{\mu }}_{\mathrm{Sat}}={\mathsf{\mu }}_{\mathrm{Dry}}$$

(13)

where, ${\mathrm{K}}_{\mathrm{Dry}}$ and ${\mathsf{\mu }}_{\mathrm{Dry}}$ are derived from the rock physics model, the K and G of the dry modulus (because fluid G is nil, the G of the fluid-saturated rock is equal to ${\mathsf{\mu }}_{\mathrm{Dry}}$); and $\varnothing$ is the porosity. The saturated water rocks density is

$$\mathsf{\rho }={ \mathsf{\rho }}_{\mathrm{b}}\left(1-\varnothing \right)+{\mathsf{\rho }}_{\mathrm{fluid}}\varnothing$$

(14)

The solid-phase modulus can be computed as the “effective” fluid phase K. The elements can incorporate water, oil, and gas. The immiscible system of the effective K (${\mathrm{K}}_{\mathrm{fluid}-\mathrm{av}}$) is expressed by the authors in [62] if all the individual liquid phases are perfectly connected hydraulically, i.e., the gas pressure is similar to that in the oil and in the water.

$$\frac{1}{{\mathrm{K}}_{\mathrm{fluid}-\mathrm{av}}}=\frac{{\mathrm{S}}_{\mathrm{br}}}{{\mathrm{K}}_{\mathrm{br}}}+\frac{{\mathrm{S}}_{\mathrm{oil}}}{{\mathrm{K}}_{\mathrm{oil}}}+\frac{{\mathrm{S}}_{\mathrm{gas}}}{{\mathrm{K}}_{\mathrm{gas}}}$$

(15)

$${\mathsf{\rho }}_{\mathrm{fluid}}={\mathrm{S}}_{\mathrm{w}}{\mathsf{\rho }}_{\mathrm{w}}+\left(1-{\mathrm{S}}_{\mathrm{w}}\right){\mathsf{\rho }}_{\mathrm{hc}}$$

(16)

where ${\mathrm{K}}_{\mathrm{fluid}-\mathrm{av}}$= effective fluid bulk modulus, ${\mathrm{K}}_{\mathrm{br}}$ = brine bulk modulus, ${\mathrm{S}}_{\mathrm{br}}$= brine saturation, ${\mathrm{K}}_{\mathrm{oil}}$= oil bulk modulus, ${\mathrm{S}}_{\mathrm{oil}}$ = oil saturation, ${\mathrm{K}}_{\mathrm{gas}}$ = gas bulk modulus, ${\mathrm{S}}_{\mathrm{gas}}$ = gas saturation, ${\mathrm{S}}_{\mathrm{w}}$= water saturation, ${\mathsf{\rho }}_{\mathrm{w}}$ = water density, and ${\mathsf{\rho }}_{\mathrm{hc}}$ = hydrocarbon density.

The authors in [63] reported the fluid modulus (${\mathrm{K}}_{\mathrm{fluid}}$)

$${\mathrm{K}}_{\mathrm{fluid}}=\left(\frac{{\mathrm{S}}_{\mathrm{w}}}{{\mathrm{K}}_{\mathrm{w}}}+\frac{\left(1-{\mathrm{S}}_{\mathrm{w}}\right)}{{\mathrm{K}}_{\mathrm{hc}}}\right)$$

(17)

where ${\mathrm{K}}_{\mathrm{hc}}$ = hydrocarbon modulus and ${\mathrm{K}}_{\mathrm{w}}=$ water bulk modulus.

${\mathrm{K}}_{\mathrm{Sat}}$ is needed to inquire about fluid substitution, and the V_p, and V_s are then lastly calculated by using these equations:

$${\mathrm{V}}_{\mathrm{p}}=\sqrt{\left(\frac{{\mathrm{K}}_{\mathrm{w}}+\mathsf{\mu }}{{\mathsf{\rho }}_{\mathrm{b}}}\right)}$$

(18)

$${\mathrm{V}}_{\mathrm{s}}=\sqrt{\left(\frac{\mathsf{\mu }}{{\mathsf{\rho }}_{\mathrm{b}}}\right)}$$

(19)

4. Results and Discussion

In the initial phase, we estimated the RPM and RPT parameters, which were utilized to discriminate lithologies and fluid behavior using multiple cross-plots. In the last phase of the result and discussion, we utilized the results of RMP and RPT with mechanical parameters to discriminate the rock behavior by log plot and cross-plot analyses.

4.1. Rock Physics Template (RPT) Analysis

The RPT technology has grown into a lithology and fluid prediction instrument. A cross-plot between V_p/V_s ratio vs. AI, V_p/V_s ratio vs. PR, and SI vs. AI are the most popular types of RPT. RPT analysis of well log data can be an effective, independent exercise for analyzing and monitoring the quality of well data and assessing seismic detectability of various fluid and lithology scenarios [57]. RPT models must consider regional geological factors. In RPMs, the geological restrictions contain lithology, mineralogy, burial depth, diagenesis, pressure, and temperature. Table 2 shows reservoir parameters and computational fluid properties of the Miano-07 well under the specified pressure and temperature. Gassmann fluid substitution modeling employs the density and K of the various fluids, as described in

Table 3. Reservoir parameters for RPT model.

Parameters	Symbols	Numerical Values	Unit
Bulk modulus of sand	$k_s$	$3.7\times {10}^{10}$	Pa
Shear modulus of sand	${\mu }_s$	$4.4\times {10}^{10}$	Pa
Density of sand	${\rho }_s$	2.65	g/cc
Bulk modulus of clay	$k_c$	$2.15\times {10}^{10}$	Pa
Shear modulus of clay	${\mu }_c$	$1.18\times {10}^{10}$	Pa
Density of clay	${\rho }_c$	2.59	g/cc
Reservoir temperature	T	174.85	°C
Reservoir pressure	P	37.14	MPa
Salinity of brine	$S_b$	20,000	Ppm
Specific gravity of the gas	$G_{gas}$	0.6	API
Bulk modulus of gas	$k_{gas}$	$7.78\times {10}^7$	Pa
Density of gas	${\rho }_{gas}$	1.7123	g/cc
Bulk modulus of brine	$k_b$	$3.18\times {10}^9$	Pa
Density of brine	${\rho }_b$	1	g/cc
An aspect ratio of sand	${\alpha }_s$	0.12	Unitless
An aspect ratio of shale	${\alpha }_c$	0.035	Unitless

.

The outcome of RPD for sands and shales, with log data from 3335–3339 m from Miano-07 and 3533–3597 m from Rehmat-02, were overlaid in the V_p/V_s vs. AI, V_p/V_s vs. PR, and SI vs. AI, in shading by GR values (Figure 11, Figure 12 and Figure 13). The shale and sand plainly showed up in the line of two trends. Each line of porosity sand had its saturation line starting with 100% brine-saturated sands and ending with gas-saturated. Figure 11a,b provides the context of the upper shale-trend line demonstrating pure shale, although the line below displays clean compacted quartz sand filled with brine. In the V_p/V_s vs. AI cross-plot field, there is a line reflecting that gas saturations are nearly perpendicular to the sand trend line on a RPT. This template classifies seal rock (i.e., reservoirs of gas sand), brine-saturated sand, and cap shale for the Miano-07 (Figure 11a) and Rehmat-02 (Figure 11b) wells.

The RPT is shown in Figure 12 as the cross-plots of AI in competition with S-impedance (SI) measured from the well log data of the Miano-07 (Figure 12a) and Rehmat-02 (Figure 12b) wells. A clear distinction between gas reservoir, brine sand, and shale cluster points shows lithology discrimination, with full concentration of the sand clusters at an average porosity of 24% to 26% for both Miano-07 and Rehmat-02 wells. The points with minor AI exhibited high porosity and gas saturation for the gas reservoir plots. The points between sand and shale porosity patterns demonstrate the shaly sand and brine sand. The slight increase in AI at these porosities (24–26%) for the two wells shows that the reservoirs are composed of shaly sand.

Figure 13a,b shows the RPT deployed as the cross-plot of V_p/V_s with PR computed from well log data, focused on available V_p and V_s velocities and the RHOB logging data at the Miano-07 and Rehmat-02 wells. Moreover, the ratio of compressional V_p and V_s velocities and PR could also be calculated. The PR is a diagnostic lithological indicator. It measures the change of a body of material when it is affected by compression or tension [25]. Any change in PR could indicate a change in pore fluid. The gas-saturated zones have a PR (0.23 to 0.4) comparable to the surrounding zones. The shale lithology was plotted at a PR greater than the gas-saturated and brine sand zones.

The identification of gas zones can be established by examining the relationship between V_p/V_s and PR. The gas-saturated (red box), brine sand (green box), and cap shale (blue box) in the research area can be classified using the guideline plot based on extra log information using the V_p/V_s versus PR cross-plot. The area on the cross-plot shows a good fit with our preceding observations of cross-plots about the location of the gas reservoir.

4.2. RPM and RPT-Based Discrimination Characteristics of Reservoir

Both the Miano-07 and Rehmat-02 wells showed a dominant sequence of a sand–shale pattern. These two wells were studied in terms of their lithology and fluid forms. The high GR values represents the lithology of shale. Due to its high conductivity, shale lithology causes the resistivity to shift to the far-left side. Sand lithology with high resistivity and high AI is a highly gas-saturated area. The overall porosity of the rocks (sand–shale lithology) depends on grain size, accumulation, shape, sorting, intergranular matrix, and cement. Rock physics modeling and RPT studies have provided an approach to realize reservoir microstructures, cementation, fluid types, and lithology under specific temperature and pressure conditions. For the chosen intervals of 3325–3335 m from Miano-07 (Figure 14a), the NPHI and RHOB cross-plots showed the existence of gas-bearing Lower Goru sand which is highlighted by a red box. The cross-plot in Figure 14b demonstrates that the presence of laminated shale in the Lower Goru formation affects the Rehmat-02 well at the interval of 3533–3597 m. The log examination, along with the NPHI and RHOB cross-plot (Figure 14), identified clean sand gas-bearing zones (red colored boxes) and shaly sand gas-bearing zones that lie below the SS (sandstone) line.

Rock physics modeling and RPT were used to characterize the active seismic sands away from well regulations. For the Miano-07 (Figure 15) and Rehmat-02 (Figure 16) wells, conventional log curves along with the mechanical and elastic properties agree with each other and reliably indicate the existence of shaly sand and sand reservoirs. Overall, low GR, low V_p/V_s, low PR ratio value, higher resistivity, and a cross-over between NPHI and RHOB logs were observed. As a result, log data displayed on the developed RPT fit well with conventional wireline log data, as seen in the intersection between the NPHI and RHOB logs. The application of the RPT is a more effective tool for characterizing the reservoir in any geological setting. It can be known as reservoir identification, local geological prediction, and risk minimization in undeveloped areas. In addition to applying RPT to the whole Lower Goru tight sand reservoir in Pakistan and elsewhere with similar geological trends and reservoir distribution, the anticipated model may be used to estimate porosity from seismically determined impedance anywhere in the globe.

5. Conclusions

A model was designed for an unconsolidated reservoir in a Lower Goru gas field on the supposition that the rock’s internal structure is identical to the idealized structure anticipated by the theoretical model if well log data points drop close to a theoretical line in the velocity–porosity plane. Several models were compared: friable sand, contact cement, and constant cement with varying degrees of cementation. The gas sands water saturation in the range of 21–60%, density in the range of 2.25–2.35 g/cc, with P-velocity and S-velocity in the range of 2800–2950 m/s and 2251–2372 m/s, respectively, inferred that the reservoir rocks are unconsolidated.

The RPM established in the study would improve the understanding of characteristics for the relatively high V_p/V_s unconsolidated reservoir sands under the study. The integration of RPMs would improve the prediction of reservoir properties from the elastic properties.

Through extensive cross-plot analysis of various rock properties, lithology and fluid discrimination were achieved with shear velocity versus porosity, shear impedance, and V_p/V_s versus acoustic impedance cross-plots, in agreement with the estimated water saturation and oil sand porosity. The shear velocity versus porosity cross-plots estimated oil sand porosities ranges from 10–15%, with a 2−3% degree of diagenetic process of constant cement model, alongside clear lithology distinction between oil sands and shales. These templates could be correlated with the seismic inverted results from the same field or with deposition environment to predict the lithology of undrilled areas and their fluid content.

Future development will include close work with real seismic and well log data to further calibrate and validate the proposed methods of reservoir characterization as well as to chart the areas of their applicability.