3.2.1. Seismic Data Interpretation

Seismic interpretation in this study is focused on a particular area of interest within the two-way time (TWT) range of −1800 to −2600 (ms) (Figure 4). Depending on the seismic reflection discontinuities and terminations, manual horizon picking followed by seeded horizon auto-tracking was adopted to interpret the target horizons on 2D cross-sectional (i.e., vertical) slices of a 3D seismic volume (Figure 4a). In the middle parts of the seismic cross-sections, the reflections are mostly moderately chaotic and difficult to correlate due to the complexity of the geology and faults resulting from tectonic compression (Figure 4b). The modeling step of the seismic interpretation involves 3D TWT contour surfaces construction. Consequently, 3D TWT contour surfaces were constructed by marking the tops of each stratigraphic interface on the extended 3D seismic volume. The smooth function was then applied to the generated surfaces at three iteration levels to produce geologically reasonable stratigraphic surfaces. These 3D TWT contour surfaces were then used to interpret the prevailing structural trends in the study area via conventional methods.

**Figure 3.** Workflow highlighting various steps and methods employed in this study for reservoir characterization and prospect assessment of the Kadanwari field, MIB, Pakistan. **Figure 3.** Workflow highlighting various steps and methods employed in this study for reservoir characterization and prospect assessment of the Kadanwari field, MIB, Pakistan.

Seismic interpretation in this study is focused on a particular area of interest within the two-way time (TWT) range of −1800 to −2600 (ms) (Figure 4). Depending on the seismic reflection discontinuities and terminations, manual horizon picking followed by seeded horizon auto-tracking was adopted to interpret the target horizons on 2D cross-sectional

3.2.1. Seismic Data Interpretation

methods.

**Figure 4.** (**a**) 3D time-seeded horizon auto-tracking; (**b**) amplitude display of vertical time inline cross-sections (e.g., 1970, 2000, 2030, and 2050) in the SW to NE direction, showing horizon interpretation. Mint green, white, and yellow colors represent the lateral extent of the G sand interval, the E sand interval, and the Sembar Formation, respectively. **Figure 4.** (**a**) 3D time-seeded horizon auto-tracking; (**b**) amplitude display of vertical time inline crosssections (e.g., 1970, 2000, 2030, and 2050) in the SW to NE direction, showing horizon interpretation. Mint green, white, and yellow colors represent the lateral extent of the G sand interval, the E sand interval, and the Sembar Formation, respectively.

(i.e., vertical) slices of a 3D seismic volume (Figure 4a). In the middle parts of the seismic cross-sections, the reflections are mostly moderately chaotic and difficult to correlate due to the complexity of the geology and faults resulting from tectonic compression (Figure 4b). The modeling step of the seismic interpretation involves 3D TWT contour surfaces construction. Consequently, 3D TWT contour surfaces were constructed by marking the tops of each stratigraphic interface on the extended 3D seismic volume. The smooth function was then applied to the generated surfaces at three iteration levels to produce geologically reasonable stratigraphic surfaces. These 3D TWT contour surfaces were then used to interpret the prevailing structural trends in the study area via conventional

#### 3.2.2. Three-Dimensional Structural Modeling (3D SM) 3.2.2. Three-Dimensional Structural Modeling (3D SM)

Employing 3D fault system modeling (FSM) in seismic analysis is a crucial step to constrain horizon interpretation [21,50]. One recent significant advance in FSM is the rendering available of 3D seismic data that provide detailed images of large volumes of rock and the often-complex 3D fault networks [51]. Meanwhile, most seismic datasets have signal disturbance zones, particularly in highly faulted areas. Furthermore, the discontinuities in seismic reflectors can be poorly resolved, resulting in the approximate localization or misinterpretation of faults. However, seismic discontinuities can be more clearly defined if the detection is based on multiple attributes and suitable filters [51,52]. Therefore, in this study, before proceeding to the 3D FSM, the precision of the geologic fault boundaries was assessed first for accuracy using dip magnitude and seismic edge enhancement attributes (Figure 5). Accordingly, the dominance of a normal fault system was inferred in the dip magnitude and edge enhancement cross-sections by aligning Employing 3D fault system modeling (FSM) in seismic analysis is a crucial step to constrain horizon interpretation [21,50]. One recent significant advance in FSM is the rendering available of 3D seismic data that provide detailed images of large volumes of rock and the often-complex 3D fault networks [51]. Meanwhile, most seismic datasets have signal disturbance zones, particularly in highly faulted areas. Furthermore, the discontinuities in seismic reflectors can be poorly resolved, resulting in the approximate localization or misinterpretation of faults. However, seismic discontinuities can be more clearly defined if the detection is based on multiple attributes and suitable filters [51,52]. Therefore, in this study, before proceeding to the 3D FSM, the precision of the geologic fault boundaries was assessed first for accuracy using dip magnitude and seismic edge enhancement attributes (Figure 5). Accordingly, the dominance of a normal fault system was inferred in the dip magnitude and edge enhancement cross-sections by aligning reflector discontinuities and degree of displacement. In the modeling phase, the fault surfaces were constructed by employing the fault polygon in PetrelTM software for each type of fault with various geometrical structures. Finally, 3D FSMs and their attribute models (e.g., dip angle, fault rose diagram and histogram) were constructed to evaluate the fault system's geometric distribution and mechanics within the early Cretaceous stratigraphic sequence.

stratigraphic sequence.

**Figure 5.** Fault system tracing and interpretation in the Kadanwari field using resemble brittle deformation features (seismic reflection discontinuities, terminations, and amount of displacement) on (**a**) dip magnitude and (**b**) seismic edge enhancement attributes cross-sections. **Figure 5.** Fault system tracing and interpretation in the Kadanwari field using resemble brittle deformation features (seismic reflection discontinuities, terminations, and amount of displacement) on (**a**) dip magnitude and (**b**) seismic edge enhancement attributes cross-sections.

