A groundwater potential zone mapping approach for semi-arid environments using remote sensing (RS), geographic information system (GIS), and analytical hierarchical process (AHP) techniques: a case study of Buffalo catchment, Eastern Cape, South Africa

Theme unsuitability is noted to have inhibited the accuracy of groundwater potential zones (GWPZs) mapping approach, especially in a semi-arid environment where surface water supply is inadequate. This work, therefore presents a geoscience approach for mapping high-precision GWPZs peculiar to the semi-arid area, using Buffalo catchment, Eastern Cape, South Africa, as a case study. Maps of surficial-lithology, lineament-density, drainage-density, rainfall-distribution, normalized-difference-vegetation-index, topographic-wetness-index, land use/land cover, and land-surface-temperature were produced. These were overlaid based on analytical hierarchical process weightage prioritization at a constituency ratio of 0.087. The model categorizes GWPZs into the good (187 km2), moderate (338 km2), fair (406 km2), poor (185 km2), and very poor (121 km2) zones. The model validation using borehole yield through on the coefficient of determination (R2 = 0.901) and correlation (R = 0.949) indicates a significant replication of ground situation (p value < 0.001). The analysis corroboration shows that the groundwater is mainly hosted by a fractured aquifer where the GWPZs is either good (9.3 l/s) or moderate (5.5 l/s). The overall result indicates that the model approach is reliable and can be adopted for a reliable characterization of GWPZs in any semi-arid/arid environment.


