Estimation of Litho-Fluid Facies Distribution from Zero-Offset Acoustic and Shear Impedances

Abstract

Seismic data are considered crucial sources of data that help identify the litho-fluid facies distributions in reservoir rocks. However, different facies mostly have similar responses to seismic attributes. In addition, seismic anisotropy negatively affects the facies predictors extracted from seismic data. Accordingly, this study aims at estimating zero-offset acoustic and shear impedances based on partial-stack inversion by two methods: statistical modeling and a multilayer feed-forward neural network (MLFN). The resulting impedance volumes are compared to those obtained from isotropic simultaneous inversion by using impedance logs. The best impedance volumes are applied to Thomsen’s anisotropy equations to solve for the anisotropy parameters Epsilon and Delta. Finally, the shear and acoustic impedances are transformed into elastic properties from which the facies and fluid distributions are predicted by using the logistic regression and decision tree algorithms. The results obtained from the MLFN show better matching with the impedance and facies logs compared to those obtained from isotropic inversion and statistical modeling.

1. Introduction

Seismic data have proved to be important quantitative tools that can bring valuable information about rock properties, especially when being integrated with statistical tools. For example, multiattribute transforms were used to predict log properties from seismic data [1]. Another common approach is to invert seismic data into elastic properties and then into petrophysical properties by statistical relationships obtained from logging data [2]. Porosity, lithology, and fluid properties were also estimated with the resampling of rock physics constraints [3,4,5]. Russel used multivariate statistics and neural networks to predict reservoir parameters from seismic attributes [6]. In addition, Bachrach inverted porosity and water saturation by stochastic rock physics modeling [7].

Other methodologies consider reservoir parameter inversion based on Bayesian classification theory [8,9]. Moreover, a geostatistical inversion method was used to predict reservoir properties in a higher certainty than that of deterministic methods [10]. In addition, some studies derived elastic properties from prestack seismic inversion to define facies [11,12]. Another straightforward methodology is to combine petroelastic modeling and stochastic inversion to predict hydrocarbon zones based on Bayesian theory [13]. A rock physics model was postulated to derive the elastic rock properties from mineral parameters and structure information [14]. Recently, rock properties were estimated from angle stack seismic data by a Bayesian inversion based on the Gibbs sampling algorithm [15].

Machine learning (ML) tools have been widely used in reservoir characterization [16]. For instance, formation porosity and permeability were estimated from logging and core data by using an artificial neural network (ANN) [17,18,19]. Other studies also used neural networks to predict essential rock properties such as permeability [20,21,22] and shear wave velocity [23,24,25,26]. Neural networks can be combined with seismic attributes to predict rock physics properties from seismic data [27,28,29,30]. In addition, facies distributions were estimated by both supervised and unsupervised ML algorithms such as ANN [31,32], convolutional neural network (CNN) [33,34,35], support vector machine (SVM) [36,37,38,39,40,41], bagged tree (BT) [42,43], relevance vector machine (RVM) [44], seed region growing (SRG) [45], self-organizing map (SOM) [46,47], principal component analysis (PCA) [48,49,50], and generative topographic mapping (GTM) [51]. Some challenges are common in ML facies models such as over-fitting and over-parametrization. Moreover, such models should consider the limitations of seismic data to eliminate error as much as possible. One of the problems that contributes to the uncertainty of seismic-derived methods is seismic anisotropy.

This study aimed at forecasting zero-offset P-impedance ($Z_P$) and S-impedance ($Z_S$) by using statistical modeling and MLFN after applying a partial-log-constrained inversion to the near, mid, and far angle stacks. The outputs of this objective were zero-offset impedances that were more accurate than those obtained from isotropic simultaneous inversion. Another objective was to compare the abilities of both statistical modeling and neural networks to turn the randomness of the data into identified patterns. The next step was to apply the best $Z_P$ and $Z_S$ volumes to Thomsen’s anisotropy equations to solve for the anisotropy parameters, Epsilon and Delta, assuming a vertical transverse isotropic (VTI) medium. Next, the zero-offset $Z_P$ and $Z_S$ were used as inputs to a facies model, which predicted the distribution of three classes: gas sand, wet sand, and shale. The idea of the facies model is to forecast the sand and gas probabilities by logistic regression [52] and then use them as inputs for a decision tree to predict the distribution of each facies class.

The three offset stacks were converted into angle stacks, which were the inputs of the partial-stack inversion, and angle gathers, which were the inputs of the simultaneous isotropic inversion. The ranges of the near, mid, and far stacks were from 5° to 15°, 15° to 25°, and 25${}^{\circ }$ to 40°. The partial-stack inversion resulted in the $Z_P$ and $Z_S$ volumes at the central angles, 10°, 20°, and 32.5°, from which the zero-offset properties were forecasted. The target zone consisted of two gas-bearing shaly sand reservoirs in the Malay basin: A and B. The training data of the Bayesian and MLFN models were gathered from the wells (×1) and (×2); however, the facies log was only obtained at the well (×1) due to the availability of the petrophysical cut-offs and oil-mud resistivity imaging data. Another challenge was that the $Z_S$ log was only available at the well (×1), so it was estimated from other log properties by the MLFN neural network at the well (×2).

The resulting $Z_P$ and $Z_S$ were compared to those obtained from isotropic inversion by using the impedance logs. The MLFN resulted in the most accurate zero-offset impedances. Therefore, the $Z_P$ and $Z_S$ volumes were applied to Thomsen’s anisotropy equations to solve for the anisotropy parameters Epsilon and Delta. The $Z_P$ and $Z_S$ of MLFN and isotropic inversion were transformed into the near and far elastic impedances, which were the lithology predictors, and then into the Mu–Rho (MR) ratio, Lambda–Rho/Mu–Rho (LR/MR) ratio, and Poisson’s ratio (PR), which are the fluid predictors. The litho-fluid model was applied to the inverted elastic volumes to forecast the litho-fluid facies distribution. The facies log was used to validate the predicted facies obtained from the isotropic and anisotropic approaches.

2. Seismic Anisotropy

Seismic anisotropy is the dependence of the seismic velocity on the incident angle [53,54,55]. Wave propagation in anisotropic media has been discussed by various studies [53,56,57,58]. To understand the seismic anisotropy of the Earth’s subsurface, the relationship between the internal forces (stress) and deformation occurring in the medium (strain) should be defined. When an elastic wave propagates through the subsurface, the stress–strain relationship can be obtained according to Hooke’s law as shown below [59]:

$${\sigma }_{ij}=c_{ijkl}{e}_{kl}$$

(1)

where ${\sigma }_{ij}$ is the stress tensor, $e_{kl}$ is the strain tensor, and $c_{ijkl}$ is the elastic (stiffness) tensor, which is a fourth-order tensor with 81 independent components. Both the strain and stress tensors are symmetric [60], resulting in:

$$c_{ijkl}=c_{jikl}=c_{ijlk}$$

(2)

The 81 independent components are accordingly reduced to 36 due to the symmetry in the stress and strain tensors. Therefore, the elastic tensor ($c_{ijlk}$) is replaced with a 6 × 6 matrix ($c_{\alpha \beta }$), where [61]:

$$\alpha =ij=ji,\beta =kl=lk$$

(3)

where $\alpha$ and $\beta$ are the new indices of the elastic tensor. The resulting (6 × 6) notation matrix is expressed as shown below:

$$\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& C_{14}& C_{15}& C_{16}\\ C_{12}& C_{22}& C_{23}& C_{24}& C_{25}& C_{26}\\ C_{13}& C_{23}& C_{33}& C_{34}& C_{35}& C_{36}\\ C_{14}& C_{24}& C_{34}& C_{44}& C_{45}& C_{46}\\ C_{15}& C_{25}& C_{35}& C_{45}& C_{55}& C_{56}\\ C_{16}& C_{26}& C_{36}& C_{46}& C_{56}& C_{66}\end{array}\right]$$

(4)

There are several types of symmetry for the independent components that determine the degree of anisotropy. The maximum anisotropy degree occurs at the lowest possible symmetry, which is called “triclinic” symmetry and occurs when the ($C_{\alpha \beta }$) matrix has 21 different nonzero components, which is considered the maximum number of elements required to describe an anisotropic medium [58]. On the other hand, the isotropic case occurs at the highest symmetry case with a smaller amount of independent components.

The most realistic symmetry case of the Earth’s subsurface is polar symmetry, which has a single pole of rotational symmetry (Z-axis) that is different from the other two axes (X and Y) [61]. This anisotropy case is called transverse isotropy (TI), at which rock properties are assumed invariant about an axis of symmetry. If the symmetry axis is vertical, this means that the layers are considered horizontal and is hence called vertical transverse isotropy (VTI). Another case is called tilted transverse isotropy (TTI), which is a more realistic case because it follows the variations of the dip angles of layers. Polar symmetry can be described by a matrix of five different independent components, as shown below:

$$\left[\begin{array}{cccccc}{C}_{11}& C_{11}-2C_{66}& C_{13}& 0& 0& 0\\ C_{11}-2C_{66}& C_{11}& C_{13}& 0& 0& 0\\ C_{13}& C_{13}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right]$$

(5)

