Integrated Geophysical Assessment of Groundwater Potential in Southwestern Saudi Arabia

Saudi Arabia is seeking fresh groundwater resources to face the increase in anthropogenic activities. The groundwater storage variations and occurrence were investigated and the surface and subsurface structures in ﬂ uencing the groundwater resources in the research area were de ﬁ ned using a combined study of Gravity Recovery and Climate Experiment, aeromagnetic data, and electrical resistivity data with other relevant datasets. Results are: The groundwater storage ﬂ uctuation is calculated at − 0.34 ± 0.01 mm/yr during the period 04/2002-12/2021. The area is receiving an average annual rainfall rate of 117.6 mm during the period 2002 to 2019. Three structural trends, de ﬁ ned in the directions of NS, NNW, and NNE are cutting the sedimentary cover and the basement rocks. The sedimentary cover ranges from 0 to 1.2 km thick. Vertical electrical sounding results indicate three main geoelectric layers: the surface geoelectrical layer of higher resistivity values (428-9626 Ω . m) is made up of unconsolidated Quaternary sediments; the water-bearing layer of saturated sands with a resistivity range between 5.1 and 153 Ω . m and with depths vary from 1 to 94 m, and highly fractured basement rocks with resistivity values ranging from 813 to 6030 Ω . m. The integrated results are useful in providing a comprehensive image of the study area ’ s surface and subsurface structures, as well as groundwater potential in the southwestern part of Saudi Arabia. Our integrated approach provides a reproducible model for assessing groundwater potential in arid and semiarid areas.