reflector discontinuities and degree of displacement. In the modeling phase, the fault surfaces were constructed by employing the fault polygon in PetrelTM software for each type of fault with various geometrical structures. Finally, 3D FSMs and their attribute models (e.g., dip angle, fault rose diagram and histogram) were constructed to evaluate the fault system's geometric distribution and mechanics within the early Cretaceous

Three-dimensional SM aims to better image and understand complex reservoirs' geological structures. The local tectonics of the Kadanwari field in MIB have resulted in various structural deformations that produce many uncertainties when assessing the reservoirs' 3D structural framework [7]. Understanding such deformations can be better achieved via 3D SM using relevant algorithm models, e.g., the volume-based modeling (VBM) approach in the PetrelTM software based on inferred seismic data integrated with borehole information [50,53,54]. Traditional modeling methods generally involve the oversimplification of geological settings; however, grid-generated VBM captures a realistic reservoir architecture [2]. This can be easily transferred into the dynamic realm, providing a better understanding of the reservoir for future field management and development activities [50]. In this study, 3D SM based on the VBM approach was performed in three main steps (e.g., fault modeling, pillar gridding, and horizon creation).

3.2.3. Three-Dimensional Seismic Attribute Analysis The seismic attributes were computed via mathematical manipulation of the original seismic data to highlight specific geological, physical, or reservoir properties [50,55]. Variations in the amplitude, phase, frequency, and bandwidth of the seismic waves were subsequently used to validate the spatial forecasts of the geological structure and to evaluate the spatial distribution of the facies and the DHIs in the G sand interval, the E sand interval, and the Sember Formation time windows. The variance edge attribute can visualize seismic amplitude discontinuities related to faulting or stratigraphy. It delineated the prominent faults and seismic amplitude discontinuity in both the horizon slices and the vertical seismic profiles, thereby validating the manual interpretation of faults. The sweetness attribute was applied to both horizon Three-dimensional SM aims to better image and understand complex reservoirs' geological structures. The local tectonics of the Kadanwari field in MIB have resulted in various structural deformations that produce many uncertainties when assessing the reservoirs' 3D structural framework [7]. Understanding such deformations can be better achieved via 3D SM using relevant algorithm models, e.g., the volume-based modeling (VBM) approach in the PetrelTM software based on inferred seismic data integrated with borehole information [50,53,54]. Traditional modeling methods generally involve the oversimplification of geological settings; however, grid-generated VBM captures a realistic reservoir architecture [2]. This can be easily transferred into the dynamic realm, providing a better understanding of the reservoir for future field management and development activities [50]. In this study, 3D SM based on the VBM approach was performed in three main steps (e.g., fault modeling, pillar gridding, and horizon creation).

#### 3.2.3. Three-Dimensional Seismic Attribute Analysis

The seismic attributes were computed via mathematical manipulation of the original seismic data to highlight specific geological, physical, or reservoir properties [50,55]. Variations in the amplitude, phase, frequency, and bandwidth of the seismic waves were subsequently used to validate the spatial forecasts of the geological structure and to evaluate the spatial distribution of the facies and the DHIs in the G sand interval, the E sand interval, and the Sember Formation time windows.

The variance edge attribute can visualize seismic amplitude discontinuities related to faulting or stratigraphy. It delineated the prominent faults and seismic amplitude discontinuity in both the horizon slices and the vertical seismic profiles, thereby validating the manual interpretation of faults. The sweetness attribute was applied to both horizon slices and vertical seismic profiles to identify sweet spot zones that are hydrocarbon-prone. It can be defined as the reflection strength (instantaneous amplitude) divided by the square root of the instantaneous frequency. The high sweetness anomalies highlight a seismic signal's high amplitudes and low-frequency contents and vice versa. Therefore, combining these two physical quantities helps distinguish sand bodies from shale and predicted gasprone zones [56]. The RMS amplitude attribute was also applied to both horizon slices and

vertical seismic profiles to measure amplitude anomalies to identify the spatial distribution of facies and DHIs at the G and E sand reservoir intervals. The RMS amplitude computes the square root of the sum of squared amplitudes divided by the number of samples within the window used [19].

#### 3.2.4. Petrophysical Modeling

Petrophysical modeling is critical in a reservoir study because it represents a primary input data source for integrated reservoir characterization and resource evaluation [54,57]. This study performed petrophysical modeling based on well log data and internal geological reports to evaluate the G and E sand reservoir intervals' properties, e.g., lithology, *Vshale*, ∅*avg*, ∅*eff*, *SW*, and *Shc* (Table 1). Successful evaluation of these properties is necessary when determining the hydrocarbon potential of a reservoir system. The detailed plots of the well log curves and their depth ranges within the G and E sand reservoir intervals are shown in Figure 6. These log curves express the physical motifs of the stacked geological strata as a function of depth, which can help identify lithologies and porosities, and differentiate between porous and non-porous rocks and pay zones. *Minerals* **2022**, *12*, x 11 of 27

**Figure 6.** Input conventional log curves for the petrophysical modeling of the G sand interval in the (**a**) Kadanwari-10 (**b**) and Kadanwari-11 wells and for E sand interval in the (**c**) Kadanwari-10 (**d**) and Kadanwari-11 wells. (1) **Figure 6.** Input conventional log curves for the petrophysical modeling of the G sand interval in the (**a**) Kadanwari-10 (**b**) and Kadanwari-11 wells and for E sand interval in the (**c**) Kadanwari-10 (**d**) and Kadanwari-11 wells.

The following five key steps were followed to evaluate the reservoirs' fundamental

permeability [58]. We used the GR log technique for *Vshale* estimation by firstly estimating the gamma ray index (*IGR*). The IGR was initially adopted using Equation (1) to estimate *Vshale*, utilizing the GRlog in track 1of Figure 6. Secondly, to obtain a realistic *Vshale* estimation without overestimating the content of shale (first-order approximation: V*shale* = *IGR*), a non-linear relationship (Equation (2)) proposed by

