Application of Unconventional Seismic Attributes and Unsupervised Machine Learning for the Identification of Fault and Fracture Network

: The identification of small scale faults (SSFs) and fractures provides an improved understanding of geologic structural features and can be exploited for future drilling prospects. Conventional SSF and fracture characterization are challenging and time ‐ consuming. Thus, the current study was conducted with the following aims: (a) to provide an effective way of utilizing the seismic data in the absence of image logs and cores for characterizing SSFs and fractures; (b) to present an unconventional way of data conditioning using geostatistical and structural filtering; (c) to provide an advanced workflow through multi ‐ attributes, neural networks, and ant ‐ colony optimization (ACO) for the recognition of fracture networks; and (d) to identify the fault and fracture orientation parameters within the study area. Initially, a steering cube was generated, and a dip ‐ steered median filter (DSMF), a dip ‐ steered diffusion filter (DSDF), and a fault enhancement filter (FEF) were applied to sharpen the discontinuities. Multiple structural attributes were applied and shortlisted, including dip and curvature attributes, filtered and unfiltered similarity attributes, thinned fault likelihood (TFL), fracture density, and fracture proximity. These shortlisted attributes were computed through unsupervised vector quantization (UVQ) neural networks. The results of the UVQ revealed the orientations, locations, and extensions of fractures in the study area. The ACO proved helpful in identifying the fracture parameters such as fracture length, dip angle, azimuth, and surface area. The adopted workflow also revealed a small scale fault which had an NNW–SSE orientation with minor heave and throw. The implemented workflow of structural interpretation is helpful for the field development of the study area and can be applied worldwide in carbonate, sand, coal, and shale gas fields.