Integrated Geophysical Assessment of Groundwater Potential in Southwestern Saudi Arabia INTRODUCTION Saudi Arabia is the largest arid country in the Middle East, with an area of 2.24 × 10 6 km 2 in the Arabian Peninsula. Saudi Arabia is an arid country with limited water resources. Due to the increasing demand and unreliable sources of surface water, about 75-85% of water supplies come from groundwater, which is essentially classified as a fossil water resource. Minor natural recharge might occur in the mountainous areas. Therefore, investigating groundwater occurrence plays an important role in the country.
The launch of the Gravity Recovery and Climate Experiment (GRACE) mission opened up new possibilities for monitoring and estimating the terrestrial component of the hydrological cycle, which includes groundwater resources (Tapley et al., 2004). GRACE anomalies in terrestrial water storage (TWS) are widely employed in hydrological applications ranging from regional to global scales (e.g., Wouters et al., 2014;Frappart et al., 2019). In recent years, several studies were carried out to assess the hydrological settings and components of major basins and large aquifers. For instance, GRACE and Global Land Data Assimilation System (GLDAS) data, with other observed data were used to differentiate the Mississippi River basin's water budget and came up with respectable results (Yeh et al., 1998;Yeh et al., 2006;Rodell et al., 2007). Furthermore, some studies in the Arab world have combined GRACE with other relevant information to quantify groundwater storage variability and estimate aquifer recharge and depletion rates (Mohamed et al., 2014;Mohamed et al., 2015;Fallatah et al., 2017;Mohamed et al., 2017;Fallatah et al., 2019;Mohamed, 2019;Mohamed, 2020a;Mohamed, 2020b;Mohamed, 2020c;Mohamed and Gonçalvès, 2021;Taha et al., 2021;Mohamed et al., 2022b). GRACE can be used to fill the blanks in hydrologically monitoring data (Famiglietti et al., 2011).
It offers global monthly fluctuations in terrestrial water storage (ΔTWS) (Wahr et al., 2004;Syed et al., 2008). GRACE data has a low horizontal resolution, and there is no vertical resolution for GRACE data, which makes this use difficult (e.g., Ahmed et al., 2016). GRACE is unable to distinguish between contributions from TWS's multiple compartments (for example, surface water, groundwater, and soil moisture). The combination of land surface model outputs with GRACE data is now allowing individual elements from GRACE-derived TWS estimates to be extracted and the data's horizontal resolution to be improved.
The magnetic method has been widely applied for delineation of the structural trends of the subsurface structures and estimating the depth of the basement crystalline rocks in many regions (Al-Garni, 2004a;Al-Garni, 2004b; Al-Garni  Sultan et al., 2009;Al-Garni, 2010). The aeromagnetic data is integrated with remote sensing data to investigate groundwater resources (Meneisy and Al Deep, 2020;Mohamed and Ella, 2021). Bakheit et al. (2013) used airborne magnetic data to trace and analyze structural trends to locate the position of tectonic lines and subsurface structures. Eldosouky et al. (2022a) used the analytical signal and the horizontal gradient magnitude to map the structures of Gabal Shilman area. Saada et al. (2022) have applied the regionalresidual separation, tilt derivative, and spectral analysis techniques to better understanding the structural features for the east of the Qattara depression area using magnetic and gravity data. Ranganai and Ebinger (2008) assessed the regional groundwater resources using an integrated approach of Landsat images with aeromagnetic data at south-central Zimbabwe Craton.
The geoelectrical resistivity, electromagnetic, seismic refractive, gravity, magnetic, and radioactivity techniques are among the surface geophysical methods utilized in groundwater study. Overburden thickness, aquiferous zones, bedrock architecture, and topography can all be estimated using these methods (Adagunodo et al., 2014;Joel et al., 2016;Adagunodo et al., 2017;Oyeyemi et al., 2017;Adagunodo et al., 2018). The current research is based on the use of electrical resistivity as a geophysical tool for groundwater exploration that is considered the most promising and appropriate technology for groundwater exploration and aids in the identification of the appropriate place for groundwater investigation.
The electrical resistivity method was widely used for issues associated with groundwater exploration (Mousa, 2003;Hosny et al., 2005;Nigm et al., 2008;Mohamaden et al., 2009) in different media by detecting the surface effects produced by the flow of electric current inside the earth and correlation between the electrical properties of geological formations and their fluid content. The aquifer material, pore space volume, fluid quantity, and salinity level all have an impact on electrical resistivity values. Moreover, the water-bearing formations' geometry has been assessed. (Robain et al., 1995;Robain et al., 1996). This electrical resistivity method was employed to identify the saline-freshwater interface zone (Yechieli, 2000;Choudhury et al., 2001), water-related parameters (Kosinski andKelly, 1981;Frohlich and Kelly, 1988;Troisi et al., 2000) and water quantity (Kessels et al., 1985). Furthermore, contaminated groundwater zones influence the electrical resistivity approach (Karlik and Kaya, 2001). The resistivity technique has been effectively utilized to determine the formations' thickness and groundwater potential (Raju and Reddy, 1998) and to study the groundwater potential and the aquifer system of fractured hard rocks (Chandra et al., 2012;Ologe et al., 2014; Al Deep Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 937402 3 et al., 2021). It has also been combined with borehole data to estimate geohydraulic parameters from Abi's fractured shales and sandstone aquifers (Ebong et al., 2014). In hard rock environments, geoelectrical resistivity techniques have been widely used to delineae shallow subsurface lineamnets and fractures (Ebong et al., 2014;Almadani et al., 2019) and faults Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 937402 4 (Suzuki et al., 2000;Fazzito et al., 2009), which are the targets of groundwater investigation (Titus et al., 2009). Ebong et al. (2021) have applied vertical electrical sounding (VES) and electrical resistivity tomography techniques to evaluate the geological features and groundwater potential of the Obudu basement complex in southeastern Nigeria. Because of its availability, simplicity, and developments in computer software and other inversion approaches, VES is the most often used groundwater prospecting technique (Kana et al., 2015). Morsy and Othman (2021) have applied an integrated approach of electrical resistivity and ground-penetrating radar with topographic data to evaluate the groundwater potential zones in the southwestern part of Makkah city. Their results show that the groundwater accumulation is localized in faultedbounded depressions and wadis. Sharaf (2011) has integrated the resistivity and seismic measurements with test holes drilled in the As Suqah area, Makkah ditrcit. His findings demonstrate that groundwater is found mostly in two water-bearing zones: alluvial deposits and the Haddat Ash Sham and Ash Shumaysi formations' clastic sedimentary rocks. Taha et al. (2021) used gravity and electrical resistivity data to investigate the groundwater resources of wadi Sar in Hijaz mountains on a regional and local scale. Their findings reveal a general downward trend in groundwater storage variation computed at −2.00 ± 0.34 mm/yr from 04/2002 to 07/2017, as well as the presence of a water-bearing layer with low resistivity and variable thickness in fractured basement rocks of Wadi Sar. Despite the fact that similar studies may be conducted in various locations in Saudi Arabia's southwest, no previous studies using our proposed geophysical methodologies have been conducted for the entire study area, which is located between longitudes 41.91°and 45.61°a nd latitudes 17.26°and 19.12°.
Other satellite and airborne geophysical field datasets have been utilized to explore crustal features at the continental scale , the geometry of the magma chamber and heat flow (Mohamed et al., 2022a), and land subsidence caused by heavy groundwater withdrawal Othman and Abotalib, 2019).
GRACE and GRACE Follow-On (GRACE-FO) solutions are used in conjunction with other climatic data to track spatial and temporal changes in water storage as a result of climate change and/or anthropogenic activities. The study also aims to investigate the groundwater potential of the study area using airborne magnetic data and electrical resistivity technique. Magnetic data was used to identify structural features and map the depth to underlying crystalline rocks that affect groundwater in the study area, whereas the shallow electrical survey was carried out to investigate the depth of the groundwater and the subsurface geology in the study area. Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 937402 5

