Geoelectrical parameter-based multivariate regression borehole yield model for predicting aquifer yield in managing groundwater resource sustainability

Peer review under responsibility of Taibah University.

Abstract

This study developed a GIS-based multivariate regression (MVR) yield rate prediction model of groundwater resource sustainability in the hard-rock geology terrain of southwestern Nigeria. This model can economically manage the aquifer yield rate potential predictions that are often overlooked in groundwater resources development. The proposed model relates the borehole yield rate inventory of the area to geoelectrically derived parameters. Three sets of borehole yield rate conditioning geoelectrically derived parameters—aquifer unit resistivity (ρ), aquifer unit thickness (D) and coefficient of anisotropy (λ)—were determined from the acquired and interpreted geophysical data. The extracted borehole yield rate values and the geoelectrically derived parameter values were regressed to develop the MVR relationship model by applying linear regression and GIS techniques. The sensitivity analysis results of the MVR model evaluated at P ⩽ 0.05 for the predictors ρ, D and λ provided values of 2.68 × 10⁻⁰⁵, 2 × 10⁻⁰² and 2.09 × 10⁻⁰⁶, respectively. The accuracy and predictive power tests conducted on the MVR model using the Theil inequality coefficient measurement approach, coupled with the sensitivity analysis results, confirmed the model yield rate estimation and prediction capability. The MVR borehole yield prediction model estimates were processed in a GIS environment to model an aquifer yield potential prediction map of the area. The information on the prediction map can serve as a scientific basis for predicting aquifer yield potential rates relevant in groundwater resources sustainability management. The developed MVR borehole yield rate prediction mode provides a good alternative to other methods used for this purpose.

Keywords:

• Multivariate regression yield model
• Groundwater reservoir yield
• Geophysical parameters
• GIS
• Hydrogeological

1 Introduction

The need for successful planning and management of groundwater exploration in the field of groundwater hydrology for both local and regional groundwater productivity potential mapping cannot be underestimated [1,2]. Nonetheless, existing literature has established that the success of area groundwater productivity potential assessment is largely based on the availability of well yield rate parameters [3,4]. The uniqueness of this parameter as a good, albeit indirect, indicator of groundwater occurrence in an area has been documented in the studies of [5–7]. A provisional means of assessing borehole yield information either locally or regionally is necessary to prevent the mismanagement of this invaluable but non-renewable natural resource, as this information can safeguard the sustainability of groundwater resources by avoiding saline intrusion, encrustation, borehole failures, and the excess lowering of water tables or piezometric surfaces. Consequently, the modelling of an aquifer yield potential rating through a simple and accurate mathematical technique can complement the conventional borehole pumping test (BPT) approach. This is necessary for enhancing the optimisation of groundwater resources that can provide portable water resources to address the challenges of scarce water supplies that threaten the entire world. The growing reliance on groundwater resources stem from unique attributes such as constant temperature, excellent quality, and low vulnerability to pollution and catastrophic events compared to surface water. Moreover, driving factors in various indispensable areas of human economic activity such as irrigation services and industrial utilisation often provide the basis for a nation to increase the availability of groundwater resources [8,9]. Therefore, the concept of enhancing the sustainability of groundwater resources through regional borehole yield rate assessment is timely and necessary.

The conventional approach of assessing borehole yield rate utilises a borehole pumping test (BPT) technique. This approach also provides information such as transmissivity, storability, aquifer geometry, and hydraulic conductivity properties that is useful when hydrogeologists are required to make accurate decisions concerning aquifers [10]. Ultimately, the evaluation of an aquifer's yield potential rate for the purpose of successful groundwater resource development can be best assessed by the BPT technique. However, the BPT approach is laborious, costly, uneconomical and time consuming [11–14]. In addition, the BPT output application is highly localised, thus limiting any regional evaluation of an area aquifer potential yield rating where there are no boreholes in existence. To gain better insight into regional subsurface aquifer potential yield ratings, which can be particularly useful in hard-rock terrains where the locations of groundwater reservoirs (aquifers) are varied and discontinuous, this study proposes a simpler and regionally compliant aquifer yield potential evaluation technique that can address the deficiencies of the BPT approach.

This study develops a GIS-based multivariate regression borehole yield prediction model. The approach employs empirical modelling that can accommodate the simultaneous integration of multiple factors for estimating an aquifer yield potential rating index. The efficiency of regression-based techniques in determining borehole yields from relevant borehole yield conditioning parameters obtained from different sources was established by [15,16]; these studies found that the well yield parameter was strongly correlated with well yield conditioning factors. Other studies [16–18] identified some of those well yield conditioning factors as hole depth, drawdown, screen length, geoelectrical parameters and geological rock formations. However, the regional applicability of this proposed regression model efficiency was carried out through a GIS technique application. The potential of GIS techniques has been explored in numerous environmental decision-making studies with encouraging results, including the management of natural resources—particularly in groundwater potential prediction domains [3,14,19]. Employing the proposed GIS-based MVR yield prediction model can significantly enhance the prediction of the potential yield of underlying groundwater reservoirs (aquifer units) where drilled holes are lacking. The proposed aquifer yield rating potential model output could optimise the sustainability of underlying groundwater resources.

