**2. Seismic Anisotropy**

Seismic anisotropy is the dependence of the seismic velocity on the incident angle [53–55]. Wave propagation in anisotropic media has been discussed by various studies [53,56–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} \tag{1}
$$

where *σij* is the stress tensor, *ekl* is the strain tensor, and *cijkl* 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} \tag{2}$$

The 81 independent components are accordingly reduced to 36 due to the symmetry in the stress and strain tensors. Therefore, the elastic tensor (*cijlk*) is replaced with a 6 × 6 matrix (*cαβ*), where [61]:

$$
\alpha = i \dot{\jmath} = \ddot{\jmath}, \beta = kl = lk \tag{3}
$$

where *α* and *β* are the new indices of the elastic tensor. The resulting (6 × 6) notation matrix is expressed as shown below:

$$
\begin{bmatrix}
\text{C}\_{11} & \text{C}\_{12} & \text{C}\_{13} & \text{C}\_{14} & \text{C}\_{15} & \text{C}\_{16} \\
\text{C}\_{12} & \text{C}\_{22} & \text{C}\_{23} & \text{C}\_{24} & \text{C}\_{25} & \text{C}\_{26} \\
\text{C}\_{13} & \text{C}\_{23} & \text{C}\_{33} & \text{C}\_{34} & \text{C}\_{35} & \text{C}\_{36} \\
\text{C}\_{14} & \text{C}\_{24} & \text{C}\_{34} & \text{C}\_{44} & \text{C}\_{45} & \text{C}\_{46} \\
\text{C}\_{15} & \text{C}\_{25} & \text{C}\_{35} & \text{C}\_{45} & \text{C}\_{55} & \text{C}\_{56} \\
\text{C}\_{16} & \text{C}\_{26} & \text{C}\_{36} & \text{C}\_{46} & \text{C}\_{56} & \text{C}\_{66}
\end{bmatrix} \tag{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αβ*) 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:

$$
\begin{bmatrix}
\mathcal{C}\_{11} & \mathcal{C}\_{11} - 2\mathcal{C}\_{66} & \mathcal{C}\_{13} & 0 & 0 & 0 \\
\mathcal{C}\_{11} - 2\mathcal{C}\_{66} & \mathcal{C}\_{11} & \mathcal{C}\_{13} & 0 & 0 & 0 \\
\mathcal{C}\_{13} & \mathcal{C}\_{13} & \mathcal{C}\_{33} & 0 & 0 & 0 \\
0 & 0 & 0 & \mathcal{C}\_{44} & 0 & 0 \\
0 & 0 & 0 & 0 & \mathcal{C}\_{44} & 0 \\
0 & 0 & 0 & 0 & 0 & \mathcal{C}\_{66}
\end{bmatrix}
\tag{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 (*δ*), and 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:

$$\epsilon = \frac{V\_{PH} - V\_{PV}}{V\_{PV}} = \frac{\mathcal{C}\_{11} - \mathcal{C}\_{33}}{2\mathcal{C}\_{33}}\tag{6}$$

where is the anisotropy parameter Epsilon, *VPH* is the horizontal P-wave velocity, *VPV* is the vertical P-wave velocity, and *C*<sup>11</sup> and *C*<sup>33</sup> 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{(\mathbb{C}\_{13} + \mathbb{C}\_{44})^2 - (\mathbb{C}\_{33} - \mathbb{C}\_{44})^2}{2\mathbb{C}\_{33}(\mathbb{C}\_{33} - \mathbb{C}\_{44})} \tag{7}$$

where *δ* is the anisotropy parameter Delta, while *C*13, *C*44, and *C*<sup>33</sup> 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{\mathbb{C}\_{66} - \mathbb{C}\_{44}}{2\mathbb{C}\_{44}}\tag{8}$$

where *γ* is the anisotropy parameter Gamma, *VSH* is the horizontal S-wave velocity, *VSV* is the vertical S-wave velocity, and *C*<sup>66</sup> and *C*<sup>44</sup> 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(\theta) = V\_{P0} \left[ 1 + \delta \sin^2(\theta) \cos^2(\theta) + \epsilon \sin^4(\theta) \right] \tag{9}$$

$$V\_{SV}(\theta) = V\_{S0} \left[ 1 + \left(\frac{V\_{P0}}{V\_{S0}}\right)^2 (\epsilon - \delta) \sin^2(\theta) \cos^2(\theta) \right] \tag{10}$$

$$V\_{SH}(\theta) = V\_{\text{S0}}(1 + \gamma \sin^2(\theta))\tag{11}$$

where *VP*(*θ*) and *VSV*(*θ*) are the vertical P-wave and S-wave velocities, respectively, *VSH*(*θ*) is the horizontal S-wave velocity, *VP*<sup>0</sup> and *VS*<sup>0</sup> are the zero-offset P-wave and S-wave velocities, respectively, at *θ* = 0, and *θ* 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–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].

**Figure 1.** A plot of the reflection coefficient on the Y-axis versus the angle of incidence on the X-axis. The blue and green lines represent the reflectivities obtained based on the isotropic and anisotropic assumptions, respectively.

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–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 compactiondependent 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◦, 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 *ZP* and *ZS* 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 *ZP*, *ZS*, and density volumes and inverting the angle stacks into *ZP* and *ZS* 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 *ZP* and *ZS* 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.

**Figure 2.** A flow diagram showing the methodology of estimating the facies distribution from the zero-offset acoustic and shear impedances, which are obtained from isotropic inversion, statistical modeling, and MLFN.

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.
