Using an Open-Source Tool to Develop a Three-Dimensional Hydrogeologic Framework of the Kobo Valley, Ethiopia

Abstract

Groundwater resource management requires understanding the groundwater basin’s hydrogeology and would be improved with the development of a three-dimensional hydrogeologic framework model (HFM). A wide range of methods and software exist to quantify the extent, structure, and properties of geologic systems. However, most geologic software is proprietary and cost-prohibitive for use in developing countries. GemPy is a Python-based, open-source (no-cost) tool for generating three-dimensional geological models. This study uses available data and GemPy to develop the Kobo Valley Hydrogeologic Framework Model (KV-HFM), a three-dimensional HFM for Kobo Valley in northern Ethiopia, which is part of the East African Rift System. The KV-HFM is a conceptual model that comprises the hydrostratigraphy, structural features, and hydraulic properties of the Kobo Valley groundwater system. The limited data described the extent and altitude of the hydrostratigraphic units using the GemPy implicit potential–field interpolation. The KV-HFM showed the existence of an east-to-west, structural-based groundwater divide composed of volcanic rock and clay. This divide splits the catchment into two groundwater systems with limited interconnected flow. This study illustrates the use of open-source software for developing an HFM using sparse, existing geologic data.

1. Introduction

Groundwater is a crucial resource for urban and domestic water supply, irrigated agriculture, industry, and ecosystems. With the growing scarcity of surface water sources, developing nations are currently relying more on groundwater resources, particularly in rural parts of Africa. This paradigm shift results from groundwater becoming a strategic resource for economic growth, food security, poverty reduction, and groundwater suitability to adapt to climate change impacts on urban and rural livelihoods [1,2,3]. However, properly using and managing groundwater is difficult because of a lack of experience and knowledge [1,2]. Unlike surface water, groundwater systems are challenging to describe without a hydrogeologic framework model (HFM). HFM is essentially a conceptual model that helps develop a simplified path and framework for a numerical geologic model that embodies the concepts in a numerical representation of the geologic features and honors the geologic history of the region. The HFM helps groundwater resource management by describing the thickness and area extent of the aquifers, aquitards, and bedrock and embodies a 3D numerical regional estimation of the geologic framework. It also delineates faults, sedimentary facies, and related estimates of aquifer properties that result from different provenances.

The development of an HFM has evolved over time. Prior to the 1970s, the development of an HFM to represent the groundwater flow systems was accomplished using pencil and paper [4,5]. However, the approach did not fully represent complex terrains structurally or topographically. Detailed representations demanded maps that show details of lateral lithological contrasts for 3D modeling [6]. In order to close this gap, efforts were made to provide 3D information. For example, in the late 1970s, the Illinois State Geological Survey (ISGS) produced black-and-white maps and subsequently updated maps with color and patterns [7,8]. Yet even these maps did not sufficiently help users readily grasp the regional geologic features. As a result, efforts were again made to use the geographic information systems (GIS) approach to help speed up the production of colored maps to create cross-sections and line-and-dot patterns [9].

With the advance of computer technology, the development of HFM transitioned towards using advanced numerical-computer algorithms to create 3D geologic models that represent the subsurface structure and stratigraphy [4,10,11]. The layers that typically represent formations in the geologic model are delineated and estimated within the HFM (typically, formations above what is considered “bedrock”) that provide an estimate of water-bearing layers that can be developed for groundwater supply. The availability of computer resources and 3D geological software has led to developing an HFM as an important step in developing and managing groundwater systems [10,12].

The recent development of HFMs started in the petroleum industry and included “sequence stratigraphy” to replicate the geologic history of a region and related layering, faults, folds, and associated properties. Most of the development of petroleum framework models occurred in offshore regions. The first onshore HFMs were sequence stratigraphic HFMs of terrestrial regions developed by the U.S. Geological Survey (USGS) for the Santa Clara Valley along with the Hydrostatrigraphic framework model of the Central Valley [13], which is a hybrid approach to the sequence approach. There are three types of HFMs: Sequence Stratigraphic Models (Example, Earthvision, etc.), Hydrostratigraphic Models (Example, Rockware, Surfer, LeapFrog, etc.), and Hybrid Elevation Layer Models (Example, ISATIS, FREEWARE, TPROGS, Arc-geostatistical analyst, etc.).

Modern studies have developed HFMs with computer and software resources. For example, the Pajaro Valley groundwater basin hydrogeologic framework was developed by combining a driller’s log synthesis [14], where e-logs were used with resistance limits to delineate fractions of coarse and fine-grained material with hydrostratigraphic units [15]. A more recent use of driller’s logs was combined with an access database script to delineate coarse-grained and fine-grained thicknesses to be kriged within hydrostratigraphic units [13,16]. Then, many studies emerged employing ground-based and aerial Transient Electromagnetic (EM). However, they use proprietary software and are not linked to lithology databases or well-based transient EM logs (dual induction logs) [17].