This paper develops a GIS-based multiple variate regression (MVR) yield prediction model for predicting aquifer yield potential on a regional scale. This study introduces the use of derived parameters from geoelectric surveys obtained from an electrical resistivity (ER) prospecting method and applies linear regression and a GIS technique to develop the MVR yield prediction model. The ER method has widely been used for the quantitative estimation of the water transmitting properties of aquifers, aquifer zone delineation and the evaluation of the geophysical properties of aquifer zones in several locations [20–28]. Moreover, the ER method possesses non-invasive, low-cost attributes and offers quick data acquisition and the ability to map geological layers and determine the nature and composition of unseen subsurface formations [29,30]. The unique attributes of the ER method are well exploited in the proposed GIS-based MVR yield model and enable the model to surpass the conventional BPT technique. Applying the developed GIS-based MVR yield model in an area can provide aquifer yield potential rate prediction on a regional scale. The methodologies are illustrated by using a case study in the Crystalline Basement Complex terrain in southwestern Nigeria to establish a robust hydrological database system for continued effective groundwater resources management in the area.

2 Materials and methods

2.1 Geography, hydrology, and hydrogeology of the study area

The study area, as depicted in Fig. 1, is located in the northern part of Ondo state, southwestern Nigeria. It lies between the latitudes of 6°40′00″ and 7°45′00″ N and the longitudes of 4°20′00″ and 6°05′00″ E and contains 14 local government areas. The area is bounded to the east by Edo and Delta states, to the west by Ogun and Osun States, and to the north by Ekiti and Kogi States. The topography varies from moderate to high elevation values ranging from 200 m to 250 m. The area lies in a tropical rain forest with mean annual rainfall of approximately 1300 mm. The annual mean monthly temperature ranges from 6 °C to 33 °C. According to [31–33], the underlying regional geology generally contains rocks of the Precambrian Basement Complex of Nigeria. However, the typical petrology units recognised in the study area include granitic rocks (granite gneiss, med-coarse-grained biotite granite, fine-grained biotite granite, undifferentiated older granite, quartzo-feldspathic gneiss, coarse porphyritic biotite), migmatite-gneiss complex (migmatite, megmatite, undifferentiated schist), charnockitic meta-intrusive (charnockitic rocks) and quartzite series (fine-grained quartzite and schist quartz schist) (Fig. 1). These rocks, because of their peculiar complex crystalline basement characteristics, are generally characterised with low porosity and negligible permeability. However, where they are concealed both locally and regionally, these rocks may contain highly faulted and tightly folded areas, incipient joints and fracture systems derived from multiple tectonic events that can result in the development of secondary porosity. According to [34–36], these aforementioned geologic structures are localised and discontinuous. These geologic/structural features can significantly enhance the recharge rate of an underlying groundwater reservoir (aquifer unit).

Fig. 1 The study area map both the VES location, the geologic rock types map with inset map of Nigeria.

3 Methodology

Fig. 2 presents an overview of the approach that was applied for the regional assessment of the aquifer unit potential yield prediction mapping in the study area. The adopted methodological concept is multidisciplinary and divided into five phases. (1) The hydrogeological survey phase entails borehole yield rating extraction from the existing borehole database. (2) The hydrogeophysical phase involves geophysical data acquisition, processing and interpretations to determine both primary and secondary geoelectrical parameters. The GIS operation phase (3) of the study addresses overlaying and spatial modelling analyses. The inputs for this third phase derive from the results from the first and second phases. (4) The multiple variate regression (MVR) modelling approach utilises a linear regression technique and statistical software applications. The combined outputs of phase 2–4 constitute phase (5), which results in the modelling of regional borehole yield rate estimation and prediction in the study area. The details of the aforementioned phases (1)–(5) are highlighted below.

Fig. 2 The flowchart of the adopted methodology for the study.

3.1 The hydrogeological survey phase

An inventory of drilled boreholes was surveyed for the purpose of hydrogeological investigation. The inventory map of the boreholes drilled in the area is shown in Fig. 3; based on these borehole data, the hydrogeological condition of the area was analysed. The borehole pumping test technique was employed for the quantitative appraisement of each drilled hole and their yield ratings were determined. The adopted procedures of the borehole pumping test (BPT) technique were similar to those documented in [37]. The unpublished technical BPT report containing the details of the well drawdown level, the pumping duration, the rate, etc. is catalogued by Nigerian State water board agencies (WATSAN) and the Benin Development Authority, a parastatal Federal Ministry of water resources. The summary of the borehole yield rate records extracted from the existing BPT database in the area is presented in Table 1.

Fig. 3 The borehole distribution inventory map of the area.

Table 1 Borehole yield rating inventory records of the area.