There are three anisotropy parameters that can define any anisotropy degree in terms of the stiffness components of a TI medium [53]. The anisotropy parameters are called Epsilon ($ϵ$), Delta ($\delta$), and Gamma ($\gamma$). The weaker the anisotropy degree is, the closer these parameters are to zero. The anisotropy parameters in the weak anisotropy case have values less than 1.

Epsilon is the fractional difference between the vertical and horizontal P-wave velocities. It generally has a positive sign because rocks are less compressible along bedding than across bedding. It can be obtained by the following equation:

$$ϵ=\frac{V_{PH}-V_{PV}}{V_{PV}}=\frac{C_{11}-C_{33}}{2C_{33}}$$

(6)

where $ϵ$ is the anisotropy parameter Epsilon, $V_{PH}$ is the horizontal P-wave velocity, $V_{PV}$ is the vertical P-wave velocity, and $C_{11}$ and $C_{33}$ are the stiffness components of the elastic tensor.

The Delta parameter controls the near-vertical anisotropy and does not relate to the horizontal velocity, so it can be positive or negative. Assuming weak anisotropy, the Delta parameter refers to the stiffness components by the following equation:

$$\delta =\frac{{(C_{13}+C_{44})}^2-{(C_{33}-C_{44})}^2}{2C_{33}(C_{33}-C_{44})}$$

(7)

where $\delta$ is the anisotropy parameter Delta, while $C_{13}$, $C_{44}$, and $C_{33}$ are the stiffness components of the elastic tensor.

The Gamma parameter is the fractional difference between the vertical and horizontal S-wave velocities and can be given as follows:

$$\gamma =\frac{V_{SH}-V_{SV}}{V_{SV}}=\frac{C_{66}-C_{44}}{2C_{44}}$$

(8)

where $\gamma$ is the anisotropy parameter Gamma, $V_{SH}$ is the horizontal S-wave velocity, $V_{SV}$ is the vertical S-wave velocity, and $C_{66}$ and $C_{44}$ are the stiffness components of the elastic tensor.

The P- and S-wave velocities can be defined in terms of the anisotropy parameters based on the weak anisotropy assumption [53,61], as shown below:

$$V_P\left(\theta \right)=V_{P0}\left[1+\delta si{n}^2\left(\theta \right)co{s}^2\left(\theta \right)+ϵsi{n}^4\left(\theta \right)\right]$$

(9)

$$V_{SV}\left(\theta \right)=V_{S0}[1+{\left(\frac{V_{P0}}{V_{S0}}\right)}^2(ϵ-\delta )si{n}^2\left(\theta \right)co{s}^2\left(\theta \right)]$$

(10)

$$V_{SH}\left(\theta \right)=V_{S0}(1+\gamma si{n}^2\left(\theta \right))$$

(11)

where $V_P(\theta$) and $V_{SV}(\theta$) are the vertical P-wave and S-wave velocities, respectively, $V_{SH}(\theta$) is the horizontal S-wave velocity, $V_{P0}$ and $V_{S0}$ are the zero-offset P-wave and S-wave velocities, respectively, at $\theta$ = 0, and $\theta$ is the incident angle.

3. The Effect of Seismic Anisotropy on Seismic Data

Seismic anisotropy causes a misinterpretation of seismic data, which affects the final reservoir model. For example, if seismic anisotropy is neglected, a target may be expected at the wrong depth. Another problem is that seismic anisotropy affects the angle-dependent reflectivity calculation, which is further used for prestack seismic applications such as amplitude versus offset (AVO) and elastic impedance (EI). Figure 1 shows a comparison between the isotropic and anisotropic assumptions, which yield significantly different reflectivities at the top of the C-sand reservoir in the Sawan gas field [62]. In addition to that, seismic anisotropy adds uncertainty to prestack seismic inversion, which results in inaccurate rock properties.

A common workflow of anisotropic seismic inversion is to obtain the anisotropy parameters and then process seismic data before the inversion process. The anisotropy parameters can be obtained from various types of data. However, the most confident sources of anisotropy quantification are core data. Traditional laboratory ultrasonic measurements of intrinsic anisotropy [63] depend on a limited number of orientations, resulting in inaccurate anisotropy measurements. On the other hand, more effective methods [64,65,66,67] depend on many ray paths having various inclination angles to the bedding plane.

A numerical study was applied to ultrasonic measurements according to Markov chain Monte Carlo simulation (MCMC) [68]. The idea of this study was to obtain the probability distribution curves of the anisotropy parameters based on the layering effects of five lithofacies: dry sandstone, wet sandstone, dry carbonate, wet carbonate, and shale. Accordingly, the anisotropy parameters were simulated for each facies class. The Backus averaging method helped calculate the harmonic means of the anisotropy parameters, assuming a horizontally layered inhomogeneous medium [69].

Another source of anisotropy information is the vertical seismic profile (VSP). For instance, one study aimed at obtaining the horizontal slowness (Sx) as well as the vertical slowness (Sz) and then combining them in the vector way to yield the scalar slowness [70,71]. This method results in accurate arrival time and depth values, but it assumes that layers are homogeneous and horizontal, which is not realistic [61].

The Delta parameter was estimated by using both prestack time data and sonic logs [72]. Moreover, a method, called The Empirical Relationship relates the anisotropy parameters to the shale volume (Vsh) that causes the intrinsic anisotropy [73,74,75]. The volume fraction and degree of compaction are assumed to be the controlling factors of the orientation of clay minerals and hence the anisotropy parameters. This method is convenient for heterogeneous media because it solves for the intrinsic anisotropy components. However, the problem of obtaining an accurate zero-offset velocity away from wells still exists. Another study considered the intrinsic anisotropy in shale based on an orientation distribution function [76]. This study aimed at modeling anisotropic elastic properties based on two approaches. The first approach was to apply an anisotropic differential effective medium model (DEM) based on the textural information about the kerogen network in kerogen-rich shale. The second approach was to use the DEM model with a compaction-dependent orientation distribution function to study the anisotropy in laminated shaly sand rocks.

The anisotropic elastic constants were estimated as functions of the fluid-filled porosity and the aspect ratio of the inclusions [77]. Another study introduced an anisotropic dual porosity (ADP) model to forecast the elastic properties of shaly sand rocks based on a combination of the anisotropic formulations of self-consistent approximation (SCA), DEM, and anisotropic Gassmann theory [78].

Kelter also estimated both Epsilon and Delta from synthetic modeling through some algorithms such as neural networks and regridding inversion [56]. The idea of this method is to model synthetics by using a PP and PS survey such that layer boundaries and material properties are predefined, and then the survey is simulated by the ray tracing method to generate a synthetic model from which the anisotropy parameters are forecasted. The best parameter estimates are selected according to the comparison between the synthetic model and the real seismic data. This method gives reliable results but only at well locations.

Another method was set by Liner and Fe [79], who estimated Thomsen’s parameters from routine well logging data, following a theoretical framework which uses sonic logs to provide a fine-layered isotropic elastic model that can be used to calculate anisotropy parameters as functions of depth and seismic wavelength [80]. Thomsen’s parameters can be also estimated from the average velocity, obtained from sonic data, and the mean inclination angles of deviated wells [81]. However, the results of that study were inaccurate because most of the wells’ inclination angles were below 5${}^{\circ }$, showing no anisotropy effect on seismic data. The residual moveout analysis of diving waves has been used to estimate anisotropy parameters by applying different parametrization scenarios [82]. This method is reasonably accurate even for large values of velocity gradients, but it assumes that the velocity increases linearly with depth, while the anisotropy parameters do not. However, this assumption can only apply to homogeneous media.

The previous methods can be used to obtain the anisotropy parameters and then process seismic data to obtain zero-offset images, which can be inverted into accurate rock properties. However, a different approach is to transform the layer properties, obtained from isotropic inversion, into zero-offset parameters based on Ruger’s reflectivity function in horizontal transverse isotropy (HTI) media [83]. Another anisotropic prestack seismic inversion [84] is based on the Markov random field (MRF), which can establish dependencies between the spacial wave propagation field’s nodes based on prior constraints [85]. The weights of those constraints have been optimized by the anisotropic diffusion method to remove the anisotropy effect in three different models: the symmetrical, vertical asymmetrical, and asymmetrical models. This method has produced encouraging results for Vp, Vs, and density.

Based on the elastic wave equation, the anisotropic elastic waveform inversion was applied to invert for the P- and S-wave velocities, density, and Thomsen’s parameters [86]. The method is based on the idea of eliminating the error between the synthetic and field seismic waveform data according to a misfit function [87]. The anisotropic inversion outputs were successfully used to reduce the ground-roll noise in 2D seismic data.

Another methodology depends on multiscale phase inversion, assuming an anisotropic medium [88]. The advantage of this method is that it avoids the nonlinearity of the misfit functions [89] by adding wave-number details to the velocity model by a skeletonized inversion. This method assumes a VTI medium with Delta equal to zero, while the zero-offset P-wave velocity and Epsilon are inverted, resulting in more accurate results compared to full-wave inversion.