Three-dimensional geological models can play a vital role in understanding the stratigraphic framework of aquifer systems and are used to create hydrogeologic representations of a groundwater flow system [18,19]. Three-dimensional geological models can also help estimate the geometry and related volumes of stratigraphic units or subregions delineated as aquifers or aquitards [20,21,22,23,24,25,26]. For example, digital 3D HFMs were constructed for many parts of California [13,24,25,27,28,29]. These 3D HFMs defined the aquifer system geometry and subsurface lithologic characteristics for the subsurface hydrogeologic conceptualization of the aquifer system to be used in regional hydrologic and hydrogeologic modeling. With the increased computing power and advanced 3D numerical models, datasets such as drill logs and geo-electrical resistivity data are used as input to construct detailed 3D geological models of groundwater flow systems [20,23]. For example, drill log and vertical electrical sounding (VES) data were combined to construct 3D geological models of a groundwater basin and it was reported that the combined method improved the representation of the groundwater flow system and reduced model uncertainties [21,30,31,32,33].

There are various options for using software for building a geologic model. Typically, the choice for software is either to use open-source GIS software or commercial geological modeling software. Open-source GIS options such as QGIS and gdal [34], are examples of open-source GIS software that provide access to process earth observation raster and vector data for visualization but are 2D maps, complicating efforts to visualize fault networks, complex structures, or stratigraphic sequences. Examples of commercial software that are closed-source include EarthVision [35], Leapfrog [36], GOCAD [37], Petra [38], Rockware [39], and Hydro GeoAnalyst [40], provide advanced visualization and algorithms for developing geological models [41] but do not provide access to the source code implementation. Because commercial software is closed source, the implementation of the methods is unknown, other than from the provided documentation, and cannot be modified. This limits the utility to only what is provided by the software and prevents connecting to open-source libraries for machine learning and computational inference. As a result, the integration into other computational frameworks is limited. Often, commercial software is used to enhance groundwater management. However, the cost of such software is usually prohibitive for most areas in the world. This has led to the development of open-source and free software that constitutes a cornerstone for enhancing groundwater development and management [42,43].

To fill this gap, efforts have been made to develop a fully open-source software application to create complete 3D HFMs that are comparable to implementations in commercial packages. In this regard, the popularity of using open-source programming languages such as Python is emerging and playing a crucial role in facilitating scientific programming and script-based science. GemPy, for example, is a fully open-source and easy-to-use software presented recently to create 3D visualization of geological models [4,44]. To construct 3D geological models, GemPy relies on Python and is based on an explicit and implicit potential-field interpolation approach. This interpolation approach consists of fast and automated surface formation supported with manual framing and interpolating of a scalar function based on the cokriging of point data and structural orientations. However, there are also other common approaches to developing HFMs, such as the use of spline interpolators to estimate the geologic top and bottom layer surfaces [45] and the use of 2D ordinary kriging or cokriging to spatially estimate the vertically averaged properties [46].

Groundwater is a primary source of water for irrigation and domestic water supply in rural arid and semi-arid regions in sub-Saharan Africa. However, this resource is little studied and poorly understood, partly because of the scarcity of existing hydrogeological information in many regions of sub-Saharan Africa. As a result, existing studies on assessments of groundwater resources relied on remotely sensed data combined with modeling and missed detailed information at the catchment scale [47,48]. In this study, we developed a 3D HFM considering Kobo Valley in Northern Ethiopia as a case study to understand the groundwater basin’s hydrogeology of the area, which is basic for sustainable groundwater resource evaluation and management based on field investigations and synthesizing Vertical Electrical Sounding (VES) measurements, driller’s logs, pumping tests, groundwater levels, and land surface elevation data. This study aims to fill the existing gap of implicit modeling in an open-source tool to develop 3D HFM in geosciences for a region that has a shortage of surface water and utilizes the groundwater for similar developments. The groundwater resources have not been previously studied in this way and are potentially vulnerable to considerable future overdevelopment. Hence, the methodology developed can be used for similar purposes, contributing to sustainable groundwater development and management.

2. Materials and Methods

2.1. Description of Study Area

Kobo Valley is located in the Afar Depression, the western margin of the Main Ethiopian Rift, East African Rift. Geographically, the study area is located at 11°54′ to 12°24′ N and 39°20′ to 39°48′ E, with a total area of approximately 1544.24 km². The elevation ranges between 3975 m above sea level in the mountains and 1018 m above sea level in the Kobo Valley plain area (Figure 1) with reference to UTM zone 37N (WGS 84). The valley is bounded in the west by mountains of the western plateau and in the east by a chain of mountain terrain, the rift escarpment. The western catchment is mountainous that is covered by and composed of volcanic rocks and has about 1045 km² (68%), while the valley plain is about 499 km² (32%).

The Kobo Valley is a semi-arid catchment and has limited surface water resources [49]. The use of groundwater is expanding to meet the growing demand for irrigation [50]. The Kobo Grana Valley Irrigation (KGVI) project serves an area of 3100 hectares and was developed to support the local community as they strive to maintain food security [51]. Over the past two decades, KGVI has developed more than 100 public boreholes to provide water for irrigation and domestic water supply. A series of groundwater simulation models have been developed to understand the Kobo Valley water resources and plan for future development [52,53]. However, the lack of a detailed geologic layer formation led the studies to the assumption of homogeneous aquifer properties.