GEOLOGICAL AND HYDROGEOLOGICAL SETTING
According to Greenwood et al. (1982), there are two primary subdivisions of layered rocks on the southern Arabian Shield south of latitude 22°N. An older ensimatic-arc complex and a younger marginal-arc complex are involved. Greywacke and mafic to intermediate volcanic rocks of the essentially contemporaneous Baish, Bahah, and Jidah groups make up the older arc complex. The younger arc complex is ensimatic in nature and partly superimposed over the older arc complex. The Halaban group, which was located to the east and northeast of the older ensimatic-arc complex, formed part of the younger arc group's ensimatic portion. The Ablah, Samran, and maybe Arafat groups reflect the superimposed parts of the younger arc complex. The Halaban group consists of andesitic and dacitic volcanic rocks as well as related clastic sedimentary rocks.
Based on the distribution of layered rock units of various ages and the orientation of overlying structural features, a number of N to NE-trending belts were delineated in the southern part of the Arabian Shield. Major fault zones run the length of these belts, defining their limits. The majority of these boundary faults go north and dip severely. Folds and other lineations in and around them frequently plunge down dip, both shallowly and steeply. Many of these faults are folded in random locations along their strike (Greenwood et al., 1982).
Our study region is located in the southwestern part of Saudi Arabia ( Figure 1A). Its geology is shown in Figure 1B. The study region's exposed basement rocks are made up of two separate units: Khamis Mushayt Gneiss and pegmatites. The Khamis Mushayt Gneiss rocks, which include banded orthogneiss, migmatite with little paragneiss, and amphibolite, are the major basement unit in the area. There are various aplite and pegmatite dikes in this unit. Quartz, orthoclase, plagioclase, and biotite make up the pegmatite unit. Flat-lying Cambrian-Ordovician sandstones and spatially variable alluvial deposits of pebbles, gravels, sands, and clays overlay these two basement layers in some locations. In the western part of the area, the old Wajid Sandstone rocks appear. After the end of late Proterozoic igneous activity, the Cambrian-Ordovician Wajid Sandstone rocks were deposited.
According to field research, the majority of the aquifer in the researched area is made up of cracked and jointed hard rocks, which are then refilled by precipitation via infiltration. Shallow aquifers with varied porosity and permeability are formed by cracked and jointed hard rocks, affecting the aquifer storage coefficient. Groundwater flow is reduced by fine sediments deposited in joints, cracks, and faults in hard rocks, allowing runoff loss and aquifer recharging. By digging the rocks below, the water from these aquifers is pumped through the drilled wells (Khan, et al., 2022).