It can be concluded that most of the previous methods have quantified seismic anisotropy only at well locations. The interpolation of the anisotropy parameters between wells is challenging and should be further linked to the effect of the facies distribution. Therefore, the current study has improved a previous anisotropic inversion approach [90] which solves for the zero-offset impedances and anisotropy parameters, Epsilon and Delta, based on partial-stack inversion and Thomsen’s anisotropy equations [91]. On the other hand, our study estimates $Z_P$ and $Z_S$ based on statistical modeling and MLFN rather than calculating them empirically.

The advantage of ML algorithms is that they reduce the random error of partial-stack inversion and result in accurate impedances that best match the logs. Another contribution is that the facies model uses logistic regression to model lithologies and fluids and then combines them using the decision tree algorithm. The integration of two ML algorithms adds stability to the model and raises the chance of being applied to other fields. Moreover, the model includes the depth variable, which acts as a trend factor that highly affects the facies distribution.

4. Methodology

The study workflow, as shown in Figure 2, begins with inverting the angle gathers into $Z_P$, $Z_S$, and density volumes and inverting the angle stacks into $Z_P$ and $Z_S$ volumes at the near, mid, and far incident angles. Next, the partially inverted impedances are used to forecast the zero-offset impedances by the statistical and MLFN models. The impedance volumes are then compared to each other by using the impedance logs. The best $Z_P$ and $Z_S$ volumes are applied to Thomsen’s anisotropy equations to solve for the anisotropy parameters, Epsilon and Delta. The next step is to transform the impedance and density volumes into the elastic volumes, from which the sand and gas probabilities are estimated by using the logistic regression model. These probabilities are then inserted into a decision tree model to forecast the distribution of the facies classes: gas sand, wet sand, and shale. The inversion processes in our study were carried out by using the Hampson–Russel software provided by the CGG company.

This study is not applicable to poststack data. The impedance model’s accuracy depends on the seismic data quality as well as the angle gathers’ range. The wider the angle range, the more enhanced signal-to-noise ratio and hence the more accurate the model. Moreover, both the P- and S-impedances should be available in order to solve for the anisotropy parameters: Epsilon and Delta. The study assumes a VTI medium which neglects some anisotropy-related factors such as the layers’ dip angles and fractures; however, impedance modeling by MLFNs accurately solves for the zero-offset impedances by normalizing the error and mitigating the anisotropy effect simultaneously. The study’s facies model is only applicable to the three facies classes: gas sand, wet sand, and shale. However, the model can be extended to include more classes. Moreover, the applicability of the model to other fields can be enhanced if there are more facies data available from a large number of wells.

4.1. Isotropic Simultaneous Inversion

Simultaneous inversion is a prestack seismic inversion that simultaneously inverts the angle gathers into the reflectivities of the $Z_P$, $Z_S$, and density based on the linear relationships between the $Z_P$ and each of the $Z_S$ and density. The methodology used in this study [92] directly inverted for the $Z_P$, $Z_S$, and density [93,94].

The Aki–Richards equation [95] is an approximation of the relationship between the P-wave reflectivity between two media and the incident angle as follows:

$$R\left(\theta \right)=\frac{1}{2}(\Delta V/V+\Delta \rho /\rho )-2{(W/V)}^2(2\Delta W/W+\Delta \rho /\rho )si{n}^2\theta +\frac{1}{2}(\Delta V/V)ta{n}^2\theta$$

(12)

where $R\left(\theta \right)$ is the P-wave reflectivity at an incident angle ($\theta$), while V, W, and $\rho$ are the average P-wave velocity, S-wave velocity, and density between the two media, respectively. The idea of Fatti’s equation [93] is to re-express the Aki–Richards equation in (12) in terms of the acoustic and shear impedances as follows:

$$R\left(\theta \right)=\frac{1}{2}\frac{\Delta I}{I}(1+ta{n}^2\theta )-4{\left(\frac{W}{V}\right)}^2\left(\frac{\Delta J}{J}\right)si{n}^2\theta -\left[\frac{1}{2}\left(\frac{\Delta \rho }{\rho }\right)ta{n}^2\theta -2{\left(\frac{W}{V}\right)}^2\left(\frac{\Delta \rho }{\rho }\right)si{n}^2\theta \right]$$

(13)

where I is the acoustic impedance, J is the shear impedance, the term ($\frac{1}{2}\frac{\Delta I}{I}$) is the zero-offset P-wave reflectivity, and the term ($\frac{1}{2}\frac{\Delta J}{J}$) is the zero-offset S-wave reflectivity, such that:

$$\frac{1}{2}\frac{\Delta J}{J}=\frac{1}{2}(\Delta V/V+\Delta \rho /\rho ),\frac{1}{2}\frac{\Delta J}{J}=\frac{1}{2}(\Delta W/W+\Delta \rho /\rho )$$

(14)

Fatti’s equation could be re-expressed by replacing the reflectivities by the rock properties as follows:

$$s\left(\theta \right)=\frac{1}{2}{C}_1{W}_{\theta }Dln\left(Z_P\right)+\frac{1}{2}{C}_2{W}_{\theta }Dln\left(Z_S\right)+C_3{W}_{\theta }Dln\left(\rho \right)$$

(15)

where s($\theta$) is the angle trace, W($\theta$) is the angle-dependent wavelet, D is a derivative operator, $Z_P$ is the acoustic impedance, $Z_S$ is the shear impedance, $\rho$ is the density, and $\theta$ is the incident angle. The terms $C_1$, $C_2$, and $C_3$ are given by:

$$C_1=1+ta{n}^2\theta ,C_2=-8\left(\frac{V_S}{V_P}\right)ta{n}^2\theta ,C_3=-0.5ta{n}^2\theta +2{\left(\frac{V_S}{V_P}\right)}^{2}si{n}^2\theta$$

(16)

where $C_1$, $C_2$, and $C_3$ are constants, $\theta$ is the incident angle, $V_S$ is the S-wave velocity, and $V_P$ is the P-wave velocity.

Next, the terms $ln$($Z_P$), $ln$($Z_S$), and $ln$($\rho$) were coupled by the following equations:

$$ln\left(Z_s\right)=kln\left(Z_p\right)+k_c+\Delta ln\left(Z_s\right)$$

(17)

$$ln\left(\rho \right)=mln\left(Z_p\right)+m_c+\Delta ln\left(\rho \right)$$

(18)

where k and m are the intercepts of the linear relations and $k_c$ and $m_c$ are constants, while $\Delta$$ln\left(Z_S\right)$ and $\Delta$$ln$($\rho$) are controlled by the deviations away from the best-fit lines of the linear relationships. $\Delta$$L_D$ and $\Delta$$L_S$ are the displacements of $ln$($\rho$) and $ln$($Z_S$) due to the deviation from the best-fit lines.

The next step was to apply the derivative operator (D) to Equations (3) and (4):

$$D\mathit{ln}{Z}_S=kD\mathit{ln}{Z}_P+D\Delta \mathit{ln}{Z}_S$$

(19)

$$D\mathit{ln}\rho =mD\mathit{ln}{Z}_P+D\Delta \mathit{ln}\rho$$

(20)

where D is the derivative operator, while k and m are the intecepts of the linear relations. Eventually, the algorithm was applied to an initial guess and then iterated until it reached a solution at which the initial model matched the real data with the least-squared error.

4.2. Partial-Log-Constrained Inversion

The separate inversion of angle stacks is based on a partial inversion approach that obtains the zero-offset P- and S-wave velocities from the inverted P- and S-waves at different angular ranges [90]. The idea is based on Thomsen’s equations in (9) and (10), which relate the vertical P-wave ($V_P$) and S-wave ($V_S$) velocities to the anisotropy parameters, Epsilon and Delta, in VTI media.

Partial-log-constrained inversion requires a wavelet to be extracted from each angle stack so that the wavelet is centered at the average angle of the angle stack [90]. An initial model was created from logging data to guide the inversion. This method required S-wave velocity logs in the initial model to constrain the inversion. Each partial seismic volume was inverted into a velocity volume at the central angle of its angle range. Next, the velocity-domain Thomsen’s equations in (9) and (10) were used to solve for the zero-offset P- and S-wave velocities as well as the Epsilon and Delta parameters, as shown in the following equations:

$$\left[\begin{array}{ccc}1& si{n}^2\left({\theta }_1\right)co{s}^2\left({\theta }_1\right)& si{n}^4\left({\theta }_1\right)\\ 1& si{n}^2\left({\theta }_2\right)co{s}^2\left({\theta }_2\right)& si{n}^4\left({\theta }_2\right)\\ ⋮& ⋮& ⋮\\ 1& si{n}^2\left({\theta }_n\right)co{s}^2\left({\theta }_n\right)& si{n}^4\left({\theta }_n\right)\end{array}\right]\left[\begin{array}{c}{V}_{P0}\\ V_{P0}\delta \\ V_{P0}ϵ\end{array}\right]=\left[\begin{array}{c}{V}_{qP}\left({\theta }_1\right)\\ V_{qP}\left({\theta }_2\right)\\ ⋮\\ V_{qP}\left({\theta }_n\right)\end{array}\right]$$

(21)

where $V_{P0}$ is the zero-offset P-wave velocity, $ϵ$ and $\delta$ are the anisotropy parameters, $V_{qP}$ is the quasivertical P-wave velocity, $\theta$ is the incident angle, and n is the number of partial stacks.