2.2. Geology of the Kobo Valley

The study area is found in the Western Margin of the Afar Depression, the region in East Africa that represents a key location for studying continental breakup [54,55,56]. Zwaan and Corti Zwaan, Corti [54] revealed that the south–south striking Afar Depression is still actively deforming and is characterized by NNW–SSE normal faulting and a series of marginal grabens (Figure 2). The Afar Depression forms a triangular depression near and partially below sea level between the Ethiopian and Somalian plateaus to the west and south and the Danakil and Ashia Blocks to the northeast and east (Figure 2).

The lithology of the Western Margin of the Afar Depression consists of Miocene sedimentary and volcanic infill (the Pliocene–Quaternary Stratoid units) [58]. This forms the most recent unit and covers the Afar Depression floor but has also accumulated in depressions along the Western Afar Margin and on the Ethiopian Plateau. Faulting along the margin is generally considered dominantly antithetic (i.e., dipping away from the Afar rift basin, e.g., [59], Figure 2). Furthermore, a series of faulted basins referred to as “marginal grabens” [54,60] align along the Western Afar Margin.

The formation of the geological structure is controlled by tectonic events that led to the development of the Ethiopian Rift System on the western side of the Afar Margin. Its origin is considered to be local tectonic development forming an intermountain trough. These events are characterized by tensional movements, which gave rise to fissural volcanism followed by block-faulting and tilting to form the escarpment zone, including marginal grabens. These marginal grabens are narrow, elongated depressions bounded on both sides by normal faults facing each other. The main axis of the trough runs in a north–south direction. The trough is formed in the west by the rift escarpment and east by the horst of the mountain ridge (Figure 1c). The frame is mainly composed of Tertiary volcanic rocks (Figure 3). The eastern and western ridges bounding the plain area are characterized by opposite dipping faults parallel to the plateau escarpments. The intermountain trough (Kobo Valley) is dominantly composed of poorly compacted sedimentary basin-fill deposits [57,61].

2.3. Data Collection and Processing

Hydraulic parameter data of aquifers in the Kobo Valley well field from 63 boreholes and driller’s logs from 45 wells were collected from the office of Ethiopian Construction Design and Supervision Works Corporation [63] as well as 37-point VES data from the Ministry of Water, Irrigation and Energy (MWIE) office. The VES field data were collected by Metaferia Consulting Engineers [64]. The positions of the data points are validated using QGIS mapping and field visits made from January–February 2021. The data are also carefully checked for consistency in geological interpretation for each driller’s log and the descriptions of the geology of the area from different published reports.