BH location	Easting	Northing	Borehole yield records (l/s)
1	812,450	733,340	1.20
2	796,868	744,129	1.12
3	800,301	741,145	1.14
4	814,326	740,086	0.20
5	807,079	747,427	1.20
6	793,129	739,270	0.20
7	806,085	737,591	0.90
8	803,218	739,455	0.80
9	819,920	726,332	0.80
10	806,377	737,676	1.20
11	793,126	739,258	1.10
12	795,941	785,935	1.24
13	809,576	807,360	1.10
14	827,590	802,671	0.59
15	826,342	802,678	0.70
16	801,783	743,350	1.14
17	801,798	743,271	0.80
18	803,232	761,824	1.21
19	799,246	758,367	0.80
20	850,119	820,582	1.23
21	819,602	822,630	1.23
22	802,716	741,710	5.90
23	738,805	664,528	2.64
24	748,867	678,543	1.39
25	748,148	683,394	1.80
26	749,226	686,809	1.24
27	754,617	683,934	1.39
28	744,555	709,448	3.89
29	767,913	718,432	1.08
30	774,022	716,635	3.89
31	774,740	720,588	2.28
32	780,850	717,713	2.90
33	785,521	716,276	1.29
34	790,912	720,589	1.40
35	788,755	710,885	2.40
36	790,912	713,042	0.80
37	796,302	711,604	1.26
38	773,662	741,790	11.8
39	773,303	746,102	1.80
40	775,819	747,180	3.00
41	790,436	746,927	5.10
42	791,573	751,746	7.60
43	794,188	749,472	3.30
44	791,118	747,312	1.07
45	793,392	743,446	1.80
46	794,643	741,512	1.72
47	796,917	745,834	1.47
48	797,030	741,058	1.72
49	796,979	745,692	1.53
50	801,620	743,066	0.70
51	801,728	743,318	1.23
52	802,628	741,735	0.90
53	803,203	739,504	2.70
54	804,067	736,482	0.20
55	807,017	747,455	2.10
56	807,269	744,433	7.10
57	810,111	740,547	1.26
58	812,240	730,798	1.20
59	812,456	733,352	0.80
60	812,996	734,719	1.14
61	814,327	739,972	1.10
62	816,054	735,619	1.00
63	816,522	740,332	6.40
64	779,817	777,386	6.50
65	786,336	782,389	3.50
66	794,386	785,755	1.10
67	795,921	785,840	19.4
68	795,566	774,839	6.60
69	800,185	776,886	1.44
70	804,960	783,111	2.70
71	840,095	801,862	1.25
72	833,727	800,195	2.36
73	827,815	803,985	7.60
74	823,721	810,807	1.67
75	767,652	803,568	2.30
76	771,486	803,568	3.50
77	828,421	817,174	0.83

3.2 The geoelectrical method data acquisition approach

A one-dimensional vertical electrical resistivity sounding (VES) survey was carried out at 450 locations (Fig. 1) using a Schlumberger electrode configuration. The apparent resistivity values of the underlying soil formation were measured using an ABEM Terrameter (SAS 1000/4000 series). The electrode current spread length (AB) for the adopted Schlumberger array was varied between 2 and 200 m, whereas the spacing between the potential electrodes (MN) was intermittently varied between 0.5 m and 10 m to achieve suitable current penetration and a depth of investigation that will enable proper delineation of the subsurface lithological sequence to achieve a reliable estimate of aquifer parameters. The coordinates at each VES station were recorded with a Garmin handheld GPS device.

3.2.1 Data processing and interpretation technique

The acquired vertical electrical soundings data were processed by plotting the measured apparent resistivity values against half-current electrode spacing (AB/2) at each station on log-log graph sheets to generate the resistivity model curve. The typical resistivity model curves obtained based on the underlying geologic units in the area are shown in Fig. 4. The field curve types were further quantitatively interpreted using the conventional partial curve matching technique utilising the master curves and the corresponding auxiliary curves for the subsurface layer resistivity and thickness value determination. The estimated model parameters (layer resistivity values and thicknesses) were enhanced through computer iterations using a Win-Resist program. The outputs of both the qualitative and quantitative interpretations of the Schlumberger-VES data are called the geoelectrical parameters. The determined geoelectrical parameters were further grouped into two classes, namely the primary and secondary classes (see Fig. 2). The estimated layer resistivity (ρ) and layer thickness (D) comprise the primary class; the secondary class parameter is the coefficient of anisotropy (λ). The coefficient of anisotropy (λ) parameter is computed from the layer resistivity (ρ) and layer thickness (D) by applying the Dar–Zarrouk relationships documented in [38]. The λ is described as the ratio of the total transverse resistivity (ρ_t) to the total longitudinal resistivity (ρ_L) as expressed in Eq. (1)(1) $\lambda =\sqrt{\frac{{\rho }_t}{{\rho }_L}}$(1) , where ρ_t is defined as the ratio of the total transverse resistance (T) to the total thickness of the layers (H), i.e., T/H, whereas ρ_L is defined as the fraction of total thickness of the layers (H) with the longitudinal conductance (S), i.e., H/S [39]. The substitution of ρ_t = T/H and ρ_L = H/S in Eq. (1)(1) $\lambda =\sqrt{\frac{{\rho }_t}{{\rho }_L}}$(1) creates Eq. (2)(2) $\lambda =\frac{\sqrt{TS}}{H}$(2) .(1) $\lambda =\sqrt{\frac{{\rho }_t}{{\rho }_L}}$(1) (2) $\lambda =\frac{\sqrt{TS}}{H}$(2)

Fig. 4 Typical resistivity model curves obtained in the area.

3.2.2 The geoelectrically derived borehole yield potential conditioning parameters mapped with GIS approach

