Utilizing Remote Sensing and Satellite-Based Bouguer Gravity data to Predict Potential Sites of Hydrothermal Minerals and Gold Deposits in Central Saudi Arabia

Abstract

This article aims to aid in exploring and forecasting hydrothermal minerals and gold deposits in Central Saudi Arabia (SA), with a focus on structural contexts. Remote sensing (RS) and satellite-based Bouguer gravity (SBG) data were integrated in order to create a mineral prediction map for the researched location. Data from the Landsat Operational Land Imager (OLI) and Shuttle Radar Topography Mission (SRTM) were transformed and enhanced using a variety of approaches. The delineation of hydrothermal alteration zones (HAZs) and highlighting of structural discontinuities in the OLI data were made possible using band ratios and oriented principal component analysis (PCA). Additionally, the underlying structural features were successfully exposed by processing the SBG using a variety of edge detection techniques, like the analytical signal (AS), total horizontal derivative (THD), tilt angle (TA), horizontal tilt angle (TDX), theta map (TM), horizontal derivative of the tilt derivative (HD_TDR), horizontal gradient of the tilt angle (HGTA), tilt angle of the analytical signal (TAAS), and soft sign function (SF). As a result, more prominent lineaments were found in the NW–SE, NNW–SSE, NE–SW, and NNE–SSW directions than in the N–S and E–W directions. The GIS incorporated surface/subsurface geological structure density maps with zones of hydrothermal alteration. It was found that the lineaments derived from the analysis of the RS and SBG data were more in line with the HAZs, which demonstrated the common connection between alteration zones and deep lineaments. The findings revealed a mineral prediction map with extremely low to extremely high probabilities. Overall, combining RS and SBG data effectively identified probable mineralization sites associated with hydrothermal processes and made it easier to create this study’s final predictive mineralization map.

1. Introduction

The Arabian Shield (ASh) is considered an economically significant shield worldwide, as it has promising possibilities for exploring metallic ore deposits. The ASh of Saudi Arabia (SA) consists of stratiform, vein, and magmatic to late magmatic contact metamorphism deposits [1]. Over 800 occurrences of gold (Au) and other metal deposits are found throughout the ASh, including porphyry copper, volcanic-hosted massive sulfide (VHMS) (Al-Amar, Gabal Sayid, and Al-Nuqrah), and Au ore deposits, such as those in the Bulghah and Mahd Ad Dahab mines (in the production stage) [1,2].

Significant gold mineralization is found in the central zones of the Najd fault system (NFS), specifically in Ghadarah, Fawarah, Hamdah, and Gariat Avala. Major fault zones parallel to the Red Sea also run northward from the south, becoming more connected as they go deeper underground, and they overlook several gold mines, including the Wadi Bidah, Jadmah, Wadi Laif, and Mamilah mines (Figure 1a) [3].

The Arabian-Nubian Shield (ANS) is considered significant for gold exploration and mining in SA. It is crucial for relationships between gold genesis and orogenic tectonics. Based on the tectonic settings and gold host rocks, primary gold mineralization can be divided into three basic types [4], as described in the following paragraphs.

The first type includes volcano-sedimentary connections, such as the Siham Group, and epithermal deposits, like the Mahd Ad Dhahab and Al-Amar ore deposits, as well as volcanogenic massive sulfide (VMS) gold deposits.

Based on Al-Shanti [5], the second type of gold mineralization is found in listwaenite rocks, which were formed after ophiolitic rocks were carbonatized. These rocks have several suture zones. For the Saudi Ma’aden Mining Company (SMMC), Bi’r Tawilah (study area, e.g., Masarah and Al Mansourah) is a suitable candidate for gold-mineralization-hosted listwaenite in the Nabitah suture (NS) zone and the western section of the Afif terrane (AT) [4].

The third form of gold deposit is hosted by late- to post-tectonic (640–610 Ma) plutons of diorite–granite and/or their subvolcanic substitutes that intruded the Bani Ghayy and Murdama Groups, chiefly abundant in the western (AT) and eastern (NS) zones, such as in Bi’r Tawilah, Sukhaybarat, Bulghah, Zalm, Jabal Ghadarah, and Ad Duwayhi [6,7,8,9,10]. Mesothermal quartz-vein-type gold ore mineralization, developed in intraplate settings, controls these deposits, which are magmatically and/or structurally controlled [11].