Introduction
The water shortage is a global issue, owing to the nexus of water-food-energy and its influence on livelihood and global economics. In South Africa, the water shortage may persist for a longer time due to the rate of increase in population, urbanization, and industries, as well as the regional severity of aridity on water resources, which is further complicated by climate change (Owolabi et al. 2020b). Ad hoc effort towards water security has involved data gathering and evaluation, creation of impoundments and water infrastructures for water transfer schemes, water policies, optimization programs, and management measures (Schreiner and Hassan 2010). However, water supply has been outmatched and this is impacting other vital organs of development in the country. The exploitation of groundwater resources as an alternative has not been harnessed as a result of limited knowledge on its development (Cobbing 2014). Groundwater status and Responsible Editor: Biswajeet Pradhan development were relegated to private operations with limited restrictions (Pietersen et al. 2012). As a result, vast numbers of groundwater and surface water researches have been carried out as a separate entity until the reform of the Water Act (Act 36. of 1998) which also depolarizes their disjunctive management (Tanner and Hughes 2015). Information on groundwater distribution of groundwater is still poorly documented (Cobbing and de Wit 2018). Institutional framework and arrangement for data collection, private sector data handling, integrated developmental plans, and management strategy for groundwater resources are issues yet to be properly accounted for (Pietersen et al. 2012). Till present, few aquifers have been assessed for groundwater potential. Hence, the current research attempts to improve the awareness of groundwater availability and to proffer a regional approach to map the zones of groundwater potential.
Groundwater is a component of water resources that fills soil pore spaces, joints, and voids within geologic structures and strata. Groundwater occurrence in rocks depends on the hydraulic conductivity of lithologic materials which is a function of the porosity, permeability, and the flow of fluid through the geologic aperture or structures (Barlow and Leake 2012;Burberry et al. 2018;Owolabi et al. 2020a). Identification of the water-bearing structures and stratigraphic layer with significant hydraulic conductivity is therefore considered crucial to groundwater development. Groundwater potential zones (GWPZs) are areas enclosing the occurrence of a considerable and economically exploitable quantity of groundwater resources (Mandel 2012;Waikar and Nilawar 2014;Thompson 2017). Exploration of GWPZs is essential for water resources reserve estimation, zone budgeting, water quality protection, vulnerability mapping, and environmental management (Waikar and Nilawar 2014).
Conventionally, groundwater prospect has been explored directly through geologic, geophysical, and hydrogeologic approaches. In recent times, the advancement in geoscientific knowledge of data access and processing as well as the multifarious applications of geoinformatics has enabled the regional exploration of areas of groundwater potentials. The strong awareness gained by the geospatial technology for GWPZ mapping has facilitated drastic improvement especially due to the inculcation of multi-criteria decision-making (MCDM) statistical classifier (Machiwal et al. 2011). Some of the statistical tools include the weight of evidence, evidential belief (Tahmassebipoor et al. 2016), weighted overlay (Senthilkumar et al. 2019), multi-influencing factors (Anbarasu et al. 2019), analytical hierarchical process (Sahoo et al. 2017;Sandoval and Tiburan 2019), logistic regression , and frequency ratio . Among these, the analytical hierarchical process (AHP) has been rated as the most efficient multi-criteria decision-making tool (Jha et al. 2010). The evolution in the domain of forecast modeling and the advancement in computer programming have facilitated the discovery of machine learning techniques (MLT) for resource modeling (da Costa et al. 2019). Due to the versatility of MLT for analysis of intricate structures and stochastic data, it has offered a more reliable exploration outcome for groundwater resources (Lee et al. 2020;Pourghasemi et al. 2020;Prasad et al. 2020). On this note, several regression models of MLT were proposed for GWPZ mapping. These include multivariate adaptive regression splines (MARS), boosted regression tree Kordestani et al. 2019), support vector machine (Lee et al. 2018;Naghibi et al. 2018), artificial neural network (Lee et al. 2012;Lee et al. 2018), radial basis function, multiplelayer perception, standalone logistic regression (Pradhan 2010), and random forest (Golkarian et al. 2018;Arabameri et al. 2019). The robustness of machine learning techniques over unsupervised statistical approaches is due to its accuracy, speed, large database computation capability, and its ability to improve its self (Mueller and Massaron 2016).
Remarkable improvement in the performance of machine learning models has been intensified also. For instance, Rizeei et al. (2019) developed an ensemble multi-adoptive boosting logistic regression (MABLR) technique from the hybridization of a multi-adaptive boosting model and logistic regression. The model reports an excellent performance for groundwater aquifer potential maps for bias and variance error reduction due to insufficient sample size, oversimplification, and outliers sensitivity . Tien Bui et al. (2019) designed the hybrid computational intelligence approach called AB-ADTree from the integration of alternating decision tree classifier and adaptive boosting ensemble model for the assessment of groundwater spring potential zones. Kalantar et al. (2019) present two data mining approach based on the application of mixture discriminant analysis and linear discriminant analysis (LDA) compared with random forest for mapping the groundwater potential zone. The redundancy within the groundwater conditioning factors was controlled by normalizing the 15 conditioning factors in variance inflation factor, chi-square factor optimization, and Gini importance. Their outcome revealed the two data mining techniques are satisfactory and moderate even though the random forest approach proved to be better in performance . Also, Chen et al. (2019) developed a hybrid approach that involves the integration of Fisher's LDA, rotation forest LDA, and bagging LDA for mapping groundwater spring potential.
The computation capability of MLT is strengthened by its meta-algorithm compartment and the ensemble learning model that facilitates the multivariate analysis of a huge dataset (Prasad et al. 2020). MLT works by improving itself scientifically as sample distribution, sizes, and randomness increase while the increase in randomness and sample size increases the skewness and operational complexity of MCDM statistics (Mueller and Massaron 2016). Due to the huge dataset that machine learning works with, decision-making is based on probability, whereas in MCDM statistics, the variance and uncertainty of statistical operands are concealed by the results (Mitchell 1997). Machine learning does not require the definition of sample distribution and does not provide any room for assumption meanwhile MCDM statistics require the definition of distribution and enables assumption, thereby, providing the chance to hypothesize (Duan et al. 2016). As a result of the huge dataset, the outcome of machine learning can only be generalized since it works on probability and best fit, while the outcome of MCDM statistics can be fit to the defined data distribution (Kelleher et al. 2015). In this work, a few datasets (seven) are being used to assess their applicability for sitespecific GWPZ exploration. Hence, this study is rather research-oriented than being result-oriented, and to this end, AHP was adopted as the statistical classifier.
Pieces of literature that provide information on groundwater potential in the study area are relatively few. These include DWA (2010), Madi and Zhao (2013), Cobbing (2014), and Owolabi et al. (2020a, b). DWA (2010) characterized the regional attributes of the water management area and indicates that the study area is characterized by a groundwater yield of 2.0 l/s to 5.3 l/s range within aquifer of fractured and intergranular hydrostratigraphy domain. Madi and Zhao (2013) examined the significance of neotectonic belt to groundwater potential across Eastern Cape and noted the existence of a shallow quartzite vein with high groundwater potential in the study area through field mapping of geology structures. Cobbing (2014) assessed the factors influencing the underdevelopment of groundwater and poor management of groundwater boreholes and noted the existence of boreholes in the west of the study area that have been abandoned due to the mechanical breakdown of borehole hardware. Owolabi et al. (2020a) investigated the relationship between the environmental flow of Buffalo and the hydrostratigraphic properties of the catchment. Their work noted that the perennial attribute of the Buffalo River is due to spring and baseflow discharges from Buffalo aquifer. The work specifically indicated through flow duration curve, baseflow analysis, and hydrostratigraphy domain classification that Yellowwood, Tshoxa, and Mgqakwebe low-flow are driven by groundwater discharges, thus indicating the excellent potential for groundwater in the catchment. Owolabi et al. (2020b) assessed the degree of association and interdependence of the Buffalo streamflow on rainfall. The work indicated through bivariate statistical analysis that the hydrodynamic response of Buffalo streamflow is not only controlled by storm-flow but the existence of spring. With the numerous calls made on the need for improvement of knowledge on groundwater resources of South Africa, this work serves as a distinctive response to the call with a concentration on river basin-scale assessment (DWA 2010; Kahinda et al. 2016).
The present paper focuses on the assessment of zones of groundwater potential based on the agglomeration of statistically weighted geoscientific layers within the study area. Seven thematic layers were geospatially concerted using AHP for their weightage prioritization at a catchment scale. These include the thematic layers of surficial lithology (SL), rainfall distribution (RD), lineament density (LD), drainage density (DD), topographic wetness index (TWI), land use/ land cover (LULC), and land surface temperature (LST). Importantly, the selection of the thematic factors was based on the regional geology of the study area. Hence, this work intends to address the gap in the holistic investigation of groundwater potential in the Buffalo River catchment, Eastern Cape, South Africa. It proposes to project the interrelationship among influencing factors of groundwater development peculiar to the semi-arid environment. Based on the intent of the work, it is hoped that the selected themes would improve the application of AHP to an integrated approach for groundwater exploration.
The aim of this study is therefore to present the feasibility of developing a high-precision map of groundwater potential zone for a semi-arid environment. The rationale behind this approach is based on the integrated water resources management scheme that encourages a conjunctive development of water resources in an environmentally sensitive manner (Cosgrove and Loucks 2015).
The catchment-based geohydrological mapping is achieved based on the following objectives: (i) by producing the maps of the thematic layers (SL, RD, LD, DD, TWI, LULC, and LST maps); (ii) by quantifying their critical weights to be used for overlay analysis and their relative consistency ratio; (iii) by delineating and validating the groundwater potential zones for the assessment of the model reliability; and (iv) by corroborating and conceptualizing the geohydrology and hydrogeological variability across the environment of study. The paper presents a hybrid approach of geology, geophysics, geomorphology, and geoinformatics for the geohydrological characterization of groundwater enrichment spots. The approach developed in this study would be of tremendous benefit to the decision-makers, stakeholders, and the host community at large.