> ீோ = − ௫ −

Dolan was employed [59].

properties.

The following five key steps were followed to evaluate the reservoirs' fundamental properties.

(1) Volume of Shale (*Vshale*)—The presence of shale in the productive zone severely impacts the petrophysical properties and can cause a reduction in the ∅*eff*, ∅*avg*, and permeability [58]. We used the GR log technique for *Vshale* estimation by firstly estimating the gamma ray index (*IGR*). The *IGR* was initially adopted using Equation (1) to estimate *Vshale*, utilizing the *GRlog* in track 1 of Figure 6. Secondly, to obtain a realistic *Vshale* estimation without overestimating the content of shale (first-order approximation: *Vshale* = *IGR*), a non-linear relationship (Equation (2)) proposed by Dolan was employed [59].

$$I\_{GR} = \frac{GR\_{log} - GR\_{min}}{GR\_{max} - GR\_{min}} \tag{1}$$

$$V\_{\text{shale}} = 1.7 - [(3.38 - (I\_{GR} + 0.7)^2)]^{\frac{1}{2}} \tag{2}$$

where *IGR* stands for gamma ray index, *GRlog* shows the reading of the gamma ray log, *GRmin* and *GRmax* are the lower and upper limits of the *GRlog* value in shale, respectively.

(2) Average Porosity (∅*avg*)—Total porosity or ∅*avg* represents all the voids or pore spaces of the rock, including interconnected and isolated pores and pore spaces occupied by clay-bound water [2]. In this study, DT, RHOB, and NPHI logs that are sensitive to sedimentary micro-facies were selected to calculate ∅*avg*, by which process the conventional logging responses of the G and E sand reservoir intervals can be summarized (Figure 6). The DT log measures the sound waves' traveling times in the rock unit. The sound waves in the rock unit depend on the shape, matrix material, and cementation (Equation (3)). Accordingly, the Sonic–Raymer (SR) porosity model was used to evaluate sonic porosity (∅*S*) (Equation (4)) [40].

$$t\_{\log} = t\_m - \left[ (1 - \mathcal{Q}\_S) + t\_f \,\mathcal{Q}\_S \right] \tag{3}$$

$$\mathcal{Z}\_S = \frac{\Delta t\_{\log} - \Delta t\_m}{\Delta t\_m - \Delta t\_m} \tag{4}$$

where ∅*S* represents sonic porosity, *tlog* represents the log reading in µsec/m, *tm* represents the matrix interval transient time, ∆*tlog* represents the formation interval transient time in µsec/m, and ∆*tm* represents formation fluids' interval transient time in µsec/m. The density porosity (*φD*) was calculated using the RHOB log via Equation (5) [2].

$$\phi\_D = \left[ (\rho\_{\rm mat} - \rho\_b) / (\rho\_m - \rho\_f) \right] \tag{5}$$

where *ρmat* represents the density of the matrix, *ρf* represents the fluid density, and *ρb* is the bulk density. The NPHI log measures the neutron porosity (∅*N*) by assuming that the pores are filled with fluid. Therefore, it measures the hydrogen concentration and energy loss. The ∅*N* can be expressed via Equation (6).

$$
\phi\_N = [\mathbf{a}N + \mathbf{b}] \tag{6}
$$

where *φN* is the neutron-derived porosity, a and b are constants, and *N* is the neutron count in the formation intervals. In addition, the density–neutron cross-plot in track 4 in Figure 6 determines the cross-over gas effects. After identifying the porosities (e.g., ∅*S*, ∅*D*, and ∅*N*) from the DT, RHOB, and NPHI logs, Equation (7) was used to calculate ∅*avg*.

$$\mathcal{Z}\_{\text{avg}} = \frac{\mathcal{Q}\_S + \mathcal{Q}\_D + \mathcal{Q}\_N}{3} \tag{7}$$

(3) Effective Porosity (∅*eff*)—The ∅*eff* was calculated using Equation (8) [40,43].

$$\phi\_{eff} = \begin{bmatrix} (\phi\_{avg}) \times (1 - V\_{shale}) \end{bmatrix} \tag{8}$$

where ∅*avg* represents the average porosity, and *Vshale* is the shale content in volume units.

(4) Water Saturation (*SW*)—The Poupon–Leveaux Indonesian (PLI) model is one of the best models for estimating S*W* in a shaly sand reservoir [60]. In this study, the constraints of *Vshale* (Equation (2)) and ∅*eff* (Equation (8)), and the resistivity variation in *Vshale* and water formation, were subsequently integrated using the PLI model, as in Equation (9), to determine S*W*.

$$S\_w = \left\{ \left[ \left( \frac{V\_{shale} 2^{-V\_{shale}}}{R\_{shale}} \right)^2 + \left( \frac{\mathcal{Z}\_{eff}}{R\_w} \right)^2 \right]^2 R\_t \right\}^2 \tag{9}$$

where *Rt* is the true resistivity of the formation obtained from the LLD log response, *Rshale* is the resistivity variation in the *Vshale*, and *Rw* is the resistivity of water formation.

(5) Hydrocarbon Saturation (*Shc*)—The *Shc* was calculated by subtracting the percentage of pore volume occupied by S*w* from 1; the remaining percentage pore volume gives the S*hc* (Equation (10)).

$$S\_{\rm llc} = 1 - S\_{\rm wv} \tag{10}$$

#### **4. Results**

#### *4.1. Stratigraphic Interfaces Interpretation*