Geological data of the area were collected from the Ethiopian Geological Study [62] and the Ethiopian Construction Design & Supervision Works Corporation Office [63]. Earthquake data collected from the USGS Earthquake catalog (accessed on 9 August 2021; https://earthquake.usgs.gov/earthquakes/search/) also helped us visualize the ongoing tectonic activity in the area. Study area location preparation and mapping of the structures within the study area were carried out through analysis in QGIS 3.2 (Accessed on 27 January 2020; www.qgis.org), programming language Python, and IPI2win software also used for terrain and VES data analysis. Detailed (30 m resolution) Shuttle Radar Topography Mission (SRTM) and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) digital topography data from NASA and METI (Accessed on 10 July 2020; https://earthexplorer.usgs.gov) provided an excellent basis for our mapping and processing using QGIS.

The final aim was to develop a three-dimensional hydrogeologic framework model (HFM), a conceptual model that comprised the hydrostratigraphy, structural features, and hydraulic properties that helped to understand the groundwater basin’s hydrogeology for groundwater resource management in the area. We used the limited existing data in the study area to achieve these results and followed several steps. Figure 4 presents our conceptual workflow. The input data were VES measurements, driller’s lithologic logs, pump tests, and land-surface altitude data. First, the input data were processed to extract the surface contact points, orientation measurements, and defined topological relationships (stratigraphic sequences and fault networks) and synthesized to create input data for the GemPy model. Then, the following three steps summarize the whole simulation hierarchically:

• Create a digitized geospatial database from the input data that contains all the raw data, topological relationships, standardized projection, and spatial extent;
• Define the spatial distribution of geological structures and discretize the 3D space regular grid geometry based on a potential-field interpolation method to define the spatial distribution of geological structures, such as layers, interfaces, and faults (computations of lithologic stratigraphic unit (LSU));
• Discretize and visualize an interactive 3D geological model using Python fundamental plotting library;
• Then, pre-process and analyze the driller’s log data to check whether they are consistent with the defined geometry and to identify the information that the contacts bring about the possible positions of the surface deviations.

The GemPy-derived 3D HFM model performance was evaluated using the built-in functionality to compute forward gravity conserving the package’s automatic differentiation and the concept of topology, a useful tool to describe adjacency relations in geological models, such as stratigraphic contacts or across-fault connectivity. As a final check, the generated aquifer profiles were mapped and evaluated with the drillhole profile data. Goodness of fit was qualitatively assessed by visual inspection and quantitatively assessed by statistical analysis of residuals (differences between the measured and calculated elevations).

2.4. VES Data

The VES point data were profiled into six lines from west to east (Figure 5). The raw data were collected using Terameter SAS 4000, which is the product of the ABEM Company according to Schlumberger electrode configuration. For each VES measurement point, the spreading was performed in the north–south direction in a way that the potential electrodes remained fixed for the defined current electrode, separating while the current electrodes were spread apart for each measurement. This arrangement is known as a Schlumberger Array. Thus, the VES points in each profile line were distanced at 1 km intervals in the west–east direction (Figure 5). The resistivity data were converted to a one-dimensional (1D) resistivity earth structure using IPI2win 3.1 software [65]. Then, using the IPI2win software, the apparent resistivity for each VES was plotted against half-electrode separation on double logarithmic paper [66,67]. The geophysical data were analyzed and interpreted for each VES point measurement. The specific resistivity values and the corresponding layers noted in each VES point were correlated along the west–east profile lines. A set of geoelectric sections was produced using the interpreted results of each VES point and correlating the values along the profile line. The inverted resistivity values were interpreted using the adopted resistivity values for earth materials [68,69,70].

2.5. GemPy Modeling Approach

GemPy is an open-source Python library for generating full 3D structural geological models based on an implicit potential-field interpolation approach [11,71]. The interpolation algorithm constructs 3D geological models, including fault networks, fault–surface interactions, and unconformities. This algorithm is applied in Python’s programming language, using a Theano library for efficient code generation that directly executes on graphical processing units (GPUs).

The method was first introduced by [71] and is grounded on the mathematical principles of universal cokriging. Later, the method was updated by integrating it into Bayesian inference frameworks and advanced machine-learning [72] for stochastic geomodelling and Bayesian inversion, making efficient implementations of automatic differentiation in novel machine-learning frameworks [4]. To efficiently compute gradients and provide optimized compiled code, GemPy 1.0 was built on top of Theano libraries [73]. In addition, Pandas for data storage and manipulation [74], Visualization Toolkit (vtk) for interactive 3D visualization [75], Matplotlib [76], and NumPy for efficient numerical computations [77] were used.

In GemPy, the main method to generate the 3D geological models is the potential-field method developed by [71], which has been successfully deployed in the modeling software GeoModeller 3D [11]. The basic idea is to construct an interpolation function $Z\left(x_o\right)$ where $x$ is any point in the continuous 3D space $\left(x,y,z\right)\in {ꓣ}^3$ that describes the domain ꓓ as a scalar field. The gradient will follow the planar orientation of the stratigraphic structure throughout the volume. It means that every possible isosurface of the scalar field will represent every synchronal deposition of the layer (Figure 6). After creating the stratigraphic layers, the fault series is considered as layer formations with a “Fault_Series” representation as the key entry in the GemPy set_series dictionary.

In the final interpolation function $Z\left(x_o\right)$, $x_o$ refers to the estimated quantity for some integrable measure $p_o$. To characterize the scalar field interpolation, two types of parameters were used: (i) layer interface points $x_{\alpha }$ for the respective isosurfaces of interest; and (ii) the gradients of the scalar field, $x_{\beta }$, which are poles of the layer or normal vectors to the dip plane in geological terms. Accordingly, gradients are oriented perpendicular to the isosurfaces and located anywhere in space. The gradient of the scalar field is referred to as $\partial Z/\partial u$ with $u$ defined as any unit vector and its samples as $x_{\beta }$. A complete description of the core functionality of the GemPy model can be found in Varga and Schaaf Varga, Schaaf [4] and additional references (See [71]). Appendix A contains the Python code that sets up the GemPy model presented in this paper, which utilizes the Kobo Valley datasets to generate a 3D visualization of the unknown hydrogeologic framework.

2.6. Model Performance Evaluation

Obtaining field geological measurements is expensive and often cost-prohibitive for developing nations. Consequently, most field measurements are sparse in time and space. Geological models must address the uncertainty that results from sparse data to reach a reasonable level of confidence in the model [4,78,79]. The advantage of GemPy is that the software model is fully designed to be coupled with probabilistic frameworks. GemPy supports stochastic geological modeling for uncertainty analysis (e.g., Monte Carlo simulations, Bayesian inference), which helps consider uncertainties in the model input data and use additional secondary information in a Bayesian inference framework. For example, GemPy can be coupled with pymc3 [72] to build probabilistic graphical models [4,80]. GemPy uses the latest developments in uncertainty visualization for 3D structural geological modeling and geological inversion [81,82].

GemPy includes a built-in functionality to compute forward gravity conserving the package’s automatic differentiation (AD). Topology helps to describe adjacency relations in geomodels, such as stratigraphic contacts or across-fault connectivity. GemPy can analyze the adjacency topology of its generated models using the topology compute method (See Appendix A (Figure A2)).

The generation of the stratigraphic geological layers using GemPy is formed by the number of points distributed in layers and orientations. The stratigraphic profiles were evaluated for goodness of fit using two approaches: visual inspection and statistical analysis of residuals. Visual inspection is a method of visually examining the map to see if the generated grid points with the model are a good representation of the original data. Residuals (difference between the measured value and an interpolated value) play an important role in interpolation characteristics; analysis is complete with a thorough examination of residuals [83]. Quantitative measures that can be used as goodness-of-fit statistics are performed, and a residual map is prepared to illustrate where the generated surface points are nearer or further away from the actual data of stratigraphic boundary elevations from drill logs. The mean absolute deviation and the standard deviation of the cross-validation residuals and the rank correlation between the measurements and the estimates are calculated for analysis. The coefficient of determination, R², is calculated by

$${\mathrm{R}}^2=1-({\mathrm{S}\mathrm{S}}_{\mathrm{r}\mathrm{e}\mathrm{s}}/{\mathrm{S}\mathrm{S}}_{\mathrm{t}\mathrm{o}\mathrm{t}})$$

(1)

where ${\mathrm{S}\mathrm{S}}_{\mathrm{r}\mathrm{e}\mathrm{s}}=\mathrm{sum}\mathrm{of}\mathrm{the}\mathrm{squares}\mathrm{of}\mathrm{the}\mathrm{residuals}$, ${\mathrm{SS}}_{\mathrm{tot}}=\mathrm{sum}\mathrm{of}\mathrm{squares}\mathrm{of}\mathrm{the}\mathrm{differences}\mathrm{from}\mathrm{the}\mathrm{mean},\mathrm{Sum}{({\mathrm{E}\mathrm{l}}_{\mathrm{i}}-{\mathrm{E}\mathrm{l}}_{\mathrm{mean}})}^2$, and El = the original elevation point data.

In addition to the built-in functionalities to the GemPy model, the model is checked by overlapping the generated 3D layer profiles with the well logger’s data to check the discrepancy of the layer marked in the well logs with the model.

3. Results and Discussion

3.1. Driller’s Log Lithology and VES Analysis

The geophysical data were interpreted for each VES measurement point. The specific resistivity values and the corresponding layers noted in each point were correlated with the west–east profile lines. The resistivity-survey data analysis showed that high resistivity values indicated greater sand and gravel content, whereas low resistivity values represented greater clay and silt content in the valley’s deposits, which were crosschecked with the drillers’ logs to confirm that low resistivity can also be caused by saline waters or poor water quality for coarse-grained layers. Six geoelectrical sections were produced, and four main layers were recognized along with the selected profiles (Figure 7) characterized as topsoil, clay and silty sand, sand and/or gravel, and weathered volcanic rock layers.

Profile lines with low resistivity segments due to clay deposition or salty water interpretation are compared to drillers’ log data to dismiss any effects from saline waters that would give a false signal that could be identified as clay deposition towards the east. Generally, the geoelectric correlation sections showed vertical and lateral variations in the profile layers and thickness because of Kobo Valley’s geomorphology.

The produced longitudinal sections (Figure 7) were constructed using the data collected at 1 km intervals. The created profile layers are not smooth because the layers were drawn by only connecting the points with a straight line. In transferring the raw data points to the GemPy model, this rough sketch was adjusted and smoothed with each layer orientation ($x_{\beta j}$) because an isosurface of the points connected the interface points and the scalar field perpendicular to the gradient. As explained in Section 2.2, the formation of the geological structure in the study area was controlled by tectonic events that affected the geophysical features of the valley plain and led to the development of uplands/remanent hills in the middle of Kobo Valley, which influenced the identified aquifer layer thickness to vary, as can be observed in Figure 7. For example, profile 3 of the generated geoelectric section showed that the layers are different in profiles and thickness to the west and east sides of the center, where the geomorphology shows upland areas/hills. This upland area is also coincident with faults, micro-seismic events, and a possible transform fault that separates the basin into two sub-basins.

To the east side of this section, either a layer of gravel is absent (profiles 2, 3, and 4 in Figure 7) or a thin layer of gravel is present (profiles 1, 5, and 6 in Figure 7). In the model, this thin layer is characterized as mainly clay. As a result, this area has a relatively thin layer of aquifer thickness. Generally, the produced longitudinal geoelectric sections show a clay layer beneath the thin topsoil layer. This clay layer has varying thicknesses throughout the whole valley. Furthermore, the sand and/or gravel and weathered volcanic-rock units follow from top to bottom as a third and fourth layer with varying thickness. Finally, the lowest layer below the weathered zone is the bedrock. The interpretation of each profile is given in

Table A1. Calculated volume of each individual layer of the study area.

Profile	VES ID	Profile Layers and Lithology Identified
Profile1	VESW1 to VESW8	Top soil (1 to 33 m), clay layer 112 m at VESW1 to 210 m thick at VESW5, sandy/gravel layer (5 m at VESW5 to 59 m thick at VESW8), weathered volcanic 50 m at VESW1 and 48 m at VESW5, bed rock
Profile2	VESW9 to VESW11	Top soil (1 to 8 m), clayey layer (174 m at VESW9 and 140 at VESW11), thin sand/gravel layer, weathered volcanic 24 m at VESW10 and 68 m thick at VESW11, bed rock
Profile3	VESW12 to VESW16	Very thin top soil, clay layer ranges from 76 m at VESW16 to 176 m at VESW12, gravel layer of 13 m at VESW13 to 19 m at VESW14, weathered rock of 19 m at VESW15 and 56 m at VESW14, bed rock
Profile4	VESK1 to VESK7	Top soil (2 to 6 m), sand/gravel layer at the western half at VESK 2 and 3 thickness of 59 m and 20.9 and clay at the eastern half and weathered volcanic at the center. Clay layer on western half has thickness of 33.5 m at VESK2 and 104 m at VESK3. The eastern clay layer is 160 m at VESK6 and 183 m at VESK7, weathered zone 24 m at VESK2 and 88.8 m at VESk7, bed rock
Profile5	VESK8 to VESK12	The sandy/gravel layer thickness varies from 105 m at VESK8 to 45 m at VESK12, the clay layer filling the central and eastern part of the profile is 105 m to 149 m at VESK12, weathered zone has thickness of 20 to 30 m, bed rock
Profile6	VESHG1 to VESHG9	Thick clay layer of max 150 m at VESHG7 and 114 m at VESHG4, sand/gravel layer with maximum thickness of 195 m at VESHG2 and minimum thickness at VESHG4 (10 m), third layer above the fresh bed rock is the weathered zone of 30 to 40 m.

in Appendix B.

This study showed that the central part of the valley is thicker in sediment deposits and fine-grained material. At the same time, the sediment thickness around the eastern side from the center of the valley is mainly clay. The profile data are digitized to obtain the x, y, and z surface contact input data for the GemPy model and defined topological relationships (stratigraphic sequences and fault networks).

3.2. Three-Dimensional Hydrogeological Framework

The KV-HFM was developed by constructing a GemPy model using a GemPy python model object (object here is in reference to object-oriented programming). The model object uses a regular grid to interpolate the 3D geological model at any point in a 3D space. GemPy relies on Theano library for efficiency [73]. Theano_optimizer is used with the fast compile option before computing the model to generate Figure 8 (which uses PyVista, a 3D visualization library and mesh analysis in Python for 3D visualization).

The generated 3D GemPy model significantly enhances the visual interpretation and understanding of the valley’s complex subsurface formations and geological profiles. This type of representation also helps non-professionals understand the subsurface profiles and complexity. The model enables each layer to be visualized individually using its spatial properties. Using this feature, the volume of each aquifer layer is calculated to estimate the water availability in each layer. This can be used to help quantify the groundwater resources in the region. Previous studies in Kobo Valley have been conducted with many generalized assumptions considering the homogenous thickness and the same structural settings throughout the catchment [50,53]. However, a GemPy model fills this knowledge gap by providing (1) the layer volume for the given geological settings; (2) visualization of the formation thickness distribution of each unit; and (3) the relations among the units (Figure 8, Figure 9, and Figure 10, respectively).

The 3D model shows that the graben is bounded on both sides, from West and East normal faults facing each other. This fault bounding is also observable from topographical maps. A cross-sectional view of the model (see Figure 10a–c) for the profiles indicated in Figure 9 shows that the West and East frames are mountain ridges bounding the plain area (marginal grabens) that are mainly composed of volcanic rocks. The valley plain forms the main aquifer system.

The 3D HFM also provides a map of the distribution and thickness of the layers in the main aquifer system. The spatial distribution of the thickness of a layer determines its volume. Hence the model facilitates the estimation of the volume of a layer by considering the spatial distribution of its thickness (

Table 1. Calculated volume of each individual layer of the study area.

Hydrogeologic Stratigraphic Unit (HSU)	HSU Order	Volume of Unit in m3
Clay, silt, and silty sand	First (top) layer	38.21 × 109
Fine sand, coarse sand, and gravel	Second layer	26.79 × 109
Weathered volcanic rock (basalts)	Third layer 1	17.22 × 109

¹ The layer below the third is hard rock.

). This is the volume of the layers in the hydrogeologic stratigraphic unit (HSU), not the volume of water or the water that could be extracted from the HSU. The aquifer volume is estimated using the valley plain area of the main aquifer (alluvial part) of the valley using the difference between the isopach of the potentiometric surface and the top of the aquifer.

3.3. Uncertainty in the GemPy Model

The first evaluation made was a visual examination. The map was visually examined to see if the generated grid points with the GemPy model closely represented the original data. Next, the top surface of the 3D structural geological map generated using GemPy was compared with the digital elevation map (DEM) of the study area. Then, the minimum, maximum, mean, standard error, standard, and coefficient of determination (R²) were calculated for the residuals from the differences between the known values of VES and the driller’s logs with the estimates (

Table 2. Summary of residual statistics.

Minimum (m)	Maximum (m)	Mean (m)	SE (m)	SD (m)	R2	HSU	Description
−5.68	7.71	0.01	0.08	0.91	0.93	HSU_1	top of clay, silt, and silty sand layer
−7.52	7.84	-0.58	0.60	6.76	0.95	HSU_2	bottom of clay, silt, and silty sand layer, and top of fine sand, coarse sand, and gravel layer
−9.68	7.71	0.01	0.10	1.01	0.90	HSU_3	bottom of fine sand, coarse sand, and gravel layer and top of basalts layer
−6.14	8.32	0.04	0.10	1.12	0.91	Base	top of hard rock

HSU = Hydrogeologic Stratigraphic Unit, SD = Standard Deviation, SE = Standard Error.

).

The residual’s minimum and maximum values indicate the magnitudes of differences between the generated grid values from GemPy and the actual measured data points. For example, the absolute maximum value is 9.68 m for HSU_3, where the elevation of the HSU varies from 1094 m to 2621 m elevation above sea level with reference to UTM zone 37N (WGS 84). Hence, 9.68 m is less than a 1% variation. The mean differences between GemPy model-generated grid point elevation values and measured elevation values were less than 1% for all HSUs, demonstrating that the gridded values were reasonably close to the original data values for most of the grid. The standard deviation was 6.76 m. These residual values reflect the accuracy of the original data points compared to the gridded values. The large values of the coefficient of correlation (R² > 0.90) for all HSUs indicate that the model explains most of the variation observed in measured data points, but some patterns were detected in the residuals (Figure 11). Figure 11 demonstrates that the model closely matches the original data points except along the main axis of the trough frame, the rift escarpment, and the horst of the mountain ridges, where measured data were scarce. Overall, we concluded that the model performed well with the available data.

To visualize the subsurface geologic unit extents and fault locations, fractures, and aquifer formations in GemPy, stratigraphic and unconformity connections were computed for each point in the grid (Figure 12a) to interpret the sparse field measurements. In the validation of the interpolation, it is assumed that the intrinsic error is zero and the validation is only for the error of the interpolation estimator. Well log data were used to evaluate the overlap variation with the layers generated by the GemPy model to minimize uncertainty in the input parameters and therefore in the model outcomes (Figure 12b).

3.4. Kobo Valley Aquifer and Sediment Layer Visualization

An aquifer is a saturated permeable geologic unit that can contain and transmit considerable quantities of water to wells and springs [84,85]. The thickness of Kobo Valley aquifer was determined from VES and drilling data. For this GemPy model, the Kobo Valley aquifer is assumed to be composed of sand, gravel, pebbles, and fractured volcanic rocks to obtain the spatial distribution of the aquifer thickness map. In addition, the driller’s log and the geophysical survey data of the sub-surface material below the water table in the catchment were analyzed to assess the thickness of the aquifer. Generally, the aquifer thickness varied throughout Kobo Valley and increased from the center to the north and south (Figure 13a). The average aquifer thickness of the study area was 104 m. All layers were considered as penetrated by boreholes, including volcanic rocks, for the sediment thickness spatial distribution (Figure 13b). The sediment thickness to the eastern side of the center of Kobo Valley is mainly clay, and as a result, the aquifer is relatively thin in this area. The sediment thickness varied from a minimum of 120 m in the eastern central part to a maximum of 220 m in the western part, with an average sediment thickness of 150 m for the entire study area.

3.5. Hydraulic Properties of the Valley

The pertinent hydraulic properties of the aquifer were hydraulic conductivity, transmissivity, and specific yield. Hydraulic conductivity is the capacity of an aquifer to transmit water and is expressed as the volume of groundwater at the existing kinematic viscosity that will move in unit time under a unit hydraulic gradient through a unit area at a right angle to the direction of flow. Transmissivity is defined as the rate at which water of prevailing kinematic viscosity is transmitted through unit width of the aquifer under a unit hydraulic gradient and can be calculated by multiplying the hydraulic conductivity by the saturated thickness of the aquifer. Specific yield is defined as the ratio of the volume of water that drains because of gravity to the total volume of saturated aquifer [85,86,87].

Hydraulic conductivity, transmissivity, and specific yield are required inputs for most numerical groundwater simulation models and proper management of groundwater resources. A pumping test is the most common technique for estimating these parameter values. For this study, pumping test data from irrigation and water-supply wells were obtained from the MWIE office and were analyzed. The first hydraulic parameter analyzed was transmissivity as a product of hydraulic conductivity and saturated thickness. Irrigation and water-supply wells were used to determine the transmissivity distribution of the alluvial aquifer. Transmissivity data were determined from performed constant pumping tests, which were carried out uninterrupted for 3 days. The spatial distribution of the transmissivity of the Kobo Valley was plotted using QGIS and varied from 7.9 m²/day to 2500 m²/day with an average value of 467 m²/day (Figure 14b). The hydraulic conductivity of the Kobo Valley was also mapped, and values ranged from 0.1 m/day to 35 m/day (Figure 14a). Generally, as noted in Figure 14a, the western side of the central part of Kobo Valley has a low hydraulic conductivity zone, which indicates a fine deposit, and the western part of the alluvial aquifer is a high hydraulic conductivity zone ranging from 7.1 m/day to 33 m/day. The third parameter is specific yield, a material boundary that represents the effective (drainable) porosity of the unconfined sediments and the specific retention of the volumetric fraction of water that remains during a unit decline in the water table [88]. Specific yields obtained from the 63 wells were mapped, and observed values ranged from 0.06 to 0.3 with an average value of 0.22 (Figure 14c). As can be seen from the figure, the specific yield decreased from west to eastward. The minimum specific yield was observed on the eastern side of the central and northern part of Kobo Valley and was likely due to thick silt and clay aquifer materials. In contrast, the high specific yields observed at the western margin of the study area were likely due to the coarse aquifer material, as observed in Figure 14c.

3.6. Groundwater Flow System

The groundwater-level contour map of Kobo Valley was generated from the GemPy model as shown in Figure 15a. Overlaying the groundwater-level contour map with the geomorphology of Kobo Valley (Figure 15b) and the hydrogeological longitudinal cross-section (Figure 10c), the data indicate a structural-based groundwater divide that divides the catchment into northern and southern groundwater systems that are prominently shown in the 3D HFM. The northern groundwater (north of the groundwater divide) flows in the northeast direction to the Selen Wuha outlet, whereas the southern groundwater (south of the groundwater divide) flows in the southeast direction to the Golina outlet (Figure 15b). Furthermore, the groundwater contour and the stream flow direction indicate that the groundwater systems have a hydraulic connection with the streams and rivers in Kobo Valley, and the flow direction is also influenced by the geomorphology and surface drainage of the area.

3.7. Groundwater Storage in the Valley

The developed 3D HFM and hydrogeological cross-section using GemPy (Figure 9 and Figure 10) were analyzed and Kobo Valley’s alluvial aquifer was identified as an unconfined aquifer with an impermeable hard rock as a base layer. Groundwater storage (GWS) is the groundwater in the pores of the alluvial aquifer and can be computed from the saturated aquifer thickness and specific yield of the aquifer [88,89,90,91,92].

$$\mathrm{G}\mathrm{W}\mathrm{S}=\left(\mathrm{A}\mathrm{\ast }\mathrm{H}\right){\mathrm{S}}_{\mathrm{y}}$$

(2)

where A is the total area of the aquifer in meter square (m²), H is saturated aquifer thickness (m), and S_y is a specific yield of the aquifer. The aquifer thickness varies throughout Kobo Valley (detail shown in Section 3.2). Using the saturated aquifer thicknesses and the average specific yield of 0.22, groundwater storage in Kobo Valley was estimated to be 4132 MCM. This estimate is higher than the groundwater storage estimates of 3081 MCM reported by [93] and 2548 MCM reported by [94]. We believe that our estimation is more accurate because unlike the previous studies, we developed a three-dimensional hydrogeologic framework model (HFM) that accurately defined the hydrostratigraphy, structural features, and hydraulic properties of the Kobo Valley groundwater basin hydrogeology.

4. Conclusions

In this study, the GemPy model, a Python-based open-source tool, was used to develop a 3D HFM for Kobo Valley, part of the East African Rift System in northern Ethiopia. This work developed a simplified representation of a groundwater system to better understand the basin’s groundwater hydrogeology and provide a tool that can be used for sustainable management of the groundwater system in the catchment. The developed model could help optimize the management of stratigraphic information; 3D visualization speeds up the process of stratigraphic setting evaluation, allowing for the identification of the existing geological aquifer layers in vertical and horizontal sections. The developed 3D HFM can help detect geological contacts, assess volumes and thicknesses, and evaluate geometric relationships between hydrostratigraphic units and their effects on the groundwater flow system of the basin. The 3D HFM, together with the geophysical and driller’s lithologic logs and literature review, facilitated the development of a conceptual model describing the dynamics of the groundwater flow system in Kobo Valley.

For Kobo Valley, GemPy provided a new perspective for understanding the groundwater basin’s hydrogeology and the development of a 3D HFM that could benefit sustainable groundwater management. Kobo Valley is an important source of fresh groundwater and has complex geological and structural settings. Estimation of the groundwater resource and modeling based on the generalized aquifer geometry may mislead the sustainable development of this aquifer. However, GemPy provides an easy, flexible, and interactive platform for incorporating the natural settings and complexity of the aquifer for the development and visualization of the 3D HFM of Kobo Valley. The groundwater storage in Kobo Valley was estimated as 4132 MCM; more accurate estimates can be used to improve groundwater resource management in Kobo Valley. The developed 3D HFM of Kobo Valley provides information to complement the development and management policy for sustainable groundwater extraction from Kobo Valley’s alluvial aquifer.

The GemPy model also demonstrated the existence of remnant volcanic hills in the middle of Kobo Valley that act as a structural-based groundwater divide that splits the catchment into northern and southern groundwater systems. The acquired visualization and understanding of the subsurface structure were essential for the quantitative modeling of groundwater flow, and we recommend modeling and evaluating the two groundwater systems as two separate sub-basins.

We believe that GemPy fills the existing gap of implicit modeling in the open-source ecosystem in geosciences and offers a reliable and easy-to-use technology to generate complex models with only a few lines of code. The advancements described in this study maximize the computational capacity of present computing systems and have the potential to help improve groundwater management in Kobo Valley.