Study area
The study area, the Buffalo River catchment, is situated in the Buffalo Metropolitan District Municipality, Eastern Cape, South Africa, within the geographic coordinates of latitudes 32°40′ 07″ S to 32°58′ 50″ S and longitudes 27°7′ 54″ E to 27°33′ 16″ E (Fig. 1). It spans over an estimated area of 1237 km 2 , at an elevation of 258-1370 m above the sea level. Its geomorphological settings can be classified into three landform patterns: the plain (South), the dissected plain (trending from West to the Northeast), and the medium gradient mountain (Northwest). The main hydrologic feature in the catchment is the Buffalo River. It runs across a 54-km channel in Southeastern direction to the mouth of the Laing dam. It is fed by six major tributaries: Ngqokweni, Tshoxa, Mgqakwebe, Quencwe, North Zwelitsha, and Yellowwoods rivers. The catchment is characterized by a highly varying climatic pattern, with an average maximum temperature of 22.3°C in summer and an average minimum temperature of 13.5°C in winter (Slaughter et al. 2014). It has a mean annual rainfall of 590 mm per annum. The area is covered by two dominant soil types: the undulating and non-stratified sandyclayey soil, occupying the entire North to the center, and the unconsolidated clayey assemblages from the center to the South.
The area is covered by Tatarian arenaceous mudstone of the Balfour Formation, while the extreme lower section of the area exposed the Kazarian argillaceous shaly-sandstone of the Middleton Formation, both belonging to the Adelaide Subgroup within the Karoo Supergroup (Owolabi et al. 2020a). More than half of the areal geology is interspersed by Jurassic dolerite dykes and sills. The formation develops from the episodic history of the Adelaide subgroup, with Koonap and Middleton formations, which also belongs to the Beaufort Group, and Karoo Supergroup (Johnson et al. 2006). As a finding from the Neotectonic study of the region, Madi and Zhao (2013) noted that there are potentials for the development of groundwater resources in the catchment. The paleoenvironment setting of the Formation was reported to be characterized by a humid-to-temperate climatic type (Catuneanu and Elango 2001). Buffalo aquifer was classified under the Ciskean water management area with fractured aquifer and groundwater yield of 2-5 l/s (Vegter 2006).

Materials and Methods
Indices of geomorphology and environment such as topographic wetness index (TWI), land use/land cover (LULC), lineament density (LD), and drainage density (DD) have been extensively annotated in researches to characterize groundwater recharge index (Srivastava and Bhattacharya 2006;Hojati and Mokarram 2016;Aquilué et al. 2017).
AHP enables an unbiased and consistent pair-wise prioritization approach for the computation of weightage of features. The criterion for scoring the weight of features is based on the eigenvector of the square reciprocal matrix of paired features (Saaty 1999). AHP has been successfully adopted in environmental management (Althuwaynee et al. 2014), water resource Fig. 1 Location and regional geological map of the Buffalo catchment (modified from the Council of Geoscience regional map sheet (Johnson et al. 2006)) management (Saaty 1992), and numerous groundwater potential zone mapping (Rahmati et al. 2014;Mohammadi-Behzad et al. 2019;Rajasekhar et al. 2019). The interface of GIS tools presents the opportunity to integrate spatially encoded and statistically weighed groundwater influencing factors in a single map (Lillesand et al. 2014).

