Geohydraulic and vulnerability assessment of tropically weathered and fractured gneissic aquifers using combined electrical resistivity and geostatistical methods

Sustainable potable groundwater supplied by aquifers depends on the protective capacity of the strata overlying the aquifer zones and their thicknesses, as well as the nature of the aquifers and the conduit systems. The poor overburden development of the Araromi area of Akungba-Akoko, in the crystalline basement of southwestern Nigeria, restricts most aquifers to shallow depths. Hence, there is a need to investigate the groundwater quality of the tropically weathered and fractured gneissic aquifers in the area. A combined electrical resistivity tomography (ERT) and Schlumberger vertical electrical sounding (VES) technique were employed to assess the groundwater-yielding potential and vulnerability of the aquifer units. The measured geoelectric parameters (i.e., resistivity and thickness values) at the respective VES surveyed stations were used to compute the geohydraulic parameters, such as aquifer resistivity (ρρoo), hydraulic conductivity (K), transmissivity (T), porosity (φφ), permeability (Ψ), hydraulic resistance (KRR), and longitudinal conductance (S). In addition, regression analysis was employed to establish the correlations between the K and other geohydraulic parameters to achieve the objectives of this study. The subsurface lithostratigraphic units of the studied site were delineated as the motley topsoil, weathered layers, partially weathered/fractured bedrock units, and the fresh bedrock, based on the ERT and the A, H, AK, HA, and KQ curve models. The K model regression-assisted analysis showed that the ρρoo, T, φφ, Ψ, and S contributed about 81.7%, 3.31%. 96.6%, 100%, and 11.63%, respectively, of the determined K values for the study area. The results, except T and S, have strong high positive correlations with the K of the aquifer units; hence, accounted for the recorded high percentages. The aquifer units in the area were classified as low to moderate groundwater-yielding potential due to the thin overburden, with an average depth of <4 m. However, the deep-weathered and fractured aquifer zones with depths ranging from about 39–55 m could supply high groundwater yield for sustainable exploitation. The estimated S values, i.e., 0.0226–0.1926 mho, for aquifer protective capacity ratings rated the aquifer units in the area as poor/weak to moderately high with extremely high to high aquifer vulnerability index, based on the estimated low Log KRR of about 0.01–1.77 years. Hence, intended wells/boreholes in the study area and its environs, as well as any environments with similar geohydraulic and vulnerability characteristics, should be properly constructed to adequately prevent surface and subsurface infiltrating contaminants.