$$\left[\begin{array}{cc}1& si{n}^2\left({\theta }_1\right)co{s}^2\left({\theta }_1\right)\\ 1& si{n}^2\left({\theta }_2\right)co{s}^2\left({\theta }_2\right)\\ ⋮& ⋮\\ 1& si{n}^2\left({\theta }_n\right)co{s}^2\left({\theta }_n\right)\end{array}\right]\left[\begin{array}{c}{V}_{s0}\\ \frac{V_{P0}^2(ϵ-\delta )}{V_{s0}}\end{array}\right]=\left[\begin{array}{c}{V}_{qSV}\left({\theta }_1\right)\\ V_{qSV}\left({\theta }_2\right)\\ ⋮\\ V_{qSV}\left({\theta }_n\right)\end{array}\right]$$

(22)

where $V_{S0}$ is the zero-offset P-wave velocity, $ϵ$ and $\delta$ are the anisotropy parameters, $V_{qSV}$ is the quasivertical S-wave velocity, $\theta$ is the incident angle, and n is the number of the partial stacks.

The log-constrained inversion of partial stacks was adopted from the nonlinear optimization inversion [90], such that the objective function is obtained by:

$$f\left(V\right)=\Vert S-D\Vert \to min$$

(23)

where V is the inverted P- or S-wave velocity, D is the angle-domain seismic data, and S is the corresponding prediction, which can be expressed according to the convolution model [96]. The reflectivity function is based on the Aki and Richard approximation of the angle-dependent reflectivity [97], as shown below:

$$R\left(\theta \right)=R\left(0\right)+Asi{n}^2\left(\theta \right)+Bta{n}^2\left(\theta \right)$$

(24)

where $R\left(\theta \right)$ is the reflectivity at an incident angle $\theta$, while $R\left(0\right)$ is the zero-offset reflectivity, which is given by:

$$R\left(0\right)=\frac{1}{2}\left(\frac{\Delta V_P}{V_P}+\frac{\Delta \rho }{\rho }\right)$$

(25)

$$A=-2{\left(\frac{V_S}{V_P}\right)}^2\left(\frac{2\Delta V_S}{V_S}+\frac{\Delta \rho }{\rho }\right)$$

(26)

$$B=\frac{1}{2}\left(\frac{\Delta V_P}{V_P}\right)$$

(27)

where $V_P$, $V_S$, and $\rho$ are, respectively, the average P-wave velocity, S-wave velocity, and density of the underlying and overlying layers, while $\Delta V_P$, $\Delta V_S$, and $\Delta \rho$ are the differences between these properties at the underlying and overlying strata.

The term $\frac{\Delta \rho }{\rho }$ can be obtained by Gardner’s equation as follows [98]:

$$\frac{\Delta \rho }{\rho }=\frac{1}{4}\left(\frac{\Delta V_P}{V_P}\right)$$

(28)

According to Equations (24)–(28), the reflectivity will be as follows:

$$R\left(\theta \right)=\left(\frac{5}{8}-\frac{1}{2}\frac{V_S^2}{V_P^2}si{n}^2\theta +\frac{1}{2}ta{n}^2\theta \right)\frac{\Delta V_P}{V_P}-4\frac{V_S^2}{V_P^2}si{n}^2\theta \frac{\Delta V_S}{V_S}$$

(29)

After inverting the partial stacks separately into $V_P$ and $V_S$ volumes at the near, mid, and far angles, the zero-offset velocities Epsilon and Delta were solved by Equations (21) and (22).

The current study applied the same methodology to $Z_P$ and $Z_S$ by replacing the velocities of Equations (9) and (10) by impedances, as shown below:

$$Z_P\left(\theta \right)=Z_{P0}\left[1+\delta si{n}^2\left(\theta \right)co{s}^2\left(\theta \right)+ϵsi{n}^4\left(\theta \right)\right]$$

(30)

$$Z_S\left(\theta \right)=Z_{S0}[1+{\left(\frac{Z_{P0}}{Z_{S0}}\right)}^2(ϵ-\delta )si{n}^2\left(\theta \right)co{s}^2\left(\theta \right)]$$

(31)

Partial inversion is sensitive to the quality of seismic data because the magnitude of the noise may be more than that of the anisotropy effect. Hence, the zero-offset impedances were estimated from the impedances at the near, mid, and far incident angles by using statistical modeling and MLFN. Seismic data consisted of three 3D volumes of angle ranges, from 5${}^{\circ }$ to 15°, 15${}^{\circ }$ to 25°, and 25${}^{\circ }$ to 40°. Accordingly, three partial stacks were created at the near, mid, and far central incident angles: 10°, 20°, and 32.5°, respectively. The advantage of impedance modeling is that it normalizes the error and reduces the anisotropy effect simultaneously.

An important step prior to the modeling process is to calculate Epsilon and Delta at the well locations in order to have an idea about the degree of anisotropy within the study area and to specify the model’s assumptions accordingly. The anisotropy parameters were calculated by the methodology of Liner and Fei [79], who calculated the layer-induced stiffness parameters from the averaged properties of the isotropic elastic layers and then obtained Epsilon and Delta as functions of those parameters. The stiffness parameters were obtained from the following equations:

$$a={\langle \frac{\lambda }{\lambda +2\mu }\rangle }^2{\langle \frac{1}{\lambda +2\mu }\rangle }^{-1}+4\langle \frac{\mu (\lambda +\mu )}{\lambda +2\mu }\rangle$$

(32)

$$c={\langle \frac{1}{\lambda +2\mu }\rangle }^{-1}$$

(33)

$$f=\langle \frac{1}{\lambda +2\mu }\rangle {\langle \frac{1}{\lambda +2\mu }\rangle }^{-1}$$

(34)

$$l={\langle \frac{1}{\mu }\rangle }^{-1}$$

(35)

where $\langle \rangle$ is the Backus averaging symbol, $\mu$ is the shear modulus, and $\lambda$ is Lame’s constant, while $a,c,f,$ and l are the stiffness parameters. Next, Epsilon and Delta were obtained from the following equations:

$$ϵ=\frac{a-c}{2c}$$

(36)

$$\delta =\frac{{\left(f+l\right)}^2-{\left(c-l\right)}^2}{2c(c-l)}$$

(37)

Another premodeling step is to explore the relations between the zero-offset impedances, which are the $Z_P$ and $Z_S$ logs, and the angle impedances, which are the composite traces obtained from the partial-stack inversion. Figure 3 shows a fair linear relationship based on the data of wells (×1) and (×2), which consist of 1292 observations for $Z_P$ 1582 observations for $Z_S$.

The linearity between variables is proven visually through the best-fit lines (in blue). In addition, the Pearson correlation coefficient (R-value) was calculated in Python using NumPy. The R-values show strong positive correlations at the near and mid angles, while the correlation is moderate at the far angle. This is matching with Thomsen’s equations in (9) and (10), which consider a linear relationship between the zero-offset velocity and the velocity at an angle theta.

Accordingly, a linear regression model was firstly created to forecast the zero-offset impedance from the partially inverted impedances based on the MCMC simulation. As the R-value between each of the P- and S-impedances and their far-angle impedances was moderate due to the anisotropy effect, the MLFN was created to better estimate the zero-offset impedances and turn the randomness in the data points into recognized patterns.

4.3. Zero-Offset Impedance Modeling

4.3.1. Statistical Modeling

The idea of the model is to fit a linear relationship between the zero-offset impedance, which is the response variable, and the partially inverted impedances, which are the explanatory variables. The output of the model consists of the posterior means of the coefficients of the three partially inverted impedances. The general form of the model is:

$$Im{p}_0=b_{1}Im{p}_{\theta 1}+b_{2}Im{p}_{\theta 2}+b_{3}Im{p}_{\theta 3}$$

(38)

where $Im{p}_0$ is the zero-offset $Z_P$ or $Z_S$, while $b_1$, $b_2$, and $b_3$ are the model’s coefficients and $Im{p}_{\theta 1}$, $Im{p}_{\theta 2}$, and $Im{p}_{\theta 3}$ are the impedances at the near, mid, and far angles, respectively. The model’s parameters are estimated based on Bayesian theory and MCMC simulation. The Bayesian model consists of three components: the likelihood, priors, and posterior. The priors of the model were the parameters of the proposed distribution of the model coefficients. The posterior component was considered the probability distribution function of all realizations of the simulated coefficients.

The linear expression of the $Z_{P0}$ model was obtained from Thomsen’s equation in (9), which can be rewritten as follows:

$$Z_{P0}=b{Z}_P\left(\theta \right)$$

(39)

where:

$$b=\frac{1}{1+\delta si{n}^2\left(\theta \right)co{s}^2\left(\theta \right)+ϵsi{n}^4\left(\theta \right)}$$

(40)

By substituting the near, mid, and far P-impedances, $Z_{P1}$, $Z_{P2}$, and $Z_{P3}$, respectively, to Equation (39), three equations were obtained as follows:

$$Z_{P0}=b_1{Z}_{P1}$$

(41)

$$Z_{P0}=b_2{Z}_{P2}$$

(42)

$$Z_{P0}=b_3{Z}_{P3}$$