The majority of gold is found as loose microscopic specks along arsenopyrite, pyrite, and occasionally chalcopyrite fractures, which were produced by CO2-rich ore fluids within these deposits, which are connected through quartz veins and stringers [4,9].

Traditional field geology surveys are extremely costly, time-consuming, and difficult for mineral exploration in hydrothermal alteration zones (HAZs). Furthermore, interpolations are widely used in the outcomes achieved. In this situation, the combination of RS and SBG data may be considered as a substitute method for carrying out regional geological mapping with tolerable precision in a manageable amount of time.

The major goal of the current research is to recognize the best locations for mineral resources and gold deposits, using a combination of surface and subsurface structure components and HAZ areas. RS data (Landsat 8 OLI and SRTM data) and SBG data were used in conjunction in order to identify HAZs and structural lineaments in the area under investigation.

Reconnaissance investigations frequently utilize remote sensing technology. Remote sensing is a crucial method for understanding several aspects, including geological, geomorphic, structural, and mineral deposits. Multiple generations of unmanned satellites, called Landsat, ranging from their initial generation to generation 9, have been developed and released by the National Aeronautics and Space Administration (NASA). These satellites collect vital RS data regarding mineral resources, structural aspects, and land use/land cover. Geomorphic and geological components are described in [9,12,13,14,15,16,17,18,19,20,21,22,23]. In order to emphasize the spectral contrasts between the bands, ratio images have been employed. This method has been applied in remote sensing in order to map lithology and HAZs, and to mark spectral differences [24], so that a high contrast range of spectral properties—not possible with separate bands—may be displayed.

Based on fluctuations in rock density, SBG data play a vital role in finding subsurface structures. These geological formations serve as conduits for mineral deposits connected to hydrothermal systems. An essential problem to address when examining potential field data concerns locating and relating subsurface density patterns to subsurface structures [25,26,27,28,29]. Mineral exploration and structural studies can use edge detection methods to delineate the boundaries of subsurface structures [29,30,31,32]. Many approaches can be applied in order to detect structural edges, with most methods depending on the source field’s derivatives [33,34]. Previously, the limits of sources buried at varying depths on an observational plane were explored using various enhancing filters [35,36]. Several studies have increased the accuracy with which filters detect the boundaries of buried HAZs [3,37,38]. Many edge detection strategies can be applied in order to achieve the goals of such investigations. The first-order vertical derivative (FVD), total horizontal derivative (THD) [23], analytical signal (AS), tilt derivative (TD) [39], tilt angle (TA), and theta filtering [40] are some of these methods. Recently, the horizontal derivative of the tilt derivative (HD_TDR), horizontal gradient of the tilt angle (HGTA), tilt angle of the analytical signal (TAAS), and soft sign function (SF), filters that are based on second-order derivatives of potential field data, were developed. These techniques effectively improve the resolution of the findings of boundary detection while avoiding the secondary boundary problem.

As a result, combining RS images and SBG data should aid in improving lithologic mapping, enhancing geologic structures, and identifying gold mineralization zones.

Figure 1. (a) Location map of the Arabian Shield (Ash) showing major ophiolite belts, tectonostratigraphic terranes, fault zones, sutures, and post-accretionary basins in Saudi Arabia. The locations of 14 significant gold deposits are shown and include our investigated site in Central Saudi Arabia and the Bi’r Tawilah gold mining area (Jabal Masarah, Al Mansourah, and Ghadarah). (b) A geological map of the Ash showing the study area location [41], with modifications [42,43].

Figure 1. (a) Location map of the Arabian Shield (Ash) showing major ophiolite belts, tectonostratigraphic terranes, fault zones, sutures, and post-accretionary basins in Saudi Arabia. The locations of 14 significant gold deposits are shown and include our investigated site in Central Saudi Arabia and the Bi’r Tawilah gold mining area (Jabal Masarah, Al Mansourah, and Ghadarah). (b) A geological map of the Ash showing the study area location [41], with modifications [42,43].

2. Geological Outlines

2.1. Geological Setting