Geology
The geology of the study area was produced from the integration of information drawn from an enhanced aeromagnetic data, borehole lithology cross-section profile, and field geologic survey on a scanned geology map as the base map. The enhancement of aeromagnetic data enables the accurate mapping of basement rocks as well as the distinct detrital of sedimentary rocks due to the differential magnetic property of a geologic environment. The effectiveness of analytical signal for mapping basement configuration and lithology variability from an aeromagnetic data have been acknowledged in numerous literature (Matter et al. 2006;Baiyegunhi and Gwavava 2017;Owolabi et al. 2020a). The aeromagnetic data gathered by Fugro Airborne Surveys, South Africa, was supplied by the Centre for Geological Survey, South Africa. The geophysical survey for the aeromagnetic data was carried out with the use of a proton procession magnetometer with a resolution of 0.01 nT, at a constant flight height of 60 m in the North-South direction within a sampling line of 250 m and line spacing of 200 m. The data was processed in Geosoft Oasis montaj by defining the domain grid to reduce the long wavelength and high amplitude effect. The reduced-to-pole filtering was carried out by adapting the data to the average magnetic inclination of 63.47 and declination of −28.67 after removing the International Geomagnetic Reference Field (Peddie 1982). This was then convolved and filtered in the wavenumber domain by applying the 2D forward and inverse fast Fourier transform algorithm. The final computation on the processed aeromagnetic data involved the calculation of the analytical signal using Eq. (1): where AS represent the analytical signal feature, x, y, and z represent the edges of the magnetic structure, and T represents the total magnetic intensity. The resulting signal layer was classified into three: the extremely low, the intermediate, and the high amplitude signal for validation of relative surficial lithology types.

Rainfall
The data of rainfall from rain gauging stations around the study area were acquired from South Africa Weather Service. These include Stutterheim, Cata, Dimbaza, Berlin, East London, Kidd Beach, and Peddie stations. These were summated into an annual average for at least thirty years record, 1987 to 2016. In doing so, missing data are averaged out. The annual average rainfall was computed and projected to generate a rainfall thematic layer using the ordinary kriging interpolation method, and a linear semi-variogram model. The resulting spatially varying rainfall map was classified for overlaying purposes.
Land use/land cover Landsat 8 OLI with less than 10 % cloud was downloaded from the Earth Explorer USGS website and processed in ArcGIS 10.5.1 to produce the map of land surface temperature and the land use/land cover. The downloaded imagery was rectified for spectral distortion and reflectance quality. This was used to produce the LULC, LST, and LD maps. Mapping of LULC was carried out by generating a composite band of bands 1, 2, 3, 4, 5, 6, and 7. The buffering of LULC features was guided by the band combination validated by Butler (2013). This enables the easy identification of the constituent cells associated with the LULC features of interest presented in Table 1.
Hence, training samples were digitized and the signature file was developed for supervised mapping through a maximum

Land surface temperature
Landsat 8 Operational Land Imager of July 26, 2019, was downloaded from the USGS Earth Explorer website to produce the LST theme. The month represents the driest month of the year while the year was linked with the period of extremely low soil moisture (Owolabi et al. 2020b). The bands 4, 5, and 10 of the satellite data are processed in the following order: 1. Top of atmospheric (TOA) spectral radiance was calculated in raster algebra using Eq. (2): where RF M = radiometric multiplicative rescaling factor for TIRS 1 which is 0.0003342, TIRS = thermal infrared sensor 1 which is Band 10, and RF A = radiometric addictive rescaling factor for TIRS 1 which is 0.1.
2. TOA was converted to brightness temperature (BT) using Eq. (3): where K 1 and K 2 = Thermal conversion constant for Band 10 extracted from metadata, 774.885 and 1321.0789 respectively and 273.15 = constant for temperature conversion from Kelvin to Celsius.
4. Vegetation proportion (P v ) was calculated from NDVI, using Eq. (5): 5. Emissivity was calculated from P v , using Eq. (6): where 0.004 is the downscaling constant for P v and 0.986 is the correction value of the equation.
6. LST is deduced from the integration of BT and Emissivity, ℇ , according to Eq. (7): The approach adopted here for LST has been extensively validated in the previous land surface temperature researches (Orimoloye et al. 2018;Suresh et al. 2016). Land surface temperature indicates zones of poor saturation thickness, high evaporation, and evapotranspiration considering the severity of semi-aridity in the study area (Urqueta et al. 2018). Hence, its spatial information is considered significant to groundwater exploration in the study environment.

Lineament density
The surficial lineaments were produced using the panchromatic band (Band 8) of Landsat 8 OLI due to its high resolution (16 m by 16 m) in PCI Geomatica, 2017 version. This was achieved through an algorithm of multi-stage line detection of Canny edge and contour detection. The contour detection enables the filtering of curves and edges. The line detection enables a four-stage transformation which includes speciation of a minimum length of the curves, maximum error, maximum angle between polylines segments, and the minimum distance between two polylines (Hashim et al. 2013). These processes produce the polylines defined as lineaments. Its density was processed in ArcGIS 10.5.1 using grid cells method based on Eq. (8) (Rahmati et al. 2014): where LD is the lineament density, L i is the sum of the length of all the lineaments (km), i represent each linear feature in the study area, and A is the effective area of lineament cell grids (km 2 ). The lineament plot was validated by plotting its rose diagram and comparing its orientation to that of the Neotectonic structure in the study area. LD's significance to groundwater exploration is on account of its spatial analyst to indicates the zones of the tectonic macro-structures that serve as conduits for groundwater influx into hydrostratigraphic units (Fenta et al. 2015).