Hydrogeologically speaking, the measures of an area's well yield rate largely depend on the degree of aquifer productivity. By exploiting the aforementioned uniqueness of the ER geophysical prospecting method in Section 1, the delineation of the aquifer layer in the subsurface can be effectively mapped [40,41]. Thus, in situ aquifer physical factors such as aquifer layer resistivity (ρ), aquifer layer thickness (D) and coefficient of anisotropy (λ) can be determined based on the delineated aquifer layer by using the interpreted geophysical results (Fig. 4 and Table 2). According to [3,36,42], the aquifer resistivity and its thickness parameters are vital factors for assessing the groundwater potential in an area. However, the coefficient of anisotropy (λ) has been established to be linearly correlated with groundwater yield in an area [43]. Contrary to [19], the determined aquifer resistivity, its thickness and the coefficient of anisotropy parameters are considered as the borehole yield conditioning factors in the study area. The determined geoelectrical parameter (primary and secondary) values characterising the delineated aquifer layer at each VES location and the measured coordinates are presented in Table 2. Based on overlaying and spatial analysis, the determined aquifer parameter values (ρ, D and λ) referenced to each VES location are finally processed in a GIS environment in order to map the geoelectrically derived borehole yield potential conditioning parameter maps for the area (Fig. 5a–c).

Fig. 5 Borehole yield conditioning factors used in developing the multi-variate regression yield model where (a) aquifer resistivity map, (b) aquifer thickness map and (c) coefficient of anisotropy map.

Table 2 Summary of the interpreted geoelectric parameters.

VES NO	Northing	Easting	Ar	At	λ	CT
1	815,896.8	807,147.3	200	4.30	0.09	H
2	815,872.5	807,196.6	561	12.6	0.29	HA
3	801,509.2	746,343.9	296	4.10	0.06	H
4	801,527.6	746,334.6	179	3.90	0.04	H
5	797,873.3	740,475.3	136	6.00	0.04	H
6	823,907.7	800,671.3	147	8.80	0.03	H
7	823,895.6	800,702	239	14.9	0.07	H
8	823,892.6	800,711.2	111	5.40	0.02	H
9	803,404.7	743,741.5	245	2.40	0.02	H
10	803,401.6	743,741.5	958	7.10	0.74	KH
11	803,312.6	743,769.6	347	2.00	0.02	HA
12	803,525.7	739,530.7	63	13.70	0.04	H
13	803,516.5	739,543.1	150	13.90	0.04	HA
14	802,386.9	743,611.4	947	32.40	1.57	H
15	803,176.6	743,607.6	121	3.80	0.02	HA
16	803,136.5	743,571	128	5.60	0.02	H
17	803,146.4	743,712.1	311	8.30	0.09	H
18	803,167.3	743,589.2	157	5.00	0.03	H
19	803,177	743,687.4	326	7.60	0.03	HK
20	803,177	743,684.3	148	3.00	0.03	H
21	816,233.7	801,973.6	312	5.00	0.01	HA
22	816,187.6	801,983.1	194	4.20	0.03	H
23	800,485.3	777,998.8	373	3.70	0.08	H
24	800,516.4	778,063.1	299	11.10	0.05	H
25	800,444	742,147.8	362	23.60	0.12	H
26	796,914.4	743,634.9	51	4.80	0.02	H
27	807,934.3	743,808.3	584	14.10	0.15	HKHK
28	798,722.3	743,859.4	23	12.20	0.05	QH
29	819,450	727,582.5	31	0.90	0.02	H
30	796,715.4	784,288.2	13	1.70	0.02	H
31	796,718.5	784,294.4	80	7.90	0.03	QH
32	796,854.3	786,565.4	50	12.00	0.04	KH
33	793,601.7	746,753.9	72	12.00	0.07	KH
34	793,493.6	746,644	63	7.00	0.14	H
35	799,085.2	743,289.9	65	9.30	0.02	QH
36	799,050.8	743,167.3	105	15.60	0.05	H
37	797,093.6	741,927.6	832	11.10	4.20	A
38	799,003	740,239.8	160	5.00	0.09	KH
39	799,000	740,255.1	22	4.70	0.05	KH
40	798,214	742,944.2	89	14.20	0.04	HKH
?.	–	–	–	–	–	–
?	–	–	–	–	–	–
450			79	10.40	0.04	H

A_r = aquifer layer resistivity, A_t = aquifer layer thickness, λ = coefficient of anisotropy and C_T = curve types.

3.2.3 The geoelectric parameters–borehole yield relationship

To establish the influence of the above determined in situ physical aquifer parameters/factors on the extracted borehole yield rate records (Table 1), the relationships between the ρ, D and λ parameters and the records in Table 1 must be established. In examining the spatial analysis and the overlay functionality of the GIS technique, the actual borehole yield values (Table 1) and the corresponding interpreted and determined ρ, D and λ parameter values can be obtained by using the spatial model map of the determined geoelectrical parameters of the delineated aquifer units for each VES location in the study area (Table 2 and Fig. 5a–c). The results of the estimated borehole yield values and the corresponding values of the ρ, D and λ parameter values are presented in Table 3. The obtained results in Table 3 were used to generate linear graphs showing the relationship of the actual borehole yield values to the aquifer layer resistivity values, aquifer layer thickness, and the coefficient of anisotropy values. This was performed to establish the linear correlation of the individual parameters to serve as the basis for their use in the proposed multivariate regression model. Fig. 6a–c presents linear graphs showing the relationship between the actual borehole yield rate values and the corresponding ρ, D and λ parameter values. According to Fig. 6a–c, the high regression coefficients of 0.82, 0.73 and 0.87, respectively, establish the high degree of direct relationship between the borehole yield rate and the determined in situ aquifer geoelectrical parameters (ρ, D and λ).