The region under investigation is situated in the central (ASh) area neighborhood, including the Bi’r Tawilah gold mining area and its surroundings. Along the Ad Dafinah faults that divide the Afif terrane, the area is located east of the Jeddah terrane [44,45]. Latitudes 22°00′ to 23°00′ N and longitudes 42°00′ to 43°30′ E serve as its boundaries (Figure 1b). The Hulyfah volcano-sedimentary pattern (Bani Ghayy and Siham Groups, 690 Ma and 620 Ma, respectively; [46]) make up the majority of the Bi’r Tawilah area (Figure 2). Younger syn-orogenic rocks (granodiorite and porphyritic rhyodacite) that are 710–685 Ma old intrude into this region. The Bani Ghayy Group comprises siltstone, greywacke, and conglomerate and is intercalated with metabasalts and pyroclastics. In contrast, the Siham Group comprises severely deformed and foliated chlorite, hornblende, and sericite ortho-schists. According to Al Jahadli [9], both groups underwent metamorphism into greenschist–lower amphibolite facies. Ortho-schists and molasse sediments are the rocks that are visible in Jabal Bani Ghayy [47]. Gold and tungsten ore resources can be found at Bi’r Tawilah (Figure 2). The Bani Ghayy Group–Siham Group–N-striking thrust fault divides them [48].

A Pan-African example of gold mineralization, the Bi’r Tawilah gold mine in Central Saudi Arabia, includes hypogene and supergene minerals. The so-called “Nabitah orogenic zone” along the thrust fault (N–S) area hosts sulfidic gold deposits in the intermediate to felsic intrusions. Serpentines and associated listwaenites, porphyries, granites, and dioritic rocks are the four types of rock units that are present and are arranged from oldest to youngest. All kinds of rocks are subject to hydrothermal alteration, which includes chloritization, sericitization, carbonatization, and silicification [49].

2.2. Gold Mineralization and Structure Setting

Three gold mines (in the Bi’r Tawilah, Al Mansurah, and Masarah regions) are located within the research area. Mafic volcaniclastic rock units (Siham Group), conglomerate, and sandstone comprise most of the rock units in Bi’r Tawilah (the first area), with diorite and granodiorite only occasionally occurring at the surface. Jabal Ghadarh is a notable example of this. The N–S thrust fault (Bi’r Tawilah fault zone) deforms these rocks. Only one good barren quartz vein exposure, which is 55–120 cm thick and oriented NW, NE, and E–W, is present in the region. The Najd fault systems, which encompass the Bi’r Tawilah area, have structural data that are congruent with this [47].

Massive diamond drilling is mainly focused on listwaenite or carbonatized ultramafics, which have gold occurrences. The rocks made of listwaenite contain excellent gold occurrences. The boreholes in the Bi’r Tawilah region to the south reveal that minor ophiolitic serpentinite, diorite, deformed granodiorite, monzogranite, and porphyry intrusions constitute the principal lithologies present.

According to Surour et al. [49], the Bi’r Tawilah supergene zone exhibited the following characteristics: (1) the diorite was highly altered and had barren quartz stringers; (2) the granodiorite was affected by varying degrees of weathering; (3) the porphyry was weathered as indicated by the oxidized pyrite; (4) gold was found as inclusions in pyrite and silica in the deep horizons (i.e., hydrothermal alteration zones); and (5) in the weathered layers, supergene gold was discovered as the consequence of the oxidization of pyrite to goethite.

3. Materials and Methods

The present work utilized several analysis techniques to integrate the available datasets, including RS and SBG data [3,14,50], for mapping the structures and their related hydrothermal alteration zones (HAZs). These methods were used to connect the presence of mineral resources with the surface and subsurface structures delineated from the integrated approach. Figure 3 shows the integrated approach’s workflow and each method’s details.

3.1. Remote Sensing Data