(43)

By adding the three Equations (41)–(43), the linear expression will be as follows:

$$Z_{P0}=0.33b_1{Z}_{P1}+0.33b_2{Z}_{P2}+0.33b_3{Z}_{P3}$$

(44)

Similarly, the linear expression of the $Z_{S0}$ model was obtained based on Thomsen’s equation in (10), which can be rewritten as shown below:

$$Z_{S0}=b.Z_S\left(\theta \right)$$

(45)

where

$$b=\frac{1}{1+{\left(\frac{Z_{P0}}{Z_{S0}}\right)}^2(ϵ-\delta )si{n}^2\left(\theta \right)co{s}^2\left(\theta \right)}$$

(46)

After adding the three equations of the near, mid, and far angles, the linear expression will be as follows:

$$Z_{S0}=0.33b_1{Z}_{S1}+0.33b_2{Z}_{S2}+0.33b_3{Z}_{S3}$$

(47)

The likelihood function of the $Z_P$ or $Z_S$ models can be written as shown below:

$$Im{p}_0|b_j,{\left(Im{p}_{\theta }\right)}_i,{\sigma }^2\sim N\left(\left[b_1{\left(Im{p}_1\right)}_i+b_2{\left(Im{p}_2\right)}_i+b_3{\left(Im{p}_3\right)}_i\right],{\sigma }^2\right)$$

(48)

Given the coefficients vector b, the explanatory variables vector $Im{p}_{\theta }$, and the variance ${\sigma }^2$, the zero-offset impedance $Im{p}_0$ follows a normal distribution with a mean equal to the linear expression of the covariates and coefficients, where (i) is the number of observations.

The next layer of the model is the coefficients vector $b_j$, which consists of $b_1$, $b_2$, and $b_3$ such that the three coefficients follow a common normal distribution. The mean and variance of this distribution were obtained from Equation (40) for the $Z_P$ model and Equation (46) for the $Z_S$ model. Based on the weak anisotropy assumption, the mean of the anisotropy parameters was considered zero and the standard deviation was 0.2. Moreover, the squared ratio of $Z_P$ and $Z_S$ was set, according to the wells’ data, to have a mean equal to 6 and a standard deviation (SD) equal to 3. The prior function for the $Z_P$ and $Z_S$ models can be written as follows:

$$bj\sim N(\mu ,\tau ),j=\{1,2,3\}$$

(49)

where $\mu$ and $\tau$ are, respectively, the mean and variance of the normal distribution of the coefficients vector b, which has an index (j) from 1 to 3.

The variance of the model was assumed to follow an inverse gamma distribution, which depends on the shape parameter $\alpha$ and the scale parameter $\beta$, which act as priors of the parameter ${\sigma }^2$. The distribution function of the ${\sigma }^2$ is shown below:

$${\sigma }^2\sim IG(\alpha ,\beta )$$

(50)

The idea of the model is to calculate the the zero-offset impedance by substituting the posterior means of the coefficients to Equation (38). The MCMC was used to draw multiple random realizations from the proposed normal distribution. Based on the law of large numbers, the average of a random sample from a distribution converges to the true mean of that distribution [99]. In other words, the Markov chains would converge to the posterior means of the model’s coefficients. The advantage of the MCMC is that it can solve for complicated posterior distributions which are difficult to be calculated mathematically.

Finally, the posterior means of the $Z_P$ and $Z_S$ models were substituted into Equation (38) to solve for the zero-offset impedances. The final graphical representation of the impedance model is shown in Figure 4, where $Im{p}_0$ is the zero-offset impedance, $Im{p}_1$, $Im{p}_2$, and $Im{p}_3$ are the impedances at the near, mid, and far angles, respectively, i is the observation index that ranges from 1 to n, bj is the common distribution of the coefficients which depends on the mean $\mu$ and variance $\tau$, and ${\sigma }^2$ is the variance of the model which depends on shape parameter $\alpha$ and scale parameter $\beta$.

4.3.2. Multilayer Feed-Forward Neural Network (MLFN)

Due to the availability of vast amounts of data in recent years, ANNs have become the most common prediction methods [100]. One of the widely used ANNs is the multilayer feed-forward neural network (MLFN) [101,102]. The MLFN consists of an input layer, one or more hidden layers, and an output layer [103]. An activation function is used to map the summation of the weighted inputs into the output layer [104]. The number of hidden layers and number of neurons in each of them should be carefully determined because this step controls the model’s accuracy [105].

The selection of appropriate activation functions has a high impact on a model’s performance. There are various types of activation functions such as the sigmoid, or logistic, function, which is commonly used in classification models. The current study considers a linear activation function to solve for the zero-offset impedance in a regression model. The choice of linear function is based on the moderate-to-strong correlation between the input and output variables as discussed earlier.

The connection between any two neurons is controlled by a random weight. A model’s random weights are adjusted until reaching the least mean square error and the best match between the predicted and actual output values. The backpropagation method is commonly used to obtain the optimal set of weights and calculate the error in a repetitive way until reaching the best training accuracy [106,107]. The backpropagation process involves an optimization algorithm such as stochastic gradient descent (SGD) [108] and adaptive moment estimation (Adam) [109]. The algorithm used in our model is the Levenberg–Marquardt (LMA) [110] due to its fast performance and ability to train MLFN models [111].

A model’s error can be represented by various loss functions such as cross-entropy [112], weighted binary cross-entropy (WCE) [113], and dice loss [114]. The loss function used in our study is the mean squared error (MSE), which is one of the best unbiased estimators used in regression models [115]. The MSE can be defined by the following equation:

$$\mathrm{MSE}=\frac{1}{N}\sum _{i=1}^N{(y_i-\widehat{y_i})}^2$$

(51)

where N is the number of samples and $y_i$ is the output variable at an observation (i).

The structure of the MLFN, in this study, consists of three layers: an input, a hidden, and an output layer. Figure 5 shows the graphical representation of the MLFN models for both $Z_P$ and $Z_S$, where the input layer consists of three neurons, which are the near-, mid-, and far-angle impedances, while the output layer consists of one neuron, which is the zero-offset impedance. Several trials were performed in MATLAB to determine the best number of neurons for the hidden layer. The best performance was observed at 20 neurons for both the $Z_P$ and $Z_S$ models.

4.4. Calculation of the Anisotropy Parameters: Epsilon and Delta

After comparing and validating $Z_P$ and $Z_S$, obtained from the isotropic inversion, statistical modeling, and MLFN, the best P- and S-impedance volumes were applied to Thomsen’s anisotropy equations to solve for the parameters Epsilon and Delta. Both the partially inverted impedances, zero-offset impedances, and central angles of the angle stacks (10, 20, and 32.5 degrees) were applied to the impedance-domain Thomsen’s equations in (30) and (31), yielding six equations, as shown below:

$$Z_{P1}=Z_{P0}\left[1+\delta si{n}^2\left(10\right)co{s}^2\left(10\right)+ϵsi{n}^4\left(10\right)\right]$$

(52)

$$Z_{P2}=Z_{P0}\left[1+\delta si{n}^2\left(20\right)co{s}^2\left(20\right)+ϵsi{n}^4\left(20\right)\right]$$

(53)

$$Z_{P3}=Z_{P0}\left[1+\delta si{n}^2\left(32.5\right)co{s}^2\left(32.5\right)+ϵsi{n}^4\left(32.5\right)\right]$$

(54)

$$Z_{S1}=Z_{S0}[1+{\left(\frac{Z_{P0}}{Z_{S0}}\right)}^2(ϵ-\delta )si{n}^2\left(10\right)co{s}^2\left(10\right)]$$

(55)

$$Z_{S2}=Z_{S0}[1+{\left(\frac{Z_{P0}}{Z_{S0}}\right)}^2(ϵ-\delta )si{n}^2\left(20\right)co{s}^2\left(20\right)]$$

(56)

$$Z_{S3}=Z_{S0}[1+{\left(\frac{Z_{P0}}{Z_{S0}}\right)}^2(ϵ-\delta )si{n}^2\left(32.5\right)co{s}^2\left(32.5\right)]$$

(57)

where $Z_{P1}$, $Z_{P2}$, and $Z_{P3}$ are the P-impedance volumes at the near, mid, and far angles, respectively, while $Z_{S1}$, $Z_{S2}$, and $Z_{S3}$ are the S-impedance volumes at the near, mid, and far angles, respectively. $Z_{P0}$ and $Z_{S0}$ are the zero-offset P- and S-impedances, respectively, and $ϵ$ and $\delta$ are the anisotropy parameters.

By adding Equations (52)–(54) to each other and Equations (55)–(57) to each other, two equations were obtained such that one of them related Epsilon and Delta to the $Z_P$ volumes and the other one related Epsilon and Delta to the $Z_S$ volumes. Finally, the two resulting equations were solved by MATLAB for the anisotropy parameters Epsilon and Delta, as shown below:

$$\mathrm{Epsilon}=2.3(Z_{P0}{Z}_{P1}+Z_{P0}{Z}_{P2}+Z_{P0}{Z}_{P3}+Z_{S0}{Z}_{S1}+Z_{S0}{Z}_{S2}+Z_{S0}{Z}_{S3}-3Z_{P0}^2-3Z_{S0}^2)/\left(Z_{P0}^2\right)$$

