Constructing a 3D geological model from geophysical data for groundwater modelling and management in the Kribi-Campo sedimentary sub-basin, Cameroon

Geological exploration for water resources, geothermal energy, oil, gas and minerals, and infrastructure projects have not been particularly successful in most African countries. Since such resources are needed for some degree of self-sufficiency and prosperity (growth), it is necessary to understand the geometry of the subsurface using 3D geological modelling (Hassen et al. 2016). 3D modelling is useful and has been used increasingly to integrate available datasets, leading to more realistic representations of given geological settings (Wu et al. 2005).

Water is vital and usually in limited supply for its innumerable uses, the most important of which are domestic purposes. Of earth's water resources, 97.2% are saline (mainly in oceans), and only 2.8% fresh. Groundwater accounts for about 26% of global renewable water resources (Elbeih 2015). Groundwater is the principal source of domestic water for many communities and the most extracted raw material globally. Withdrawal rates are currently estimated at about 982 km³/a (Senthilkumar & Rajakumar 2014). For groundwater resources in any area to be managed and monitored effectively, a good understanding of the geological setting is required. This will enable quantification of the resource and determination of its quality, which is important because of the potential effects of drought, floods, agriculture, urbanization, saline intrusion, and climate change.

A geological model, a spatial representation of the subsurface distribution of sediments and rocks, is traditionally presented in 2D cross-sections; for example, geological maps and profiles, and isopach maps showing the horizontal extent and thickness of subsurface structures (Andersen et al. 2018). 3D models are used increasingly (Sausse et al. 2010; Høyer et al. 2015; Hassen et al. 2016; Malquaire et al. 2017), and are constructed by interpolation from sparse geological and/or geophysical data. Geological and hydrogeological studies in urban areas are also compounded by the presence of anthropogenic deposits, making the subsurface difficult to predict (Andersen et al. 2018).

Collecting the geophysical datasets applied to characterize urban areas is time consuming, and expensive in terms of the survey team and resources mobilized. The methods used include GPR (ground penetrating radar), EM (electromagnetism), DC (direct current) and seismic reflection profiles (Andersen et al. 2018). Simple and robust technologies are now available that save time and money in data collection, while yielding a high-density geophysical dataset. One of these is the electro-seismic (ES) method (Mikhailov et al. 1997). Electro-kinetic prospecting methods are especially promising for detecting materials whose electrical resistivity varies markedly from the background. EM methods may also be used for this, but their spatial resolution is usually less than desirable because of EM's long wavelengths in the earth. Seismic methods have better spatial resolution but respond to smaller contrasts in material properties (White 2005). ES combines the virtues of both approaches – electro-kinetic and seismic – in the context of this study. The datasets used to develop the geological model in this study were collected using the ATS GeoSuite Application installed on a mobile handset and operating under the premise of the ES method (ATSGeoConsultants 2019).

The Kribi-Campo Sedimentary Sub-Basin aquifer system presents many challenges to geological interpretation, mainly because of the limited availability of data in most areas. No geological model has ever been developed for the area, geared towards groundwater resource management. Extensive geological data exist on the area, however, but are not available for academic research (Ntamak-Nida et al. 2010), because they relate to both on- and off- shore petroleum exploration. The lack of sufficient knowledge about the geological and hydrological settings introduces significant uncertainties in projections of the fates of contaminants and public health (Andersen et al. 2018).

Some 250 ES sounding points were surveyed, with an average inter-point profile grid spacing of 200 m inland and variable along the shore. This low-resolution grid spacing was used to provide representative coverage of the area, to extract the hydrological structures and parameters. Some survey points were difficult to access due to the terrain and, in general, it was not possible to use regular, straight grid lines because soundings had to be made along existing tracks and open fields.

The aim of this work was to develop a straightforward procedure for collecting geophysical datasets, and developing geological and hydrogeological models, to enable assessment of aquifer geometries and volumes, and defining their connectivity. This should enable answers to be found for questions concerning water availability, management and monitoring, especially in an urban setting.

Kribi town is in Kribi 2 Sub-Division, of the Ocean Division of the South Region of Cameroon. It is between latitudes 2.940 and 3.007° N, and longitudes 9.90 and 9.95° E, and covers about 13 km². The Atlantic Ocean bounds it to the west – Figure 1. The area is urban and has several beautiful beaches along approximately 187 km of coastline, which attract many tourists (Nfomou et al. 2004). The local aquifer system in the coastal area is very shallow, as indicated by the depth to the water table in some hand-dug wells. The area has an equatorial climate, so annual rainfall is sufficient. The top soil varies little.