The Landsat-8 OLI (L8) dataset covering the investigated area originated from a mosaic of four scenes. The dataset was acquired on 4 July 2019 and was downloaded for free from the USGS Earth Explorer website (https://earthexplorer.usgs.gov, accessed on 20 July 2023). There are 11 bands in this dataset, each with various wavelengths and resolutions (Table S1).

For the L8 raw data to be prepared for further processing and to be in an adequate condition for investigation, pretreatment had to be carried out. Radiometric calibration, atmospheric correction, and minimum noise fraction (MNF) were all included in the preliminary processing.

Radiometric calibration and the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction model were applied using the ENVI software version 5.3 by L3Harris Geospatial USA (Broomfield, CO, USA), as described by Jahvar et al. [51], to calibrate the L8 data with the spectral radiance, subset them, and convert the digital number (DN) to reflectance.

MNF segregates noise in the data from the spectral information [52,53]. L8 bands in the initial space are converted to eigenspaces via MNF. A data space can be divided into two distinct regions: one region is characterized by large eigenvalues and coherent eigenimages, while near-unity eigenvalues and noise-dominated images represent the other region. Eliminating unwanted information from the data and only using coherent portions will improve the outcomes of spectral processing [53]. The MNF bands were converted into their initial data space minus the noise band(s) for bands containing eigenvalues less than 2, which were then regarded as noise [54]. There were only seven bands for which the MNF approach for the L8 subset covering the research area was employed.

3.1.1. Band Ratios (BRs)

The band ratio technique is applied mainly in geological applications [55,56,57]. The band ratio approach divides the digital data value of one band with high reflectance by the digital data value of another band with low reflectance [58,59,60]. This method reduces the negative impact of environmental details, such as the Earth’s surface, while enhancing contrast and compositional data. Some characteristics that cannot be seen with raw data are highlighted by the shadow cast by the terrain’s surface.

3.1.2. Crosta Analysis

To distinguish between hydrothermally altered minerals, Crosta analysis is used to decrease the number of spectral band inputs for principal component analysis (PCA). This is carried out to enhance the likelihood that the unwanted spectral bands will be identified, which can then be simply mapped onto just one of the principal component images. The main component images that directly emphasize the spectral characteristics of particular targets were chosen using the Feature-Orientated Principal Component Selection (FPCS) technique. This technique is used to examine the loadings of PCA eigenvectors to identify the components that capture information relevant to the target spectral signatures. The features that are most relevant to the theoretical spectral signatures of specific targets are then used to detect the targets in the image. The pixels in those selected principal component images that correspond to the targets are assigned as either dark or bright pixels, indicating the presence of the targets [61].

The USGS spectral database was used to determine the spectral characteristics of iron oxides and minerals that contain hydroxyl [62]. Bands 2, 5, 6, and 7 of the Landsat 8 OLI were chosen based on this spectral database to emphasize the spectral responses of iron-oxide-bearing and hydroxyl-bearing (clay) minerals.

3.1.3. Digital Elevation Model (DEM) of SRTM

The Shuttle Radar Topography Mission (SRTM) data were obtained using the radar C-band at 3 arcsec (90 m at the equator) during a collaborative effort by the National Aeronautics and Space Administration (NASA) and the National Geospatial-Intelligence Agency (NGA—previously known as the National Imagery and Mapping Agency, or NIMA), along with the participation of the German and Italian space agencies [63]. Digital elevation model (DEM) data extracted from the SRTM data [64] were used to recognize and extract geological lineaments by calculating and interpreting DEM derivatives, including shaded relief maps, slope maps, and traverse profiles [65].

Fill Sink

The resolution, or the separation among sample points, is a key factor in determining the correctness of DEM data [66]. Sink (and peak) refers to a standard error in DEM images brought on by resolution issues or by rounding elevation values to the closest integer. In the ArcGIS 10.9 environment, a method known as fill sink was used for the SRTM DEM data of the research region to remove sink errors. A hill shadow thematic map was created using the rectified image.

Shaded Relief

The topography is depicted visually on shaded relief thematic maps using gray values from raster imagery [59,60,61]. The borders separating regions that are darkened and those that are not shaded could be lineaments. Examination of the shade map obtained from the DEM also helps to characterize the lineaments. Therefore, a minimum–maximum range along a color ramp was added to a thematic shaded relief map that was created using the DEM. An azimuth of 135° (southeast), which represents the sun’s position in relation to the horizon, was selected after examining the sun’s various positions and relative angles since it may reveal the shaded and unshaded areas more clearly than other values. The sun’s angle of elevation over the horizon was set to 45 degrees. The spatial analysis hill shading tool of the ArcMap mapping program was applied to complete the work.

Lineament Extraction

The PCI Geomatica 2014 program was used to automatically extract lineaments from the shaded relief image created from the SRTM DEM. The two primary processing stages used for automatic lineament extraction are described below. The first step is the recognition of edges that reveal regions of sharp variations in the values of nearby pixels, and the second step is the identification of the existence of lines [67,68]. These procedures mostly employ the PCI Geomatica software’s LINE module, a popular package for automatically extracting lineaments [55].

An experimental test of automatic lineament extraction was conducted using various value ratios for each LINE module parameter. The test was run on an optical dataset to determine the best set of parameters that returned accurate results for the lines that were extracted [67]. These parameters were applied to the used data in accordance with the literature, and their validity was confirmed by visual assessment. Among this study’s LINE module settings, the default, suggested, and defined values are all shown in Table S2.

3.2. Satellite-Based Bouguer Gravity Anomaly

This study used a high-resolution Earth gravity field model, SGG-UGM-2, as a satellite-based Bouguer gravity anomaly. The SGG-UGM-2 model is a combination of marine gravity data derived from satellite altimetry data and the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE), the normal equation of the Gravity Recovery and Climate Experiment (GRACE), and EGM2008-derived continental gravity data up to degree 2190 and order 2159 [69]. The SGG-UGM-2 data are open-access and available at (http://icgem.gfz-potsdam.de/tom_longtime, accessed on 24 July 2023).

Edge Detection Operators

• Previous Edge Detection Methods

Boundary analysis operators are utilized to find the edge of the subsurface source. These methods are used for potential data, such as gravity and magnetic data. The formulae for these methods are presented below.

The analytical signal (AS) method is presented as the square root of the sum of the squares of the x, y, and z derivatives, which is known as the 3D AS [70]:

$$\left|\mathrm{A}\mathrm{S}\right|=\sqrt{{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{y}}\right)}^2{+\left(\frac{\partial \mathrm{T}}{\partial \mathrm{z}}\right)}^2}$$