The interpretation of the stratigraphic interfaces of the early Cretaceous sequence is integrated into an internally consistent 3D workflow, which is easy to incorporate and use to form conclusions that exceed known points by constraining 3D structural models. The interpreted cross-sectional (i.e., vertical) slices of a 3D seismic volume show the displacement and deformation of the stratigraphic interface (Figure 4b). The structural variations in the stratigraphic interfaces are predicated primarily on the TWT contour surfaces because contour lines connect the same elevation; this is why they are essential tools for analyzing and interpreting seismic data. The time surfaces plot the TWT of seismic signals from the surface to the horizon and reflect the interpreted stratigraphic interface's distribution. In Figure 7, the TWT surfaces indicate the extension and propagation of the regional stratigraphic structure in the subsurface. The TWT values of the interpreted stratigraphic interfaces decrease towards the central part of the field, which gives rise to structural highs in the southwestern and southeastern portions. The southwestern and southeastern portions of the interpreted stratigraphic interfaces are structurally high; hence, this is a region of interest for hydrocarbon exploration. The minimum, mean and maximum TWT variations of the G sand interval, the E sand interval, and the Sembar Formation in the southwestern and southeastern transect are presented in Table 2.

**Table 2.** Minimum, mean, and maximum TWT variations of the stratigraphic interfaces in the study area.

**Figure 7.** Three-dimensional TWT contour surfaces of the stratigraphic interfaces show different structural distributions: (**a**) TWT top surface of the G sand interval; (**b**) TWT top surface of the E sand interval; (**c**) TWT top surface of the Sembar formation in 3D space. **Figure 7.** Three-dimensional TWT contour surfaces of the stratigraphic interfaces show different structural distributions: (**a**) TWT top surface of the G sand interval; (**b**) TWT top surface of the E sand interval; (**c**) TWT top surface of the Sembar formation in 3D space.

#### *4.2. Three-Dimensional Fault System Models (3D FSMs) 4.2. Three-Dimensional Fault System Models (3D FSMs)*

The study area's complex and composite subsurface morphology resulted from multiple episodes of tectonic activity. Multi-stage tectonic movements contribute to the intersection of faults formed simultaneously or at different stages of tectonic movement (Figure 8a). These faults occurred during various phases of adjustment and deformation, along with the structure inversion and loss of reservoir quality in the Kadanwari field, The study area's complex and composite subsurface morphology resulted from multiple episodes of tectonic activity. Multi-stage tectonic movements contribute to the intersection of faults formed simultaneously or at different stages of tectonic movement (Figure 8a). These faults occurred during various phases of adjustment and deformation, along with the structure inversion and loss of reservoir quality in the Kadanwari field, MIB. In order to effectively represent the spatial distribution of the fault system and its influence on the frag-

mentation of the stratigraphic interfaces, individual 3D fault surfaces have been illustrated, along with the 3D seismic volume (Figure 8a). The orientation and spatial distribution of the fault surfaces reveal that the intact depositional environment of the Kadanwari field was influenced by the tectonic regime of the surrounding plate boundary, which has continuously influenced the stratigraphic structure (Figure 8a). The 3D seismic volume interpretation at the field scale shows that the middle part of the field has undergone greater deformation than both sides. This is why the number of faults and their complexity decrease from the center to either field side. These faults have been interpreted as normal faults, generally showing NW to SE directions, controlling the distribution of depositional facies, reservoir compartmentalization, and hydrocarbon up-dip migration in the study area. These NW–SE normal faults have governed the complex structural configuration of the stratigraphic sequence (e.g., G sand interval, E sand interval, and the Sembar Formation). The overall pattern of these faults can be regarded as a negative flower structure, which significantly increases the likelihood of the successful positioning of hydrocarbon traps. The degree of completion of these negative flower structures is positively correlated with improved hydrocarbon migration and abundance in the reservoir window (e.g., G and E sand intervals) (Figure 8b). In addition, the presence of a negative flower structure indicates the combined effects of extensional and strike-slip motions. The fault orientation results derived from the 3D dip angle models, the stereonet, the rose diagram, and the histogram show that most of the faults are oriented S30◦–45◦ E and N25◦–35◦ W, with an azimuth of 148◦–170◦ and 318◦–345◦ , and exhibiting minimum, mean and maximum dip angles of 28◦ , 62◦ , and 90◦ respectively (Figure 9a). The histogram in Figure 9b displays the relationship between fault frequency in percentage and dip angle in degrees, which shows that most of the faults have dips ranging between 35◦ and 75◦ . In comparison, 20% of the total fault planes have dips in the range 80◦–90◦ , and the lowest fault dip observed is approximately 28◦ (Figures 8a and 9a). influence on the fragmentation of the stratigraphic interfaces, individual 3D fault surfaces have been illustrated, along with the 3D seismic volume (Figure 8a). The orientation and spatial distribution of the fault surfaces reveal that the intact depositional environment of the Kadanwari field was influenced by the tectonic regime of the surrounding plate boundary, which has continuously influenced the stratigraphic structure (Figure 8a). The 3D seismic volume interpretation at the field scale shows that the middle part of the field has undergone greater deformation than both sides. This is why the number of faults and their complexity decrease from the center to either field side. These faults have been interpreted as normal faults, generally showing NW to SE directions, controlling the distribution of depositional facies, reservoir compartmentalization, and hydrocarbon updip migration in the study area. These NW–SE normal faults have governed the complex structural configuration of the stratigraphic sequence (e.g., G sand interval, E sand interval, and the Sembar Formation). The overall pattern of these faults can be regarded as a negative flower structure, which significantly increases the likelihood of the successful positioning of hydrocarbon traps. The degree of completion of these negative flower structures is positively correlated with improved hydrocarbon migration and abundance in the reservoir window (e.g., G and E sand intervals) (Figure 8b). In addition, the presence of a negative flower structure indicates the combined effects of extensional and strike-slip motions. The fault orientation results derived from the 3D dip angle models, the stereonet, the rose diagram, and the histogram show that most of the faults are oriented S30°–45° E and N25°–35° W, with an azimuth of 148°–170° and 318°–345°, and exhibiting minimum, mean and maximum dip angles of 28°, 62°, and 90° respectively (Figure 9a). The histogram in Figure 9b displays the relationship between fault frequency in percentage and dip angle in degrees, which shows that most of the faults have dips ranging between 35° and 75°. In comparison, 20% of the total fault planes have dips in the range 80°–90°, and the lowest fault dip observed is approximately 28° (Figures 8a and 9a).