Fig. 6 Linear relationship between the borehole yield rate records and the geoelectrical derived parameters where graphs (a)–(c) are resistivity, thickness and coefficient of anisotropy parameters models, respectively.

Table 3 Results of the geoelectric parameters and the corresponding actual borehole yield values.

BH Location	Easting	Northing	Geoelectrical parameters			Borehole pumping test analysis
			Aquifer resistivity (Ωm) [ρ]	Aquifer thickness (m)[T]	Coefficient of anisotrophy (λ)	The actual borehole yield results (l/s)
1	812,450	733,340	105	10.49	0.25	1.20
2	796,868	744,129	130	8.80	0.10	1.12
3	800,301	741,145	115	12.42	0.11	1.14
4	814,326	740,086	115	4.58	0.01	0.20
5	807,079	747,427	145	6.85	0.08	1.20
6	793,129	739,270	120	4.58	0.01	0.20
7	806,085	737,591	128	5.00	0.12	0.90
8	803,218	739,455	140	4.75	0.06	0.80
9	819,920	726,332	125	4.75	0.04	0.80
10	806,377	737,676	105	12.25	0.13	1.20
11	793,126	739,258	85	9.50	0.37	1.10
12	795,941	785,935	140	13.54	0.22	1.24
13	809,576	807,360	115	10.42	0.10	1.10
14	827,590	802,671	100	4.65	0.14	0.59
15	826,342	802,678	110	9.88	0.08	0.70
16	801,783	743,350	75	6.21	0.51	1.14
17	801,798	743,271	125	8.60	0.06	0.80
18	803,232	761,824	128	7.06	0.20	1.21
19	799,246	758,367	122	7.96	0.07	0.80
20	850,119	820,582	135	7.89	0.13	1.23
21	819,602	822,630	136	9.70	0.11	1.23
22	802,716	741,710	225	13.90	0.78	5.90
23	738,805	664,528	135	12.65	0.44	2.64
24	748,867	678,543	65	10.00	0.84	1.39
25	748,148	683,394	125	12.65	0.31	1.80
26	749,226	686,809	65	12.07	0.49	1.24
27	754,617	683,934	95	9.80	0.38	1.39
28	744,555	709,448	170	13.54	0.59	3.89
29	767,913	718,432	105	9.82	0.22	1.08
30	774,022	716,635	168	13.54	0.60	3.89
31	774,740	720,588	125	12.42	0.42	2.28
32	780,850	717,713	156	12.73	0.39	2.90
33	785,521	716,276	132	12.22	0.11	1.29
34	790,912	720,589	75	12.25	0.42	1.40
35	788,755	710,885	148	12.32	0.42	2.40
36	790,912	713,042	122	8.75	0.06	0.80
37	796,302	711,604	106	12.00	0.32	1.26
38	773,662	741,790	380	17.50	1.43	11.8
39	773,303	746,102	118	11.88	0.30	1.80
40	775,819	747,180	148	13.13	0.45	3.00
41	790,436	746,927	240	13.56	0.65	5.10
42	791,573	751,746	320	15.44	1.23	7.60
43	794,188	749,472	152	13.22	0.49	3.30
44	791,118	747,312	131	7.28	0.13	1.07
45	793,392	743,446	115	11.70	0.31	1.80
46	794,643	741,512	120	11.70	0.28	1.72
47	796,917	745,834	112	11.60	0.25	1.47
48	797,030	741,058	121	10.90	0.28	1.72
49	796,979	745,692	110	11.84	0.27	1.53
50	801,620	743,066	109	8.00	0.10	0.70
51	801,728	743,318	116	10.20	0.18	1.23
52	802,628	741,735	105	5.77	0.18	0.90
53	803,203	739,504	135	13.90	0.46	2.70
54	804,067	736,482	100	4.58	0.08	0.20
55	807,017	747,455	127	12.07	0.37	2.10
56	807,269	744,433	342	15.44	1.20	7.10
57	810,111	740,547	230	10.20	0.11	1.26
58	812,240	730,798	165	11.00	0.18	1.20
59	812,456	733,352	185	8.50	0.24	0.80
60	812,996	734,719	259	10.99	0.14	1.14
61	814,327	739,972	127	10.39	0.14	1.10
62	816,054	735,619	247	9.90	0.12	1.00
63	816,522	740,332	340	13.94	0.84	6.40
64	779,817	777,386	380	15.21	1.18	6.50
65	786,336	782,389	161	13.56	0.51	3.50
66	794,386	785,755	143	9.00	0.01	1.10
67	795,921	785,840	456	20.47	1.70	19.40
68	795,566	774,839	390	15.22	1.90	6.60
69	800,185	776,886	198	11.50	0.24	1.44
70	804,960	783,111	225	13.90	0.5	2.70
71	840,095	801,862	106	13.22	0.16	1.25
72	833,727	800,195	207	10.48	0.43	2.36
73	827,815	803,985	320	15.48	1.22	7.60
74	823,721	810,807	134	9.80	0.26	1.67
75	767,652	803,568	245	12.21	0.40	2.30
76	771,486	803,568	270	13.56	0.52	3.50
77	828,421	817,174	125	7.04	0.04	0.83