(58)

$$\mathrm{Delta}=2.3(Z_{P0}{Z}_{P1}+Z_{P0}{Z}_{P2}+Z_{P0}{Z}_{P3}-0.3Z_{S0}{Z}_{S1}-0.3Z_{S0}{Z}_{S2}-0.3Z_{S0}{Z}_{S3}-3Z_{P0}^2+0.9Z_{S0}^2)/\left(Z_{P0}^2\right)$$

(59)

4.5. Estimation of Facies Distribution

After obtaining the zero-offset P- and S-impedances, they were transformed into elastic properties in order to forecast the lithology and fluid probabilities separately by two logistic regression models. The lithology and fluid models were merged together to estimate the distribution of gas sand, wet sand, and shale by using the decision tree algorithm. $Z_P$ and $Z_S$ were firstly transformed into the lithology predictors: near and far elastic impedances [116]. Next, the fluid predictors, Mu–Rho, Lambda–Mu–Rho, and Poisson’s ratios, were calculated based on the LMR theory [117]. The idea of the facies model is based on the combination of the logistic regression [52] and decision tree algorithms.

Logistic regression is one of the linear regression models that deal with binary response variables. The idea of logistic regression is to provide model estimates which lie between zero and one [118]. These estimates can be considered the probabilities of occurrence of both zero and one. Therefore, the likelihood of a logistic regression model follows a Bernoulli distribution which takes on the value 1 with a probability $\varphi$ and 0 with a probability (1 – $\varphi$), as shown in the following equation:

$$Y_i|{\varphi }_i\sim Bernoulli\left({\varphi }_i\right)$$

(60)

where $Y_i$ is a binary event at an observation (i), while ${\varphi }_i$ is the probability of the success of this event at the observation (i). Instead of directly relating the expected value of $Y_i$ to the model’s variables, it can be assigned to the value of ${\varphi }_i$, which can be related to the model’s variables and coefficients by using a link function, as shown in the following equation:

$$logit\left({\varphi }_i\right)=log\left(\frac{{\varphi }_i}{1-{\varphi }_i}\right)={\beta }_0+{\beta }_1{\left(X_1\right)}_i+\dots {\beta }_n{\left(X_n\right)}_i$$

(61)

where $logit$(${\varphi }_i$) is the logarithmic function of the odds of ${\varphi }_i$, ${\beta }_0$ is the intercept, ${\beta }_1$ to ${\beta }_n$ are the model’s coefficients, and ${\left(X_1\right)}_i$ to ${\left(X_n\right)}_i$ are the model’s explanatory variables at an observation (i). After simplifying the previous equation in (61), the ${\varphi }_i$ can be obtained as shown below:

$$E\left(Y_i\right)={\varphi }_i=(1/1+e^{-({\beta }_0+{\beta }_1{\left(X_1\right)}_i+\dots {\beta }_n{\left(X_n\right)}_i)})$$

(62)

where $E\left(Y_i\right)$ is the expected value of $Y_i$. Assuming a three-class medium, which consists of gas sand, wet sand, and shale, the lithology model aims at predicting the sand probability from the lithology predictors, while the fluid model aims at estimating the gas probability from the fluid predictors. Given the observations of the elastic attributes and the facies log, the model’s coefficients can be obtained based on Bayes’s theorem [119].

The lithology model aims at predicting the sand probability from the near-EI and far-EI attributes. On the other hand, the fluid model aims at estimating the gas probability from the Mu–Rho, Lambda–Rho/Mu–Rho, and Poisson’s ratios. The posterior means of the coefficients of the two models were obtained by MCMC simulation. The facies and fluid models were merged together by using the sand and gas probabilities, along with depth values, as inputs to a decision tree, which turned the probabilities into facies estimates.

A decision tree is a classification method that generates questions over training samples and answers them using a set of statistical criteria for data classification [42]. There are two types of decision trees, which are regression and classification trees. As the response variable in this study was categorical, a decision tree was considered.

A decision tree begins with a root node where the samples are classified into two branches such that if the answer to the question is “Yes”, the sample goes under the “Yes” branch, while if the answer is “NO”, the sample goes under the “No” branch. A distinct series of questions and branches are set in the internal nodes until satisfying the stopping rule in the last terminal of the tree, which is called the leaf node. According to the classification rule [120], each leaf in a tree represents a class (A or B).

A popular method used to construct decision trees is called classification and regression trees (CART) [121]. The CART method splits observations into two parts by minimizing the relative sum of squared errors [122] between the two partitions, according to the principle of binary recursive partitioning (BRP) [123]. The advantages of this method are that it is nonparametric, free of significance tests, and has low sensitivity to outliers [124].

The process of tree shortening is called pruning. Avoiding the large size of a tree is crucial to eschew over-fitting. The size of a decision tree is controlled by many factors, such as the minimum number of samples that should be left in a node to achieve splitting [125]. A pruned tree mostly leads to results which are close to those of the original tree. However, the pruned tree may be more accurate than the initial one in some cases [126].

This study used the “rpart” package in R to model the three-class facies variable based on the CART method. The inputs of the tree were the sand and gas probabilities that had been estimated by logistic regression as well as the depth variable as trend factors. The output of the tree was a categorical variable having three values, 1, 2, and 3, which corresponded to the three litho-fluid classes: gas sand, wet sand, and shale, respectively. The tree was assessed by calculating the correct classification rate (CCR), which is the number of correctly classified observations over the total number of samples.

5. Results and Discussion

The wells were tied to seismic data by check shots of wells ×1 and ×2. A constant wavelet was extracted from the wells to achieve the tie between the synthetic logs and seismic traces. Next, three horizons were picked: A, B, and C, which appear as troughs in seismic data, as shown in Figure 6.

Both the conventional well logs and oil-mud resistivity imager (OMRI) were interpreted at the well (×1) such that the zone of interest extended from 1300 m to 1480 m. Firstly, the sand and shale could be preliminarily discriminated by using the shale volume at the cut-off 0.3. The gas sand and wet sand were differentiated by $V_P/V_S$ and water saturation at the cut-offs 2.8 and 0.45, respectively. The thicknesses of the gas sand zones were validated by the net pay thickness mentioned in the well report.

Figure 7 shows the OMRI-derived facies against the well-log-derived facies obtained from the petrophysical cut-offs at the well (×1). Both facies logs confirm that the dominant facies is a massive shale that corresponds to the high shale volume and water saturation zones. The massive and planar-laminated sandstone layers correspond to the gas intervals of the two reservoirs: A and B. Those gas zones appear as fining-upward cycles in the shale volume (Vsh) log where $V_P/V_S$ and water saturation are low. Minor siltstone zones appear in track (a) and coincide with some of the wet sand zones in track (b). The layers become thinner as the depth increases, coinciding with abrupt changes in the well logs.

The top of the target zone (1300–1316 m) consists of massive shale, which can be interpreted as a delta plain. The next interval extends from depth 1316 m to 1330 m, where the A reservoir appears as a distributary channel that consists of fining-upward sediments. The top of this interval is considered a fine-grained sandstone, which appears as wet sand in track (b). The rest of the A reservoir consists of equal amounts of the planar and cross-bedded gas-bearing sandstones, where $V_P/V_S$, Vsh, and water saturation are low. A sharp contact is clear at the base of the channel at depth 1330 m due to the change from sandstone to shale.

Another depositional cycle begins at depth 1330 m, where the massive and laminated shales are interbedded with minor siltstone. The OMRI images show broad and spiral bands in darker clay-rich portions. The shale extends until depth 1422 m, where another distributary channel appears within the B reservoir. This unit is dominated by massive-to-weakly stratified porous sands. A sharp contact is defined at depth 1438 m, below which the lower interval consists of laminated sandstone with low porosity. The depth interval from 1452 m to 1481 m consists of interbedded sandstone, siltstone, and shale. The reservoir quality in this zone is poor, with high water-saturation and shale-volume. This is matching with the facies log in track (b), where the sandstone is wet.

The OMRI images are crucial for data filtering due to their high resolution relative to the conventional well logs. For example, there is a possible streak of coal at depth 1316.8 m which could not be detected by the conventional logs, as shown in Figure 7. Accordingly, all thin layers and coal streaks have been excluded from the facies model. This data preparation step is essential to reduce misclassification. In other words, if the OMRI images are not available, the collected training data cannot be accurate due to the presence of coal streaks and thin layers which are below logging data resolution.

A group wavelet was extracted from the angle gathers such that it consisted of three wavelets at the three central angles: 10°, 20°, and 32.5°. The initial model of the simultaneous inversion was created from wells ×1 and ×2 such that the target zone ranged from two-way time 1000 to 1800 ms with a two-millisecond (ms) sample rate. A smoother was applied to the modeled trace in the output domain after lateral interpolation. The inverted properties of the inversion were the P-impedance, S-impedance, and density. Similarly, an initial model was created for each angle stack to invert for P- and S-impedances at the near, mid, and far incident angles.