Introduction
Small scale faults (SSFs) are responsible for a minor portion of all seismic activity in an active region and have an offset of about 25-50 m [1] that can extend to about 2-4 km with small throws of about 40 m [2]. The study on the assessment of the small scale faults and fractures in heterogeneous sandstone reservoirs is a global issue [3]. However, it is challenging to accurately categorize fractures due to their heterogeneous nature and complicated features [4,5]. SSFs and fractures control the storage and impact distribution of the natural gas storage in reservoirs [6] and provide necessary information related to tectonics, overpressure, burial history, and diagenesis [7]. SSFs and fractures tend to make gas reservoirs very complicated, particularly in those fields that show a heterogeneity and poor continuity of lateral facies distribution [8].Research of SSFs and fractures plays a vital role in establishing the porosity and permeability of these heterogeneous and tight sandstone reservoirs. The evaluation of SSFs and fractures provides controls and guidance for the prospect evaluation and development of heterogeneous sandstone reservoirs [9]. Interpreters usually manually mark faults on seismic sections to carry out fault interpretation, which can be problematic and time consuming. The uncertainty in manually interpreting the faults and fractures is very high, particularly in basins with sparse seismic data [10]. Several authors have illustrated the importance of classifying faults through various techniques [11][12][13] and fractures through an ant-tracking algorithm [14][15][16][17]. Consequently, an effective and reliable method is required to characterize these faults and fractures.
Contemplating past studies, the application of seismic attributes has played a vital role in aiding structural interpretation. Numerous authors have successfully applied structural attributes for a comprehensive analysis of the structural interpretation of seismic data [18][19][20][21][22]. Seismic edgedetection technology is the most commonly used method in the oil and gas industry to deduce the structural and stratigraphic features of seismic data [23]. However, this method is very challenging and makes it difficult to trace discontinuities that are concealed in the seismic sections, which may appear as slight changes in seismic waveforms, which makes it very challenging to interpret SSFs and fractures [24]. Recently, several authors have used curvature attributes to identify the orientations of faults and fractures [11,23,[25][26][27][28]. The coherency-based similarity attribute enhances the seismic edge of the seismic data [29]. The use of conventional attributes such as dip [11], curvature [23], coherency [30], dip and azimuth [31], chaos [32], and variance [33,34] can be used to successfully detect the subtle faults and fractures of the seismic amplitude volume in sedimentary basins. In our study, we also adopted unconventional geometrical attributes that are not routinely used for characterizing faults and fractures. These include thinned-fault likelihood (TFL) [35,36], fracture density [37][38][39][40][41][42], fracture proximity, and fault-enhancement filter (FEF) similarity [43], the last of which can be further observed through a voxel connectivity filter (VCF) to yield a 3D fracture network; interpretation was performed to rank these fracture geobodies. Afterwards, artificial neural networks (ANNs) were applied to these attributes to characterize the fracture network. These seismic attributes provided the basis for detecting small scale reflection continuities in the seismic data, which is difficult to do when using conventional interpretation techniques. Contemporary studies of multiattribute analysis via artificial neural networks have been widely acknowledged because of their effective approaches of identifying structural features [44][45][46][47][48][49]. ANNs are proficient in the recognition of complicated characteristics using numerous parameters [41] and are frequently applied in the structural and stratigraphic identification and interpretation of seismic features [50].
Fault attributes are accompanied and followed by an ant-tracking attribute to interpret the discontinuities of the seismic amplitude [11,32,51,52]. Ant-colony optimization (ACO) is an efficient method for detecting small scale faults and fractures within the reservoirs, and it is more reliable, accurate, and intuitive than other methods [12,[53][54][55]. Ant-tracking technology is a proficient tool in interpreting and augmenting the production of deployed fields [53]. The identification of small scale faults plays a vital role in detecting promising well-locations and explaining injection contradictions [8]. The small scale fault structure helps in identifying remaining reserves by analyzing fracture development using the ant-tracking method [56]. The application of an ant-colony optimization algorithm conveys the basis for the future development and exploration of remaining hydrocarbons, and it also gives a reference for analyzing the micro-scale faults and fractures of similar reservoirs.
Multiple studies have been carried out to explore and classify reservoir hydrocarbons. Still, the vast majority remains unexplored due to insufficient attention to comprehend the fracture network of reservoir sands. Several authors have studied thin sections to illustrate the presence of fractures in their study areas [57,58].
In some scenarios of these exploration ventures, the geophysicists and geologists have had solitary seismic data for interpretation. Due to the absence of well logs, image logs, and core data in these situations, unconventional seismic attributes and unsupervised neural networks play a vital role in characterizing fault and fracture networks by using conventional seismic data. Therefore, the present research aimed to provide an innovative workflow of data conditioning by using integrated geostatistical filtering for the dip and azimuth calculation, as well as structural filtering for the enhancement of structural features such as faults and fractures. Moreover, the study also provides an advanced workflow using the multi-attribute analysis, neural networks, and artificial ant-tracking for the identification of fracture parameters. Additionally, the current study also aimed to identify the small scale fault and network of fractures in reservoir sands using 3D seismic data within the Sawan gas field of the Lower Indus Basin, which has been previously overlooked. The methods applied in the study for the delineation of the small scale faults and fractures incorporated the combined use of data conditioning, multi-attribute analysis (using dip, curvature, variance, similarity, thinned-fault likelihood, fracture density, and fracture proximity), ANNs, and the ACO approach. The applied ant-colony optimization algorithm was followed by the automatic fault extraction technology to reveal the dip azimuth, dip angle, fracture length, and surface area of the fractures.

Geological Settings
The study area (Sawan gas field) is located at the junction of the Middle Indus Basin and Southern (lower) Indus Basin of Pakistan. The lower Indus Basin is located in such a way that Indian Shield is located in the East; Sulaiman Fold, Thrust Belt, and Kirther Ranges are located in the West; Jacobabad-Khairpur High is located in the South; and Sargodha High is located in the North [59]. During the early Cretaceous time, the Indian plate drifted northward, where marine shales and limestones of the lower Cretaceous section of the Sembar and Goru formations were deposited on the northwest margin, whereas the carbonates were deposited on the northern shelf of the late Cretaceous [60]. The Precambrian basement of the basin and the majority of the Cretaceous rocks within the study area are tectonically stable and controlled [61][62][63][64]. Several structural lows and highs have been identified from the Precambrian, and these played important roles in the lateral and vertical thicknesses of the reservoir. Amongst these lows and highs, Jacobabad-Khairpur High, illustrated as a horst, is situated in our research area ( Figure 1a). The Paleocene Ranikot clastics pinch out along the Khairpur High, which gives an indication of a preliminary Paleocene uplift episode that affects the Khairpur High. The enlargement of different sub-basins and complicated basinal topography during the Eocene suggests the occurrence of carbonate accumulations along with deep marls that were accumulated with Ghazij shales in the later stage. The tectonostratigraphic of the study area resulted from three tectonic episodes: erosion and uplifting in the late Cretaceous, basement rooted NNW-SSE wrench faults, and the recent uplifting of Jacobabad-Khairpur Highs of late Tertiary [65]. These faults were formed as a result of transtensional tectonics that were linked with Indian and Eurasian plate's collision [66]. The stratigraphy of the study area incorporates Chiltan limestone (Jurassic age) at the base, which is overlain by different rocks of Sembar and lower Gou Formation (Cretaceous age) to alluvial rocks (Quaternary age). The primary lithologies of the lower Indus Basin are shelf carbonates and clastics, marine marls and shales, and turbidities [67] ( Figure 1b). Several authors have researched the existence of wide Cretaceous deltaic settings and related shelf margins [68,69]

Methodology
The data were acquired from the Directorate General Concession Petroleum (DGPC), Ministry of Petroleum, Islamabad Pakistan. The seismic data were composed of a high-resolution 3D poststack seismic volume that covered an area of approximately 25 km 2 of the Sawan gas field. The seismic grid contained an inline that ranged from 643 to 867, and the crossline ranged from 893 to 1031. A total of 225 inlines and 139 crosslines were used. The length of the inline was 3459.09, and the interval of the inlines was 25.03, whereas the length of the crossline was 5606.15, and the interval of the crosslines was 25.07. The sample value format was the International Business Machines (IBM) floating point. The sample interval was 4 ms, and the number of traces per sample was 1125.
The workflow employed for the current study incorporated the analysis of data conditioning, evaluation, and selection of seismic attributes, neural network computation, and artificial anttracking ( Figure 2). In the first step, the data were conditioned, a process that incorporated seismic conditioning. The seismic conditioning included geostatistical and structural filtering to enhance and sharpen the structural features of the seismic data. In the second step, multiple structural attributes were analyzed to execute the selection process. The seismic attributes, which were extracting the best possible results of discontinuities, were shortlisted. In the third process, the shortlisted seismic attributes were tested and trained for the neural network analysis. In the fourth and last step of the study, the optimum attributes were selected and conditioned to process the ACO algorithm. The ACO was applied to evaluate the fracture network parameters. Figure 2. The workflow implemented for the current study comprised four steps; the first was dataconditioning, the second was the selection of attribute computation, the third was artificial neural network (ANN) computation, and the fourth was ant colony optimization (ACO) computation.

Artificial Neural Network
ANNs are computer-based mathematical models that are inspired by human brains that use continuous and parallel layers. These layers consist of nodes that are interconnected through weights [70]. These interconnected weights create associations amid the input layer and output layer of every node [50]. In the current study, unsupervised vector quantization (UVQ) was used for characterizing the faults and fractures. The UVQ used non-linear filtering and consisted of three steps: initialization, training, and application. The initialization evaluated the set of the classes, training improves these set of classes, and the application was used for the classification of the data. The number of classes was hard to determine in advance. Too many classes would have given redundant class centers, but too few classes would not have matched the seismic. Therefore, different quality checks were performed using different classes to proceed with the training.
In unsupervised training, network performance is tracked in a graph that shows the average match (confidence) of clustered input. Typically, the average match increases in a step-function. Each step indicates that the network has found a new cluster. Training can be stopped as soon as the average match has reached a stable situation (a straight line is obtained). Training that results in maximum match is accepted and is used to make the output maps.
Unsupervised classification maps are an advanced method that represents vectors that utilize a set of classes that lack user control [71]. The representing vector uses integrated, seismic attribute values that act as test data sets. Since the structural features are associated with erosional surfaces, reflection terminations, the variation of dip orientations, and changes in the amplitude signal, a set of seismic attributes are chosen as the input, which is very useful in identifying the network of discontinuities as an output.

Ant-Colony Optimization
ACO is a probabilistic method that is capable of resolving computational altercations that are constrained to finding paths in several ways [72]. It monitors the behavior of ants in the pursuit of food and fetching back to their nest by analyzing their paths [73]. The assemblage of the food is a unified effort of a group of ants for fetching the food, and they diffuse pheromones while making their way back; this diffusion appears as a fracture in a seismic volume. The diffusing pheromone can be recognized as an odor, and the scouting ants hunt the food by smelling this pheromone. The diffused pheromone is very volatile and, hence, evaporates. Thus, the solitary route remains, which is of shorter length and is followed by the majority of the ants to detect the noisy trends in the seismic data [11]. The scouting group of ants shows intelligence to extract the discontinuities that lead to colony optimization. The scouting ants detect their ideal track to food based on their intelligence and pheromone strategy [74], which is termed as swarm intelligence [75]. Swarm intelligence augments the discontinuous features in an edge-detection seismic volume since it tacitly tracks those trends that are continuous and are likely to be fault or fracture.
The algorithm supposes that the number of ants is proportional to the amount of diffused pheromone on a branch [15]. The pheromone on a branch at an intensity of path (i, j) at time t is denoted by , while the branch tracking (0 ≤ ≤ 1) and the pheromone amount per unit length of ant k at a path (i. j) is denoted by ∆ .
The branch length of ant k in a circle is represented by , and ∆ / , where Q is a constant. The factor η shows the visibility of the edge path (i, j), 0 shows the relative importance of the path visibility, 0 shows the relative importance of the branch, U represents the possible vertex set, and gives the transition probability of the ant k at time t. The probability of electing the branch was used in this study by using the 3D seismic volume to identify and sharpen the tracked faults and fractures [76].
The application of ant tracking entails several factors that extract fault information in a seismic amplitude volume. These factors incorporate an initial ant boundary, ant track deviation, ant step size, permissible illegal steps, obligatory legal steps, and stopping criteria [77]. The initial ant boundary governs that how narrowly the preliminary ant agents may be kept within a seismic volume. The ant track deviation illustrates the examining tolerance of the ants along the side of its tracking direction. A higher value of ant track deviation gives higher connections of the ant agents. The step size shows how far an ant can progress for every increment of its search. The permissible illegal steps showcase the ability of ants to image discontinuities further than the existing position. The obligatory legal steps factor works with the illegal step feature by acquiring the designated number of legal steps after an illegal step. The stop criteria also work with the illegal steps factor and do not permit the ants' movement when plenty of illegal steps are involved.

Fault and Fracture Extraction
Interpreters usually manually mark the fault on the seismic section to carry out fault interpretations, which is a very difficult and time-consuming process to accurately interpret a fault. The uncertainty in manually interpreting faults and fractures is very high. The automatic extraction of faults and fractures can minimize human interference by providing the maximum proficiency and accuracy in extracting the faults [78]. Automated fault extraction technology allows interpreters to comprehend a fault and mark selections amid the extracted surfaces instead of generating surfaces. However, the outcome of the automatic fault extraction technology is constrained by various factors. These factors have a direct influence on the extraction of faults and fractures, and the factors are patch downsampling, minimum patch size, connectivity constraint, extraction background threshold, extraction sampling threshold, extraction sampling distance, and deviation from a plane.
Downsampling is effective in controlling the density of the various points of every patch. The minimum patch size sets the minimum number of points for fracture patches. The extraction background threshold and extraction sampling threshold control the least threshold signal value for the extraction and estimation of fracture patches. The extraction sampling distance controls the least distance by defining the radius amid extraction points. The deviation from a plane (DFAP) controls the orientation and deviation of a fracture from the surface of the plane. The DFAP plays an essential role in distinguishing fractures. A higher value of the DFAP decreases the possibility of extracting more fractures and vice versa. A DFAP of 9 was used in the study because the fractures merged into each other when a higher value of the DFAP was used.

Seismic Conditioning
It is necessary to apply various filters to obsolete the redundant noise while keeping geometrical features of seismic sections, because seismic attributes are susceptible to the noise within the seismic data [79]. For this purpose, multiple filters can be used to remove noise in seismic data; these filters are mean, median, and diffusion. The mean filter improves long-wavelength trends and reduces short-wavelength features. Moreover, it measures the mean amplitudes of provided data at chosen intervals to eradicate noise. The median filter measures the dip and azimuth of the targeted window and substitutes the dominant sample position with the median position of the amplitude. In addition, it also improves the spatial features, which are continuous, and reduces the noise of the seismic volume. Together, the duo of the mean and median filters reduces the noise while improving the signal-to-noise ratio that may deplete the fault and fracture features that can be solved using a diffusion filter [80]. The diffusion filter enhances the features that are continuous by keeping discontinuous features such as faults and fractures. If discontinuous features are detected, then smoothing does not take place, and filtering only applies to incoherent noise and insignificant stratigraphic features.
Initially, a geostatistical filter was applied here, and a steering cube was generated using the fast Fourier transformation (FFT) of the background steering. The steering cube was prepared by using the dip and azimuth values [48]. The background steering cube was comprised of the structural dip and contained multiple filters to study the dip-associated faults, fractures, and sedimentary structures. Hence, FFT was applied to delineate the stratigraphic and structural information in the seismic volume. The dip was computed within a cube of 3 × 3 × 3 samples around each sample of the seismic window.
The generated background steering cube was further incorporated while applying structural filtering. The dip-steered median filter (DSMF) was applied to remove random noise and to improve the dip and azimuth information of the generated steering cube. The noise removed was analyzed by the DSMF using the mathematics attribute. A diffusion filter was then applied on the dip-steered median filter to make the smoothing influence obsolete on the discontinuous features. The dipsteered diffusion filter (DSDF) enhanced the low-quality seismic data nearby faults and fractures. The two DSMF and DSDF filters were joined together to design an FEF. The FEF was based on a similarity threshold, and the data were smoothed (DSMF) away from the discontinuities and sharpened (DSDF) at the fault locations. Both similarity and steering cubes were incorporated in the study to establish a filter that had the same-sized window as the median filter. The results of the fault enhancement filter enhanced the sharpness of the faults and fractures, reduced the noise, and increased the smoothing of the reflectors. The FEF was then applied to the seismic volume to delineate the discontinuous features. The results provided a comparison of the original seismic (Figure 3a), dip-steered median filter (Figure 3b), dip-steered diffusion filter (Figure 3c), and the fault enhancement filter (Figure 3d) at inline 665. The results showed that the FEF effectively showcased the discontinuities marked by black arrows, which were not easily identified on the original seismic, DSMF, and DSDF data. The FEF and detailed steering cube were further used in the multi-attribute analysis to characterize the SSF and fracture network.

Shortlisting of Seismic Attributes
Multiple structural attributes were applied to delineate the maximum structural variations. Initially, the attributes were performed on the specific inlines or crosslines using various parameters. Those attributes that extracted the maximum information of structural discontinuities were selected. The shortlisted attributes were then applied to the whole seismic cube. A short overview and the results of the shortlisted geometric attributes are presented in this section.

Dip Attributes
Dip attributes were useful for evaluating the curvature and the dip-steered attributes, and they displayed the faults with displacement offsets, which were considerably less than the width of the seismic wavelet [11]. Dip attributes can be calculated using various components of the dip of a volumetric seismic cube. These include gradient structure tensor [81], discrete semblance-based dip searches [82], and wavenumbers and instantaneous frequencies [71].
Several dip attributes were applied using the detailed steering cube as the input that incorporated the azimuth, inline-dip, crossline-dip, and polar-dip to study the effect of dip on the structural features. The polar-dip was the true dip, which was measured from the horizontal, and its range was always positive. The dip in the direction of inlines was inline-dip, while the dip in the direction of the crossline was crossline-dip. The results of the dip attributes were analyzed on the time slices at 1780 ms, which delineated a fault that was present near the southwest of the study area shown by a red arrow and that was located around crossline 984 and inline 665 (Figure 4a-d). The results showed that the polar-dip proved to more beneficial in recognizing the discontinuities (marked by yellow arrows) compared to the crossline-dip, inline-dip, and azimuth results.

Curvature Attributes
Here, curvature measured the extent of the bending of a surface at a specific point using volumetric methods and analyzed the geometry distortion of the seismic volume instead of analyzing the variations in the seismic amplitude. Curvature is an effective tool to detect small scale faults and fractures with low displacement that would otherwise not be possible to identify with other fault attributes such as coherency. Compared to coherency, the curvature attribute draws the discontinuous features using minimum and maximum anomalies. Thus, it does not depict the precise position of a fault [83]. The surface of a reflector is supposed to be a quadratic surface to compute the curvature attribute [28]. , Whereas, where p is the volumetric estimation of the inline-dip, q is the crossline-dip, and the coefficients become x = y = 0 [25].
Multiple curvature attributes were applied in the current study; these were maximum curvature, minimum curvature, the most positive and the most negative curvature attributes. Each of these curvature attributes had different parameters regardless of their signs to detect the fault and fracture discontinuities [11]. The results showed that the maximum curvature attribute was the most beneficial in terms of delineating the fractures in the study area (Figure 5a-d).

Similarity
Similarity attributes are a type of coherency that has a scale of 0-1 and that displays the similarity of the continuous features that are related to continuity of structural and stratigraphic trends [22]. Thus, similarity attributes are very effective in terms of delineating discontinuities. The algorithm of the similarity is semblance-based, which is evaluated for dips that have the maximum semblance selected for that dip [84].
The triple ( , p, q) represents the local planar event at specific time. , p, and q are the apparent dips at inline and crossline directions. H shows the Hilbert transformation of the real seismic trace (u).
A similarity algorithm was incorporated in the study to recognize the coherency in the inline and crossline directions, which allowed us to interpret the major discontinuities in specific directions. The full steered and non-steered similarity attributes were applied by using different windows to study the differences of the dip influence towards the discontinuities. Dip-steering eliminated the structural features of the similarity attribute, and precise results of faults were produced, whereas the non-steering similarity excluded the steering cube. The full-steering attribute proved to be more helpful than the non-steering attribute to delineate the discontinuous features (Figure 6a,b). Since the similarity attribute was very sensitive to the noise, the fault enhancement filter was also applied to the similarity attribute. The FEF improved the precision and sharpened the continuity of the faults. The discontinuities were compared through full-steered similarity, non-steered similarity, and FEFsimilarity (Figure 6a-c). A comparatively large window was used for increased continuity, highresolution results, and to extract the maximum information of the discontinuous features. Since the study area was relatively stable and the small scale faults and fractures are present. Thus, the effects of discontinuities were sharply analyzed through similarity attributes. The results of FEF-similarity were more continuous and sharp than the dip-steering similarity attributes. The results were also analyzed on the FEF similarity cube, which showed the significant discontinuities (yellow arrows) in the studied seismic volume (Figure 6d).
The FEF-similarity highlighted the major and thick discontinuities in the study cube. Still, these low similarity values were not able to detect the minor discontinuities and the association between these major and minor discontinuities. To study the connectivity of the discontinuities, TFL was applied.

Thinned Fault Likelihood
TFL is a novel innovative fault attribute that yields precise and exact discontinuities. It is defined as a power of semblance, and it has a range between 0 and 1. It generates seismic volumes with razorsharp edges that are very well suited for structural interpretation. The algorithm detects the range of fault dips to find the maximum likelihood to delineate the fractures in the area of interest and gives razor-sharp fault images on vertical sections and horizontal slices. A comprehensive study by [35] proposed a method of fault likelihood to recognize locally planar discontinuity. The TFL technique measures the likely combinations of dips and strikes to obtain a solitary direction, which gives the maximum fault likelihood for every sample image [35].
The TFL for every sample was documented in TFL images across the whole seismic cube, as shown in Figure 7. The inserted slice on the resultant TFL cube at 2176 ms showed that the area was affected by small scale fractures. The results showed that maximum fractures were present at the top from 600 to 900 ms, which were large scale fractures of higher lengths. Additionally, the TFL values were at minimum below 3000 ms, which showed the presence of relatively fewer fractures of a very small scale. The discontinuity features were further observed using the fault discontinuity attribute volume to create and rank the fracture bodies. The resultant TFL volume was used as input for estimating the fracture bodies using the VCF. The VCF is an effective tool to create continuous bodies based on the amplitudes in a stored volume. A voxel is defined as the volume around one sample. The sample's interconnection was computed based on amplitude criteria and geometrical spreading settings. The voxel connectivity filter was applied using the envelope of the amplitudes, which were higher than the 0.6 value; the envelope was defined for the fracture bodies to be computed. The propagation was done in all directions (26) using the six faces, twelve edges, and eight corners adjacent to the current seed. The bodies that were larger than 100 voxels were processed as relatively large scale fractures, and the output volume was processed using the body-size. The results of the generated connected fracture bodies were displayed using the voxel connectivity filter, as shown in Figure 8a. The results showed the orientations of the connected small scale and large scale fracture bodies. The facture bodies were oriented mostly in the east-west directions, and large-scale fracture bodies were present at the top, whereas the small scale fracture bodies were present at the bottom of the study area.
The fracture planes and sticks were also extracted using the discontinuous TFL volume as input. The fracture planes proved to be quite beneficial in recognizing the sweet spots for prospect evaluation, field development, and future drilling of exploration and appraisal wells [85]. The fractures were interpreted and visualized on the cube as a stick, and all sticks that belonged to one fracture were merged into a single stick set. Thus, a single fracture stick set was comprised of an unordered collection of the interpreted sticks. The minimum TFL value of 0.01 was used in the current study. The minimum value of sticks per fracture was applied, which helped to identify the small scale fractures that had small throws. The minimum vertical overlap rate of 0.6 was used to extract maximum tilted fractures, because a value greater than 0.6 may have extracted maximum fractures with a risk of picking up more noise. The minimum fracture length per z-slice of 250 was applied to extract the maximum fractures with smaller heaves. The feature of merging fractures along z-slices was also applied, and this increased the possibility of capturing the lateral continuity of fractures. The results of the extracted fracture planes and sticks generated for the seismic study volume are shown in Figure 8b.

Fracture Density
Fracture density highlights regions of high-density fractures. It requires a radius for scanning and computing the density of fracture anomalies. Thus, it is defined as the ratio of the number of traces classified as fractures to the total number of traces present in a circle of given radius along zslices. Here, the fracture density showed the improved visualization of potential fracture anomalies. The maximum curvature attribute was used as an input for the generation of fracture density. The results of the fracture density are shown in Figure 9a. The results of the fracture density showed maximum values along the erosional surfaces and in the proximity of the fracture activities. It was also essential to consider that structural discontinuities, such as fractures, did not represent the flow channels because channels and fractures had different morphology. Moreover, the regions of maximum fractures were more conductive than regions of minimum fracture activities.

Fracture Proximity
Fracture proximity helps to identify the locations of maximum fracture activity in a specified radius. It computes the lateral distance (i.e., along z-slice) from a trace location classified as a fracture. Fracture proximity is also helpful in tracking radius for drilling purposes. The maximum curvature attribute volume cube was used as an input here for the generation of fracture proximity. The fracture proximity highlighted the distance from the fracture center. The result of the fracture proximity is shown in Figure 9b and suggested that the studied seismic volume was mostly affected by fracture activities because there were multiple regions across the whole studied area that showed the distance from the fractures.
These fracture regions were also displayed as generated 3D volume that showed high-density and maximum fracture activity regions (marked by yellow arrows) (Figure 10a,b).

Neural Network Computation
After the selection process of seismic attributes, both supervised and unsupervised neural networks can be used to evaluate discontinuities [11]. In the current study, UVQ was used because the TFL and fracture density identified the most likely fault and fracture features. Therefore, the ANN was computed using the unsupervised algorithm.
The fracture positions were connected with bed terminations that were categorized by maximum variation of dip values, maximum curvature attributes, low values of similarity, maximum values of TFL, and maximum values of fracture density and fracture proximity. Henceforth, those seismic attributes that were already shortlisted in the current study (such as polar-dip, FEF-similarity, TFL, maximum-curvature, fracture density, and fracture proximity) were chosen for neural network computation.
The shortlisted seismic attributes were tested for input, and, after various quality checks (QC's), they went through a training program for ANN computation by using the UVQ of ten classes to yield an optimal output. The neural network tried to cluster the set of vectors, and similar vectors went into the same segment. Segmentation was performed using the shortlisted multi-attributes, and class centers were analyzed. The ANN-UVQ process comprised three distinct layers: the input, hidden layer, and output layers. These layers were connected with one another to make an entirely connected UVQ network that was comprised of ten nodes. Among them, the input layer consisted of six nodes (three related to the hidden layer), and the output layer was represented by one node. (Figure 11). The number of vectors assigned to each node and their average match to the corresponding class center were also analyzed, appearing as a robust classification. Figure 11. The ANN-UVQ (unsupervised vector quantization) network that was comprised of three layers of 10 nodes. The input layer had 6 nodes, the hidden layer had 3 nodes, and the output layer had 1 node. The colored legend shows the relative importance of every shortlisted attribute as an input to compute an output.
The input test data and the training data were applied in training, and this was continued until a straight line was obtained (that resulted above 75%), which suggested minimum error. The training was then ceased, and the results of the network were analyzed regarding the volume, section and time slices. The relative contribution of the seismic attributes UVQ classification is shown in Table 1. The results of the ANN-UVQ highlighted the fracture network using the weighted average of all the selected attributes. The seismic cube and the seismic section of the ANN-UVQ showed the discontinuities with maximum probability, as shown by the red rectangles in Figure 12a,b. The maximum likelihood of fractures was also computed at different slices, which indicated that maximum fractures were present at slice 700 ms. In contrast, the probability of fractures at 2180 ms was minimum (Figure 13a-d). Since the reservoir in the study area lies from 2100 to 2400 ms. the likelihood of the fractures was also analyzed at 2250, 2300, 2350, and 2400 ms, thus giving the idea of the trend of the fractures in the study area (Figure 13e,f). These slices could be further exploited for future drilling prospects, if further evaluated.

Attribute Conditioning
Before applying the ant-tracking algorithm to the attributes, the conditioning of multiple attributes was done to calibrate the suitability of finding the best outcomes for fault and fracture delineation. Therefore, attribute conditioning was used to improve the signal-to-noise ratio to show continuous fault and fracture features.

Structural Smoothing
The obtained amplitude volume of the Sawan gas field interrupted amplitude continuity at several locations. Henceforward, structural smoothing based on the Gaussian smoothing algorithm (GSA) was incorporated to reconstruct the discontinuities. The structural smoothing measured the guided signal of the structural features by enhancing the continuity of the reflectors [69]. The azimuth and dip parameters were used to measure the structural features. The Gaussian smoothing algorithm was then used adjacent to these structures. GSA is an efficient tool to use before auto-tracking to steady results [77]. The estimation of structural smoothing use the value of the Gaussian filter, which governs the number of traces horizontally and samples vertically. The value symbolizes the standard deviation of the applied Gaussian filter. The maximum value was used in the study so that the maximum number of traces and samples could be smoothened. A standard deviation of 1.5 was applied to get the smoothed volume.

Chaos
The term chaoticness is used for highlighting the reflector discontinuities using differences in azimuth and dip [77]. The maximum value of chaoticness indicates the regions of discontinuities such as fractures, faults, or unconformities [86]. The signal pattern of the chaotic texture also indicates the migration of gas pathways. Here, the 3D map of the chaos attribute showed many discontinuities in the upper portion of the seismic section, while there were fewer discontinuities in the lower area. (Figure 14a)

Variance
Variance is an edge detection technique, and it calibrates the dissimilarities from a mean value to produce computationally proficient results that are much sharper than coherency [86]. The variance was incorporated to reduce the noise and to yield possible vertical smoothing. The variance was applied using the dip to highlight discontinuities. The value of the 0.6 plane confidence threshold was used for evaluating the discontinuities.
The variance attribute was able to more precisely and sharply detect discontinuities than the chaos attribute (Figure 14b). Henceforth, the variance was further used in the ant-tracking analysis for fault and fracture study.

Ant-Tracking Result
Ant-tracking was used to extract the faults from the pre-processed variance volume. Initially, structural smoothing was applied to enhance the continuity of the structural features (faults and fractures) using the Gaussian smoothing algorithm, and then the variance attribute was used as an input on the structural smoothing attribute to get the ant-tracking attribute; the results were displayed at various locations of inline 665 and crossline 972. The results showed that the original seismic volume interrupted amplitude continuity (Figure 15a,e), which was rectified with the structural smoothing effect to enhance and sharpen the visibility of discontinuous features ( Figure  15b,f). The variance attribute was applied to the structural smoothing seismic volume, which showed the discontinuous features (Figure 15c,g). The results of the ant-tracking applied to the variance attribute precisely characterized the faults and fractures (Figure 15d,h). The direction of fractures was evaluated by applying various azimuthal filters to the ant-tracking algorithm. The aim was to restrict the ants in specific directions to get optimum results. The anttracking using east-west (E-W) and north-south (N-S) azimuthal filters were applied to recognize the orientation of the fractures. The results showed that the north-south azimuthal filter helped to identify the north-south oriented fractures. In contrast, the east-west azimuthal filter was beneficial in recognizing the east-west directed fractures. Moreover, the results of applied filters revealed that the study area has more east-west oriented fractures than north-south oriented fractures (Figure 16ac).

Automatic Fault and Fracture Extraction Using ACO
Fracture detective attributes generate responses along the surfaces of a fracture and help to detect automatic fracture interpretations. These attributes may incorporate noises and residual responses, which makes it very challenging to interpret fractures. To overcome these challenges to map the fractures, attributes conditioning was applied through an innovative approach of the ACO algorithm.
The 2D (bird's-eye view) (Figure 17a) and the 3D view (Figure 17b) of the extracted fractures show the orientation of fractures. Fracture surfaces generated in the ACO method are commonly fractured segments instead of fracture surfaces [87]. The results of the extracted fracture network were displayed on the stereonet, which showed that 607 visible patches (fractures) were present, and these were scattered in different orientations. The results of the dip azimuth showed that most fractures were clustered in the NW-SE direction (Figure 18a). The histogram of the fracture dip angle showed that the most fracture patches lie between 16° and 32° (Figure 18b). The histogram of the fracture length suggested that most fractures have fracture lengths between 200 and 500 m (Figure 18c). The histogram of the fracture surface area indicated that many fracture patches are present between surface areas of about 10,000-30,000 m 2 (Figure 18d). All the results of the histograms showed that the frequency of the number of patches decreased with increasing fracture angle, fracture length, and fracture surface area. The results showcased the presence of small-to-medium scale fracture segments with an average length of about 625 m and that were gently dipping to about average 38°. The corollary is that the overall results of the workflow showed that these numerous fractures played an important role in the development of the reservoir potential in the Sawan gas field.
The current study also revealed a fault in the study area that had been previously overlooked. The fault is a small scale fault, and the orientation of the fault is NNW-SSE (Figure 19a-c). The throw and heave of the fault are minor and are present in the Upper Goru Formation of Upper Cretaceous age. This NNW-SSE oriented fault distributes the study area into east and west parts and is named the Sawan fault [65]. The generated steering cube for the dip and azimuth estimation, as well as the application of structural filtering to remove the noise, considerably enhanced the sharpness and continuity of the structural features such as SSFs and fractures. Hence, the current study provided a workflow that put forward an effective data conditioning method for fault and fracture identification. Moreover, the adopted workflow also presented an advanced and novel approach of attribute computation using ANN-UVQ for the identification of the fracture network.
The current study proved helpful in identifying the fault and fracture probabilities of the study area. However, various technical aspects can be upgraded for future works. Instead of using the unsupervised ANN, other unconventional techniques, like the supervised ANN such as multi-layer perceptron network (MLP), can be used. The limitation of the presented work is that our ANN model was intended and trained to trace the fault and fracture probabilities only. In future research, the ANN model should be designed in such a way to extract information regarding fault slip and fault displacement. Another limitation of our work is that the current research lacked the usage of image log and core data, which were unfortunately not available. For future study, the use of image logs, conventional logs, and their correlation with core data along with 3D seismic is suggested to give a more comprehensive picture of large and small scale fracture development in the study area.

Conclusions
The current study provided an effective and advanced workflow for the identification of SSFs and fractures of a heterogeneous sand reservoir. The conclusions are as follows; 1. The structural features of the seismic volume were enhanced and sharpened using a DSMF, a DSDF, and a FEF. The FEF proved to be most effective in recognizing discontinuities. 2. The novel technique of TFL was used for the generation of the likelihood, dip and the strike of the seismic cube. The generated TFL cube highlighted the maximum likelihood of the dips and the strikes of the fractures. The interconnected VCF ranked fracture bodies, automated fracture surfaces, and fracture sticks using TFL, and it delineated the orientations of fractures. 3. The use of fracture density and fracture proximity are powerful tools in visualizing the highdensity and maximum fracture activities regions that can be exploited for future drilling. 4. The ANN-UVQ using multi-attribute computation identified the maximum and subtle fault, as well as the fractured locations and orientations, which will be beneficial for the field development of the study area. 5. The ACO algorithm was effectively applied to study the dip, length, azimuth, and surface area of the fractures. The results of the ACO showed that there are more E-W oriented fractures than N-S fractures. The automatic extraction of fractures using ACO identified 607 subtle fracture patches. The dip azimuth of these fractures are clustered at SE-NW direction, the dip of most of the fractures lies between 16° to 32°, the fracture length lies between 200 and 500 m, and the surface area lies between 10,000 and 30,000 m 2 . 6. The ANN-UVQ and ACO revealed an NNW-SSE oriented fault that has minor heave and throw. 7. The adopted workflow in the study for the recognition of SSFs and fractures is automated, adaptive, time-saving, cost-saving, effective, innovative, and novel. The applied workflow is advanced and can be further utilized for the delineation of SSFs, large-scale faults, and fractures in any reservoir worldwide.

Conflicts of Interest:
The authors declare no conflict of interest.