3.3 Multiple variate regression (MVR) yield model

Consider the following generalised multiple variate regression model:(3) $Y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+\cdots +{\beta }_n{x}_n+{\ell }_i$(3) where β₀ is the intercept; β₁,β₂ and β₃ are the slopes of the regression line for x₁, x₂ and x₃ which are the independent (predictor) variables, respectively; ℓ_i is the error term; and Y is the dependent variable (response) as reported in [44,45]. Eq. (3)(3) $Y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+\cdots +{\beta }_n{x}_n+{\ell }_i$(3) is an example of a linear regression model where the dependent variable Y is linearly related to the predictor variables x₁, x₂ and x₃; however, the relationship between these variables is not exact and is subject to individual variation.

In this study, the borehole yield rate is the dependent variable (Y) and the geoelectrical parameters (ρ, D and λ) are the independent (predictor) variables (x₁, x₂ and x₃). The linear relationship between the borehole yield rate (Y) and the individual dependent variables ρ, D and λ has been established as shown in Fig. 6a–c. However, because the relationship between the dependent and independent variables is not exact, an error term ℓ_i is introduced to accommodate any individual variation. Eq. (3)(3) $Y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+\cdots +{\beta }_n{x}_n+{\ell }_i$(3) can be written as Eq. (4)(4) $Y={\beta }_1(\rho )+{\beta }_2(D)+{\beta }_3(\lambda )+{\beta }_0+{\ell }_i$(4) when x₁, x₂ and x₃ are substituted with ρ, D and λ, respectively;(4) $Y={\beta }_1(\rho )+{\beta }_2(D)+{\beta }_3(\lambda )+{\beta }_0+{\ell }_i$(4) The estimation of coefficients β₀, β₁, β₂ and β₃ were determined through integrative model regression analysis of the records in Table 3 using R software. Based on the integrative model regression analysis, β₁ = 0.02, β₂ = 0.15, β₃ = 3.52, β₀ = −2.6 and the R-square value = 0.84. By substituting the estimated results, Eq. (4)(4) $Y={\beta }_1(\rho )+{\beta }_2(D)+{\beta }_3(\lambda )+{\beta }_0+{\ell }_i$(4) can be modified as Eq. (5)(5) $\hat{Y}=0.02(\rho )+0.15(D)+3.52(\lambda )-2.6$(5) ; the carat over Y indicates that it is an estimate and the error term ℓ_i is thus automatically resolved.(5) $\hat{Y}=0.02(\rho )+0.15(D)+3.52(\lambda )-2.6$(5) Eq. (5)(5) $\hat{Y}=0.02(\rho )+0.15(D)+3.52(\lambda )-2.6$(5) is a multiple variate regression (MVR) equation with the borehole yield rate (Y) as the dependent variable and (ρ), (D) and (λ) as the multiple independent variables. Based on the submission of [46] reported in [27], Eq. (5)(5) $\hat{Y}=0.02(\rho )+0.15(D)+3.52(\lambda )-2.6$(5) is referred to as a Model. Thus, Eq. (5)(5) $\hat{Y}=0.02(\rho )+0.15(D)+3.52(\lambda )-2.6$(5) is the proposed multivariate regressions (MVR) yield prediction model developed for the study area.

4 Results and discussion

4.1 Sensitivity analysis of the developed MVR prediction yield model

According to [47] and as cited by [48], a sensitivity analysis is required that involves an analysis of the contribution of individual variables and input parameters on the resultant output of an analytical model. Adopting this approach, the developed MVR yield prediction model was evaluated for a parameter significance assessment. This analysis was carried out using R statistical software. Table 4 presents the parameter significance evaluation results of the developed MVR yield prediction model. The results in Table 4 show that the evaluated aquifer layer resistivity, the aquifer layer thickness, and the coefficient of anisotropy parameters have a significant relationship to the response variable borehole yield rate (Y) at P ≤ 0.05 (5%) in the examined area. This implies that the considered borehole yield rate conditioning factors (ρ, D and λ) have high significance at probability (P) ≥ 95% for estimating borehole yield rate (Y) using the MVR yield rate prediction model (Eq. (5)(5) $\hat{Y}=0.02(\rho )+0.15(D)+3.52(\lambda )-2.6$(5) ). It can thus be implied that the integrative modelling of these aforementioned geoelectrical derived borehole yield rate conditioning parameters with their varying influences are appropriate for estimating and predicting borehole yield rate in the study area. The used borehole yield rate conditioning parameters (ρ, D and λ) in the developed MVR prediction yield model are easily derivable from simple geophysical measurement and thus make Eq. (5)(5) $\hat{Y}=0.02(\rho )+0.15(D)+3.52(\lambda )-2.6$(5) a legitimate model for aquifer yield potential determination. Therefore, the estimation and prediction of an area groundwater reservoir potential yield can be reliably determined from knowing these geoelectrical parameters, even in areas without boreholes.

Table 4 Parameters estimations analysis of the developed multiple variate regression (MVR) yield model developed in the area.

Developed yield model	Parameters	Standard error and parameters significance testing using standard value of α at 5% significance level		Remark: parameters significance OK at Pr-value <5%
		t-Value	Pr (>〉|t|)-value