MIB. In order to effectively represent the spatial distribution of the fault system and its

*Minerals* **2022**, *12*, x 15 of 27

**Figure 8.** (**a**) The individual and combined distribution of tectonic extensional fault surfaces along the early Cretaceous stratigraphic sequence in the 3D seismic volume (SW-SE) and (**b**) 3D dip angle models (SW-SE) of the individual and combined tectonic extensional fault surfaces.

11b).

the early Cretaceous stratigraphic sequence in the 3D seismic volume (SW-SE) and (**b**) 3D dip angle

models (SW-SE) of the individual and combined tectonic extensional fault surfaces.

**Figure 9.** (**a**) Steroenet showing the distribution of faults' dips and orientations identified, along with the entire seismic survey in the Kadanwari field, MIB, Pakistan; (**b**) a histogram showing the dip angles of the interpreted faults along with their frequency. **Figure 9.** (**a**) Steroenet showing the distribution of faults' dips and orientations identified, along with the entire seismic survey in the Kadanwari field, MIB, Pakistan; (**b**) a histogram showing the dip angles of the interpreted faults along with their frequency.

#### *4.3. Three-Dimensional Structural Models (3D SMs) 4.3. Three-Dimensional Structural Models (3D SMs)*

The 3D SMs and several 2D time-domain structural cross-sections have been derived to illustrate the detailed structural and stratigraphic setting of the sequence in the study area. Figure 10a displays the southwestern and southeastern transects of the 3D TWT model, while Figure 10b represents the 2D TWT cross-sections derived from the 3D TWT model result. Figure 11a displays the southwestern and southeastern transect of the 3D SMs derived from VBM. Moreover, the fault system is incorporated into the 2D crosssections during model computation to constrain the horizons at each fault offset (Figure The 3D SMs and several 2D time-domain structural cross-sections have been derived to illustrate the detailed structural and stratigraphic setting of the sequence in the study area. Figure 10a displays the southwestern and southeastern transects of the 3D TWT model, while Figure 10b represents the 2D TWT cross-sections derived from the 3D TWT model result. Figure 11a displays the southwestern and southeastern transect of the 3D SMs derived from VBM. Moreover, the fault system is incorporated into the 2D cross-sections during model computation to constrain the horizons at each fault offset (Figure 11b). *Minerals* **2022**, *12*, x 17 of 27

**Figure 10.** TWT domain structural model of the early Cretaceous stratigraphic sequence (e.g., G sand interval, E sand interval, and the Sembar Formation): (**a**) SW–SE view; (**b**) extracted inline 2D TWT cross-sections. **Figure 10.** TWT domain structural model of the early Cretaceous stratigraphic sequence (e.g., G sand interval, E sand interval, and the Sembar Formation): (**a**) SW–SE view; (**b**) extracted inline 2D TWT cross-sections.

**Figure 11.** Final result of the 3D SM: (**a**) SW–SE view and (**b**) SW–NE oriented 2D cross-sections

drawn from the 3D SM results.

*4.4. Seismic Attributes Interpretation*  4.4.1. Variance Edge Attribute

TWT cross-sections.

**Figure 11.** Final result of the 3D SM: (**a**) SW–SE view and (**b**) SW–NE oriented 2D cross-sections drawn from the 3D SM results. **Figure 11.** Final result of the 3D SM: (**a**) SW–SE view and (**b**) SW–NE oriented 2D cross-sections drawn from the 3D SM results.

**Figure 10.** TWT domain structural model of the early Cretaceous stratigraphic sequence (e.g., G sand interval, E sand interval, and the Sembar Formation): (**a**) SW–SE view; (**b**) extracted inline 2D

*4.4. Seismic Attributes Interpretation*  4.4.1. Variance Edge Attribute Detailed structural analysis reveals that the geological structural complexity is a consequence of different tectonic phases of deformation (e.g., extensional and strike-slip deformation from the early Cretaceous to Quaternary). The complex structural mechanics of both extensional and strike-slip movements have been observed in the 3D SMs. These complex structural mechanics were controlled by the normal fault system's NW–SE dipping. The above-explained normal fault patterns (e.g., negative flower structure) show brittle deformation features, where the interpreted sequences have been displaced relative to each other (amount of displacement) (Figure 11b). These significantly control the structural domains, consisting mainly of half-graben, horst, half-graben, graben, half-graben, and horst, from SW to NE. The thickness of the sequence increases towards the NE. The thickness decreases in the central portion of the region from the top E sand interval to the top Sembar Formation. The hydrocarbon migration direction can be determined from the spatial distribution and composition of the fault system within the early Cretaceous sequences. The geometrical trend of these fault systems creates a pathway that is crucial for hydrocarbon migration in the vertical direction. This up-dip migration of the hydrocarbons from the Sembar formation towards the G and E sand intervals has resulted in hydrocarbon accumulation, validated through seismic attribute analysis and petrophysical modeling.

#### *4.4. Seismic Attributes Interpretation*

#### 4.4.1. Variance Edge Attribute

The variance edge attribute indicates discontinuities related to faulting or stratigraphy in the vertical seismic cross-sections and is proved to help depict significant fault zones and fractures, thereby validating fault interpretation. Figure 12 shows the results of the variance edge attribute calculated from the 3D seismic volume, appropriate cross-sections (e.g., A–A0 , B–B0 , C–C0 , D–D0 , E–E0 and F–F0 ), and horizon slices. The horizon slices include the G sand interval, the E sand interval, and the Sembar Formation. The variance edge attribute values in the 2D variance cross-sections and horizon slices range from 0.0 to 1.0. A variance value equal to 1 indicates discontinuity (fault), while a value of 0 variances

represents continuous seismic events. The darkest regions (e.g., values ranging from 0.8 to 1) in vertical strips can be interpreted as faults or fractures (Figure 12c). These faults and fractures create an essential pathway for vertical hydrocarbon migration. ranging from 0.8 to 1) in vertical strips can be interpreted as faults or fractures (Figure 12c). These faults and fractures create an essential pathway for vertical hydrocarbon migration.