The Kribi-Campo sub-basin comprises an early to mid-Cretaceous series from West Africa's Atlantic coast and is in the Central African equatorial rain forest in Southern Cameroon. It is the smallest coastal basin in Cameroon and forms the southern part of the Douala/Kribi-Campo basin (Ntamak-Nida et al. 2010). The region is tectonically active, and several earthquakes of moderate magnitude (3.4–5.7) have been felt and reported, indicating that the region could be subject to seismic risk and also that the Kribi-Campo Fault (KCF) that crosses it might be active (Owona Angue et al. 2011). The KCF seems to trend in continuation of the offshore Kribi Fracture Zone, which has been inferred as the eastern end of the Ascension Fracture Zone (Ntamak-Nida et al. 2010; Malquaire et al. 2017; Nfor et al. 2018).

The Kribi-Campo sub-basin is sedimentary (Figure 2), and lies across the Cameroon coast between latitudes 2°20′ and 3°20′N, and longitudes 9°15′ and 10°00′E. It covers about 6,195 km² (Malquaire et al. 2017). The Ntem Complex (Archaean) basement is overlain by Paleoproterozoic rocks of the Nyong Unit. The basement is composed mostly of greenstone belt rocks, Congo Craton charnockites and potassic granitoids (Owona Angue et al. 2011). The sub-basin is bordered to the south by the Campo high, which separates it and the Rio Muni Basin (Iboum Kissaaka et al. 2011). The sub-basin shows three classical stages in its structure: syn-rift (late Jurassic – early Cretaceous); transition (Aptian-Albian) and post-rift (Albian-Recent) (Ntamak-Nida et al. 2010; Iboum Kissaaka et al. 2011).

Electromagnetic (EM) and mechanical waves are coupled through electro-kinetics, in a porous medium such as the earth's subsurface (White 2005). The pore fluid ions are attracted to ions of the opposite sign in the solid at the pore walls, so that there is an electrical double layer (a thin charged layer coating the mineral grain surfaces) at the pore boundary – the Debye Layer. An electric field acting on this layer moves the ions relative to one another, causing movement of both the fluid and the solid. Conversely, a mechanical wave that moves the fluid and solid relative to each other will create an EM wave (White 2005). Techniques based on electro-kinetics have been proposed for use in sub-surface exploration and field experiments conducted. The experiments were variously seismo-electric and electro-seismic. In seismo-electric experiments, a mechanical (seismic) source is used and an EM wave detected, and vice versa in an electro-seismic experiment (White 2005). While such phenomena can be caused by various different physical mechanisms, they are referred to collectively as ‘electro-seismic’ (Mikhailov et al. 1997).

ES phenomena in porous media have the potential to enable detection of subsurface zones of high fluid mobility and fluid chemistry contrasts, supported by observation, with explicit comparison to full waveform modelling results (Mikhailov et al. 1997). A seismic wave propagating in a medium can induce an electric field or generate an electromagnetic wave.

The ES effects employed (ATSGeoConsultants 2019) look at results of a seismic wave crossing or travelling along the interface between two porous media. When a seismic source crosses an interface, there is ES conversion (Figure 3(a)). When a spherical P-wave crosses an interface, however, the imbalance of the streaming currents induced by the seismic wave on opposite sides of the interface creates a dipole charge separation. The induced dipole radiates an electromagnetic wave, which can be detected by remote antennae (Mikhailov et al. 1997; ATSGeoConsultants 2019).

When a seismic head wave travels along an interface between media (Figure 3(b)), it causes charge separation across the interface, which induces an electric field that moves along the interface with the head wave and can be detected by ground-dipole antennae when the head wave passes under them (Mikhailov et al. 1997; ATSGeoConsultants 2019). The electric fields are transformed into a set of time-varying potential differences, which are transmitted through the signal cables to the mobile device running the GeoSuite Application, and recorded and stored for processing (ATSGeoConsultants 2019).

ES equipment consists of a single mobile seismic source – any form of seismic impulse generator, including a low energy hammer and plate, a medium energy drop-weight, or a high-power explosive source such as a buffalo gun. The seismic source is placed over the field point and an operating area marked around it, and two stainless steel pins inserted into the ground as the ground–dipole antenna. A signal cable connected to the pins enables recording and is calibrated to the mobile device's audio port. The seismic waves are produced by striking the ground, a process repeated from 5 to 20 times at each location (referred to as data stacking).

Data processing includes filtering, extraction, correlation, stacking, and quality evaluation, all of which are completed automatically to produce up to 25 different datasets. The datasets include both hydraulic conductivity tomography (ESKT) and interface tomography (ESIT), which enable aquifer characterization, as well as others. Amongst the parameters used to determine the quality of the datasets are the correlation indexes. The spatial distribution of data points with their correlation indexes, which describe their reliability, is shown in Figure 4. In a study involving drilling, points with correlation index values below 50% would probably be re-drilled. As there was no drilling in this study, no repeated data collections were carried out at data points with correlation index below 50%. Instead, only those with correlation index values exceeding 50% were analyzed further. Water table depths were obtained from 61 hand-dug wells and shallow boreholes in the survey area, and the resulting dataset used to construct a ‘layer’ comprising the probable phreatic aquifer.