Introduction
Globally, the percentage of people who use potable water has increased twice as fast as the global population [1,2,3]. Sustainable groundwater yield in aquifer zones depends on the detailed characterization of subsurface strata, the water-retaining capacity of the strata due to porosity and permeability, water-rock interactions, subsurface conduits, and storage zones, as well as the hydrodynamics of the aquifer units [3,[4][5][6][7][8]. The sustainability of groundwater supplies also depends on the quality of the aquifers' yield, which is a function of the protective capacity of the strata overlying the aquifer zones and the depths of the aquifer zones. The occurrences of aquifer zones at shallow depths, especially within the crystalline basement terrain, give easy access for the percolation of surface runoffs and pollutants from dumpsites' leachate flows, surface and buried oil tank spillages, dissolved chemicals from mining activities, sewage from sanitation systems, etc., to degrade the stored groundwater [5,[9][10][11][12]. In addition, over-stretching of aquifers caused by over-abstraction of groundwater and silt/clay intrusion from improperly cased boreholes degrades the quality of groundwater [3, [13][14][15].
Geoelectrical resistivity methods, i.e., resistivity profiling and vertical electrical sounding (VES), and the estimation of geohydraulic parameters (e.g., hydraulic conductivity, transmissivity, porosity, permeability, transverse resistivity, longitudinal conductance, hydraulic conductance, etc.) from georesistivity datasets and/or pumping tests have been employed in the determination of groundwater-yielding potential and vulnerability of aquifer units in several geological terrains [3,5,13,16,17]. However, the pumping test method is time-consuming and expensive; hence, it has not often been utilized recently in geohydraulic evaluation since georesistivity methods are cost-efficient, rapid, and produce quality results with a higher success rate, e.g., [3, 10,11]. The advantage of the georesistivity methods is that the measured resistivity values have strong correlations with groundwater hydraulic characteristics; hence, the data offer the determination of geohydrodynamics of aquifer units, protective capacity of the near-surface strata, and selection of suitable points for sustainable potable groundwater development [5,10,11,18].
The study area covered a part of the Araromi area of Akungba-Akoko, southwestern Nigeria, which is situated between the main Akungba-Akoko town and Etioro-Akoko and is characterized by complex subsurface geology, e.g., [3,8,[19][20][21]. The study area has become the choice location for many local settlers, staff, and students of Adekunle Ajasin University, Akungba-Akoko (AAUA) due to its serenity and prospect for rapid urbanization. According to previous workers in the study area and the surrounding towns, e.g., [3,8,[19][20][21], the overburden in the areas is poorly developed and hence causes incessant water shortages to meet the growing population due to perennial failures of hand-dug wells and boreholes. In the report of Mohammed et al. [19], the groundwater potential of the northern section of the present study area was evaluated using the vertical electrical sounding (VES) technique.
Akingboye et al. [8] investigated the near-surface crustal architecture and geohydrodynamics of the Araromi area, Akungba-Akoko (i.e., the study area) using the integrated coplanar loop electromagnetic conductivity method, electrical resistivity tomography (ERT), and Schlumberger VES technique, to ameliorate the difficulties of sufficient groundwater availability and the failure of engineering structural foundations in this area. Furthermore, the subsurface geological, hydrogeophysical, and engineering characterization, as well as the vulnerability of the southern part, (i.e., Etioro-Akoko), some few meters away, has been carried out and reported, e.g., [3,8,21]. From these studies, it was reported that most of the aquifer zones occurred at shallow depths, except for the localized ones that are characterized by deep-weathered troughs and fractures.
Despite the several detailed studies mentioned above, the vulnerability of the tropically weathered and fractured aquifer zones in the Araromi area of Akungba-Akoko has not been evaluated. Hence, it becomes necessary to determine the vulnerability of the aquifer zones to contamination due to the occurrence of most of the aquifers in the study area at shallow depths. The fracture densities and groundwater-yielding potential of some identified aquifers can also offer clues on the migration rate of possible contaminants in the area. Given the above, geoelectrohydraulic method, involving combined ERT and Schlumberger VES from which geohydraulic characteristics were determined through various equations and regression-assisted analysis to evaluate the vulnerability of the aquifer units in the study area. The study area is characterized by the Nigerian rainforest belt climate, with average yearly rainfall between 1000 mm and 1500 mm, and temperature is around 33 o C. The area has a topographic relief consisting of hills, 4 low-lying outcrops, plains, and valleys, between 280 m and 400 m above the mean sea level. The dendritic drainage system follows these topographic features, with trellis drainage patterns in a few places, e.g., [3,8,21].

Geological setting of the study area
Geologically, the Araromi area of Akungba-Akoko is characterized by the Nigerian Southwestern Precambrian Basement Complex, which is part of the reactivated Pan-African mobile belt, occupying the east of the West African Craton and northwest of the Congo-Gabon Craton [22][23][24][25], as shown in Fig. 1a. The Southwestern Basement Complex of Nigeria is made up of three major rock suites, namely the Migmatite-Gneiss Complex, ranging in age from Neoproterozoic to Paleoproterozoic and Archean, e.g., [22][23][24]; the Neoproterozoic Schist Belts, consisting of low-grade, younger metasedimentary, and metavolcanic rocks with ages ranging between 690 and 489 Ma, e.g., [23,24], and the Pan-African Older Granites, which intruded the two earlier lithologies, have ages ranging between 650 and 580 Ma, e.g., [26,27]. Early magmatic phases dating from 790 to 709 Ma have also been reported in some of the Older Granites rocks. The Younger Granites, i.e., the Mesozoic anorogenic calcalkaline ring complexes, as shown in Fig. 1b, intruded the Nigerian Eastern Basement Complex terrain before the formation of any of the sedimentary basins in the eastern and western terranes of Nigeria [26][27][28].
The entire Akungba-Akoko is underlain by the Migmatite-Gneiss Complex rocks of southwestern Nigeria, which were intruded by the Pan-African Granitoids as shown in Fig.   1c. The Migmatite-Gneiss Complex rocks in the area are typically migmatite, granite gneiss, and biotite gneiss, as well as granitoids consisting of charnockites and granites. Granite gneisses are the most abundant rock type in the area, and this particular rock underlies the Araromi part of the area. This particular rock type has a blastoporphyritic to porphyroblastic fabric and is light grey, medium to coarse-grained, and moderately foliated, for example (light and dark-colored bands). Far to the west, the rock is extensively deformed and migmatized, forming migmatite with an ENE-WSW trend. In addition, important intrusives in the granite gneisses in the area include quartz veins, pegmatite, aplite, basic dykes, and sills [3,8,23,29]. The tropical climatic conditions combined with the metamorphic activities in the Akungba-Akoko Basement Complex terrain have assisted in the weathering and fracturing of the subsurface strata.
The subsurface hydrogeological features of the Araromi area of Akungba-Akoko are similar to those of the surrounding towns, e.g., the main Akungba-Akoko town and Etioro-Akoko community, e.g., [3,[19][20][21], as well as some places in the crystalline basement of southwestern Nigeria, e.g., [11,12,30]. The groundwater in the study area occurs in weathered and/or fractured aquifer zones. Groundwater occurrence in these hydrogeologic units is unevenly distributed, just like other parts of the Precambrian basement terrain.
Generally, the aquifers in the area are characterized by shallow depths with low porosity and permeability, and hence depend on the secondary porosity resulting from deep weathering and fracturing of the rock units to conduit and store fluids in the subsurface strata sufficiently.
It has been reported that the aquifer zones in the study area and surrounding communities have an average depth of about 12 m, with depths exceeding 25 m for deeply weathered and fractured aquifers, e.g., [3,8,21]. Some of the hand-dug wells and boreholes sited in the study area take advantage of the former and latter aquifer depths, respectively, for the required groundwater supplies.

Methods of Study
Due to the urgent need of the inhabitants of the Araromi area for steady and sufficient drinkable groundwater supplies to suit their everyday activities, a detailed subsurface geologic condition is highly necessary. The field data collection began with the identification of prominent areas with records of failed hand-dug wells and low yield boreholes in the study area. Predictions from earlier studies in the northern and southern areas of the study region, such as [3,8,19,21], enabled the selection of locations for establishing the geophysical traverses for data collection quicker and simpler. Six traverses (TRs) were occupied in the study area as shown in Figs Because of the population increase, structural building blockage caused a reduction in geophysical survey spread lengths on some traverses.
The ABEM Resistivity Imaging System was used for the ERT field data acquisition, utilizing the dipole-dipole electrode configuration protocol array because of its high sensitivity to vertical and lateral subsurface structural variations and low electromagnetic coupling effects, e.g., [8,31,32]. A station interval of 5 m was used for the detailed subsurface imaging of the anomalous features of interest for this study. Although the adopted n-level of 5, i.e., (n = 5), for dipole-dipole resistivity surveys could limit the depths of probing. However, the station interval is considered suitable to derive more cluster near-surface information and to avoid nuisance surface artifacts arising from the complex geological condition of the study area terrain. The Schlumberger electrode configuration for the VES technique, on the other hand, was carried out at the selected conductive or relatively conductive survey station points to address depth limitations. The approach was aimed at constraining the modeled ERT results, and to image deep-weathered bodies and the penetrative fractures. Figure 2b depicts the spatial distribution of the investigated VES station sites, whereas Figs. 2b and 2c depict the elevation of the surface topography. The current electrodes AB/2 varied from 60 to 160 m, whereas the potential electrodes spread MN were varied from 0.5 to 15 m. The depth of penetration in a homogenous subsurface geologic structure is proportional to the distance between the current electrodes, whereas, varying the electrodes distance offers information regarding the subsurface lithostratigraphic units, e.g., [3,8,12]. When a remarkable resistivity of fresh bedrock, typically with values >1200 Ωm, was attained more than twice at each VES station, the survey was stopped, indicating that there was no possibility of obtaining a fracture at deeper depths even for further probing. However, because of the barriers encountered due to primarily buildings, the surveys were sometimes halted before reaching depth to the fresh basement.
The results from the ERT surveys were processed and inverted using RES2DINV software.
Forward modeling and data inversion, utilizing the least-squares inversion approach, which depends on a mathematical inverse problem to derive the subsurface resistivity distribution from apparent resistivity data sets, were used in the inversion process. Many works have reported on the adopted inversion processes used in the RES2D data inversion, including Akingboye et al. [8], Loke [31], DeGroot-Hedlin & Constable [33], Dahlin & Loke [34], Akingboye & Bery [14,35], and other. The ERT field data sets with topography were inverted using the finite-element method of 4 nodes with L2-norm as the least-squares constraint parameter to minimize the difference between the measured and calculated apparent resistivities. A damping factor of 0.05 with a minimum value of 0.01 was employed to increase the accuracy of the calculated apparent resistivities and the resolution of the generated apparent resistivities. The root-mean-square (RMS) error limit for inverse model convergence was set to less than 10% for a maximum of 7 iterations. The desired cut-off error was set at 40% with a maximum error of 200% to achieve the desired results. Figure 3 shows the composite ERT inversion results, which include measured and calculated pseudosections as well as inverted resistivity sections for the analyzed TR1. The VES field results were successively inverted using the IPI2win software to curve match the field data to generate the model resistivity curve, which included the thicknesses and depths of the geoelectric layers, as well as the resistivity values of the delineated layers. High anomalous peaks in a curve compared to the surrounding stations were reduced in comparison to the surrounding data, owing to poor electrode grounding, circuit relay, or current transmission issues due to dry ground. Following these corrections, the iteration of such VES field data was repeated. The RMS error of the iteration convergence limit was set at less than 10%. The results were used to compute the geohydraulic characteristics of aquifer zones in the studied area. Other software, involving Oasis Montaj TM and Geosoft Surfer TM , was used to produce two-and three-dimensional (2-D and 3-D) maps for this study.
The hydraulic conductivity (K) and transmissivity (T) of the aquifers, expressed in m/day and m 2 /day, were computed using Equations 1 and 2 given by Heigold et al. [36] and Niwas & Singhal [37], respectively. Transmissivity gives the areal extent of pore-water flow per day in the saturated hydrogeologic units. In addition to these parameters, the porosity ( ) and permeability (Ψ) of the hydrogeologic units at each VES point were also computed using Equations 3 and 4, respectively. The longitudinal conductance (S) of the overburden at each VES station point was also estimated using Equation 5, as suggested by [38,39]. The hydraulic resistance (K ), expressed in years, was estimated to ascertain the AVI of the hydrogeologic unit in the studied site using Equation 7 given by Van Stempvoort et al. [40].
The logarithm of the hydraulic resistance, i.e., K , was estimated to measure the AVI of the APC of the overburden unit to the vertical flow of fluid.
Where 0 is the resistivity (Ωm) of the aquifer, , , and h are the conductivity, longitudinal conductance (mho), and the thickness (m) of the aquifer, respectively. , , and are the water dynamic viscosity adopted as 0.0014 / / according to Fetters [41], the density of water (1000 / 3 ), and acceleration due to gravity, respectively. and ℎ are the resistivity and thickness of the i th layer, respectively.
To further substantiate the analyses of the geohydraulic parameters for the Araromi area of Akungba-Akoko, regression analysis was performed using geohydraulic conductivity (K), as the independent variable, to predict the values of the dependent variables, i.e., transmissivity (T), porosity ( ), permeability (Ψ), transverse resistance ( ), longitudinal conductance (S), and hydraulic resistance (K ), as well as aquifer resistivity ( ). These geostatistical analyses illuminate the relationship between the predicted parameters (i.e., dependent variables) and the independent variables, as well as the percentages of surface contaminants that pose a vulnerability risk.

Subsurface lithostratigraphic and structural characterization
The results of ERT for the surveyed traverses are presented in Figs. 4a-f, while Table 1 presents the results of the VES survey stations points, including the curve types, thicknesses, depths, and descriptions of the delineated subsurface layers. The subsurface layers in Fig. 4a are distinguished by four distinct subsurface layers: motley topsoil, weathered layer, partially weathered/fractured bedrock, and fresh gneissic bedrock, with resistivities ranging from 10-600 Ωm, 600-1000 Ωm, 1200 Ωm, and >1200 Ωm, respectively. The deep-weathered and  Table 1. This affirms that the delineated deep-weathered/fractured zones between 80 and 100 m are deeper than the imaged depths, as shown in Fig. 4c.   Table 1.
However, the depth of the weathered column, as presented in Fig field and data inversion) approaches adopted for generating the resistivity models for the studied site.

Geohydraulic characteristics of tropically weathered and fractured gneissic aquifers and empirical relationships: insights into groundwater yielding potential
The identified weathered and fractured aquifer zones in the study area are characterized by varied aquifer resistivity (i.e., ), and thickness values ranging from 138-838.10 Ωm, and 3.5-50.10 m, respectively, as presented in Table 2. The estimation of aquifer parameters, e.g., K, T, , and Ψ, is important for evaluating the groundwater potential of tropically weathered and fractured bedrock aquifers. The estimated K and T values for the aquifer zones in the study area ranged from about 0.7246-3.8985 m/day and about 53079-60.8597 m 2 /day, respectively, as listed in Table 2. The value of T above 20 m 2 /day was estimated for VES surveyed points 1, 2, 7, and 8, while the rest of the VES points recorded values far below this estimated value. This variation may have been due to the higher thickness and aperture of the aquifer zones, which provided room for a higher rate of water-rock interaction from water flows, e.g., [3,8,16,17]. As presented in Table 2, the estimated values of and Ψ for aquifers in the study area ranged from 24.05-31.62%, and 0.103-0.556 µm 2 , respectively.
The measured range of values for the recorded porosity for the aquifer depicts the dominancy of consolidated weathered materials typical of clay, sand, and lateritic clay, e.g., [3,7]. The reduction in porosity at some VES surveyed stations may be due to a decrease in fluid contents, bulk conductivities, and/or a decrease in the hydraulic conductivities, e.g., [7]. The geostatistical results derived from regression analysis as presented in Table 3 100%, respectively, based on the coefficient of determination ( 2 ) results. However, the percentage contribution of T to the K of aquifers in the study area is about 3.31%, as listed in Table 3. The and T decline with the aquifer K at a rate of 0.0046 and 0.0112, while , and Ψ increase at a rate of 0.4342 and 7.0089, respectively. The low percentage contribution of T, which may reduce aquifer transmissivity, could be attributed to the high resistivity values of some aquifers and the occlusion of the deep-weathered and fractured zones by soils produced secondary weathering. The observed 2 values suggest that the model fitted reasonably well with the used variables for the determination of K [3, 42,43], except for T. In addition, the high percentage 2 contributions of both porosity and permeability to the K model suggest a significant contribution of both parameters to water-rock interactions, and the ease of fluids transmissibility in water-bearing aquifer units, e.g., [3,7,8]. The estimated standard errors offered the variability determination of the coefficients and a significant non-zero slope. The above-evaluated parameters have low coefficient standard errors; hence, suggesting very low statistical variation. The T-stat values with corresponding p-values were used to determine the accuracy and robustness of the analysis. The above dependent variables (i.e., , , and Ψ) yielded p≤0.05, except for T; hence, this significantly validates the accuracy of the geostatistical model. The lower and upper 95% confidence limits ranging from about -0.0067 to -0.0024, -0.0714 to 0.0491, 0.3525-0.5158, and 7.0021-7.0157 for , T, , and Ψ, respectively in Table 3, account for the unknown K values [3, [42][43][44]. The empirical relationships between K and the parameters , T, , and Ψ, were derived from their listed values in Table 3, and are given in Equations 7-10, respectively.  Table 4, classify the aquifer potential types in the study area into the low and moderate groundwater-yielding capacity aquifer zones. These aquifer types are efficient for water withdrawal at smaller quantities, (i.e., local groundwater supply), the typical rate for private/personal consumption, and also for smaller communities based on water transmissible rates. However, delineated aquifer zones with deeper depths fractures exceeding 39 m and higher porosity and permeability values above 30% and 0.4 µm 2 respectively (see Table 2), have higher groundwater-yielding capacity. The results further affirm that aquifer transmissibility depends on the physical characteristics of the subsurface geologic units.

Evaluation of the vulnerability of aquifer zones in the study area
The evaluation of aquifer vulnerability in the study area depends on the estimation of the longitudinal conductance S and hydraulic resistance K , to classify the aquifer protective capacity (APC) and aquifer vulnerability index (AVI). The K is an essential geological formation factor to determine the resistance of an aquifer to vertical fluids flow through the protective subsurface strata because it depends on K but varied inversely proportional to it.
The relationship between the AVI and the logarithm of the hydraulic resistance, i.e., Log K , for the VES surveyed station points are as shown in Table 5. In addition to the parameters, the thicknesses of aquifers in the study area were taken into consideration while evaluating the aquifer vulnerability. The S of the tropically weathered and fractured gneissic aquifers units in the study area ranged from about 0.0226-0.1926 mhos, as presented in Table 2. The  Table 3. The K yielded a moderate positive correlation of about 0.5769 with K. This parameter contributed about 33.29% to the estimated K model. The weak positive correlation between K and S implies that both parameters are independent of one another just as explained above, the lower and upper 95% confidence limits derived for the S and K account for the unknown K values that are not presented in Table 3. The K is statistical related to S and K through Equations 11 and 12, respectively, which were derived from model results in Table 3. It is worth noting that the K values for the tropically weathered and fractured gneissic aquifer zones in the study area can be determined using any of the empirical relations provided in Equations 7-12 from the regression analysis.  Table 5. These characterized ranges of values coincide with the range of values for Log K , as suggested by Van Stempvoort et al. [40], which measures AVI by hydraulic resistance.
Hence, the aquifer zones were classified by the AVI into poor/weak and moderate APC, with extremely high vulnerability and high vulnerability as presented in Table 5. These results suggest that the Araromi area of Akungba-Akoko has poor/weak to moderate APC with extremely high AVI and high AVI, respectively, due to the generally thin overburden soil material cover and low hydraulic resistance, e.g., [3,46].

Conclusions
The geohydraulic characteristics and vulnerability of tropically weathered and fractured gneissic aquifers in the Araromi area of Akungba-Akoko, southwestern Nigeria, have been assessed using a combination of ERT, VES, and geostatistical (i.e., regression analysis) methods. The results of the resistivity models revealed four distinct layers. The motley topsoil is generally thin (i.e., <1.7 m) in most sections, but thickness values ranging from about 2.5-3.5 m were recorded at VES 4 and VES 5 along TR3. Just like the areas within the Akungba-Akoko and Etioro-Akoko, the thickness of the overburden is also generally <4 m, except for deep-weathered and fractured bedrock sections that extend to depths above 39 m, and such characteristic features are localized in the study area.
The regression analysis revealed that the , T, , and Ψ, contributed significantly, about 81.7%, 3.31%. 96.6%, and 100%, respectively, to the determination of K for aquifers in the study area. Based on the delineated lithostratigraphic units, subsurface geologic structures, and their varying depths, as well as the geohydraulic characteristics, the aquifer units in the study area are classified as having low to moderate groundwater-yielding potential. However, the enhanced fractured aquifer zones to depths above 55 m can produce adequate groundwater yield in some places in the area. The reported S values (i.e., 0.0226-0.1926 mho) for APC ratings in the Araromi area of Akungba-Akoko were too low, and thus, the APC ratings are rated poor/weak to moderately high with an extremely high to high AVI. It is, therefore, important that intended wells and/or boreholes in the study area and any environments with similar geohydraulic characteristics and vulnerability ratings should ensure proper construction designs for adequate protection against both surface and subsurface infiltrating contaminants. This study has provided significant insights into the assessment of sustainable potable groundwater development in crystalline basement geologic environments using integrated geophysical resistivity and regression analytical approaches.

Acknowledgments
The Department of Earth Sciences, Adekunle Ajasin University, is appreciated for providing the field equipment used for this study. Ayanfe Moses Asulewon (formerly of the Department of Earth Sciences, Adekunle Ajasin University) is acknowledged for his assistance during the field data acquisition. I also thank the Geophysics Unit, School of Physics, Universiti Sains Malaysia, for providing adequate facilities and a congenial environment to conduct this research.

Data Availability
All data generated or analyzed during this study are included in this published article. Other supporting analyzed data can be made available by the corresponding author upon reasonable request.  showing the study area in the Nigerian Southwestern Basement Complex (modified after [28]). (c) Geological map of Akungba-Akoko and its surroundings in Ondo State, southwestern Nigeria (modified from [21]).