The variance edge attribute indicates discontinuities related to faulting or stratigraphy in the vertical seismic cross-sections and is proved to help depict significant fault zones and fractures, thereby validating fault interpretation. Figure 12 shows the results of the variance edge attribute calculated from the 3D seismic volume, appropriate cross-sections (e.g., A–A′, B–B′, C–C′, D–D′, E–E′ and F–F′), and horizon slices. The horizon slices include the G sand interval, the E sand interval, and the Sembar Formation. The variance edge attribute values in the 2D variance cross-sections and horizon slices range from 0.0 to 1.0. A variance value equal to 1 indicates discontinuity (fault), while a value of 0 variances represents continuous seismic events. The darkest regions (e.g., values

*Minerals* **2022**, *12*, x 18 of 27

**Figure 12.** (**a**) 3D variance edge attribute volume, (**b**) 2D extracted variance cross-sections, and (**c**) 3D variance horizon slices of the G sand interval, the E sand interval, and the Sembar Formation showing discontinuities. **Figure 12.** (**a**) 3D variance edge attribute volume, (**b**) 2D extracted variance cross-sections, and (**c**) 3D variance horizon slices of the G sand interval, the E sand interval, and the Sembar Formation showing discontinuities.

#### 4.4.2. Sweetness Attribute 4.4.2. Sweetness Attribute

Figure 13 shows the sweetness attribute computed from the 3D seismic volume, corresponding cross-sections (e.g., A–A′, B–B′, C–C′, D–D′, E–E′ and F–F′), and sweetness horizon slices. The high sweetness anomalies (at −1900 to −2300 ms) on the 2D sweetness cross-sections (Figure 13b) and the 3D horizon slices (e.g., G sand and E sand intervals) (Figure 13c) contributed to the high amplitude and low frequency. In contrast, the low sweetness anomalies (at −2300 to −2500 ms) on the 2D sweetness cross-sections and the SW—NE parts of the 3D horizon slices (e.g., the Sembar Formation) due to the seismic reflection low amplitude and high frequency. The high amplitude (high acoustic Figure 13 shows the sweetness attribute computed from the 3D seismic volume, corresponding cross-sections (e.g., A–A0 , B–B0 , C–C0 , D–D0 , E–E0 and F–F0 ), and sweetness horizon slices. The high sweetness anomalies (at −1900 to −2300 ms) on the 2D sweetness cross-sections (Figure 13b) and the 3D horizon slices (e.g., G sand and E sand intervals) (Figure 13c) contributed to the high amplitude and low frequency. In contrast, the low sweetness anomalies (at −2300 to −2500 ms) on the 2D sweetness cross-sections and the SW—NE parts of the 3D horizon slices (e.g., the Sembar Formation) due to the seismic reflection low amplitude and high frequency. The high amplitude (high acoustic impedance as opposed to shale) and low-frequency anomalies on the 2D sweetness cross-sections and the 3D sweetness horizon slices represent cleaner and more payable sand zones. These sweet spots suggest the presence of a high proportion of porous sand and seem to hold potential for producing gas in the SW and NE parts of the G and E sand reservoir intervals (Figure 13c). In contrast, areas with low amplitude and high-frequency anomalies within the −2100 to −2300 (ms) in the 2D sweetness cross-sections non-reservoir window (e.g., the Sembar Formation) are shale-prone; the sands here may be interbedded with shale. Although the sweetness attribute effectively distinguishes sand bodies from shale, it does so via the high acoustic impedance contrast between sand and shale.

**Figure 13.** (**a**) 3D sweetness attribute volume, (**b**) 2D extracted sweetness cross-sections, and (**c**) 3D sweetness horizon slices of the G sand interval, the E sand interval, and the Sembar Formation. **Figure 13.** (**a**) 3D sweetness attribute volume, (**b**) 2D extracted sweetness cross-sections, and (**c**) 3D sweetness horizon slices of the G sand interval, the E sand interval, and the Sembar Formation.

impedance as opposed to shale) and low-frequency anomalies on the 2D sweetness crosssections and the 3D sweetness horizon slices represent cleaner and more payable sand zones. These sweet spots suggest the presence of a high proportion of porous sand and seem to hold potential for producing gas in the SW and NE parts of the G and E sand reservoir intervals (Figure 13c). In contrast, areas with low amplitude and high-frequency anomalies within the −2100 to −2300 (ms) in the 2D sweetness cross-sections non-reservoir window (e.g., the Sembar Formation) are shale-prone; the sands here may be interbedded with shale. Although the sweetness attribute effectively distinguishes sand bodies from shale, it does so via the high acoustic impedance contrast between sand and shale.

#### 4.4.3. RMS Amplitude Attribute 4.4.3. RMS Amplitude Attribute

Figure 14 shows the results of RMS amplitude inferred from the 3D seismic volume, appropriate 2D cross-sections (e.g., A–A′, B–B′, C–C′, D–D′, E–E′ and F–F′), and 3D horizon slices (e.g., G sand interval, E sand interval, and the Sembar Formation). The extracted 2D RMS amplitude cross-section anomalies range from 0 to 4000 ms. These amplitude variations evaluated the structure's influence on depositional facies. The moderate to high-amplitude anomalies, e.g., between 2500 and 4000 (ms), are often associated with channel sand bodies, high porosity (porous sands), and sand-rich sand shoreward facies, especially gas saturated sand zones. The lateral and horizontal spatial distributions of facies show that faults and fractures significantly reduce the reservoir quality (G and E sand intervals). The gas-saturated sand in the SW and NE parts have high reflectivity, indicating high porosity within the G and E sand intervals (Figure 14c). These sand-rich gas-saturated zones can be considered for future gas exploration in the study area. In Figure 14 shows the results of RMS amplitude inferred from the 3D seismic volume, appropriate 2D cross-sections (e.g., A–A0 , B–B0 , C–C0 , D–D0 , E–E0 and F–F0 ), and 3D horizon slices (e.g., G sand interval, E sand interval, and the Sembar Formation). The extracted 2D RMS amplitude cross-section anomalies range from 0 to 4000 ms. These amplitude variations evaluated the structure's influence on depositional facies. The moderate to high-amplitude anomalies, e.g., between 2500 and 4000 (ms), are often associated with channel sand bodies, high porosity (porous sands), and sand-rich sand shoreward facies, especially gas saturated sand zones. The lateral and horizontal spatial distributions of facies show that faults and fractures significantly reduce the reservoir quality (G and E sand intervals). The gas-saturated sand in the SW and NE parts have high reflectivity, indicating high porosity within the G and E sand intervals (Figure 14c). These sand-rich gas-saturated zones can be considered for future gas exploration in the study area. In comparison, low amplitude anomalies, e.g., between 0 and 2000 (ms), may indicate that these zones contain sandy shale, making them unfavorable for gas production. Such unfavorable zones mainly include the Cretaceous organic-rich shales of the Sembar Formation (Figure 14c).

#### *4.5. Petrophysical Modeling*

Petrophysical modeling unveils the reservoir traits and offers suggestions of hydrocarbonbearing zones. The average petrophysical properties such as *Vshale*, ∅*avg*, ∅*eff*, and *SW* of the G sand interval in both wells are 36.11%, 12.5%, 7.5%, and 45%, respectively (Table 3). Similarly, the derived average petrophysical properties for the E sand interval in both wells are 30.5%, 17.4%, 12.2%, and 33.2%, respectively (Table 4). A graphical representation of these properties is

presented in Figure 15. The overall description of the G and E sand reservoir intervals depends on these petrophysical properties. This may significantly influence the decision-making process in all phases of the planning and executing hydrocarbon activities in the Kadanwari field, MIIB, Pakistan. comparison, low amplitude anomalies, e.g., between 0 and 2000 (ms), may indicate that these zones contain sandy shale, making them unfavorable for gas production. Such unfavorable zones mainly include the Cretaceous organic-rich shales of the Sembar Formation (Figure 14c).

**Figure 14.** (**a**) 3D RMS amplitude volume, (**b**) corresponding 2D extracted cross-sections, and (**c**) 3D RMS amplitude horizon slices of G sand, E sand, and the Sembar Formation. **Figure 14.** (**a**) 3D RMS amplitude volume, (**b**) corresponding 2D extracted cross-sections, and (**c**) 3D RMS amplitude horizon slices of G sand, E sand, and the Sembar Formation.

*4.5. Petrophysical Modeling*  Petrophysical modeling unveils the reservoir traits and offers suggestions of **Table 3.** Petrophysical properties of the Cretaceous G sand interval in the Kadanwari-10 and Kadanwari-11 wells.

significantly influence the decision-making process in all phases of the planning and executing hydrocarbon activities in the Kadanwari field, MIIB, Pakistan. **Table 4.** Petrophysical properties of the Cretaceous E sand interval in the Kadanwari-10 and Kadanwari-11 wells.

Kadanwari-10 G sand interval 36.11 7.8 12.2 45.4 Kadanwari-11 G sand interval 36.11 8 12.56 45.4 Kadanwari-11 wells.

**(***Vshale***) %**

Kadanwari-11 E sand interval 34.05 11.3 16.7 36.42

**Well No. Intervals Volume of shale**

**Table 4.** Petrophysical properties of the Cretaceous E sand interval in the Kadanwari-10 and

**Average Porosity (**∅*avg***) %**

**Water Saturation (***Sw***) %**

**Effective Porosity (**∅*eff***) %**

**Figure 15.** Relationship between the *Vshale* and volume of sand (track 1), distribution of ∅*avg* (track 2), distribution of ∅*eff* (track 3), and the relationship between *SW* and *Shc* with respect to the depth at the G sand interval in (**a**) Kadanwari 10 and (**b**) Kadanwari-11 wells, and at the E sand interval in (**c**) Kadanwari 10 and (**d**) Kadanwari-11 wells.

The GR log response is sensitive to the radioactive emissions predominantly concentrated in the clay minerals of shale and clean sand (feldspar-rich) [43]. The GR response in track 1 of each understudy well confirms the reservoir lithology to be sandstone (Figure 15a). Figure 15a represents the relationship between *Vshale* and sand content (track 1), ∅*avg* distribution (track 2), ∅*eff* distribution (track 3), and the relation between S*W* and S*hc* (track 4) at the G sand interval in the Kadanwari-10 and 11 wells. Similarly, Figure 15b represents the relationship between *Vshale* and sand content (track 1), ∅*avg* distribution (track 2), ∅*eff* distribution (track 3), and the relation between *SW* and *Shc* (track 4) at the E sand interval in the Kadanwari-10 and 11 wells.

#### **5. Discussion**

#### *5.1. Integration of 3D SM and JGC for Hydrocarbon Evaluation*

Knowledge of structural, lithological, and petrophysical characteristics is essential to reducing the uncertainties associated with reservoir description. Herein, 3D seismic data and well logs data have been critically analyzed to evaluate the complex and heterogeneous depositional environment of the early Cretaceous stratigraphic sequence, along with structural and stratigraphic characteristics, distribution of associated petrophysical properties, and spatial facies for reliable reservoir characterization. The result indicates that the early Cretaceous stratigraphic sequence appears irregularly, with faulted structural highs bounded by extension-related normal faults, heterogeneous nature of reservoir intervals, and potential gas-saturated zones.

The structural interpretation of the early Cretaceous stratigraphic sequence indicates that the interpreted horizons (e.g., G sand, E sand, and the Sembar Formation) appear normal faulted extension-related horsts, half-graben, and graben structures (Figure 11). The normal faulted inversion-related horsts, half-graben, and graben structures show how the sedimentary layer's structural features have contributed to the formation of traps and conduit mechanisms. In general, due to the exact kinematic mechanisms involved in the sequential deformation of the Kadanwari field, the structural deformation of each horizon is equal in orientation and extent (Figure 11b). However, the E sand interval seems to be more deformed than the G sand interval and the Sembar Formation. The geometric tendency of the extension-related normal fault in the early Cretaceous stratigraphic sequence is essential because it has created channels for hydrocarbon migration in the vertical direction. Hydrocarbon enrichment in the G and E sand reservoir intervals was positively related to the complexity of the internal structure of normal faults. The hydrocarbon concentration and dominance of these normal faults are more intense towards the central part of the seismic volume, where the faults are closely spaced. The normal fault plane profiles show the hanging wall and footwall cutoffs (Figure 11b). These profiles are useful for understanding the seal behavior for hydrocarbon potential and prospect evaluation.

Using 3D SM, we determined that the extension–related normal faults' geometrical trend created an essential pathway for the up-dip migration of hydrocarbons from the source rock (the Sembar Formation) towards the G and E sand reservoir intervals, resulting in hydrocarbon accumulation. These statements are validated here, as the moderate to high RMS amplitude anomalies are often associated with gas–saturated sand zones (Figure 14c). The hydrocarbon prospects have high reflectivity, indicating high porosity, which can be seen on both the sweetness (Figure 13b) and RMS amplitude (Figure 14b) attributes within the G and E sand intervals. In addition, the RMS amplitude attribute is advantageous compared with the sweetness attribute because the RMS amplitude attribute has a higher resolution when depicting the porous zones and DHIs. These results can significantly reduce the risk associated with hydrocarbon exploration and development in the Kadanwari field, MIB.

The calculated values of the *Shc* form the basis of future production forecasts and the determination of the economic viability of a discovered reservoir. Therefore, high accuracy is needed when determining *SW*, as it is used to calculate the estimated *Shc* reserves. The low values found for *Vshale* content in the Kadanwari 10 and 11 wells indicate the cleanliness of the sandstone. Accordingly, the shale and sandstone facies were separated by a 40% cutoff value in the targeted reservoir intervals, i.e., G and E sands. The *Vshale* contents in the G and E sand intervals are influenced by clay minerals, which reduce the ∅*avg* and ∅*eff*. The high values of ∅*eff* refer to better volume estimation, and thus, a theoretically good reservoir, and vice versa. These high values of ∅*eff* indicate the amounts of connected pore spaces in the reservoir intervals [43]. The G and E sand reservoir interval properties, as derived from our petrophysical modeling, are in good agreement with the results of the regional study conducted by [40,43].

#### *5.2. Analysis of Gas Reserve of the Kadanwari Field*

Graphical representations of the *Shc* at the G and E sand reservoir intervals are presented in track 4 in Figure 15. The petrophysical modeling results show that these intervals have good *Shc*, as the *SW* is below 60% (Table 5). A comparison between the G and E sand interval petrophysical properties shows that the E sand interval has good ∅*eff* and *Shc*, with an apparent gas effect revealed by the cross-overs of the DT and NPHI logs curves (Figure 6). Therefore, it can be considered as an economically viable reservoir interval. The *Shc* percentage of the E sand interval is satisfactory for exploration purposes.

**Table 5.** Depth range, thickness, and *Shc* (%) of the Cretaceous G and E sand intervals in Kadanwari-10 and 11 wells.

The annual report (2010–2011) of petroleum exploration and production activities in Pakistan analyzed the reserves of the Kadanwari field (Table 6) [7]. The Kadanwari field exhibits a significant amount of original recoverable gas reserves, i.e., 1110 billion cubic feet (Bcf) equivalent to 190 million barrels of oil equivalent (Mboe), respectively. From these original recoverable gas reserves, 420 (Bcf) has been extracted from the Kadanwari field. As of 30 June 2011, the balance recoverable gas reserves were 690 Bcf, equivalent to 110 Mboe. The collective original recoverable reserves of the Kadanwari field—280 Mboe—are now limited to 129 Mboe [7]. Thus, many hydrocarbon reserves are present in the Kadanwari field. The high level of gas reserves in the Kadanwari field compared to other fields (such as the Miano field) may be attributed to its complex structural configuration (e.g., negative flower structure of fault system), reservoir compartmentalization, and the up-dip migration of the hydrocarbons into the reservoir intervals (e.g., G and E sands).

**Table 6.** The gas reserves of the Kadanwari and Miano fields as of 30 June 2011.


**2. Materials and Methods**

Participants were recruited from social media posts (e.g., Facebook, Twitter) and email list serves from 26 August to 17 September 2020. An image, with a link imbedded, describing the study was posted to social media. (Example post, "If you live in a rural area, please take our survey for a chance to win an Amazon gift card!") Various local, regional, and national advocacy groups distributed social media posts and emails. Emails contained a brief description of the study with a link to the survey's consent form. Social media posts contained a link embedded in an image describing the study—when the potential participant clicked the image, they were routed to the consent form. If potential participants agreed, they were then taken to the survey. All responses were anonymous unless they wanted to be entered into the raffle for a gift card in which the participant would provide contact information. Inclusion criteria for participation were ≥18 years of age and the ability to read English. All data were collected in Qualtrics software [28]. The Washington University in St. Louis Institutional Review Board approved the study with an exempt status.