(1)

In the result of the AS, the edge of a source is shown at the boundaries around the maximum values. The advantage of the AS method is that it is less dependent on the direction of the magnetization vector when using magnetic data. However, the AS method poorly balances the signals from shallow and deep anomalies [71].

The total horizontal gradient (THG) method is another common method, and the formula for the THG is given as [72]:

$$\mathrm{T}\mathrm{H}\mathrm{G}=\sqrt{{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{y}}\right)}^2}$$

(2)

In the result of the THG, the edge of the source is shown at the maximum values. The THG method is less noise-sensitive but poorly balances the signals from shallow and deep anomalies [71].

The tilt angle (TA) shows the edge of the source at zero values. The TA method can generate a balanced image for sources at different depths. The limitation of the TA method is that it produces secondary edges around the true edge [71]. The formula for the TA is given as [39]:

$$\mathrm{T}\mathrm{A}={\mathrm{t}\mathrm{a}\mathrm{n}}^{-1}\left(\frac{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{z}}\right)}{\sqrt{{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{y}}\right)}^2}}\right)$$

(3)

The TA method normalizes the vertical derivative (VDR) to the total horizontal derivative [73].

Another common technique in geophysical interpretation is the horizontal tilt angle (TDX), which is normalized to the absolute value of the vertical derivative [74]. The TDX method is represented as follows [75]:

$$\mathrm{T}\mathrm{D}\mathrm{X}={\mathrm{t}\mathrm{a}\mathrm{n}}^{-1}\left(\frac{\sqrt{{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{y}}\right)}^2}}{\left|\left(\frac{\partial \mathrm{T}}{\partial \mathrm{z}}\right)\right|}\right)$$

(4)

Wijns et al. [40] proposed the theta map (TM) method, which produces balances between edges located at different source depths. The TM method shows the edge of a source at the minimum values. In addition, the TM method has secondary edges around the true edge, like the TA method. The TM method is given as follows [40]:

$$\mathrm{T}\mathrm{M}={\mathrm{c}\mathrm{o}\mathrm{s}}^{-1}\left(\frac{\sqrt{{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{y}}\right)}^2}}{\sqrt{{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{y}}\right)}^2{+\left(\frac{\partial \mathrm{T}}{\partial \mathrm{z}}\right)}^2}}\right)$$

(5)

• Recent Edge Detection Methods

In recent years, many new boundary analysis methods have been proposed, such as the horizontal derivative of tilt derivative (HD_TDR) [76], horizontal gradient of tilt angle (HGTA) [77], soft sign function (SF) [78], and tilt angle of analytical signal (TAAS) [79] methods.

The HD_TDR method is like the analytical signal method but is sharper and generates better-defined maxima centered over the body edges of the source [76].

$$\mathrm{H}\mathrm{D}\_\mathrm{T}\mathrm{D}\mathrm{R}=\sqrt{{\left(\frac{\partial \mathrm{T}\mathrm{D}\mathrm{R}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}\mathrm{D}\mathrm{R}}{\partial \mathrm{y}}\right)}^2}$$