GRACE/FO Data
GRACE is a pair of satellites that measure changes in the Earth's gravitational field in both spatial and temporal variations. It was launched in 2002 as a collaborative US-German project (Tapley et al., 2004). Changes in water content are the primary cause of changes in the Earth's gravity field. There are two sources of the GRACE/FO datasets, represented by the spherical harmonic and mass concentration solutions (mascon). In this study, the mascon products were used in the current study. With improved spatial resolution and minimal inaccuracy, these mascons capture all of the signals detected by GRACE within the GRACE noise limits.  (Watkins et al., 2015;Wiese, et al., 2016;Wiese, et al., 2018;Landerer et al., 2020). A subset of these JPL mascons crosses coasts and includes signals from both land and sea mass changes. As a post-processing step, the Coastal Resolution Improvement filter was used for the entire mascon solution to detect land and ocean mass from individual mascons that cross coastlines (Watkins et al., 2015). Therefore, the scaling factor must apply to the JPL-M to recover the leakage signals. In comparison to the RL05 version, the CSR mascon products of release 06, and with 0.25°× 0.5°grids (Save et al., 2016;Save 2020) use a freshly designed grid. The hexagonal tiles that cover the shoreline in this new grid are split into two tiles along the coast to reduce signal loss between land and sea. The cubic-spine approach was used to interpolate the missing monthly data. The time series and secular trends of TWS and groundwater storage variations (GWS) were constructed using the two potential solutions, and the TWS trend was produced. Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 937402 6 GLDAS, TRMM, and SRTM Data GLDAS is an observational data assimilation system with satellite and ground-based components. It produces optimal fields of land surface states and fluxes using advanced land surface modelling and data assimilation techniques. For soil moisture storage (ΔSMS), outputs of GLDAS (Rodell et al., 2004; https://disc. gsfc.nasa.gov/datasets) are employed. In this work, the averaging of three GLDAS versions (CLSM, VIC and NOAH) Mascon products CSR-M, and JPL-M; GWS: groundwater storage change; SMS: soil moisture change; Annual Average Rainfall is abbreviated as AAR. Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 937402 7 was used. To calculate the ΔGWS, the ΔSMS changes are subtracted from the ΔTWS using Eq 1

ΔTWS ΔGWS+ΔSMS
(1) The Tropical Rainfall Measuring Mission (TRMM; https:// disc.gsfc.nasa.gov/datasets) is a collaborative space mission designed to monitor and analyze tropical rainfall (Kummerow, 1998). TRMM and Shuttle radar topography mission (SRTM) data are used. The monthly TRMM data, covering the January 2002-December 2019 period with a spatial resolution of 0.25°× 0. 25°, are used to create the average annual rainfall (AAR) and to generate the monthly rainfall time series and the AAR rate over the study area. A Digital Elevation Model (DEM) was created using SRTM data with a 90 m resolution that was used for the delineation of the stream networks in the area. Landsat eight datasets were used to extract the surface lineaments affecting the area.

Aeromagnetic Data
The total magnetic intensity (TMI) map used in the current study was provided from the Saudi Geological Survey. The TMI map was subjected to different filtering techniques such as the reduction to the pole, tilt derivative, and analytical signal. These filters were used to enhance the image's signal intensity, evaluate the magnetic anomalies qualitatively, and calculate the depth to the magnetic sources. The Geosoft program (Geosoft oasis montaje, V.8.2.4, 2015) was used to handle and analyze the aeromagnetic data. Figure 2 shows the reduced-to-pole (RTP) magnetic map.
The RTP filter is used to eliminate the skewness of magnetic anomalies induced by non-vertical magnetization directions. The RTP anomalies' maxima are directly over the causal bodies, making explanations easier (Lu et al., 2003). The TMI data were corrected for the RTP in this investigation using the Oasis Montaj Fourier transformation approach (Geosoft oasis montaje, V.8.2.4, 2015). The parameters of the International Geomagnetic Reference Field (IGRF) are as follows: The total magnetic field intensity was 39,582 nT, and the inclination and declination were 25.5 and 1.6, respectively. These characteristics were utilized to produce the RTP map (Figure 2) from the TMI field.
The magnetic anomalies' edges and locations were sharpened using the Tilt angle (TDR) derivative (Mushayandebvu et al., 2001). It is applicable in mapping mineral investigations and shallow basement structures (Geosoft oasis montaje, V.8.2.4, 2015). It was defined by Miller and Singh (1994), as the ratio of a vertical to a combined horizontal derivative as in Eq.2: TDR tan −1 [(zf/zz)/sqrt zf/zx 2 + zf/zy 2 (2) The potential field is represented by f, and the vertical derivative of the field is (›f/›z). The tilt angle is positive when above the source, at the edge where the vertical derivative is zero and the horizontal derivative is a maximum, and negative when beyond the source region. The tilt angle can be used to trace the pattern of the edges because it produces a zero value over the source edges (Miller and Singh, 1994). The tilt amplitudes are in the range of -π/2 to +π/2.
To detect the borders of the magnetic structures, the total horizontal derivative (THDR) technique was used extensively (Pilkington and Keating, 2004;Cowan et al., 2005;Cooper, 2009). It calculates the potential field's rate of change in horizontal directions (Cordell, 1979). The faulted boundaries may have high gradient zones, which can be identified using horizontal derivative maps. The horizontal gradient approach has the advantage of being less vulnerable to data noise because it only involves the computation of the field's two first-order horizontal derivatives (Phillips, 1998). The horizontal gradient's amplitude (Cordell and Grauch, 1985) is represented as : The horizontal derivatives of the field in the x and y directions are (zf/zx) and (zf/zy), respectively. When applied to field observations, TDR and THDR (Figure 3) both operate as an automatic gain control filter.
To detect the magnetic and gravity source edges, Pham et al. (2020) designed a new edge detection filter. The real part (R) of an arcsine function (asin) of the ratio of the vertical gradient to the total gradient of the amplitude of the field's horizontal gradient is used to create the new Enhanced Horizontal Gradient Amplitude (EHGA) filter. They tested its usefulness in defining source edges using synthetic and real data, and compared the results to those produced by other methods. Edges produced by the filter are more accurate and have a higher resolution. Eq 4 gives the EGHA edge detector: Where HGA is the potential field's entire horizontal gradient amplitude. HGA x , HGA y, and HGA z are the x, y, and z gradients of the HGA. K is a positive number decided by the interpreter.

Geoelectrical Data
On the local scale, the use of geophysical techniques for groundwater investigation has increased considerably during Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 937402 8 the previous decade due to the high advances and developments in electronic technology and numerical modeling solutions (Metwaly et al., 2009;Ndlovu et al., 2010). The electrical resistivity is a common method for groundwater exploration due to its low cost, simple operation, and efficiency.
Electrical resistivities of subsurface materials can be measured using geoelectrical techniques. These measurements provide information on underlying geoelectrical layers, structures, and groundwater occurrence (Van Overmeeren, 1989;Dahlin et al., 1999a;Muchingami et al., 2012). Archie's Law (Archie, 1942) describes the occurrence of fluid in the fluid-bearing formations and suggests that the resistance decreases with increased water content according to an inverse power law as follows: Where ρ o is the fully saturated rock's resistivity, ρ f is the water filling the pores' resistivity, φ is the porosity, and a and m are empirical parameters (Keller and Frischknecht, 1966).
The cementation exponent is defined by the m parameter, which has values ranging from one to 5. (Glover, 2010). Regarding the soil moisture content, Grellier et al. (2005) demonstrated that Archie's law can be reduced to Eq 6: Where θ w is the medium's gravimetric moisture content, and α is a quantity that includes the medium's bulk density, as in Eq 7: Passing an electrical current (I) between two metallic electrodes and measuring the potential difference (V) between them is the electrical resistivity method. The apparent resistivity is calculated using the following equation (ρ a , Eq 8) (Dahlin, 2001): ρ a KΔV/I K, a geometrical factor, depends only on the position of the electrodes. Twelve vertical electrical resistivity soundings (VESes), with current electrode (AB) separation, ranging between 200 and 500 m, were conducted along eight profiles ( Figure 2) using a Schlumberger array. The apparent resistivity value depends on the geometrical configuration of the electrode array, as defined by the geometric factor.
The IP2WIN (2005) application was used to process the VESes. The experimental data curve has been divided into smaller curves. Each one is a geoelectrical unit with a known resistivity Ω.m) and thickness (m), and it could represent a geologic layer with different physical attributes than the layers above and below it.

RESULTS
The results of the processed TRMM data are shown in Figure 4. Analyses of the rainfall data show a general positive AAR trend over the study area estimated at 117.6 mm/yr. The results of the spatial distribution of the secular trend in GRACE-derived ΔTWS data and the monthly ΔTWS time series over the study area are shown in Figures 5, 6. Both the ΔTWS and ΔGWS from the two mascon datasets show a negative trend estimated at −0.34 ± 0.01 cm/yr over the entire period. The RTP, TDR, THDR, and EHGA maps are shown in Figures 2, 3. These maps show that the study area is affected by subsurface structural trends in the directions of NS, NNW, and NNE. The ground surface relief is forming streams (Figure 8) taking the surface water away Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 937402 11 towards the E, N, and W. The surface structures were extracted from the Landsat images (Figure 9). Surface and subsurface main trends in the direction of NNW are delineated from the Landsat images ( Figure 9) and the RTP map (Figure 10), respectively. The subsurface geology is represented by two layers (Figure 11) of different magnetic susceptibility values reflecting the top sedimentary cover and the lower basement crystalline rocks. The inversion of the resistivity data resulted in three geoelectrical layers (Figure 12) with varied physical attributes. Higher resistant unconsolidated Quaternary sediments make up the first geolectrical layer, low resistant sands saturated with water make up the second layer, and low to high resistant fractuted and/ or massive basement rocks make up the third one ( Table 2).

DISCUSSION
In this work, we used a combined method using GRACE satellite data with other airborne and ground-based geophysical data to investigate the groundwater potentialities of the southwestern part of Saudi Arabia. The monthly rainfall rate ( Figure 4A) is higher in the March to September months and lower in the October to February months. The coastline region had a higher AAR rate of about 300 mm, whilst the eastern parts had lower values of up to 50 mm ( Figure 4B). The AAR time series ( Figure 4C) shows a generally positive trend with the AAR of 117.6 mm throughout the period 2002-2019. The AAR rate is one of the most important characteristics of a region's climate, as well as its associated ecosystem types and habitats. The amount of rainfall a region receives has an impact on stream density, water availability, land cover type, and agricultural output. Based on the linear regression analysis of the AAR; three time periods are distinguished ( Figure 4C). Period I (2002)(2003)(2004)(2005)(2006) shows an AAR of 124.8 mm/yr; Period II (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) shows a minimal AAR of 97.4 mm/yr; Period III (2016-2019) shows the highest AAR of 161.3 mm/yr. Period II is consistent with the onset of the 2007 drought and the dry climatic conditions that affected the Middle East area (Trigo et al., 2010). As a result, the region has suffered from water scarcity and limited water resources since the beginning of that drought (Wolf and Newton 2007;IRIN FIGURE 12 | The VES inverted models. Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 937402 12 2010; Voss et al., 2013;Mohamed, 2020b). However, the higher AAR rate during Period III indicates that the study area is receiving a substantial precipitation rate.
The research area had an overall negative TWS trend over the entire period (04/2002-12/2021), as shown in Figure 5. The southwestern part of the research area has a little positive TWS rate (+0.08 cm/yr), but the northeastern part has a negative TWS rate (−0.5 cm/yr). Figure 6 shows the monthly time series of TWS fluctuations from the two mascon products, as well as their mean. Witnessing of this figure shows a general depletion in ΔTWS rate estimated at −0.30 ± 0.02, −0.38 ± 0.01, and −0.34 ± 0.01 cm/yr from the CSR-M, JPL-M, and their mean, respectively during the entire period. The ΔTWS time series has three trends throughout the entire study period based on the linear regression analysis of the average (Figure 6). Period I shows a slightly negative trend, calculated at −0.09 ± 0.09 cm/yr ( Table 1). A negative TWS Δtrend is estimated at −0.31 ± 0.03 cm/yr for Period II, whereas Period III shows a highly negative trend, calculated at −0.67 ± 0.09 cm/yr. The GLDAS-derived ΔSMS for Periods I, II, III, and the entire period are +0.12 ± 0.19, +0.04 ± 0.04, −0.18 ± 0.14 and +0.03 ± 0.02 mm/yr, respectively assuming that large crystalline rocks occupy the majority of the surface area. To estimate the fluctuations in the ΔGWS (Figure 7), the non-groundwater component, represented by the ΔSMS, was subtracted from the ΔTWS. During Periods I, II, III, and the entire period, average depletion rates of −0.10 ± 0.08, −0.32 ± 0.03, −0.66 ± 0.09 and −0.34 ± 0.01 cm/yr were obtained for the ΔGWS trend values, respectively. The ΔTWS and ΔGWS trends both follow the same pattern.
Inspection of Figure 4C indicates that there is a general increase in the AAP rate during Period III; however, the highest ΔGWS depletion rate (Figure 7) is observed in that Period. This is largely attributed to the drainage pattern and the surface material of the study area. Inspection of Figure 8 shows the stream networks overlay the DEM. It shows higher relief of about 2000-2900 m in the Mountainous area close to the Red Sea coastal area. It decreases north and eastward to low values of about 1,000 m in its eastern part. It shows a steep slope toward Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 937402 13 the west at the coastal area. Surface water that may be formed from the substantial rainfall is being drained away from these stream networks toward the east and north. The rest of the surface water drains into the Red Sea's coastline area. This may explain the negative trend in ΔGWS over the study area. Moreover, the eastern and central parts of the area are occupied by massive crystalline rocks that help the surface water drain away north and eastwards. Part of the surface water is drained toward the coastal area which shows a positive TWS value ( Figure 5).
Witnessing Figure 2 shows that the magnetic intensity values vary from −1,377 to +524 nT. The lithological changes, basement relief, and/or faulting and folding may be the cause of the variations in the magnetic intensity values. Lower and higher magnetic anomalies are almost arranged in N-S and NNW directions. The enhanced edges of linear features in the TDR map ( Figure 3A) could be attributable to deep and shallow magnetic sources (Miller and Singh, 1994). The most important advantage of the TDR is its zero contour line that suggests the locations of the subsurface structures caused by sharp changes. The magnetic anomalies are practically aligned in definite orientations, generating structural trends in the NS, NNW, and NNE directions, as seen in the TDR map ( Figure  3A) and its zero contour line. The edges of the maxima are welldefined in the THDR map ( Figure 3B), assuming that amplitude enhancement is preserved. It calculates the potential field's rate of change in horizontal directions. The structural trend in the NS direction can be detected from the magnetic anomalies that are oriented in a slightly NS direction, according to the THDR ( Figure 3B). The EHGA map ( Figure 3C) reveals that the structural trends are almost in the NS direction, with higher resolution and sharper responses over the edges of the study area's magnetic sources.
The surface geologic features (lineaments) of the study area were mapped using Landsat eight images ( Figure 9A). The extraction of these lineaments from remote sensing data can be done through automatic extraction (Anwar et al., 2013;Rayan, 2013). Several computer-assisted lineament extraction methods have been proposed. The majority of the methods rely on edge filtering. The PCI Geomatica's LINE module is the most extensively used software for automatic lineament extraction. The predominant trend derived from the analysis of geological lineaments is in the NNW direction ( Figure 9B). The subsurface trends of fault lines were also analyzed using the trend analysis technique according to their lengths, abundance, and magnitude concerning the azimuth. The magnetic structures have a significant relation with the intensity, direction, and geometry of the magnetic anomaly trends (Hall, 1964). Subsurface tectonics are delineated using magnetic trend patterns (Affleck, 1963). The trend analysis technique was applied to delineate the subsurface trends controlling the groundwater flow pathways in the area. The trend analysis shows that the main structural trends are in NNW, NNE, and NE directions ( Figure 10). Minor trends are represented in the NS direction. The lineaments could also be interpreted as the edges of geological bodies and directions of structures. The integration of the delineated surface and subsurface structural trends may indicate the continuity of the NNW trend within the basement rocks and the overlying sedimentary deposits. This leads to a connection between the shallow Quaternary and the deeply fractured basement rocks aquifers that may feed the deeper aquifer. The subsurface structural trend and related fractures, on the other hand, are draining groundwater and recharging water away to the north, west, and east. This causes also a depletion in the groundwater storage of the study area, as estimated from the GRACE data at −0.34 ± 0.01 mm/yr during the period 04/2002-12/2021. However, groundwater can be accumulated at the intersection zones of the faults.
The underlying geology was investigated and the thickness of the sediments overlying the crystalline basement rocks was estimated using 2D modeling of magnetic data along profiles. To create 2D models from magnetic data, we used the Geosoft GM-SYS software package. We can optimize the model using the GM-SYS inversion option. Modeling programs require adequate initial estimations of model parameters such as body form, susceptibility, topography, depth, and magnetization of suspected sources to minimise errors of non-unique solution between observed and calculated magnetic fields. The profiles' models are illustrated in Figure 11. The geological succession is mostly represented by two rock types, according to these models: (a) the upper sedimentary layer that is composed of Quaternary and/or Paleozoic sediments, and (b) the lower Pre-Cambrian crystalline rocks that are composed of granitic and metamorphic rocks. A two-layer model was applied in the modeling, given that the Quaternary cover is the non-magnetized layer, and the basement is the magnetized layer. The magnetic susceptibility values of 0.00001-0.0005 CGS and 0.001-0.038 CGS units were chosen for the non-magnetized sedimentary and the more magnetized basement rocks, respectively. The thickness of the Quaternary sediments shows wide variations, it varies from 0 tõ 1.2 km in the study area. The fewer thickness of the sedimentary cover in that area is suggesting the presence of lower groundwater potential in the upper sedimentary aquifer.
Three geoelectrical units ( Figure 12; Table 2) of known resistivity and thickness values are recognized, which may reflect geologic layers of physical properties that differ from those located above and below. The details of the geoelectrical units are described in the following: The surface geoelectrical layer has different resistivity and thickness values varying from one place to another. The resistivity values of the first layer are varying from 428 Ω m at V3 (Race 4) to 9626 Ω m at V4 (Najran 37), whereas the thickness values are varying from 1 m at VES V1 (bish-4) and VES Vt (Najran 10) to about 94 m at VES 3 (Najran 48). The second geoelectrical layer represents the water-bearing unit in the study area. This layer has resistivity values varying from 5.1 Ω m at V1 (bish-4) to 153 Ω m at V4 (Najran 37), and thickness values varying from 8 m at V7 (Edaby-2) to 107 m at V1 (bish-5). The third geoelectrical layer is represented by low to high resistant basement rocks with resitivity values varying from 20 Ω m at V1 (P1, bish-5) and 30 Ω m at V2 (P1, bish-4) to 6125 Ω m at Vt (P3, Najran 48) and 6030 Ω m at V3 (P6, Najran Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 937402 15 9). The lower resistivity values for that layer at particular VES points could result from the presence of fractured basement rocks saturated with water that lowers the resetivity values. In contrast, the higher resistivity values for that layer at certain VES points could indicate the presence of more resistant massive basement rocks. Our current VES interpretation and results are in good agreement with the results of resistivity data at Wadi Sar , which is located in the northwestern part of the research area. Piezometric wells should be drilled in the wadis of the western portion of the area in the future to test and verify the study's findings, however, this is dependent on budget availability.

CONCLUSION
For the study region, an integrated strategy combining gravity, magnetic, and electrical data with other remote sensing datasets was used. The study area is receiving an average rainfall rate of 117.6 mm/yr forming surface stream networks. Part of the surface water is drained to the Red Sea by the streams, while the rest is drained to the eastern and northern parts of the land. The groundwater shows a depletion trend estimated at −0.34 ± 0.01 mm/yr during the entire period 04/2002-12/2021. The surface trends are in NNW, and the deep trends are in NS, NNW, and NNE directions according to an examination of the structural trends retrieved from Landsat images and magnetic data. The thickness of the sediments varies greatly, ranging from a low of 0 in the mountainous areas to a high of 1.2 km in the valleys. The water-bearing unit has a thickness ranging from 8 to 107 m with low resistivity values of 5.1-153 Ω m. The depth of that unit fluctuates between 1 and 94 m. The average yearly rainfall over the research region is ranging between~300 mm in the coastal area and up to 50 mm in its eastern parts. The eastern part, as well as the lowlands of the wadies in the west, could be potential areas for agricultural growth and the construction of new villages.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.