Drainage density
ASTER Digital Elevation Model (DEM) was downloaded from the USGS earth explorer website for the development of drainage density (DD) and topographic wetness index (TWI). The production of the drainage density map was achieved in the spatial analyst of ArcGIS. This begins with the transformation of the raw DEM into Fil > > Flow direction > > accumulation raster. Calculation of the effective area of drainage then enables the delineation of the watershed hydrological pattern. The drainage density index of the watershed was calculated using Eq. (9): where DD means drainage density and A is the area. Drainage density (DD) provides about the degree of ponding within a hydrologic unit (Srivastava and Bhattacharya 2006). Its computation is an essential inference for deciphering important hydrogeologic properties such as infiltration and permeability. Hence, it is DD is significant to groundwater exploration.

Topographic wetness index
TWI was mapped from the slope map, prepared from DEM in ArcGIS 10.5.1. This was achieved by calculating the rate of change of an aspect of a cell grid within its neighborhoodbased on Eqs. (10)-(11) (Hojati and Mokarram 2016): where a is the specific catchment area, A is the catchment area, L is the contour length, and tan(β) is the slope. TWI map was employed as the integrator of slope, elevation, and landform impacts on groundwater development (Pourali et al. 2016). Its computation sums up the influence of topographic roughness, hillslope, and foothill on lateral groundwater flow. Areas of high TWI enables the identification of areas of soil moisture accumulation and infiltration potential peculiar to foothills (Hojati and Mokarram 2016;Neilson et al. 2018). A significant highlight in the selection of the themes is the application of a surficial lithology map as a replacement for a conventionally scanned geologic map.

Computation of priority scores for the thematic layers and their consistency ratio
To improve the accuracy of the groundwater potential zone model, the relative degree of influence of each theme to groundwater development was computed as a priority score, based on the field experience of groundwater experts. AHP was employed to generate the priority score, and this was computed in a Microsoft Excel sheet in the following order: 1. Opinions of groundwater experts were sought on the degree of relevance of the groundwater factors on each other and to groundwater development. These were used to formulate the criterion rating. 2. The ratings of seven groundwater factors were collated as input for the groundwater criterion pairwise comparison based on Saaty's 1-9 scale (Saaty 1980). 3. The overall sums of scores of each element under the pairwise comparison were calculated. Each criterion cell was normalized by calculating the ratio of each pairwise cell to its sum in the columns. 4. Each row of the normalized matrix was computed for its average to derive the priority scores.
The consistency of criterion ratings across the pairwise comparison cells and the derived priority scores was assessed as it further enables the optimization of priority scale and subjectivity among the GWPZ factors (Jha et al. 2010). This was achieved in three stages as follow: 1. The consistency measure, λ max , was deduced as a principal eigenvalue based on the eigenvector technique. This is done by computing the matrix multiplication function of criterion ratings (in the pairwise comparison matrix row) and the normalized average of all the factors (within the normalized matrix column) divided by the criterion normalized average. 2. The consistency index was calculated using the consistency measure as presented in Eq. (12): where λ max is the consistency measure and n is the number of GWPZ factors.
3. The consistency ratio is calculated using Eq. (13) (Saaty 1980): where CI is the consistency index based on Eq. 13 and RCI is the random consistency index obtained from Saaty's 1-9 scale (Saaty 1980). The resulting final criteria weights are considered as the normalized values as long as the consistency ratio lies within the expected limit. For a consistent normalization, the value of CR is expected to fall within 0.01 to 0.09, otherwise, the priority scores have to be adjusted.

Overlay analysis for the delineation of groundwater potential zone
The overlay analysis of the seven thematic layers was carried out in ArcGIS 10.5.1 as a summation of the effective influence of the factors based on their criterion weight to generate the ultimate potential zone. To arrive at these, all the layers were converted into an integral raster while their units were ranked. The groundwater potential zone was calculated based on Eq. 11: where GWPI is the groundwater potential index, W i is the criteria weight, R i is the ranking of parameter factors, and i denotes each of the seven influencing factors with serial number from 1 to 7.