(6)

The horizontal gradient of tilt angle (HGTA) method is given as follows [77]:

$$\mathrm{H}\mathrm{G}\mathrm{T}\mathrm{A}=\sqrt{{\left(\frac{\partial \mathrm{T}\mathrm{A}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}\mathrm{A}}{\partial \mathrm{y}}\right)}^2}$$

(7)

The soft sign function (SF) method was proposed by [78] as follows:

$$\mathrm{S}\mathrm{F}=\frac{\mathrm{k}\times \frac{\partial \mathrm{T}\mathrm{H}\mathrm{G}}{\partial \mathrm{z}}-(\mathrm{k}+2)\sqrt{{\left(\frac{\partial \mathrm{T}\mathrm{H}\mathrm{G}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}\mathrm{H}\mathrm{G}}{\partial \mathrm{y}}\right)}^2}}{\sqrt{{\left(\frac{\partial \mathrm{T}\mathrm{H}\mathrm{G}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}\mathrm{H}\mathrm{G}}{\partial \mathrm{y}}\right)}^2}+\left|\mathrm{k}\times \frac{\partial \mathrm{T}\mathrm{H}\mathrm{G}}{\partial \mathrm{z}}-(\mathrm{k}+1)\sqrt{{\left(\frac{\partial \mathrm{T}\mathrm{H}\mathrm{G}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}\mathrm{H}\mathrm{G}}{\partial \mathrm{y}}\right)}^2}\right|}$$

(8)

where $\mathrm{T}\mathrm{H}\mathrm{G}=\sqrt{{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{T}}{\partial \mathrm{y}}\right)}^2}$, and $1\le \mathrm{k}\le 10$ will yield the best results [72].

The TAAS method shows the edges of the sources at the maximum values; however, the detected edges have low resolution. The formula for the TAAS is represented as [79]:

$$\mathrm{T}\mathrm{A}\mathrm{A}\mathrm{S}={\mathrm{t}\mathrm{a}\mathrm{n}}^{-1}\left(\frac{\frac{\partial \mathrm{A}\mathrm{S}}{\partial \mathrm{z}}}{{\left(\frac{\partial \mathrm{A}\mathrm{S}}{\partial \mathrm{x}}\right)}^2+{\left(\frac{\partial \mathrm{A}\mathrm{S}}{\partial \mathrm{y}}\right)}^2}\right)$$

(9)

4. Results and Discussion

4.1. Remote Sensing (RS) Data

4.1.1. Landsat 8 OLI Data

The investigation of the Landsat data allowed for the outlining of the HAZ areas. The spectral Sabin ratios (4/2, 6/5, 6/7) are represented in the RGB images and shown in Figure 4a [80]. The OH-bearing minerals (clay minerals, such as kaolinite, montmorillonite, and alunite) mainly manifested in a 6/7 ratio, appearing in the RGB images as blue-colored regions. The crimson hue of ferric minerals, such as limonite, hematite, and goethite, was mirrored in the band ratio 4/2. Areas with a high concentration of iron oxides were detected as reddish spots in the band ratio 6/5. A false-color composite map representing an alteration map for the research area was created using band ratios 6/7, 6/5, and 5 [81,82,83], with the HAZs denoted in yellow (Figure 4b).

The corresponding eigenvectors and eigenvalues of SPCA for iron oxide layouts are shown in Table S3, with band 4 (reflection) and band 2 (absorption) as distinctive bands for iron oxides and the typical PC4 for iron oxides. Both the eigenvectors and individual values of the hydroxyl minerals were calculated using the Crosta analysis of the Landsat 8 OLI. Table S4 shows that in PC4, there was a difference in eigenvector importing among bands 6 (−0.774) and 7 (0.504). As a result, the darkest pixels represented band 7 (0.59) and reflected the OH-bearing minerals. The hydrothermally changed portions emerged in a brilliant tone following the negation of PC4.

During the generation of the Crosta alteration image (Figure 4c,d), the iron oxide minerals, hydroxyl minerals, iron oxide minerals (F, H + F, H), and RGB images were shown in grayscale. The white pixels with bright reflections that identify alteration areas in the grayscale images are likely to be zones with the most future potential since they display both iron-stained and clay alterations. While hydroxyl minerals are placed in yellow to orange reflection color in the RGB images, areas rich in iron oxide minerals are identified in cyan to bluish coloration [83].