Y = 0.02(ρ) + 0.15(D) + 3.52(λ) − 2.6	ρ	4.49	2.68 × 10^−05	OK
	D	2.38	0.02	OK
	λ	5.15	2.09 × 10^−06	OK

(ρ): aquifer layer resistivity (Ωm); (D): aquifer layer thickness (m) and λ: coefficient of anisotropy.

4.2 Appraisal of model prediction accuracy

Furthermore, the predictive power of the developed MVR yield prediction model was appraised to determine the feasibility of using the model to predict and estimate the area's groundwater reservoir potential yield. [44,49] suggested a systematic measure of accuracy for any forecast obtained from a model. This measure is called the Theil inequality coefficient, which is given by(6) $K={\sum _{l=1}^n\left(\frac{y_i-\hat{y}}{\hat{y}}\right)}^2$(6) where y_i is the actual determined borehole yield rate observed in the area, $\hat{y}$ is the corresponding estimated borehole yield rate of y_i from the MVR yield prediction model (see Table 5) and K is the Theil inequality coefficient. Eq. (6)(6) $K={\sum _{l=1}^n\left(\frac{y_i-\hat{y}}{\hat{y}}\right)}^2$(6) was applied to assess the estimation and prediction accuracy of the MVR prediction yield model. The determined Theil inequality coefficient K value was then gauged by the critical value of ${\chi }_p^2$, α, where P = n − 1, n is the number of occupied borehole locations, and α is the 5% significance level. A smaller value of K compared with the χ²-tabulated value indicates better prediction accuracy of the model under investigation [50]. The accuracy appraisal result of the MVR yield model is shown in Table 6. In addition, the regression line fitted to the actual borehole yield rate and the estimated borehole yield rate based on the MVR prediction yield model data was generated (Fig. 7), and the result shows a high regression coefficient of 0.88 for the linear relationship. These results confirmed the reliability and accuracy of using the MVR yield rate prediction (Y) model for predicting and estimating borehole yield rate in the non-investigated part of the area. Therefore, the output of this yield rate estimate using this developed MVR yield rate prediction model can be harnessed for groundwater resource evaluation and management in the study area.

Fig. 7 Linear relation between the actual borehole yield rate and the estimated borehole yield rate based on the MVR prediction yield model.

Table 5 The records of the actual borehole yield determined and the predicted borehole yield in the area.

BH location nos	Actual borehole yield determined in the area (yi)	Estimated borehole yield rate from the MVR prediction yield model $(\hat{y})$
1	1.20	1.23
2	1.12	1.06
3	1.14	1.08
4	0.20	0.10
5	1.20	1.13
6	0.20	0.20
7	0.90	0.78
8	0.80	0.79
9	0.80	0.80
10	1.20	0.94
11	1.10	1.16
12	1.24	1.25
13	1.10	0.89
14	0.59	0.26
15	0.70	0.67
16	1.14	1.19
17	0.80	0.80
18	1.21	1.21
19	0.80	0.72
20	1.23	1.19
21	1.23	1.23
22	5.90	5.76
23	2.64	2.64
24	1.39	1.40
25	1.80	2.00
26	1.24	1.25
27	1.39	1.42
28	3.89	3.90
29	1.08	1.06
30	3.89	3.96
31	2.28	2.37
32	2.90	2.91
33	1.29	1.30
34	1.40	1.36
35	2.40	2.82
36	0.80	0.75
37	1.26	1.61
38	11.8	11.43
39	1.80	1.77
40	3.00	2.99
41	5.10	5.10
42	7.60	7.40
43	3.30	3.22
44	1.07	1.06
45	1.80	1.73
46	1.72	1.65
47	1.47	1.45
48	1.72	1.68
49	1.53	1.50
50	0.70	0.57
51	1.23	1.17
52	0.90	0.60
53	2.70	2.83
54	0.20	0.05
55	2.10	2.21
56	7.10	9.70
57	1.26	3.20
58	1.20	2.21
59	0.80	2.62
60	1.14	3.95
61	1.10	1.26
62	1.00	3.55
63	6.40	8.27
64	6.50	10.37
65	3.50	3.50
66	1.10	1.02
67	19.40	19.40
68	6.60	13.11
69	1.44	3.12
70	2.70	4.77
71	1.25	1.14
72	2.36	3.89
73	7.60	9.33
74	1.67	1.78
75	2.30	4.68
76	3.50	5.72
77	0.83	0.60

Table 6 MVR yield model prediction accuracy analysis.

S/N	Proposed MVR yield model	Nos of borehole locations	${\chi }_p^2,\alpha =5\%$	K-value
1	Y = 0.02(ρ) + 0.15(D) + 3.52(λ) − 2.6	77	57.79	0.15

(ρ): aquifer layer resistivity (Ωm); (D): aquifer layer thickness (m) and λ: coefficient of anisotropy.

4.3 Spatial modelling of the borehole yield prediction map