Figure 8 and Figure 9 show the change in the P- and S-impedance values by changing the incident angle for each facies class at the well (×1). The anisotropy effect appears in the shale-dominated zones, as shown in Figure 8d. The P-impedance seems directly proportional to the incident angle, especially in the black circles, where the values gradually increase from the near to the far impedance sections. On the other hand, the gas sands of the A and B horizons appear in green color and almost have the same impedance values at the near, mid, and far angles.

These results seem consistent with a previous study that calculates the P-wave velocity for the Floyd shale by three different models, as shown in Figure 10 [127]. The first and second methods, T_45 and T_60, are based on Thomsen’s model [53], where the Delta parameter is obtained using the $C_{13}$-component on the 45${}^{\circ }$ and 60${}^{\circ }$ plugs, respectively. The third model, $Vp\_b$, is based on the strong anisotropy assumption [128]. The three models show that the P-wave velocity increases by increasing the incident angle due to the intrinsic anisotropy of shale. The results are confirmed by the ultrasonic measurements applied to the core sample.

A previous study shows the measurements of the P- and S-wave velocities in shale core samples which indicate moderate intrinsic anisotropy, where the Epsilon ($ϵ$) is 0.16, Delta ($\delta$) is 0.08, and Gamma ($\gamma$) is 0.29 [129]. Figure 11 shows a plot of velocity versus angle, where the quasivertical P-wave velocity (qP), quasivertical S-wave velocity (qSV), and horizontal S-wave velocity (SH) increase by increasing the incident angle at the angle range up to approximately 35°. These results are consistent with the current study that claims that the P- and S-impedances increase by increasing angles, as shown in the black circles of Figure 8 and Figure 9.

An important note about seismic data is that the velocity change with angle may be due to other reasons rather than seismic anisotropy, such as tuning effect and noise. Accordingly, the next step is to estimate the zero-offset impedance from the partially inverted volumes to normalize the random error and anisotropy effect simultaneously.

The total data points collected for the $Z_{P0}$ and $Z_{S0}$ models are 1292 and 1582, respectively. These data points were gathered from the Backus-averaged impedance logs of wells ×1 and ×2 and the depth-converted composite traces of the partially inverted impedances at the wells’ locations. The data points were randomly divided into three parts such that 70% of them were used for training, 15% for testing, and 15% for validation.

Figure 12 and Figure 13 show the trace plots of the Markov chains, which aim at estimating the coefficients for both the $Z_P$ and $Z_S$ models, respectively. The number of iterations is shown on the X-axis and the value of each coefficient is shown on the Y-axis. Each coefficient was simulated in R by three chains, which appear in green, black, and red colors. For both the $Z_P$ and $Z_S$ models, the chains reached their stationary distribution in less than 5000 iterations. However, the chains of the $Z_P$ model seem more converged than those of the $Z_S$ model due to the high random error of the $Z_S$ inversion.

The posterior means and standard deviations of the coefficients were finally obtained for each model, as shown in

Table 1. The statistical properties, mean (Mu) and standard deviation (SD), of the coefficients of the zero-offset $Z_P$ and $Z_S$ statistical models.

Coefficient	${\mathit{b}}_1$		${\mathit{b}}_2$		${\mathit{b}}_3$
Property	Mu	SD	Mu	SD	Mu	SD
$Z_P$ Model	0.3532	0.01	0.3555	0.01	0.286	0.01
$Z_S$ Model	0.3672	0.01	0.364	0.01	0.266	0.01

. It is clear from the posterior means of the coefficients that the rate of change of the zero-offset impedance with the angle of incidence is almost constant at the near and mid incident angles (10${}^{\circ }$ and 20°), while the rate of change decreases at the far angle (32.5°). In other words, the zero-offset $Z_P$ changes with the far $Z_P$ at a lower rate compared to the near- and mid-angle impedances.

After applying the posterior means of the coefficients to Equation (38), the $Z_P$ and $Z_S$ equations will be:

$$Z_{P0}=0.3532Z_{P1}+0.3555Z_{P2}+0.286Z_{P3}$$

(63)

$$Z_{s}0=0.3672Z_{S1}+0.364Z_{S2}+0.266Z_{S3}$$

(64)

The next step is to estimate the zero-offset impedances by the MLFN. The training, testing, and validation data were the same as those used for the statistical model. Two networks were trained for the $Z_P$ and $Z_S$ models by the Levenberg–Marquardt algorithm such that the best performance was reached at 21 iterations for the $Z_P$ model and 58 for the $Z_S$ model, as shown in Figure 14, where the number of iterations is on the X-axis and the mean square error is on the Y-axis. The error reduction trends of the training, validation, and testing data appear in blue, green, and red colors, respectively. The three trends are close to each other, which means that there is no over-fitting.

The networks reached the least mean square error (MSE) error values at 92,840 for the $Z_P$ model and 36,893 for the $Z_S$ model. In addition, the error values were plotted in a histogram, as shown in Figure 15, where the errors of the training, validation, and testing data appear in blue, green, and red colors, respectively. The error distribution for the $Z_P$ and $Z_S$ models is centered around the zero-error (yellow) line.

The quality check of the resulting zero-offset $Z_P$ and $Z_S$ starts with color matching with the impedance logs, as shown in Figure 16 and Figure 17, respectively. It is clear that the MLFN’s impedances (Figure 16c and Figure 17c) better match the impedance logs at the wells (×1) and (×2) than those of simultaneous inversion (Figure 16a and Figure 17a) and statistical modeling (Figure 16b and Figure 17b). In addition, there are extreme impedance values shown in isotropic inversion compared to the anisotropic methods that have more normalized inversion error values.

Another visual comparison is made by the histogram plots of the inverted $Z_P$ for the three models, as shown in Figure 18, where the $Z_P$, obtained from the simultaneous inversion, has extreme low and high values of impedance with relatively high standard deviation (1445) and low mean (6195). The extreme values of simultaneous inversion appear in purple and green colors in Figure 16a. On the other hand, the least standard deviation (1068) belongs to the $Z_P$ of the MLFN. Moreover, the means of the anisotropic models are higher than those of simultaneous inversion, which show extremely low $Z_P$ values down to 2000 ((m/s)·(g/cc)).

Similarly, Figure 19 shows the histogram plots extracted from the the S-impedance sections, where the isotropic method has extremely low $Z_S$ values, down to 1000 ((m/s)·(g/cc)), high values, up to 9000 ((m/s)·(g/cc)), and relatively high standard deviation (1671). On the other hand, the statistical model and MLFN have lower standard deviations, which are 821 and 886, respectively.

Another comparison is to evaluate the residuals’ plots of all models by calculating the difference between the impedance logs and the predicted impedance for each model. Figure 20 shows the plots of the calculated residuals versus the observation indexes for the $Z_P$ models (Figure 20a–c) and $Z_S$ models (Figure 20d–f). It can be noticed that the residuals of the MLFN model are centered around zero and have low magnitude compared to those of the simultaneous inversion and statistical modeling. In order to determine the best $Z_P$ model, the R-value was calculated for each model, as shown in the

Table 2. The correlation coefficients of the predicted $Z_P$ and $Z_S$ obtained from isotropic inversion, statistical modeling, and MLFN at the wells (×1) and (×2).

Property	Isotropic Inversion	Statistical Modeling	MLFN
$Z_P$	84%	88%	94%
$Z_S$	73%	85%	92%

.

The simultaneous (isotropic) inversion shows low R-values for $Z_P$ and $Z_S$ (84% and 73%, respectively), which indicates high inversion errors. The statistical model has enhanced the correlations to be 88% and 85% for $Z_P$ and $Z_S$, respectively. However, the statistical models’ results are still close to the isotropic case because the random error has been reduced by using only one value of each coefficient for the entire volume. Therefore, the error normalization is not as much as desired. On the other hand, the MLFN method has the highest R-values, 94% and 92% for $Z_P$ and $Z_S$, respectively. This shows how ML can reduce the randomness in the data and turn it into recognized patterns.

It can be concluded that the statistical model has resulted in more accurate impedance volumes than the simultaneous inversion; however, the statistical model has only one value for each coefficient, which is not enough to normalize the random error as desired. The model can be more hierarchical if there are more facies data in several wells. On the other hand, the MLFN model could estimate the zero-offset impedance efficiently by turning the randomness in data into identified patterns.

A previous method [90] aimed at calculating the zero-offset velocities and anisotropy parameters empirically by inverting the partial stacks into velocities at different incident angles and solving Thomsen’s anisotropy equations [61]. However, this method is applicable to high-resolution seismic data. On the other hand, the current study uses statistical modeling and neural networks (MLFN) to reduce the anisotropy effect and normalize the error simultaneously.

Ma estimated the P- and S-impedances from the P- and S-wave reflectivities by using a joint inversion [130]. The author reduced the nonuniqueness of the inversion output by adding the $V_S$/$V_P$ ratio as a lithology constraint. Compared to simultaneous inversion, joint inversion follows the background lithology trend and hence yields more accurate impedances. The current study also follows a lithology constraint by estimating the zero-offset impedances from the impedances at different incident angles. In other words, the impedances at non-normal angles add facies information to the model and yield a unique solution for the zero-offset impedance.