The grayscale image produced by ratioing bands 6 and 7 (6/7) traced clay and hydroxyl-mineral-bearing areas as white-colored areas or linear patterns scattered along possible tectonic elements and associated stream channels, implying that the distribution of these minerals may be structurally governed (Figure 4e,f) [84].

The images incorporate data from the band ratios and PCA-defined rock alteration caused by hydrothermal activity. These zones can be used to differentiate between locations with a high number of hydroxyl and iron oxide changes. As a result, the images were divided into five ranks. Figure 4g shows that the higher the rank, the better the mineral’s favorable prospective grade. The richness of HAZs often indicates mineral alteration regions, particularly those rich in gold mineralization, which are frequently associated with highly sheared areas [85].

4.1.2. STRM DEM Data

The topographic variances varied from 857 to 1251 m (a.s.l.) according to the analysis of the SRTM DEM (Figure 5a). Additionally, the DEM offered geometry-related information not made available by previous RS approaches. The geological analysis of the Landsat photos and SRTM data revealed structures with lineaments (Figure 5b). From the findings of the lineament analysis, a total of 552 extracted lineaments from the STRM DEM are displayed in Figure 5c, along with a density map (Figure 5d) and a rose diagram showing the spatial distribution of the lineaments. The rose diagram is used to indicate the observed lineament trends, with patterns in the NW–SE, NNW–SSE, NE–SW, and NNE–SSW directions being considerably more apparent than those in the N–S and E–W directions (Figure 5c).

Lineament density analysis was performed to determine the frequency of the lineaments per unit area. Also referred to as the lineament frequency, the lineament density is a measure of a collection of linear features in a location. This analysis led to the creation of a map (Figure 5d) that depicts the locations of the lineaments throughout the investigated area. According to [86], areas with a high lineament density are suitable for mineral prospecting.

4.2. Satellite-Based Bouguer Gravity (SBG) Data

In this study, the SGG-UGM-2 model was selected to obtain Bouguer gravity anomaly data between longitudes 42° and 43.3° and latitudes 22° and 23°, with a grid interval of 4 km. The Bouguer anomaly values changed from −62.8 to −99.6 mGal (Figure 6).

Several threshold analysis operators were used, including the AS, THG, TA, TDX, and TM. Furthermore, modern edge recognition filters such as HD_TDR, HGTA, TAAS, and SF were used, and the findings were correlated with the research area’s tectonic structure, according to the Hulyfah volcano-sedimentary pattern (Siham and Bani Ghayy Groups), which makes up most of the Bi’r Tawilah area and its surroundings. The Bouguer gravity map revealed higher values in the research region’s central, southern, and northeastern parts. At the low–high-anomaly interface, −75.1 m Gal values were reported near the Bir Tawilah gold mining area (Figure 6).

The regional high-gravity anomalies are associated with crustal thickening in the research region’s central, southern, and northeastern parts [87], implying that lithological crust changes cause this regional gravity anomaly gradient.

In the results of the AS method, the discontinuities cannot be seen, and the greatest values near the map’s north reflect a northwest–southeast trend. This result could be explained by the AS’s poor performance in recognizing shallow and deeper sources, and the discontinuities could be shallow sources in this location (Figure 7a).

The NW–SE and NE–SW trends in the THG, TA, TDX, and TM data can, however, be seen clearly. The THD results are consistent with the nearby active and potential faults (Figure 7b–e). The maximum values depict the boundaries of tectonic plates and discontinuities. Additionally, although the TDX approach provides the edges at the maximum values (Figure 7d), the TA result displays the fault edges at the zero values rather than the maximum values (Figure 7c). In the TM filter result, when the TM filter provides zero values, but the TDX filter yields maximum values at the exact locations, this may be connected to how the TM method produces secondary edges (Figure 7e). The TM and TDX maps have interconnected lineaments, which can make the interpretation more challenging.

The high-resolution SBG map underwent enhanced filtering using HD_TDR, HGTA, TAAS, and SF (Figure 8).

The HGTA generated from the tilt angle filter is shown in Figure 8b, along with its effects and outputs. The lineaments brought about by lithological variation were visible along the tilt angle’s zero-contour line, but this filter did not pick up on sharp edges [88].