The ESKT dataset effectively describes the aquifer systems under a site in 1-, 2-, or 3- D. The ESKT data are normalized to give relative hydraulic conductivity values ranging from 0 to 100%. In other words, the data are not calibrated to an absolute field reference point, at which the hydraulic conductivity with depth is known. The maximum normalized value is obtained by assigning 100% to the maximum recorded hydraulic conductivity value and normalizing all other values against this (ATSGeoConsultants 2019), thus, identifying aquifers' locations. The resulting normalized hydraulic conductivity values were classified and each class was then assigned a given lithology that could possibly correspond to the hydraulic conductivity values of that class. This classification is shown in Table 1, together with the probable lithologies used to construct a lithologic model. The classified dataset was then imported into Leapfrog Geo software (ARANZ 2018), and the aquifer model developed.

The ESIT dataset is obtained by processing ES data, and indicates the positions of interfaces between formations with differing electrical and elastic properties. The interfaces are detected on the basis of the amplitude of the seismic reflection coefficient, and all interfaces are detected irrespective of its amplitude. This highlights geological features that might be missed if viewed in the context of reflector strength. The data show the reflector's total gradient polarity (Table 2). Blue responses indicate negative polarity and red positive polarity reflectors (ATSGeoConsultants 2019). Visual inspection and the resulting combination of some interfaces/media led to development of the geological model.

The datasets analyzed were converted to resemble drilling data and imported into the Voxler software (Golden Software 2019b). The water table surface grid was prepared in Surfer 16 software (Golden Software 2019a) and also imported into Voxler. The resulting model shows a system that fits together well, comprising a shallow unconfined aquifer along that length of coast; that is, the datasets reflect the geological and hydrogeological formations (see Figure 5).

The geological model component of Leapfrog Geo software incorporates five sub-components – (model) boundary, fault system, lithologies, surface chronology and output volume. The boundary inherits the clipped boundary set at topography level, which defines the upper boundary automatically. The fault system and lithologies sub-components are self-explanatory, although it is noted that the lithology field can be over-ridden manually if the dataset contains no lithologic information. The surface chronology component deals with all contact surfaces, which are arranged in chronological order, youngest first. It is these surfaces and their chronology that determine how the volume inside the model is divided into lithologic units. The output volumes comprise all those generated in building the model, in chronological order from youngest to oldest.

Leapfrog Geo can create a number of different surface types. The stratigraphic sequence was chosen for two reasons. Firstly, it should be used when there is a series of continuous layers that can be modeled separately. Secondly, it is known that this is a sedimentary basin, so that the lithologies can be modelled as layered horizontal blocks. Figure 6 shows the datasets in Leapfrog Geo and Figure 7 the ensuing geological model.

An aquifer system was developed from the geological model. The unconfined aquifer modelled is shown in Figure 8 as the Upper Sandy Aquifer. The latter was compared with that generated from the water table grid – see Figure 9 – and the two layers generally coincide, with minor discrepancies in regions of extreme extrapolation. The final aquifer system determined for the study area is shown in Figure 10.

Two fence sections were drawn in Leapfrog Geo to show the lithologic variations within the geological model. They run along the landward side of the study area (Figure 11(a) and 11(b)) and the other across the area towards the ocean (Figure 12(a) and 12(b)).

The conceptual system developed using Leapfrog Geo was converted into surfaces and exported as Surfer grid files to build the hydrogeological model in Visual MODFLOW Flex software (Hydrogeologic Waterloo 2015). The latter, a basic high-level representation of the hydrogeological system, will form the basis of other numerical models developed in Visual MODFLOW Flex, since it is grid and simulator independent.

The flow condition chosen for the groundwater model was saturated (constant density) flow, with no transport objective defined. After the model structure had been defined, the model's property zones were defined using the synthetic hydraulic conductivity values given in Table 3. The hydrogeological model is shown in Figure 13.

The idea behind this study was conceptualization of the procedure for generating hydrogeological models, based on ES dataset collection. Construction of a hydrogeological model starts with site and data acquisition point definition, data collection and processing, and using software to simulate the 3D hydrogeological system; for example, Leapfrog Geo and Visual MODFLOW Flex. This method of geophysical exploration is easily implemented in mobile devices and the groundwater exploration easily handled, so good understanding of the procedure will limit the time taken for groundwater studies. Although this model was not calibrated using down-hole log data, the high-density geophysical sets gave coherence in interpolating the surfaces, which could be helpful in groundwater modelling.

In this study the attempt was made to show how easy it can be to carry out groundwater studies. Due to the maiden nature of the work in general and the paucity of data in the study area, the work suffered from shortcomings that included, particularly:

• The model domain required further refinement by re-evaluating all sounding points with correlation index values (Figure 4) below 50% – they were excluded from analyses – to further reduce interpolation errors;

• Borehole data should be used to validate the ES dataset in the area. In other words, the lithology obtained from the datasets should conform with that determined from drilling logs and hydraulic conductivities. If this is done, the method can be reliably exported to other areas.

The datasets generated and/or analyzed during this study are available in the 4TU.Center for Research Data repository, http://researchdata.4tu.nl/.