The $Z_P$ and $Z_S$ volumes, obtained from the MLFN method, were applied to Equations (58) and (59) in order to solve for the anisotropy parameters Epsilon and Delta. Figure 21 shows the resulting Epsilon and Delta sections plotted against those obtained from the methodology of Liner and Fei [79] at the well (×1).

Figure 22 shows the histograms of Epsilon and Delta sections, where the values of Epsilon range from −0.2 to 0.2 with a mean value of 0.002, while the Delta ranges from −0.075 to +0.075 with a mean value of 0.001. These values are within the weak-anisotropy ranges, according to Thomsen’s model [53].

The composite traces of the resulting Epsilon and Delta were extracted from their volumes and plotted against the facies log at the well (×1), as shown in Figure 23. The tops of the gas horizons, A and B, show high positive Epsilon values (blue rectangles) due to the reduction in the P-wave velocity. The positive sign is reasonable because the vertical velocity decreases relative to the horizontal velocity, which is in the same direction as the bedding.

Another observation is that the anisotropy degree is directly proportional to the contrast of the elastic properties between layers. This phenomenon appears in the abrupt change in facies, as shown in Figure 23, where the black rectangle shows high contrast boundaries between layers coinciding with significant fluctuations in the Epsilon parameter. The same conclusion is mentioned by Bandyopadhyay [76], who plotted the anisotropy parameters, Epsilon and Delta, versus the depth for a laminated shaly sand succession, as shown in Figure 24. The magnitude of the anisotropy is low for the water-saturated sand (Figure 24a), which has Lame’s parameters that are similar to those of shale. On the other hand, the Lame’s parameters of the gas-saturated sand (Figure 24b) and shale are significantly different, resulting in a high magnitude of the anisotropy parameters.

Most of the previous studies stated that Epsilon is always positive or zero. However, the resulting Epsilon, as shown in Figure 21a, has a negative sign in some zones. This is matching with Bandyopadhyay’s study [76], which stated that Epsilon can be negative in laminated shaly gas sand layers, as shown in the anisotropy–depth plot (Figure 24b), where the depth interval is divided into three zones, (A), (B), and (C), which are the shallow, compaction, and diagenesis regimes, respectively. Epsilon is significantly fluctuating around zero in the gas-saturated sand such that it becomes negative at the shallow and diagenesis zones.

The Epsilon and Delta sections are plotted in Figure 25, showing a strong positive correlation between them. The color legend represents the two-way time value for each observation. The best-fit line gives the following linear relationship between Epsilon and Delta:

$$\delta =0.31ϵ+0.004$$

(65)

Interestingly, the observed relationship is mostly similar to that obtained by Li [73], who implemented the following equation:

$$\delta =0.32ϵ$$

(66)

The next step is to transform the $Z_P$ and $Z_S$ into the elastic properties from which the lithology and fluid distributions should be forecasted. The facies integrated model was firstly applied to logging data at the well (×1) for training and then used to predict the facies distribution from the elastic volumes. The lithology model was postulated by using R, based on logistic regression, to forecast the sand probability from the near and far elastic impedances. Similarly, the fluid model was created to predict the gas probability from the Mu–Rho, Lambda–Rho/Mu–Rho, and Poisson’s ratios [52]. Finally, the two models were merged by a decision tree to obtain the distribution of each facies class.

The predicted sand and gas probabilities (Phi-sand and Phi-gas-sand, respectively) of the training data were plotted versus the lithology and fluid variables, respectively, as shown in Figure 26. As expected, the sand probability values are high where the lithology variable is closer to 1, indicating sandstone, and low where the lithology variable is closer to 0, indicating shale (Figure 26a). Like the sand probability, the gas probability is high where the fluid variable is 1, indicating gas sand, and generally low where the fluid variable equals 0, indicating wet sand (Figure 26b). The best cut-off values were selected to differentiate between the sand and shale observations in the lithology model as well as the gas sand and wet sand in the fluid model.

A decision tree was created to estimate the litho-fluid classes from the sand probability, gas probability, and depth. The total data points are 1084, starting from depth 1300 m to 1851 m. According to the facies log of the well (×1), the data consist of 75 gas sand, 286 wet sand, and 723 shale data points. The data were divided into three parts for training (80%), testing, (10%), and validation (10%). The best minimum split value has been observed to be 10, which leads to the best performance of the tree. Figure 27 shows the graphical representation of the decision tree, which classifies the training data into gas sand in red, wet sand in blue, and shale in green. The input variables are the gas probability (Pg), sand probability (Ps), and depth (d).

The accuracy of the tree was calculated for the training, testing, and validation data to be 91%, 89%, and 90%, respectively, and the average accuracy is 91%. These results show that the tree is accurate enough to be applied to the depth-converted seismic data. A function was created in R-Studio to transform the zero-offset $Z_P$ and $Z_S$ volumes, obtained from the isotropic inversion and MLFN, into elastic properties and then into sand and gas probabilities based on the facies model.

Eventually, the decision tree was applied to predict the facies classes from the probabilities and depth. The composite traces of the facies volumes were extracted at the well (×1) to compare the isotropic simultaneous inversion and the anisotropic MLFN outputs, as shown in Figure 28, where the gas sand, wet sand, and shale are represented in the red, blue, and green colors, respectively. The CCR was calculated for each litho-fluid class for both the isotropic and anisotropic cases, as shown in

Table 3. The correct classification rates of the litho-fluid classes for isotropic inversion and MLFN.

Method	Gas Sand	Wet Sand	Shale	All Facies
Isotropic Inversion	51%	82%	55%	56%
MLFN	97%	83%	73%	77%

.

The CCR values of the gas sand are 51% and 97% for the isotropic and anisotropic cases, respectively. This improvement is evident in Figure 28, where the gas sand is correctly detected at the anisotropic case. Similarly, the CCR values of the shale are 55% and 73% for the isotropic and anisotropic cases, respectively. The CCR values of the wet sand are almost the same for the two cases, 82% and 83% for the isotropic and anisotropic cases, respectively.

The average CCR value is 56% for the isotropic facies and 77% for the anisotropic case. This means that the accuracy of litho-fluid detection has been enhanced by about 21%. This improvement is not only due to the reduction in the anisotropy effect but also because the MLFN has normalized the inversion errors due to the tuning effect.

These results seem comparable to those of another case study, from the Okoli field in the Sava depression, which aimed at predicting clastic facies using neural networks [107]. Although the lithology was estimated based on various logs, the results showed correlation levels ranging from 78.3% to 82.1%, which are close to the results of this study.

A recent study set various ML algorithms to estimate four facies classes from seismic attributes [131]. The logistic regression showed encouraging success rates, up to 0.98, compared to other algorithms, such as the Gaussian process, linear discriminant analysis, and stochastic gradient descent. However, a high success rate does not indicate the model’s stability, especially in heterogeneous fields. The advantage of the current study’s facies model is that the litho-fluid probabilities as well as the depth variable are inserted into a decision tree which turns the probabilities into facies estimates. The integration of two different algorithms guarantees the model’s stability and reduces the chance of over-fitting.

6. Conclusions

The zero-offset $Z_P$ and $Z_S$ volumes were obtained by an isotropic simultaneous inversion and then by two anisotropic methods: statistical modeling and MLFN. The validation of the P-impedance volumes shows that the MLFN model had near-zero residuals and a high R-value (0.94) compared to simultaneous inversion and statistical modeling, which had low R-squared values (0.84 and 0.88, respectively). Similarly, the MLFN resulted in a high correlation for the S-impedance (0.92), while simultaneous inversion and statistical modeling had lesser R-values (0.73 and 0.85, respectively).

The $Z_P$ and $Z_S$ volumes, obtained from the MLFN models, were applied to the impedance-domain Thomsen’s anisotropy equations to solve for the anisotropy parameters: Epsilon and Delta. The anisotropy magnitude was observed to be low in the shaly wet sand successions, where the contrast of the elastic properties between layers was low, while it showed higher magnitude in the shaly gas sand successions, which had high-contrast fluctuations in the elastic properties.

The $Z_P$ and $Z_S$ volumes, obtained from the isotropic inversion and MLFN, were converted into elastic properties from which the sand and gas probabilities were forecasted by logistic regression. The predicted probabilities were inserted into a decision tree to estimate the distribution of gas sand, wet sand, and shale classes within the target zone. The predicted facies distributions were validated by the facies log at the well (×1), where the MLFN showed a high average success rate (0.77) compared to the simultaneous inversion, which had a low success rate (0.56).

Conventional inversion algorithms have proved to be misleading if seismic anisotropy is neglected. On the other hand, statistical methods and ML are more efficient tools for obtaining more accurate rock properties, especially when considering the dependence of these properties on incident angles.

The decision tree of the facies model shows encouraging results due to the contribution of the depth constraint. This proves that the trend factor is crucial for accurate reservoir characterization. Therefore, the availability of geochemical data would be an added value that links other properties, such as elastic properties and anisotropy parameters, to the pressure regime and compaction profile.

The study considers that the subsurface follows the VTI assumption, which applies to horizontal layers. However, the methodology can be extended to include more realistic cases, such as the TTI and HTI assumptions, which can enhance the accuracy of zero-offset parameter modeling. In addition, core data are needed to validate anisotropy measurement and yield confident models.