The NE–SW, NW–SE, and N–S faults had prominent signatures on the HD_TDR and HGTA maps (Figure 8a,b). Therefore, they could not be used to identify the structural patterns in this situation. Strong filters called the TAAS and SF enable simultaneous equalization of signals with low and high amplitudes.

The geologic characteristics were considerably simpler to comprehend visually and qualitatively because the TAAS and SF maps (Figure 8c,d) provided a more detailed characterization of the features. A more exact length and width delineation can be seen on the SF map (Figure 8d). The visual comparison makes it simpler to understand the SF-derived lineaments. Subsurface linear structures detected in the research region according to the SF map (Figure 8e) exhibited lineaments, with the trends in the NW–SE, NNW–SSE, NE–SW, and NNE–SSW directions being more evident than those in the N–S and E–W directions. The lineament map derived from the subsurface gravity data highlighted the bulk of elongated NW–SE zones and their intersecting NE–SW trends (Figure 8e).

4.3. Integration of RS and SBG Results

Integrating multiple data sources in a spatial evaluation approach is essential in mineral prospect mapping [89]. To estimate the likely location of gold resources, a combined strategy has been used to characterize the areas of mineral occurrences [90]. This section combines the lineaments and HAZs produced from the RS and SBG data. This is helpful for investigation strategies in searching for minerals and serves as adequate evidence of mineralization.

By integrating the lineaments derived from the RS and SBG data, an integrative lineament map (Figure 9a) was constructed, which displays the critical surface and subsurface features (faults/fractures). Analysis of the RS- and SBG-derived lineaments revealed that the lineaments in the NW–SE, NNW–SSE, NE–SW, and NNE–SSW directions were more pronounced than those in the N–S and E–W directions. A density map was created using the data’s recovered lineaments (Figure 9b), and we classified it into four classes. This map showed a high density of fracture/fault zones in the southeast of the research region (south of the Bir Tawilah gold mine) and a medium density in some sections of the western and eastern parts. The density map also showed a prominent NW–SE tendency that crossed the NE–SW trend.

The geodynamic evolution of the studied area has been significantly controlled by the linear trends that have been found. These patterns, which govern the rock units and represent ideal locations for hydrothermal solutions, were confirmed in the Landsat data. The integrated and extracted lineaments from the RS and SBG data and the geological map were superimposed on the derived hygrothermal zone map to indicate the favorable and predicted sites of gold and hydrothermal mineralization (Figure 10). High-class sites indicate a desirable location and a greater density of mineral resources, considering that many mineral and gold deposits are linked to HAZs, including fluids that rise through fault lines and fissures.

As a consequence of this, regions that have a greater lineament density, which is associated with tectonic/structural lithosphere deformation [91], are more likely to have deposits of gold mineralization. Therefore, the areas with the greatest change and lineament density are compatible with the most appropriate predicted mineralization regions. In addition, the location of the existing gold mines also correlates well with these locations. As a result, the predicted Au mineralization zones are found in the southeastern, eastern, and northern regions (Figure 10).

5. Conclusions

The research region is situated in the center of the Arabian Shield (Ash) and includes the Bi’r Tawilah gold mine area and its surrounding area, which is considered to be one of Saudi Arabia’s largest gold mining areas. The current study’s primary objective was to identify the optimal locations of mineral resources and gold deposits by combining surface and subsurface structure components and hydrothermal alteration zone locations. The detection of HAZs and structural lineaments in the investigated area was achieved using an integrated approach based on remote sensing (RS) data (Landsat 8 OLI and SRTM data) and SBG data.

The merging of the RS and SBG data allowed for the detection of structural discontinuities and HAZs in the research area. In order to locate HAZs in the Landsat OLI data, image processing techniques including band ratios and PCA were applied. From the final image, four zones representing different degrees of alteration were identified. Lineaments were taken from the SRTM data analysis. Using the SBG data, several edge detection methods were utilized to build subsurface lineaments. The lineaments in the NW–SE, NNW–SSE, NE–SW, and NNE–SSW directions were more significant than those in the N–S and E–W directions, according to the analysis of the RS and SBG data’s lineaments. The HAZs and lineament density were combined to create a mineral potential map that maximized the area of possible mineral resources.

The lineaments derived from the RS and SBG data analysis were found to be more compatible with the HAZs. This demonstrates that deep lineaments are frequently linked to alteration zones. This combination reflects the efficacy of using RS and SBG maps to explore gold occurrences and the associated hydrothermal metallic mineralization.