Multi-criteria analysis using satellite images and GIS for mapping groundwater recharge zones in Saudi Arabia: A case study of the Riyadh Region

In arid regions, such as the Arabian Peninsula, precipitation is often lower than 200 mm/year while the potential evapotranspiration exceeds 2,000 mm (Baban 2022). The precipitation patterns in the entire Arabian Peninsula have become torrential, and huge amounts of rainwater precipitate in a short time period evidencing the anticipated climatic variability in the region (IPCC 2007). The intensive patterns of precipitated water resulted in increasing the number of flash floods in several parts of Saudi Arabia and this is often associated with severe damage to the infrastructure and the environment (Al Saud 2015a).

The investment of groundwater in the Kingdom of Saudi Arabia is mainly from deep aquifers (i.e., fossil water), which are usually at a depth exceeding 1,200 m reaching the Nubian Sandstone rock formation, Minjur and Saq Sandstones aquifers; and ranging between 300 and 500 m in the carbonate rocks (i.e., Wasia-Biyadh and Wajid rock, Umm Er Radhuma, and Dammam Limestone formations) (Al Saud 2010). The lowering in the piezometric level and the decreased discharge rate are well pronounced. For instance, the depletion of groundwater in agricultural lands of the northern part of Saudi Arabia (e.g., Ha'al, Al-Jouf, and Tabouk) was calculated between −6.9 × 10⁻² and −8.6 × 10⁻² cm/month as it was investigated from GRACE TWS (Wehbe 2022).

The deficit in water supply in Saudi Arabia exists in both the urban and rural areas. Besides, the water availability of less than 70 m³/capita/year, there is an estimated water demand exceeding 700 m³/capita/year, which is equivalent to 10 times (UNESCO 2021). With such exceeded water demand, the state-of-the-art collection of rainwater for non-potable or irrigation uses is now being adopted in many urban areas (Leggett & Shaffer 2002; Peters 2006; Mahmoud & Alazba 2016). This is also associated with the groundwater recharge (GWR) in valley systems in order to reduce the velocity of streamflow and the related flood occurrence. It has become a known geo-environmental problem in urban areas in Saudi Arabia which are witnessing serious water deficit besides recurrent occurrence of floods and torrents. The Riyadh Region is a typical example where the precipitation rate does not exceed 150 mm/year, while renewable freshwater is few enough (i.e., less than 50 mm/year). This region has been struck lately by a series of damaging floods. It is a paradox that the region is suffering from a water deficit while flood water immerses large geographic patches. In this respect, studies to calculate hydrological parameters are still few to diagnose water flow/storage regime as well as to calculate the water budget, and there are only studies on flood assessment, vulnerability, and management issues (Rahman et al. 2016; Sharif et al. 2016; Ledraa & Al-Ghamdi 2019).

This study investigates a principal hydrologic parameter that interlinks between the movement of surface water, and the feeding mechanism of groundwater reservoirs. Thus, surface water percolation into subsurface rock layers must be well identified, and this can be illustrated in a map form. The produced GWRP map can assist in identifying the most suitable localities to adopt artificial GWR; therefore, flood water can be invested in this man-made hydrologic process instead of immersing urban and agricultural lands and causing damages and unfavourable environmental impact to the state-of-the-art in the Riyadh Region.

In order to produce a GWRP map, geospatial data will be retrieved from the available sources mainly from the thematic maps, and from satellite images with different temporal and spectral resolutions, in particular Landsat 7 ETM⁺, Spot-5, and Aster images which are characterized by spatial resolution that enables detecting a miscellany of hydrologic (and hydro-geologic) features which are not feasible to be identified on-ground. This must be integrated with the applications of the geographic information system (GIS) which helps in digital data production, visualization, analysis, and overlapping. The novelty of the applied methodology in this study includes the manipulation of factors of impact, and this was represented by the adoption of weighting and rating for each factor on the water flow process from surface to subsurface rock layers. Also, it must be made clear that this study aims to identify potential lands for groundwater feeding, which can be used for positioning suitable sites for groundwater artificial recharge; however, these sites are not necessary to identify lands with potential groundwater storage.

There are several studies applied to assess groundwater potential zones (Shaban et al. 2006; Nolan et al. 2007; Yeh et al. 2008; Boughariou et al. 2015; Souissi et al. 2018; Arshad et al. 2020). These studies followed different tools and methodologies; and therefore, the results were almost contradictory.

The Riyadh Region has a characteristically continental desert climate, dominated by hot summers and relatively cold winters with low precipitation. Average temperatures in the Riyadh Region are 25 °C (77°F), ranging from more than 50 °C (122°F) in summer, to 0 °C (32°F) in winter. Humidity is 33% with annual precipitation of 85 mm. In addition, the region is prone to occasional sandstorms.

The topography of the Riyadh Region is almost mild and includes gently sloping plateaus with slope gradients not exceeding 12 m/km%. It occupies 30 major drainage systems with a surface area ranging between 5,590 and 70,477 km² for the existing watersheds. The geology of the Riyadh Region shows that the exposed rocks are from the Precambrian age where metamorphic, igneous, and sedimentary rocks exist with a variety of rock types including largely granite, schist, sandstone, and limestone with volcanic eruptions. This has been well reflected in its hydrogeology where groundwater is mainly stored in deep sandstone aquifers and few aquiferous layers are located in the carbonate rocks.

The Riyadh Region occupies several agricultural lands which are irrigated either from the deep boreholes in the non-renewable sandstone aquifers or from flood water along valleys and collected in reservoirs behind constructed dams. Sometimes surface water naturally accumulates in lowlands and depressions where it lasts for a couple of months before the largest part of it is evaporated and the rest percolates into the underlying rocks.

This study aims to generate a digital cartography (i.e., mapping) for the GWRP zones in the Riyadh Region. The GWRP map represents a helpful tool and can be used as a thematic layer for further applications, notably, it shows the ability of surface water to percolate into the underlying rocks. For instance, it can help stakeholders and decision-makers to recognize suitable areas for groundwater artificial recharge, areas vulnerable to groundwater pollution, as well as it can be used for agricultural planning, notably those related to irrigation and water efficiency. The production of the GWRP map will significantly depend on the use of satellite images for geospatial data extraction as well as requires use of the GIS for manipulation of this data and then their cartography in a digital form.

Seepage of surface water, whether from runoff or from overland flow, into the beneath rocks is a significant hydrologic process, and it is governed by a number of factors that differ between regions according to their geographic setting and relative degree of impact. These factors have been illustrated by many researchers who also adopted different methods for analysis to perform GWRP zones (e.g.,Shaban et al. 2006; Yeh et al. 2008; Boughariou et al. 2015; Souissi et al. 2018; Abdekareem et al. 2023). However, the majority of these factors include mainly the characteristics of materials located on the terrain surface and also the substratum materials where water can seep from. In this study, five factors were adopted for mapping GWRP, and the selection of these factors was based on the natural setting of the Riyadh Region. These factors are rainfall distribution, lithological characteristics of the exposed rocks, density of rock fractures, density of streams, and slope. There are some other factors that might be used but with less impact on the recharge process. For example, soil characteristics can be a factor; however, it was not included because of the shallow thickness (i.e., tens of centimeters) of soil in the study region which does not impact surface water seeps, as well as due to its deformed structure which is dominated mainly by cracks and piping due to the high evaporation rate. This was also considered for the other factors such as land use/cover which was not also included with the manipulated factors, because the surficial components of the land cover/use in the region are often changing whether due to rapid human interventions and activities or due to the seasonal and even monthly diversity in weather and the resulted impact on the cover and use of lands.

Therefore, these factors were primarily converted into digital data forms and then manipulated in order to harmonize their impact by a multicriteria systematic analysis.

Generally, areas with high rainfall rates are expected to capture a considerable volume of precipitated water where a part of this water infiltrates the underlying rocks, while dry areas have the opposite. However, this factor is anomalous and is controlled by topographic characteristics, because when water from rainfall reaches the terrain surface, it may divert into different areas (i.e., areas with the least rainfall rate) due to the topography.

Normally, rainfall is significant for the natural replenishment of surface water bodies (e.g., streams, lakes, etc.) and for groundwater reservoirs (i.e., aquifers). The sustainability of water supply, whether from surface or groundwater source depends on the volume of the precipitated meteoric water. However, the concept of rainfall water recharge into groundwater remains undefined between different studies, and this is also the case in Saudi Arabia. This gap in knowledge is due to the unavailability of continuous data series and the non-uniformity in the distribution of climatic stations for evaluating surface water storage versus groundwater level.

For GWR, studies based on the piezometric observations and stable-isotope ratio indicated that rainfall largely contributes to GWR (Owor et al. 2009; Jasechko & Taylor 2015), but the role that is played by the phreatic zone in transmitting water often remains unclear (e.g., infiltration rate and the hydrologic controls). Most studies that quantify recharge rate from groundwater-level fluctuations are constrained by uncertainty in unmeasured coefficients (e.g., specific yield) of groundwater storage (Sibanda et al. 2009; Obuobie et al. 2012).

The significant contribution of rainfall in recharging groundwater is well pronounced since wet areas are mainly found in humid zones and dry areas exist in arid zones evidencing rainfall distribution and rate. For this reason, rainfall and its geographic distribution must be accounted while studying GWR (Wang et al. 2015; Amin et al. 2016; Kotchoni et al. 2019).

The characteristic of the exposed rocks is considered one of the most significant factors controlling the mechanism of water seeps (i.e., recharge) from the surface to the substratum. Hence, surface water easily flows in permeable lithologies, while it is retarded to flow in compacted rocks. Therefore, GWR is affected by the exposed rocks on the terrain surface. In this regard, there are many investigations considered lithology in GWR assessment and mapping (e.g.,Shaban et al. 2006, 2020; Nag & Ghosh 2013; Senanayake et al. 2016), while investigators replaced this factor with the relevant hydraulic characteristics of rocks, when they are available, such as porosity and permeability (El-Baz & Himida 1995).

In studying recharge potential zones, only the exposed rock materials are considered while the underlying rocks at depth are not. This is because determining GWRP describes the mechanism and rate of surface water percolation, and not groundwater flow/storage, where the latter requires integrating other factors, such as the stratigraphic sequence of rocks and the related geological structures (e.g., the inclination of rock beds).

Fracture systems are directly associated with the lithological characteristics, and rock deformations (i.e., fractures) and they represent a significant factor in enhancing rock permeability. Normally, the higher degree of rock fractures allows the percolation of large volumes of surface water to the underlying rocks; especially when these fractures occur as large-scale elongated features; therefore, water flows rapidly along these fractures. According to Kumar et al. (2007), Selvam et al. (2014) and Al Saud (2015b); in the hard rock areas, fractures represent zones that result in increased secondary porosity and permeability, and thus they compose rock domains saturated with water.

There are various forms of rock fractures including faults, fissures, joints, and cleavage, and each of them has its own hydraulic property and acts in transmitting water at different levels. For example, fissures and joints enable water to percolate locally within their fracture domains, while faults can transmit water even for hundreds of kilometers. These fracture forms are often presumed as potential groundwater zones (Al Saud 2010). Hence, higher fracture density is always accounted for while studying GWRP zones.

Streams are important geomorphological and hydrological features on the terrain surface, and they are always used as indicative features while studying water seeps from surface to subsurface rocks. The patterns and density of these surficial elements are fundamentally controlled by the underlying lithology; and hence, they affect water percolation rate. Therefore, streams are always involved in GWR studies, but the majority of their consideration includes their density (i.e., the total length of streams within a defined area). It is well known that the denser the streams the lower the recharge rate and vice versa. Thus, lower stream density is an indication of a higher recharge rate.

Slope, in many instances, is not involved while studying GWR, especially in areas with flat surface terrain. However, this factor must be included, since along steep surfaces rainfall flows away instead of infiltrating into the ground. This factor governs the velocity of surface water flow. Hence, in steep slopes, the velocity of water flow is high, and it does not allow sufficient time for water to infiltrate into the phreatic zone, but this is not the case in gentle slopes where water slowly flows, and this allows sufficient time for water to percolate. Therefore, a steep slope enhances surface water flow and vice versa (Al Saud 2010).

Many types of materials were used in this study in order to prepare the geospatial datasets for the influencing factors on GWRP, and these materials were useful for data analysis in systematic approaches. This was also associated with field surveys to ensure the reliability of the retrieved data from satellite images. The majority of this study has been initiated from a technical work performed by the author to produce the Digital Geological Map for the Riyadh Region, with a 1:1,000,000 scale (Al Saud 2015b & Available at: www.rcrc.gov.sa/en/research/). Therefore, geospatial datasets were the main required materials for analysis and manipulation approaches which can be summarized as follows:

• 1.

  Thematic maps:

• • –

    Geological maps for Saudi Arabia, 1:50,000 scale (US Geological Survey 1963).

  • –

    Digital Geological map for the Riyadh Region, 1:1,000,000 scale (Al Saud 2015b).

  • –

    Topographic maps for Saudi Arabia, 1:50,000 scale and 25-m contour interval (Ministry of Petroleum and Mineral Resources 1983).

• 2.

  Remote sensing products:

• • –

    SRTM DEM (Shuttle Radar Topography Mission)

  • –

    Landsat 7 ETM⁺ images with 30-m spatial resolution (Dates: 11 April 2018, and 3 July 2018 and selected scenes for the following Paths and Rows: 163-45, 163-44, 163-43, 163-42; 164-45, 164-44, 164-43, 164-42; 165-45, 165-44, 165-43; 166-45, 165-44, 165-43; 166-45, 166-44, 166-43; 167-44, 167-43)

  • –

    Aster (The Advanced Space-borne Thermal Emission and Reflection Radiometer) images, with 15-m spatial resolution, product: Aster L1A Reconstructed Unprocessed Instrument Data V003.

  • –

    SPOT-5 images with 2.5-m spatial resolution (Dates: 14 December 2015 and selected scenes for the following paths and Rows: 164-45, 163-44, 163-43, 163-42; 164-44, 164-43, 164-42, 165-45, 165-44, 165-43; 166-44, 165-43, 165-42; 166-44, 166-43)

  • –

    Rainfall map (1999–2019), scale of 1:1,000,000 extracted from TRMM (Tropical Rainfall Mapping Mission) (TRMM 2015), and it was also supported by Basahel & Mitri (2018).

• 3.

  Software:

• • –

    ERDAS-Imagine 11 for image processing (Leica Product).

  • –

    Arc-GIS 10.2 for the manipulation of geospatial data (ESRI Product).

All thematic maps were prepared for analysis after conversion into digital forms for further data manipulation. Satellite images were processed using mainly ERDAS-Imagine 11. The processing of satellite images can be summarized as follows:

• 1.

  Pre-processing of satellite images: This includes a series of digital operations, including mainly atmospheric correction or normalization, image registration, and geometric and radiometric corrections. In addition, ‘images subletting and mosaicking’ was applied to cut the required scenes from large images and to link these cut images together in one scene; respectively. Consequently, the selected images were primarily saved in tiff Format and the linking approach was obtained using the ERDAS-Imagine software.

• 2.

  Band combination: It is to discriminate terrain features using bands VNIR (15 m with three bands): bands ordered as 3-2-1; SWIR (30 m with six bands): bands ordered as 7-9-5, 5-7-6 and 8-4-7; TIR (90 m with five bands): bands ordered as 14-12-14.

• 3.

  Images enhancement: There are a number of image enhancement tools available in the ERDAS-Imagine, such as filtering, and different levels of linear stretching, Gaussian, equalization and square root.

• 4.

  Interactive stretching: This is another advantage of ERDAS-Imagine identifying specific features on the satellite images using histogram regulators to apply (for example) linear, Gaussian, linear, Arbitrary and User-defined LUT.

• 5.

  Density slice coloring and contouring: These methods enable classifying the color of DNs (digital numbers) that are within the same range.

Thermal interpretation from thermal bands, i.e., bands 10–14 (90 m resolution), was also applied, especially for fracture detection.

For slope generation, SRTM DEM was used as a tool for the Digital Elevation Model (DEM) which is achieved from Triangulated Irregulated Networks (TINs). Therefore, SRTM DEM was processed by Arc-GIS software, and specifically in the Toolbox system. The navigation was through the Toolbox-then Spatial Analyst Tools-Surface-Slope.

In order to reach the main objective of this study in elaborating GWR potential map, the integration of the influencing factors in a systematic approach will be the primary step. This requires conversion of all materials and datasets of the influencing factors (e.g., maps or satellite data) into digital form (i.e., GIS datasets). Then, each GIS dataset of the influencing factors will be considered a ‘GIS layer’. In this study, each GIS layer was classified into five classes, where each class represents a degree of impact on GWR property. These classes were categorized between very high and very low impact on GWR. For example, dense fracture systems will be given a very high impact, because it enhances water infiltration, while compacted rock will be given a very low level of impact because it retards the infiltration rate.

The produced rainfall map shows eight classes for rainfall rates ranging between 60 and 130 mm/year, where the largest part of the Riyadh Region is characterized by rainfall rate below 100 mm/year, while the average annual rainfall rate is calculated at about 91 mm (Figure 2 and Table 1).

Mapping the lithologies of the Riyadh Region includes the use of the available geological map for this region which was done by the author (Al Saud 2015b). This map was mainly obtained from the previous geological maps (i.e., U.S. Geological Survey 1963), where all lithologies were digitally traced from these maps. However, satellite images were also supportive in this respect. Thus, ASTER images (pixel size 15 m) and SPOT-5 (pixel size 2.5 m) images were processed using ERDAS-Imagine software, which has a number of optical and digital advantages.

Each satellite image has a number of spectral bands. For example, the ASTER image has 14 bands, while Spot-5 encompasses five bands. ERDAS-Imagine is capable of arranging each of these bands separately in order to have the most suitable and distinguished observation of the images. In addition, many other digital features in the ERDAS-Imagine software were used, such as filtering, color slicing, contrasting; thus visual tracing of many lithological units and their boundaries which were not clearly identified in the available geological maps, and was also helpful for determining various geologic structures.

The exposed rocks in the Riyadh Region are from the Precambrian, which contains ultramafic rocks with basement complex, to Quaternary deposits with dominant alluvium deposits, sand dunes and reef limestone. There are 65 rock units in the Riyadh Region belonging to 23 rock formations. These rock units were classified into five classes to integrate rock formations with similar hydrologic properties in regard to the ability to transmit water from the terrain surface into the beneath rocks. These lithologies and their areas are shown in Table 2 and Figure 3.

In this study, the extraction of lineaments was achieved by using Landsat 7 ETM⁺ images, and more certainly by analysing the thermal band (i.e., band No. 6) which has 120 m × 120 m spatial resolution. The processing of these images, using ERDAS-Imagine software, has the advantage of applying ‘edge detection’ which is a digital criterion to enable the recognition faults along wet horizons. Thus, thermal bands can determine the alignments of these lineaments.

The produced lineaments' density map was classified into six classes describing it from very high to very low lineaments' density. Such a classification depended on the maximum and minimum number of lineaments in the area of study, and it was found that low and very low lineaments' density is the most dominant in the Riyadh Region (Table 3).

The number of watercourses in an area evidences the recharge rate of the terrain surface. For the Riyadh Region, streams were drawn using the Global Digital Elevation Model (GDEM) with 30 m spatial resolution; therefore, streams up to fourth order were extracted. This was associated with the use of topographic maps in order to confirm the delineation of major streams as well as to identify their local names. For this work, Arc-View (a licensed level of Arc-GIS) was used for streams' extraction and then the representation of streams in a map showing for stream' density.

Therefore, a similar method (i.e., sliding windows) which was used to create a lineaments' density map was also applied to generate streams density map, and five streams' density classes were adopted (Table 4), starting from very high to very low streams' density. It is obvious that very low streams' density is the most dominant (about 54%) in the Riyadh Region (Figure 6).

Slope mapping for the Riyadh Region was achieved using a radar-based remote sensing product, the SRTM (Shuttle Radar Topography Mission) which is a useful tool in this regard. SRTM data are made publicly available by NASA, with a three-arc-second pixel size (1/1,200 of a degree of latitude and longitude), which has 90 m pixels' spatial resolution. Therefore, Arc-View was capable of generating slopes based on determining the flow direction of streams in each pixel.

The performance of the GWRP map for the Riyadh Region was produced following a number of criteria including the determination and classifications of the influencing factors on this hydrological process, and then these factors might be systematically manipulated in the GIS system. The obtained results were based on a number of methods, whether for geospatial data extraction and dataset manipulation, to prepare the thematic data and numeric values into digital files in the GIS system. In addition, the assignment of the degrees of influence (i.e., weights and rates) followed empirical calculations applied to reach the most accurate results.

The selection of factors for the systematic manipulation was given a concern because many factors can be elaborated but some of them are characterized by negligible impact in this hydrological process (i.e., water seeps from the terrain surface to the underlying rocks). This was considered for many slight-impact factors, such as soil, karstified carbonate rocks and land cover/use, which might be used in other regions but not for the Riyadh Region because of shallow depth, slightly karstified carbonate rocks and rapid changes in the land components; respectively.

There is another point of concern that was accounted for while assigning weights of the different factors, that was the differentiation between the terrain exposure, including the land characteristics and the materials situated above it (i.e., lithology, fracture, streams, and slope), and the inputs on this terrain, with a special emphasis to rainfall. In other words, the study aims at identifying the land characteristics and the ability of the land to percolate water regardless of the existing amount of water derived from rainfall. This was the reason why rainfall, which is the only input for water in the studies region, was given little weight of impact which equals one-third of fracture weight for example. This concept of differentiating between exposure and inputs was also viewed from the understanding that zones with relatively high rainfall rates might not capture the precipitated water in that zone, and water may flow along valleys to other places.

The Riyadh Region, with crowded human activities, still lacks a thematic map showing the characteristics of its terrain surface to seep water into the underlying rock layers (i.e., recharge potential map). This map can be used for several purposes including but not limited to, identifying localities suitable for waste disposal sites and for groundwater artificial recharge, especially along valleys (Wadis) which transport huge amounts of surface water causing severe floods in the urbanized areas.

In addition, the GWRP map identifies the natural recharge zones, where high recharge zones represent the feeding domains for groundwater; and they must be protected from chaotic human activities, notably the processes result in creating artificial barriers (e.g., concrete pavements, settlements, etc.) that prevent water percolating into the underlying rock layers.

As resulted from this study, the Riyadh Region is generally characterized by low recharge potential, and this results in an exceeded overland flow and floods, as well as it gives more time for evaporation and then water loss. This must be well considered in water management approaches to cope with the increased demand due to population growth and the changing climatic conditions.

The production of the GWRP map for the Riyadh Region has been obtained before by the author (Al Saud 2015c), but it was generated by using only two factors including the lithology of the exposed rocks and the lineaments density maps. However, in this study, three more acting factors were added including the rainfall distribution, streams' density and surface slope, notably these factors can play a role in the surface water flow regime. Many of the factors used in this study were tackled in many previous studies. Nevertheless, these studies used some of these factors, as a result of data availability and the natural setting of the studied regions. In addition, the used methodologies in these studies were different than the current study, notably in the approaches used for geospatial data extraction and manipulation.

Results show that the Riyadh Region is mostly characterized by low-to-very low GWRP rates, where more than 75% of its territory reveals this hydrologic parameter (Table 7). Besides, there is less than 10% attributed to high and very high GWRP rates. This must be anticipated especially because the rocks' lithology comprises compacted alluviums and friable sand as well as argillaceous materials.

As an outcome of this study, which is associated with the author's vision to invest in Wadis' water in the Kingdom of Saudi Arabia, an in-depth understanding of the hydrological processes must be considered, with a special emphasis on the infiltration rate and rock bed characteristics, as well as the surface water flow regime. This vision would not be reached unless Wadis' waters are harvested and then injected into groundwater reservoirs (aquifers). This in turn will provide a considerable amount of water to the exhausted aquifers, as well as reduce floods' impact. Moreover, it can reduce the problem of saltwater intrusion into aquifers in coastal areas of the Kingdom of Saudi Arabia. Hence, artificial GWAR can also address this hydrologic problem along the coastal zones.

Finally, it is highly recommended to adopt the understanding of GWRP mapping for several regions in Saudi Arabia. In doing so it must be made clear that GWRP does not mean the artificial application of GWR, but it represents the natural hydrologic process.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.