**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 *VP*/*VS* 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 *VP*/*VS* 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.

**Figure 6.** A seismic section extracted from the near-stack volume and tied to the wells ×1, ×2, and ×3, where the horizons' tops A, B, and C are in yellow, green, and cyan colors, respectively.

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 *VP*/*VS*, 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.

**Figure 7.** Well (×) facies and conventional logs: (**a**) OMRI-derived facies, (**b**) well-log-derived facies, (**c**) well logs, and (**d**) well tops.

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-toweakly 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.

Figures 8 and 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.

**Figure 8.** The partially inverted *ZP* sections at the near, mid, and far incident angles (from (**a**–**c**), respectively) against the facies log at well ×1 (**d**).

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◦ and 60◦ plugs, respectively. The third model, *V p*\_*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.

**Figure 9.** The partially inverted *ZS* sections at the near, mid, and far incident angles (from (**a**–**c**), respectively) against the facies log at well ×1 (**d**).

**Figure 10.** A plot of the P-wave velocities vs. the incident angles for the Floyd shale.

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 (*δ*) is 0.08, and 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 Figures 8 and 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.

**Figure 11.** A plot of the velocity versus incident angle, obtained from ultrasonic measurements on shale cores.

The total data points collected for the *ZP*<sup>0</sup> and *ZS*<sup>0</sup> 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.

Figures 12 and 13 show the trace plots of the Markov chains, which aim at estimating the coefficients for both the *ZP* and *ZS* 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 *ZP* and *ZS* models, the chains reached their stationary distribution in less than 5000 iterations. However, the chains of the *ZP* model seem more converged than those of the *ZS* model due to the high random error of the *ZS* inversion.

**Figure 12.** The trace plots of the Markov chains that simulate the coefficients of the variables for the *ZP* statistical model: b1, b2, and b3.

The posterior means and standard deviations of the coefficients were finally obtained for each model, as shown in Table 1. 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◦ and 20°), while the rate of change decreases at the far angle (32.5°). In other words, the zero-offset *ZP* changes with the far *ZP* at a lower rate compared to the near- and mid-angle impedances.

**Figure 13.** The trace plots of the Markov chains that simulate the coefficients of the variables for the *ZS* statistical model: b1, b2, and b3.

**Table 1.** The statistical properties, mean (Mu) and standard deviation (SD), of the coefficients of the zero-offset *ZP* and *ZS* statistical models.


After applying the posterior means of the coefficients to Equation (38), the *ZP* and *ZS* equations will be:

$$Z\_{P0} = 0.3532Z\_{P1} + 0.3555Z\_{P2} + 0.286Z\_{P3} \tag{63}$$

$$Z\_s0 = 0.3672Z\_{S1} + 0.364Z\_{S2} + 0.266Z\_{S3} \tag{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 *ZP* and *ZS* models by the Levenberg–Marquardt algorithm such that the best performance was reached at 21 iterations for the *ZP* model and 58 for the *ZS* 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 *ZP* model and 36,893 for the *ZS* 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 *ZP* and *ZS* models is centered around the zero-error (yellow) line.

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

**Figure 14.** The performance plots of the (**a**) *ZP* and (**b**) *ZS* MLFN models, 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.

**Figure 15.** The error histograms of the (**a**) *ZP* and (**b**) *ZS* MLFN models, where the error value is on the X-axis and the error frequency (instances) is on the Y-axis. The errors of the training, validation, and testing data appear in blue, green, and red colors, respectively.

Another visual comparison is made by the histogram plots of the inverted *ZP* for the three models, as shown in Figure 18, where the *ZP*, 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 *ZP* of the MLFN. Moreover, the means of the anisotropic models are higher than those of simultaneous inversion, which show extremely low *ZP* 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 *ZS* 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 *ZP* models (Figure 20a–c) and *ZS* 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 *ZP* model, the R-value was calculated for each model, as shown in the Table 2.

(**c**)

**Figure 16.** Arbitrary sections extracted from the zero-offset acoustic impedance volumes obtained from (**a**) isotropic inversion, (**b**) statistical modeling, and (**c**) MLFN. The inserted curves are the acoustic impedance logs of the wells ×1 and ×2.

(**c**)

**Figure 17.** Arbitrary sections extracted from the zero-offset shear impedance volumes obtained from (**a**) isotropic inversion, (**b**) statistical modeling, and (**c**) MLFN. The inserted curves are the shear impedance logs of the wells ×1 and ×2.

**Figure 18.** The histogram plots of the zero-offset P-impedance obtained from (**a**) isotropic inversion, (**b**) statistical modeling, and (**c**) MLFN.

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


The simultaneous (isotropic) inversion shows low R-values for *ZP* and *ZS* (84% and 73%, respectively), which indicates high inversion errors. The statistical model has enhanced the correlations to be 88% and 85% for *ZP* and *ZS*, 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 *ZP* and *ZS*, 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.

**Figure 19.** The histogram plots of the zero-offset S-impedance obtained from (**a**) isotropic inversion, (**b**) statistical modeling, and (**c**) MLFN.

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 *VS*/*VP* 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 *ZP* and *ZS* 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 21.** Cross-sections extracted from the anisotropy parameters' volumes Epsilon (**a**) and Delta (**b**) calculated by applying the MLFN-based *ZP* and *ZS* to the impedance-domain Thomsen's anisotropy equations. The inserted curves are the Epsilon and Delta calculated 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].

**Figure 22.** The histograms of the anisotropy parameters (**a**) Epsilon and (**b**) Delta sections, which were calculated by applying the MLFN-based *ZP* and *ZS* to the impedance-domain Thomsen's anisotropy equations.

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.

**Figure 23.** The anisotropy parameters, Epsilon and Delta, against the facies log at the well (×1), where A, B, C, and D are sand zones.

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.

**Figure 24.** Plots of the anisotropy parameters, Epsilon (blue curve) and Delta (red curve), on the Xaxes versus the depth on the Y-axes, for both the water- and gas-saturated sand (**a** and **b**, respectively). The depth range is divided into three zones, (A), (B), and (C), which are the shallow, compaction, and diagenesis regimes, respectively [76].

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\epsilon + 0.004\tag{65}
$$

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

$$
\delta = 0.32 \epsilon \tag{66}
$$

The next step is to transform the *ZP* and *ZS* 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.

**Figure 25.** A linear regression plot between Epsilon, on the X-axis, and Delta, on the Y-axis. The data points were obtained from the Epsilon and Delta volumes, which were calculated by applying the MLFN-based *ZP* and *ZS* to the impedance-domain Thomsen's anisotropy equations.

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).

**Figure 26.** (**a**) Plot of the sand probability (Phi-sand) versus the lithology (facies) variable. (**b**) Plot of the gas probability (Phi-gas-sand) versus the fluid variable.

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 *ZP* and *ZS* volumes, obtained from the isotropic inversion and MLFN, into elastic properties and then into sand and gas probabilities based on the facies model.

**Figure 27.** The graphical representation of the decision tree model, which classifies the training data of the well (×1) 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).

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.

**Figure 28.** The litho-fluid facies at the well (×1): (**a**) log facies, (**b**) simultaneous inversion facies, and (**c**) MLFN facies.


Inversion 51% 82% 55% 56% MLFN 97% 83% 73% 77%


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.