The application of the MVR yield prediction model was carried out using extracted parameter values from the spatial model interpreted geoelectrical parameters (ρ, D and λ) as shown in Fig. 6a–c. The results of the GIS-based spatial analysis extraction processed with Eq. (5)(5) $\hat{Y}=0.02(\rho )+0.15(D)+3.52(\lambda )-2.6$(5) is presented in Table 5. The estimated yield rates based on the MVR yield prediction model (Table 5) were processed in a GIS environment using a geostatistical interpolation Kriging technique to spatially model the aquifer yield rate potential of the area. To classify the area aquifer potential yield prediction into zones, the area yield rate classification information obtained from the state water board agency, where yield rate values range from 0.2 l/s to 1 l/s, >1 l/s, <10 l/s, and >10 l/s, are apportioned to these actual yield descriptions for the low rate, the medium rate, and the high rate, respectively. Using the above criteria, the natural jerk classification algorithm embedded in a geostatistical wizard module was employed in a GIS environment to classify the area into five aquifer yield potential predicted zones: Low (0.06–2.55/s), Low medium (2.55–4.13 l/s), Medium (4.13–5.71 l/s), Medium high (5.71–7.74 l/s) and High (7.74–19.4 l/s). The final aquifer yield spatial prediction model map of the area is shown in Fig. 8. Furthermore, the area of extended coverage by these classified zones are 5318 km², 2869 km², 1786 km², 1064 km² and 186 km², respectively.

Fig. 8 The borehole yield rate prediction model map based on the MVR borehole yield model developed.

4.4 Borehole yield model validation

The produced aquifer yield rate potential prediction map model was validated by adopting both qualitative and quantitative approaches. According to [3,51], model validation enables the reliability assessment of any proposed model and increases its usefulness in environmental decision-making studies. The indices used for the prediction map model validation in this study include the drain patterns (river/stream) and geologic lineament features. Through overlaying analysis, these surface geologic features were spatially correlated with the predicted yield rate potential zones (Fig. 8). Quantitatively, the computation of the density values for both the drain patterns and the lineament features were cross-tabulated with the predicted borehole yield classified zones to provide the results presented in Table 7. According to Table 7, across the predicted borehole yield zones, i.e., from Low to High, the mean river density varies from 0 to 0.0503 km⁻¹. The low borehole yield rate classified zone is relative, with the highest river density values reducing downward to the other classified zones (Low medium, Medium, Medium high and High). Geologically, this suggests that the lithology underlying these zones possess varying porosity and permeability degrees that could be responsible for the varying degree of river drains in the area. Because the groundwater reservoir potential yield rate largely depends on river drain percolation resulting in recharge rate enhancement [52], the possible predicted groundwater reservoir yield potential rate zones are generally structurally controlled and thus validated. However, the lineament intersection density (LID) features computed from the lineament features have mean values in the range of 0–0.009 across the classified aquifer unit potential yield zones (Table 7). The predicted high aquifer yield potential zone is characterised by the highest value of LID analysis compared to the other zones. Because high aquifer potential yield is synonymous with high groundwater potential, then the findings of [47,48] that established high LID is associated with high groundwater recharge potential thus substantiates the LID result analysis (Table 7). Therefore, the analysed LID result further validates the predicted borehole yield zones (Fig. 8). As a result, the produced borehole yield prediction model map (Fig. 8) is a viable tool for monitoring the assessment of groundwater quantity potential that can enhance groundwater resource sustainability in an area.

Table 7 The GIS based cross-tabulation analysis results for the model map validation.

Borehole yield zones predicted	Drains/rivers density		Lineament intersection density
	Range	Mean (km^−1)	Range	Mean
Low	0–0.1595	0.0503	0–0.0298	0.0040
Low medium	0–0.1329	0.0481	0–0119	0.0053
Medium	0–0.1064	0.0287	0–0.0179	0.0055
Medium high	0	0	0	0
High	0–0.774	0.0387	0.006–0.0119	0.0090

5 Conclusions and future works

This study adequately evaluated the mapping and assessment of an area with underlying aquifer unit yield potential on a regional scale through the output of a GIS-based MVR yield prediction model. This newly proposed borehole yield rate prediction model was based on relating the records of in situ borehole yield rate measurements to geoelectrically derived parameters interpreted from one-dimensional VES-Schlumberger soundings acquired in the Crystalline Basement Complex terrain of southwestern Nigeria. Sensitivity and prediction accuracy analyses of the newly proposed aquifer yield potential prediction model using χ² distribution at an α = 0.05 significance level was conducted using R statistical software. The MVR yield rate potential prediction model was used to estimate borehole yield values that were spatially processed in a GIS environment to produce a groundwater reservoir potential yield prediction model map of the area. The regional aquifer yield potential model prediction map provided an excellent insight into assessing the viability of the possible underlying aquifer unit yield to the groundwater resource sustainability in the area. The information on this prediction map can serve as a scientific basis for groundwater resource exploration and management in the area. Furthermore, the proposed MVR yield prediction model that was developed with variables from multifaceted geologic settings can be used in any area with similar geology for groundwater resource potential evaluation if the required geoelectrical parameters are known.

Compared with other borehole yield assessment methods, the approach used in this study can provide a quick, independent, and cost-effective assessment of aquifer unit yield potential rates through simple geophysical measurement. However, the efficiency of this developed aquifer yield prediction model can be greatly enhanced by deriving the required borehole yield conditioning parameters from a 2D or 3D resistivity imaging technique. This can enable an accurate modelling of subsurface geological features, from which enhanced borehole yield conditioning parameters can be derived for re-evaluating this model in the future.

Notes

Peer review under responsibility of Taibah University.