Validation of groundwater potential zone
The delineated GWPZs were validated using data of groundwater yield of exploration boreholes acquired from the National Groundwater Archive of the Department of Water Affairs, South Africa. The 84 geo-referenced yield data was overlaid on the GWPZ map to filter out the corresponding grid code identity of the borehole spot. The regression equation of the scattered diagram of the GWPZ grid codes against the groundwater yield was obtained to simulate the expected yield value. The coefficient of determination (R 2 ), coefficient of correlation (R), and the p-value of the relationship between the observed value and the simulated value were calculated in Microsoft Excel. R was obtained using Eq. (15): where O i is the observed value which is the borehole yield (l/sec), E i is the expected value drawn from the simulation of borehole yield, O ̅ n is the mean borehole yield, and n is the sample size of the borehole yield involved.
The relationship between the derived GWPZs and the actual borehole yield was further interpreted using a box-and-whisker plot. Representative borehole yield for each of the class of GWPZs was computed for their minimum value, lower quartile, median, upper quartile, and maximum value. These were used to derive the estimate for the bottom, 2Q box, 3Q box, whisker− and whisker+ (

Surficial Lithology
The surficial lithology map has been prepared using the geological map sheet of the Council of Geological Survey at a 1:250,000 scale as the base map. The surficial lithology is dominated by Mudstone (359 km 2 ) covering 29% of the area. Dolerite (281 km 2 ), sandstone (189 km 2 ), shale (56 km 2 ), and Quaternary sediment (4 km 2 ) cover 23%, 15%, 5%, and 0.3% of the river catchment, respectively (Fig. 2a). The mudstone and silty/sandy mudstone intercalation units were estimated to cover approximately 340 km 2 (27%) and 8 km 2 (0.7%) of the area respectively. Hydraulic conductivity is one of the foremost aquifer properties that determine groundwater recharge and storage potential of a lithology. Quaternary Sediment was considered as the highest importance for groundwater storage on account of its high hydraulic conductivity, followed by sandstone, silty sandstone, and silty/sandy mudstone (Barlow and Leake 2012). The least important for the groundwater storage is the dolerite rock, followed by the shale and the mudstone.

Lineament density map
The lineament density map is classified into five: the very high (4%), high (10), moderate (17), low (25), and the very low/ none lineament density area (Fig. 2b). The map shows the bias of lineament concentration to the North. This aligns with the topographic complexity in the North where the relief is steep and abrupt. Importantly, the zones of very high to moderate lineament density lie around the edges of the Karoo dolerite where the intrusion of magma creates a contact zone between the dolerite and the pre-existing sedimentary rock. Areas of extensive fracture system are replicated by high lineament density (Mostafa and Bishta 2005;Khosroshahizadeh et al. 2016;Meixner et al. 2018). Hence, the zones of effectively high/moderate lineament density, therefore, serve as a modifier for high-drained fractured dolerite that was ranked as very poor for groundwater potential under the surficial lineament. Rose diagram shows that the dominant trend of the surficial lineaments is WNW-ESE. This, therefore, conforms to the Maximum value-upper quartile direction of the Neotectonic structure shown in the study area geology map (Fig. 1). Zones of higher lineament density are expected to have higher potential for groundwater accumulation; hence, they are ranked higher.

Drainage density
The drainage density map of the Buffalo River catchment typifies dendritic drainage. Within the tributaries, only the Quencwe and Mgqakwebe rivers are characterized by dendritic drainage with well-developed stream orders while the Yellowwood, Ngqokweni, and Tshoxa rivers are associated with a trellis or subparallel drainages with underdeveloped stream orders (Figs. 1 and 2c). The dendritic drainage pattern of the Buffalo catchment and the northern sub-basins depict coarse dissection of the geomorphological layer impacted by the land cover system and climatic contribution. Meanwhile, the drainage pattern of the southern sub-basins indicates structural and lithological controls (Waikar and Nilawar 2014). The drainage density map was classified into four with areas of high concentration surrounding the river confluence where enormous groundwater discharge to the watercourse indicates a negative influence on aquifer development (Fig.  2c). The four spots of high drainage density in the South and one at the North can be associated with the existence of springs. They are shown to be proximal to contact zones of dolerite and the country rocks where lineament density is high. As a result, areas of higher drainage density are inversely ranked while the area of very low drainage density is assigned a higher rank.

Rainfall map
The mean annual rainfall for the hydrologic year records of 1989 to 2016 in the study area ranges from 544 to 610 mm (Fig. 2d). The rainfall trend shows a positive linear variation with relief in a Northwest-Southeast trend. This suggests the possible influence of relief on atmospheric circulation in such a way that favors orographic downpour at the high relief. Similarly, Owolabi et al. (2020a, b) noted that relief exhibits a linear spatial influence on a regional hydro-climatic pattern. The spatial variation in rainfall intensity influences the distribution of groundwater recharge rate across the study area. Due to the positive influence of rainfall on groundwater recharge, the areas with higher rainfall range are ranked higher.

Topographic wetness index
The result of the topographic wetness index for the study area is presented in Fig. 2e. The complexity of Buffalo topography is revealed by the lines of concentration of the low TWI, which coincides with the edges of dolerite outcrops in the North (Fig. 2). In a way, the concentrated low TWI patches depict the influence of dolerite intrusion on the initiation, evolution, and the development of rift (Madi and Zhao 2013). Consequently, the low TWI areas are more likely to develop an overland flow rather than enabling groundwater recharge on account of the influence of the hillslope factor. Meanwhile, high TWI which lies at the foothill of low TWI is more likely to enable groundwater recharge on account of the tendency for soil moisture accumulation (Fig. 2e). Soil moisture accumulation is a significant indicator of groundwater abundance, although the hydraulic properties of soil/lithologic material are the primary determinant of infiltration (Naghibi and Dashtpagerdi 2017).
TWI has been employed to decipher the average groundwater level in a watershed characterized by low permeability soils (Rinderer et al. 2014). Since the Buffalo catchment is dominated by argillaceous sedimentary material, TWI is therefore considered applicable to groundwater potential mapping here. The very low TWI that suggests the tendency for overland flow was assigned the lowest rank while the very high TWI that indicates the tendency for soil moisture accumulation zone was assigned the highest rank. The TWI showed high concentration in the north and center and coincidence with the drainage density (Fig. 2 c and e). The distribution of TWI also shows strong geologic control, influenced by the dolerite outcrop as shown in the Northwest, East, and the South (Figs. 1 and 2e).

Land surface temperature
The land surface temperature map is classified into four: the very low, low, moderate, and high as presented in Fig. 2 f. Areas of low temperature had the largest coverage (102 km 2 ), followed by a moderate temperature class (356 km 2 ), the very low-temperature class (508 km 2 ), and the least which is hot temperature class (272 km 2 ). The similitude in the geospatial attributes of LULC and LST further depicts the relative influence of urbanization in inducing urban heat index which indirectly culminates into the spatial variability in land surface temperature. This, therefore, conforms to the findings of Orimoloye et al. (2018) on the relationship between the LULC system and LST variability. Areas of high temperatures are associated with a high evaporation rate. Evaporation is a critical issue in a water-scarce country like South Africa, where it constitutes a significant soil moisture loss and diminution to shallow unconfined aquifers. As a consequence, areas with high temperatures are assigned the lowest rank while the areas with low temperatures are assigned the highest rank.

Land use/land cover
The result of the LULC mapping among the seven land cover features indicates that the cropland (320 km 2 ), has the highest coverage (Fig. 2 g). This is followed by the grassland (295 km 2 ), the built-up (247 km 2 ), the mixed forest (188 km 2 ), the scrubs (145 km 2 ), and the bare ground (37 km 2 ), while water bodies have the least coverage (4 km 2 ). The validity of the LULC plot was confirmed using Google Earth features. The Rooikrandsdam and its watercourse located by the Quencwe River were distinctly marked out just the same way as portrayed in Google Earth, The dispersed settlement across the Mgqakwebe, Tshoxa, and Ngqokweni rivers showed similar distal segmentation just as observed in Google Earth. Also, the Evergreen nature forest showed similar sectoral demarcation in the Google Earth the same way as mapped in the LULC map. Due to the varying impact of LULC changes on groundwater development, the different elements of LULC change are rated based on their significance to groundwater investigation. A mixed forest is a naturally preserved environment and an important indication of groundwater potential. Scrubs and cropland rank second to the indicators of groundwater potential; hence, they are ranked higher than others. Meanwhile, built-up area and water bodies/courses are the least important areas for groundwater investigation, followed by bare ground which is vulnerable to evaporation. Grassland is ranked higher than the bare ground because of its significance to soil moisture accumulation.

Normalization of features of GWPZ mapping
Relevant weights were assigned to the seven themes based on the influence on groundwater development during the computation of AHP as presented in Table 3. The value of CR obtained is 0.09, implying that the criteria weight was based on a reasonable level of consistency. The factors and their classes are presented in Table 4.

Delineation of groundwater potential zone
The potential zone of groundwater was classified into five: good (187 km 2 ), moderate (338 km 2 ), fair (406 km 2 ), poor (185 km 2 ), and very poor (121 km 2 ) zones as shown in Fig.  3. The study shows that the groundwater potential is high in the north and low in the south. The potential groundwater area lies in a perpendicular direction along the West-East direction to the direction of the catchment drainage. Overall, zones with high groundwater potential are dominant in the northwest, while the extremely poor groundwater potential zone is dominant in the south and southeast.

Validation and corroboration of groundwater potential zones
The scientific significance of models depends on validation reports; hence, it is considered the most important stage of scientific research. The spatial distribution of eighty-four exploration borehole yield employed is presented in Fig. 3. The coefficient of determination (R 2 = 0.901) obtained from the scattered diagram indicates that the GWPZ model is an excellent fit for the characterization of zones groundwater yield as it explains the borehole yield variability around its mean (Fig. 4).
This was further established by the coefficient of correlation and the test for significance (p value < 0.01) for the model which that the modeling procedure shows a very significant and strong positive relationship with the borehole yield (Table 5), that is, the model can provide relevant and applicable information on the availability of groundwater in place of borehole yield. A more expository relationship between borehole yield and GWP rate is presented in the box-and-whisker Fig. 2 The thematic maps of Buffalo Catchment. a Improved geological map of the Buffalo Catchment. b Lineament density map. c drainage density map. d Spatial distribution of precipitation. e Topographic wetness index map. f Land surface temperature map. g Land use/cover map  Fig. 5. Twenty boreholes each lies within the good and the moderate GWPZs with a yield range of 6.03-12.15 l/s and 3.15-9.87 l/s respectively while 15 boreholes lie within the fair GWPZs. The poor and the very poor GWPZs were occupied by 19 and 10 boreholes with the yield range of 0.19-2.05 l/s and 0.09-1.83 l/s respectively. The good GWPZ class reports a high level of agreement with borehole yield, whereby more than half of its class is 9.3 l/sec. The moderate GWPZ class is associated with the widest range and the highest yield standard deviation with an average value is 5.5 l/sec. The good and moderate GWPZs is therefore considered the most significant groundwater capture zone. The fair, poor, and very poor GWPZ classes are associated with 3.2, 0.8, and 0.4 l/s average yield which is not considered economical for groundwater development where water demand is high. The plot reveals that groundwater borehole exploitation sited on the good and moderate GWPZs area is most likely to provide considerable yield.

Discussion
Characterization of the Buffalo catchment groundwater potential The GWPZ mapping reveals that Buffalo groundwater resources are hosted by two hydrostratigraphic units: the fractured unit, dominant in the permissive contact zones of dolerite, and the inter-granular unit comprising of the arenaceous sedimentary materials (Figs. 2a and 3). This dominant borehole yield falls within the range of 2-5 l/s, the fair GWPZ (Fig. 3), covering 33% of the study area. These findings agree with the report of Vegter (2006) on the hydrogeology of the regional water management area of Buffalo Catchment. However, excellent borehole yield is regionally biased in the study, thus occupying the Northern half of the study area, where the lineament density is high. The outcome of the GWPZ model suggests that the groundwater distribution in the Buffalo River catchment is strongly controlled by the permissive contact zone of dolerite and country rocks ( Figs. 1 and 3). The distribution of the groundwater yield shows a cluster of high yield at the West compared with the East; meanwhile, considering the orientation of the lineaments, groundwater flow could be assumed to flow from the West to the East (Figs. 2b and 3). This study, therefore, confirms the assertion of Chevallier et al. (2014) on the groundwater potential of weathered and fractured dolerite as a viable aquifer. The poor potential of the dolerite outcrop in the South is due to the limited/non-existence of surficial lineaments in the area (Fig. 2b). Owolabi et al. (2020a) noted that the Buffalo streamflow is a perennial flow and that the river is possibly sustained by the existence of springs or leaky aquifer. The findings here agrees with the submission as the perennial attributes of Buffalo drainage is possibly due to replenishment from perpendicular groundwater discharges at the suggested spring zones (Fig. 2c).

The inter-relationship across the influencing factors of groundwater
The high lineament density area depicts the influenced by tectonic activities, which also possibly accounts for the uplift in the northwest and the hillslope that deter settlement residency (Chevallier et al. 2014). The structural evolution influenced by tectonic activities in the northwest area is depicted by the high concentration of TWI features in conformity to Mandal and Mondal (2019). Hence, the northern half was validated to be associated with an extensive fractured system, joints, and macro-scale geologic structures that can facilitate groundwater accumulation and flow. This corresponds to the highly vegetated area of the land use/cover map (Fig. 2 g) and the highly dissected area of the geomorphological area shown by the drainage density map (Fig. 2c). The coarseness of the high drainage density area in the northern half, therefore, corresponds to the areas associated with high permeable soil and rock materials in the area in conformity to Waikar and Nilawar (2014). Moreover, Johnson et al. (2006) and Katemaunzanga and Gunter (2009) noted the northern half is characterized by abundant siltstone and sandstones compared with the Middleton Formation that is mainly dominated by Mudstone and shale at the southern (Fig. 1). The geomorphic evolution of Buffalo drainage can also be explained by the similitude in the trend of precipitation and drainage texture (Fig. 2 c and d). A similar discovery was noted by Ngapna et al. (2018) whose investigation of climate and morpho-tectonic interaction was through morphometric analysis. Considering the land use/cover features (Fig. 2 g), the area covered by dense vegetation at the north is associated with a high degree of high ruggedness which implies high infiltration potential and lesser runoff while the urban area dominates the south and consequently culminates into a high degree of runoff. This inference can be drawn to describe the dependence of drainage density on the LULC features (Fig. 2 c and g) and its influence on groundwater discharge at the Buffalo catchment mouth in agreement with the assertion of Brody et al. (2014). Overall, the computation approach used in this work presents the GIS-based GWPZ model as an excellent tool for exploring the prospective spots of groundwater for sitespecific geologic investigation. In the previous GWPZ mapping involving AHP, Maity and Mandal (2019) noted that the validation provides 78% credibility while Roy et al. (2020) acknowledged a discouraging validation outcome. However, in this work, a significant 90% plus assessment accuracy was obtained by the R 2 and R here. This not only validated the applicability of AHP computation but also accredited the quality of influencing factors employed and the proficiency of pairwise weightage assignment.

Conclusions
The application of the geosciences approach involving the conjunction of non-linear and spatially-independent  environmental, hydrologic, and geomorphic parameters to geologic parameters for mapping the zones of groundwater potential has been demonstrated. The uniqueness in the work lies in the selection of the parameters and the conceptual steps followed to narrow down uncertainty and ensure high precision mapping. The following conclusions can be made from the study: & The hybridization of aeromagnetic data scanned geologic map, and field geologic survey offers a high-precision map of surficial lithology. & The application of the analytical hierarchical process has affirmatively revealed the feasibility of integrating the stochastic geoscience themes for environmental management. & The robustness of site-specific factors for the production of accurate GIS-based groundwater potential was proven achievable. & Strong interrelationship was shown by the correspondence between tectonically stressed zones and zones of compacted TWI features, a trend of precipitation and drainage, the zones of coarser drainage texture, and areas of permeable earth materials.
Overall, the computational approach exhibited in this work can be adopted for the characterization of groundwater potential zones and areas requiring conservative use against groundwater contamination based on the seven parameters employed in any semi-arid/arid environment with similar fractured geology.