Water Resource Assessment of a Complex Volcanic System Under Semi-Arid Climate Using Numerical Modeling: The Borena Basin in Southern Ethiopia

Abstract

The Borena basin is located in southern Ethiopia, in a semi-arid climate, on the eastern shoulder of the south Main Ethiopian Rift (MER). The study area covers 18,000 km² and is characterized by a lack of perennial surface waters that can be used for domestic and agricultural purpose. As a result, groundwater, which occurs in complex volcanic settings, is the only source for water supply in the study area. This work is focused on the basaltic aquifers, which are intensely fractured, resulting in strong connectivity within the system. All available data (geology, hydraulic head, hydraulic parameters, well inventory and discharge, etc.) were compiled in a GIS database. The overall objective of this work is the assessment of groundwater potential, its spatial distribution and factors controlling its movement using numerical groundwater modeling to enhance groundwater management and use in the Borena basin. The modeling task was conducted at two scales: (i) regional scale; (ii) wellfields scale. The regional steady state model was calibrated using the Pilot points approach, highlighting a strongly heterogeneous system. A significant result of the regional model consisted of estimating the water balance of the whole system. The total inflow to the basin amounts to 542 × 10⁶ m³/year, of which 367 × 10⁶ m³/year are provided by superficial recharge. Groundwater resources are exploited with 7 wellfields. Exploitation of the wellfields was optimized based on the Sustainable Yield concept, which reserves a fraction of natural recharge for the benefit of the environment (surface waters, ecosystems). Each wellfield was extracted from the regional model, refined and used to simulate and optimize pumping scenarios, with the objective of maximizing discharge rates and avoiding over-exploitation of the groundwater. The optimized abstraction at all wellfields amounts to 121 × 10⁶ m³/year, which represents 33% of the natural recharge and fully agrees with the sustainable yield concept.

1. Introduction

Ethiopia is an important country in East Africa, with an area of approximately 1.104 million km², and currently (2019) has an estimated population of 112 million. The population growth rate is 3.0%. The population should reach 175 million inhabitants in two decades (in 2040) [1]. Ethiopia’s economic development has been spectacular, with a double-digit annual growth rate over the past two decades [2]. Such an increase in economic activity in various sectors (agriculture, industry, etc.) has required a corollary demand for more and more water, including in the Borena area.

Ethiopia is characterized by a complex topography and geology. The Main Ethiopian Rift (MER), stretching from north-east to south-west, divides the country into two plateaus. It creates three major topographic regions: the Western Highlands, the Eastern Highlands, and the low-lying Rift Valley. The elevation ranges from 125 m below mean sea level at Danakil Depression to 4620 m above mean sea level at Ras Dejen peak. The climate is quite diversified due to its equatorial positioning and varied topography, with high inter and intra-annual rainfall variability [3,4,5]. Rainfall ranges from less than 200 mm/year in the eastern, south-eastern and north-eastern arid areas of the country, up to 2500 mm/year in the western highlands [6,7].

Rock exposures encountered throughout the country include the Precambrian basement, the Paleozoic and Mesozoic sedimentary rocks, the Tertiary sedimentary and volcanic rocks and the Quaternary sediments and volcanic rocks [8,9]. There is still retardation in the geological mapping of the country. The Ethiopia Geological Survey has so far covered 62% of Ethiopia’s surface at 1:250,000 scale, and the whole country has been mapped at 1:2,000,000 scale [10]. Hydrogeological investigations in the country started rather recently, in the 1970s, to supply water to the cities and the dispersed rural communities. The Ministry of Water and Energy (MoWE), the University of Addis Ababa, and the Ethiopia Geological Survey jointly launched the Ethiopian National Groundwater Database (ENGDA) in 2003. This groundwater database is intended to play a significant role in promoting and supporting hydrogeological studies.

The groundwater resource potential at the country scale is not clearly known. Some estimates [6,7,9,11] set the order of magnitude of the average annual recharge at around 40 × 10⁹ m³. This figure is an indicator of the importance of this resource. At present, groundwater resource is widely used in the country to supply drinking water to municipalities and in various economic sectors, in particular for irrigation. It provides, for example, roughly 40% of drinking water [7] for larges cities, including the capital Addis Ababa, with a population close to 4 million [1]. Groundwater resources are being exploited at an accelerating pace in Ethiopia [10], as elsewhere in the world, due to lower vulnerability to pollution and aquifer properties that are able to buffer the impact of external global changes (increase or decrease in rainfall, drought, flooding, etc.). Groundwater is thus often considered to be a reliable water resource compared to surface waters and is often over-exploited, leading to undesirable and harmful consequences (deterioration of groundwater quality, groundwater level lowering, environmental harm, etc.). However, groundwater is also a key component of the hydrosystem and plays a major role in preserving ecosystems and supplying surface runoffs in low water periods. Proper management of this valuable resource is necessary in order to preserve its sustainability.

In Ethiopia, the public authorities have therefore launched water resource assessment and development studies at basin scales, aiming at the sustainable exploitation of this vital resource. The distribution of aquifers in Ethiopia, at large scales, is at present fairly well known. However, quantitative information at basin scales is generally incomplete or even lacking [12]. This article focuses on water resource availability in the Borena basin. This basin, with an area of 18,000 km², is located in southern Ethiopia, within the geodynamic context of the Main Ethiopian Rift (MER), in a border zone with Kenya. The climate context is semi-arid. Rivers are not perennial, and available water resources are restricted to groundwater only. The growth of economic activities in the near future, especially in the agricultural sector, requires an assessment of available groundwater resources and an optimization of their exploitation to ensure their sustainability.

Distributed numerical modeling was used to assess the groundwater resources of the basaltic system in this basin. The modeling task was conducted at two scales: (i) regional scale; and (ii) wellfield scale. The wellfield optimization was undertaken with the objective of avoiding future over-exploitation of groundwater and excessive piezometric lowering. The overall objective of the modeling work was to propose to the end-users a sustainable framework to exploit the Borena groundwater and preserve this vital resource for the development of the region. The use of groundwater numerical modeling for this purpose is well suited, as it represents a widely applied methodology in contemporary hydrogeology for groundwater resource assessment and management. Aquifer numerical models are simplified mathematical representations of groundwater systems. They solve groundwater flow and solute transport equations using numerical methods. They are efficient tools that can simulate the response of hydrogeological systems to various stresses and exploitation alternatives [13,14,15].

In this work, the well-known software Modflow [16] is used to model the Borena basin. Modflow is a finite-difference groundwater flow modeling program, which solves the three-dimensional saturated groundwater flow equation. Modflow uses a block-centered formulation of the finite-difference equations and can simulate all types of aquifers. It has many optional packages which make it possible to simulate all the components of the hydrologic cycle. Since its first release, Modflow has been widely used in Hydrogeology to solve numerous issues related to groundwater exploitation, and is considered an international standard in groundwater modeling.

Groundwater resource assessment, management and over-exploitation issues have been extensively and successfully addressed using numerical modeling, and specifically using Modflow. Case studies can be found worldwide, including in countries like India [17,18,19], Pakistan [20], Iran [21], Greece [22], China [23], Taiwan [24], Ethiopia [25], and the USA [26]. In Ethiopia, for example, Shishaye et al. [25] developed groundwater resources sustainability analysis. They used Modflow to analyze the responses of an alluvial aquifer to pumping and the sustainability of the groundwater over 50 years. An important issue in groundwater analysis is related to stream–aquifer interaction. Modflow was used successfully to tackle this issue [27,28,29,30]. Declining groundwater levels due to over-exploitation are a major current issue, as groundwater abstraction is constantly increasing, sometimes dramatically. Groundwater level under various constraints have been predicted using Modflow [31,32]. Protection zones around wellfields have been delineated using groundwater modeling [33]. In mining hydrogeology, Modflow has been used to analyze and quantify flows components in mine environments [34]. Groundwater resources have been significantly affected by climate change for a few decades. Modflow has also been extensively used to analyze the aquifers’ responses to these changes. Steiakakis [35] simulated the drought-induced impacts and studied the combined effects of groundwater exploitation and climate variability on a karstic aquifer. Malekinezhad et al. [36] analyzed impacts of climate change and human activities on groundwater resources using Modflow. Hunter et al. [37] analyzed groundwater sustainability under climate change in an arid environment.

Over the years, groundwater models have become essential backbones of groundwater resources planning and management. Proper knowledge and understanding of the groundwater system functioning is essential prior to undertaking modeling [38,39,40]. Groundwater models require, however, varieties of data, including geology, hydrogeology, climatology, hydrology, soil, topography, etc., to be able to efficiently and accurately simulate groundwater systems. This data issue constitutes the main limitation to groundwater models applications, specifically in developing countries, where data are often scarce and sometimes affected by high uncertainty. Some attempts have been made to overcome this issue in order to use numerical modeling in such cases [41,42,43,44].

2. Materials and Methods

2.1. Study Site, the Borena Basin

2.1.1. Location of the Borena Basin

The study area is located in southern Ethiopia, western Borena Zone of Oromia Regional State (Figure 1). The UTM Zone 37 N coordinates of the study area are: X1 = 283,040 mE, X2 = 443,451 mE and Y1 = 400,000 mN, Y2 = 588,780 mN. The southern limits overlaps with the Ethio-Kenyan frontier. The study area covers a total surface of about 18,000 km². The climate context is arid to semi-arid type. The study area is part of South and Southeastern of Ethiopian moisture region that have two distinct dry periods (December to February and July to August) and two rainy seasons (March to May and September to November). The mean annual rainfall of the study area is about 593 mm. The mean annual maximum temperature of the study area varies from 23 °C to 28 °C, while the mean minimum temperature varies from 13 °C to 16 °C. The mean annual temperature of the area varies from 18 °C to 21 °C.

2.1.2. Geology and Hydrogeological Setting

The main lithostratigraphic units in the study area consist of Precambrian metamorphic and intrusive rocks, Tertiary and Quaternary volcanic rocks, and Quaternary superficial deposits (Figure 2A). These volcanic terrains include:

-

Tertiary volcanics: Pre-rift/Syn-rift volcanics including the Teltele basaltic sequence. The prominent mountainous ridges in the west and northwestern parts of the area, are made up of Teltele flood basalts. These volcanic successions have significant thickness exceeding three hundreds of meters.

-

Quaternary volcanics: Recent/Post Rift volcanics including the Bulal basalts. The Bulal basalts are the dominant horizontally layered sheets of basalt flows mapped in the area. This formation overlies the basement and its average thickness is about 200 m. The areas of Bulal basalts occurrence include vast terrains stretching from the southeastern periphery, to areas in the western-northwestern, north-eastern and central parts of the study area.

The volcanic rocks occupy the western plateau and the central plain land stretching from the Segen River in the north to Megado in the south. They occur in the form of lava flows, pyroclastic deposits, spatter cones, and scoria cones. They are predominantly basaltic in composition, associated with minor felsic pyroclastic deposits and phonolites, trachytes, and rhyolites.

Given the geodynamic context of the area, the rocks are highly fractured. Structural features in the area (Figure 2B) relate to the Precambrian orogenic f and Phanerozoic extensional tectonics. Most of the structures are in the NW, N-S and E-W directions. Two different rift systems occur in the study area: one is N-S trending and known as “Ririba Rift”, characterized by low scarp normal faults with less than 60 m vertical throw; the other is NW-SE trending and known as “Mega Rift”, characterized by high scarp (up to 1000 m) normal faults.

Geological structures (fractures, faults, joints, etc.) play an important role in the hydrogeological functioning of fractured systems. Hydraulic properties of the system (permeability, storage) are mainly related to the fractures/faults network density and connectivity. The fractures/faults survey shown in Figure 2B indicates that the network is quite dense (at the basin scale). Pumping tests achieved on the wells indicate that the fractures network connectivity is high and that the system behaves as a continuous porous media. Pumping tests data were analyzed using classical methods, as Theis or Jacob methods [46]. The available hydraulic conductivity (K, m/d) data plotted on the fractures network (Figure 2B) highlights the relation between permeability and structures. The highest K values are located along the Ririba Rift valley.

In relation to the regional lithostratigraphy, there are 3 types of aquifer/aquiclude systems in the study area: isolated alluvial deposits; extensive volcanic aquifers and the basement which forms the regional aquiclude with localized aquifers fracture and/or weathered profiles. Alluvial deposits consist of light brown to grayish colored sand and gravel, with minor silt and clay. They are found along the river channels (e.g., Ririba River). These formations do not form a continuous layer over the basalts. They occur as isolated deposits generally unsaturated. In general, in most of the areas their thickness ranges from 2 m to 40 m. These quaternary formations play a role of transmission of inflows (recharge) to the underlying volcanic aquifer systems.

Volcanic rocks in the area are of different age and may bear different hydrogeological characteristics. Among the different types of basalt units found in the study area, the main hydrostratigraphic unit, which constitutes the major aquifer, is the Bulal basalt formation, due to its geomorphological position and hydro-lithological properties.

Basement rocks or commonly known as Pre-Cambrian rocks are the oldest rocks in the region underlying all formations. The Pre-Cambrian crystalline basement consists of older metamorphic complex with various gneisses, schists and granulites, intruded locally by younger plutonic rocks. The crystalline rocks are generally considered in the area as regional aquicludes and groundwater basin boundaries. Weathered and fractured top of basement rocks can constitute under certain conditions a permeable medium which can contribute to groundwater storage and flow [47]. In the study area however, there is no data about the top part of the basement.

In arid/semi-arid context, recharge modalities are rather complex. Groundwater recharge rate during rainy periods is more important along the stream beds than diffuse recharge occurring over the whole basin area [48]. The recharge through the rivers beds is often termed as preferential or indirect recharge [49]. Not any differential measurement of rivers discharge rates was made during flood periods in the study area. In previous studies [45], an attempt was made to evaluate the diffuse recharge using the water balance method (WBM). This approach [50] provides a first estimate of the recharge, which can then be improved during the model calibration process. The WBM was applied to estimate the recharge in each drainage basin in the study area (

Table 1. Diffuse recharge estimations using the water balance method [45].

Sub-Catchment	Recharge (mm/yr)
Ririba basin	73
Teltele	67
Bulal basin	31
Megado	43

). These diffuse recharge estimations are used as initial inputs in forthcoming modeling.

A thorough analysis of isotopic data of the Borena basin (rainfall, groundwater) was previously conducted [51]. This work made an important contribution to the understanding of the functioning of Borena aquifer system. The main findings, with respect to the modeling of the system, are summarized below:

-

In all localities for which environmental isotope data were available, the dating results are returned as modern age, i.e., post 1950s. Some waters with low tritium content are on the order of a few decades. Therefore, the aquifer system of the region is associated with modern day recharge.

-

The environmental isotope analysis, associated with hydrogeochemical analyses, clearly indicates regional groundwater flow from the eastern crystalline highlands to the central part (middle Gelchet plains) and to southeast following the lowlands.

-

A flow from Mega and Megado escarpments to Megado/Biloko lowlands is also identified.

A noteworthy finding is that recharge is sustainable, as the groundwater is primarily modern. Given the above considerations, the inflows to the volcanic system occurs as follows: (i) Input from the mountainous borders in the East; (ii) Diffuse recharge from the basalts outcrops and from the alluvial cover which acts as flow transmitter to the groundwater; (iii) Indirect preferential recharge from the main course of rivers.

Groundwater flow pattern in the basaltic system was analyzed using hydraulic head data collected in available wells in the study area (Figure 3). Head data are distributed over the entire volcanic system. The piezometric map of the volcanic system drawn from these head data is shown in Figure 3. The piezometric head remains below the top of the basalts layer, indicating that the basalts layer is not fully saturated, including 2 zones, a saturated one and an unsaturated one, separated by the piezometric surface. This map highlights important features regarding the functioning of the system. Comments are given below:

-

The inflow from the eastern boundary and the outflow at the Ethio-Kenyan boundary are clearly shown. Please note that this inflow to the aquifer from the eastern highlands is also supported by the isotopes and hydrochemical analyses.

-

A groundwater divide is illustrated in the northern part, and coincides with the surface water divide between Mansagerado and Ririba sub-basins.

-

The role of the Ririba rift is well highlighted. Groundwater flow is converging from the eastern, western and northern highlands to the Rirriba rift zone.

-

This map also clarifies the relations between the northern Teltele area and the main flow area southward. Groundwater is located at much higher altitudes (compared to Bulal/Ririba basins), and flows towards East to the Mansagerado sub-basin, South to the Ririba basin and to the North and West. High gradients between the Teltele basins and the adjacents basins (Ririba and Mansagerado) are depicted.

-

Consequently, the piezometric data lead to the conclusion that groundwater is unconfined in the study area at the basin’s scale, flows through the Ririba rift valley and leaves out the study area mainly towards Kenya.

A transverse W-E cross-section (Figure 4) highlights the relations between the aquifer and the western and eastern boundaries.

To date, a total of 19 pumping tests have been completed and interpreted, providing a data base including values of transmissivity (T, m²/d) and hydraulic conductivity (K (m/d) [45]. Summary statistics of the hydraulic parameters deduced from the tests are given in

Table 2. Summary statistics of hydraulic parameters—Bulal basalts.

Parameter	Minimum	Maximum	Average	Standard Deviation	Coefficient of Variation (%)
Transmissivity (T, m2/d)	2.4	1380.0	152.0	330.2	220
Hydraulic Conductivity (K, m/d)	0.01	17.1	2.0	4.3	220

.

Transmissivity ranges from 2.4 m²/day to 1380 m²/day and hydraulic conductivity from 0.01 m/day to 17 m/day. Both parameters span almost three orders of magnitude. Their standard deviation and coefficient of variation are very high. These parameters are very variable in space highlighting a strong heterogeneity of the volcanic rocks. Hydraulic conductivity (K) is further upscaled in the modeling phase using geostatistical non-linear methods.

2.1.3. Drainage Networks

There are three major drainage basins, identified as Bulal, Ririba (including Gooso-Goray) and Arbele-Megado, which are dry except during rainy seasons (Figure 5). The former two flow from north to south. The last is composed of two sub-basins, Arbele flowing NE-SW and Megado flowing NW-SE. All drain into Kenya out of the study area. The largest drainage is the structurally controlled Ririba River which stretches along the N-S trending graben of the Ririba Rift System. There are also other smaller catchments draining the study area in different directions, like Mansagerado, Jeso and Teltele, located in the northern part, and flowing northward.

2.2. Conceptual Model

2.2.1. Aquifer Geometry and Boundary Conditions

Aquifer geometry refers to the specification of the top and bottom elevation of the aquifer system, as well as thickness and areal extent of the water bearing layer. For reasons indicated in the above sections, the volcanic terrains are considered as the main aquifer system to be modeled. Therefore, the model domain (Figure 6) is focused on the basaltic formations and matches: (i) on its eastern, south-eastern limits, the geological boundaries between the volcanic and the basement rocks; (ii) on its western limits, the volcanic rocks forming mountains; (iii) on its northern limits, the borders of Jiso, Mansagardo and Teltele sub-basins.

The south-western limit coincides with the Ethio-Kenyan border. The digital elevation model (DEM) of the model domain is specified as the top elevation of the aquifer in the model. The bottom of the aquifer is formed by the contact between the volcanic rocks and the basement.

Boundary conditions are constraints specified on the active model to characterize the interaction between the active simulation domain and the surrounding environment. In groundwater models, which are used for analyzing groundwater flow analysis, the specification of the boundary conditions usually defines the source of water to the system and its ultimate manner of discharge [5]. Thus, boundary conditions are one of the key aspects in the proper conceptualization of a groundwater system and representation of that system in a numerical computer model.

There are generally three types of boundary conditions (BC): (i) specified head (First Type), (ii) specified flow (Second Type), and (iii) head-dependent flow (Third/mixed Type). The boundary conditions on the lateral limits of the model are shown in Figure 6.

-

The northern boundary coincides with the northern limits of the Teltele sub-basins. A general head boundary (GHB) is assigned to this limit of the Borena basin.

-

The north-western boundary is a no-flow boundary and coincides with the Mermero mountain range.

-

The eastern limit, represented by the basalts and basement contact, is designated as an inflow BC.

-

The south-western limit of the flow domain, which represents the Ethio-Kenyan border, is designated as mixed BC.

2.2.2. Model Design

Geology of the model area, scale of interest and purpose of the model are important factors that should be considered when conceptual models are developed for fractured rock aquifers. These factors play an important role when selecting suitable approaches to numerically represent fractured aquifers. When fracture densities are very high, the system can be treated as one continuum where hydraulic parameters account both fractures and matrix blocks [52,53]. In addition to this, as the scale increases, the more appropriate is to employ equivalent porous media (EPM) modeling approaches [54]. In this case, standard finite difference and finite element codes used for porous media can be applied to the EPM that represents the fractured system. The volcanic rocks of the Borena basin fit those terms quite well (highly fractured medium, strong connectivity).

As a result, the EPM approach remains the most plausible one in the case of the Borena system. In this study, the numerical simulation is performed using the modular finite differences block-centered groundwater flow code, Modflow [16]. Given the available knowledge, the model is composed of a single layer representing all the volcanic formations. The analysis of groundwater levels after stabilization in the wells, vs. the ground elevations at these wells, displays a significant linear relationship between these two parameters (Figure 7). Such a relation between topography and groundwater depth is a characteristics of unconfined aquifers [55]. This result, together with the fact that the basalts are not fully saturated, implies that the volcanic aquifer system, at the scale of the Borena basin, behaves as an unconfined system. This behavior is associated with the dense network of fractures, which creates a significant hydraulic connectivity within the volcanic system.

Therefore, in this modeling work, the response of the Borena groundwater flow system is outlined with a two dimensional steady state condition considering a single layer unconfined aquifer system. The governing flow equation is written as follows [16]:

$$\frac{\partial }{\partial \mathrm{x}}\left[{\mathrm{K}}_{\mathrm{x}}\frac{\partial \mathrm{H}}{\partial \mathrm{x}}\right]+\frac{\partial }{\partial \mathrm{y}}\left[{\mathrm{K}}_{\mathrm{y}}\frac{\partial \mathrm{H}}{\partial \mathrm{y}}\right]+\mathrm{W}=0$$

(1)

where K_x, K_y are the hydraulic conductivity along the x and y coordinate axes (m·s⁻¹). As the media is assumed to be isotropic, K_x = K_y. H is the potentiometric head (m) and W is a volumetric flux per unit volume representing sources and/or sinks, with W < 0 for flow out of the groundwater system, and W > 0.0 for flow into the groundwater system (s⁻¹). This equation is derived by combining a water balance equation in a REV (Representative Elementary Volume) Equation (2) and the Darcy’s law Equation (3):

$$\frac{\partial {\mathrm{q}}_{\mathrm{x}}}{\partial \mathrm{x}}+\frac{\partial {\mathrm{q}}_{\mathrm{y}}}{\partial \mathrm{y}}=-\mathrm{W}$$

(2)

$${\mathrm{q}}_{\mathrm{x}}=-{\mathrm{K}}_{\mathrm{x}}\frac{\partial \mathrm{H}}{\partial \mathrm{x}}; {\mathrm{q}}_{\mathrm{x}}=-{\mathrm{K}}_{\mathrm{y}}\frac{\partial \mathrm{H}}{\partial \mathrm{y}}$$

(3)

where q_x and q_y are specific discharges along the x and y coordinate axes (m·s⁻¹).

A uniform spacing of 250 m × 250 m (both rows and columns) is used in discretizing the flow domain. Each cell is 0.0625 km². This resolution is acceptable for modeling the groundwater at the Borena basin’s scale.

2.3. Calibration Method

Groundwater model calibration is a process wherein some parameters of the model (recharge, hydraulic conductivity, boundary conditions, etc.) are systematically modified and the model is repeatedly run until the computed results reproduce observed data with an acceptable degree of accuracy, in association with the modeling scale. This study concerns a very heterogeneous volcanic aquifer. The conventional approach consists to divide the model domain in several zones of homogeneous properties and perform trial and error calibration [17]. This implies a great number of trial and error runs as aquifer properties are manually modified. In addition, more and more zones are added to account for the system heterogeneity. This approach often leads to unrealistic results for several reasons (arbitrary zonation, unmapped heterogeneity, nonuniqueness of the results, etc.). The non-linear parameter estimation software PEST [56] may be used in this process for assistance, but it does not actually help to overcome all these issues, e.g., the number of zones.

As an alternative, the issue of high heterogeneity in numerical models of groundwater can be addressed by the Pilot points methodology, a novel approach which is being used more and more by hydrogeologists. [57,58,59,60]. The basis of the Pilot points methodology as a method of spatial parameterization of groundwater models is that hydraulic properties (e.g., hydraulic conductivity) are assigned to a set of points distributed throughout the aquifer rather than to the cells of the numerical model. The property values are then interpolated into the model cells from the Pilot points. The property values are assigned to the Pilot points as in any normal calibration procedure, in order to minimize the discrepancies between observed head values and model calculated head values. The interpolation method from Pilot points to the model cells used in this study is the geostatistical kriging procedure which presents several advantages: smooth interpolator, respect of the known values [61].

For large aquifer models, getting the whole panel of data (hydrogeology, hydrology, geology, topography) is usually expensive, and data are generally deficient. To limit the model uncertainty, the model should well simulate the observed head values and also an overall adjustment is necessary. The model should simulate well the general groundwater flow pattern [14,62,63].

3. Results and Discussion

3.1. Calibration of the Model

Figure 8 shows the calibration diagram, which compares observed head values vs. simulated head values. The model simulates the observed values well. The correlation between observed and simulated values is significant. The determination coefficient is high: R² = 0.99. Other measures to assess model calibration are provided in

Table 3. Model calibration statistics.

Determination Coefficient	Range in Observed Heads	Mean Error	Normalized Mean Absolute Error	Normalized Root Mean Square Error
0.99	664 m	−3.3 m	0.021	0.025

.

The average error is low. The normalized mean absolute error and the normalized root mean square error are very low. Overall, the model can be favorably considered to be calibrated. Locations of Pilot points are shown in Figure 9A. The spatial distribution of hydraulic conductivity is shown in Figure 9B. The values of the hydraulic conductivity range between 8.5 10⁻⁴ m/d and 6.5 10⁺² m/d and span 5 orders of magnitude, revealing the strong heterogeneity of the system. The highest hydraulic conductivities are located along the Ririba rift valley. The lowest values of K characterize the Teltele sub-basin and explains the high hydraulic gradients between the Teltele and the surrounding basins.

3.2. Simulated Piezometry

The simulated piezometry, reported in Figure 9B, reproduces well the known characteristics of the underground flow according to the piezometric field data. The general pattern of the flow is simulated efficiently, enhancing the model calibration. The following points are noted:

-

An important groundwater flow zone appears along the Ririba rift valley towards the Kenyan border to the south; this flow is probably associated with existing N-S structures in this area;

-

The basaltic system is supplied on the eastern boundary, in contact with the basement rocks;

-

The surface water divide line between the Teltele sub-basins and the main basins of Bulal/Ririba also marks a divergence of the underground flow over the entire study area: south of this line, groundwater flows to the Kenyan border, north of this line, it flow towards the Segen river;

-

The high hydraulic head values in the Teltele basin are well simulated;

-

The very high flow gradients on the boundaries of the Teltele basin are simulated correctly;

-

The Kenyan border is the main outlet for the groundwater flow to the south.

3.3. Water Balance Simulated by the Model

The water balance of the whole model is given in

Table 4. Water balance of the Borena basin simulated by the model.

	Inflow (m3/year)	Outflow (m3/year)
Inflow from the Eastern boundary	175 × 106	---
Diffuse Recharge	328 × 106	---
Indirect recharge through rivers beds	39 × 106	---
Outflow to the North	---	25 × 106
Outflow to the Kenyan border	---	517 × 106
Total	542 × 106	542 × 106

. The simulated recharge rates are reported in

Table 5. Recharge rates simulated by the model.

	Diffuse Recharge Rate (m3/m2/day)	Indirect Recharge Rate through Rivers Beds (m3/m2/day)
Teltele basin	1.8 × 10−5	5.0 × 10−5
Ririba basin	1.3 × 10−4	3.9 × 10−4
Laga Bulal basin	5.9 × 10−5	6.1 × 10−5
Megado basin	2.9 × 10−5	5.0 × 10−5

.

Comparing the diffuse and indirect recharge rates shows that in all of the sub-basins, the recharge rate through the rivers bed is higher than the diffuse recharge rate over the basalt outcroppings. In the Ririba and Teltele basins, the indirect recharge rate is almost 3 times higher than the diffuse recharge. The recharge rates (diffuse or indirect) are the highest in the Ririba basin, probably in association with the structures. At the scale of the Borena basin, the volume of the diffuse recharge remains much higher (61% of total inflow) than that of indirect recharge (7% of total inflow) because of the areas involved in each recharge modality. The inflow from the Eastern boundary represents 32% of the total inflow in the basin. Outflow to the North towards Segen river represents only 5% of the total outflow, whereas, the Kenyan border constitutes the main outflow of the Borena system (95%).

3.4. Sensitivity Analysis

Sensitivity analysis consists of evaluating how much the model input parameters can affect model outputs, i.e., heads and/or flows [64,65]. Evaluating the relative effect of the parameters on the model results provides a better and more fundamental understanding of the functioning of the simulated groundwater system. In general, the model input parameters taken into account in a sensitivity analysis are hydraulic conductivity and recharge. Hydraulic conductivity of the Borena basin was regrouped into 54 zones. Recharge was divided according to the 4 main watersheds. In each both diffuse and preferential recharge are differentiated. Thus, 8 recharge zones are considered in the sensitivity analysis (Figure 10).

Results of the sensitivity analysis are reported in Figure 11 (sensitivity with respect to Hydraulic Conductivity) and Figure 12 (sensitivity with respect to Recharge). The multiplier coefficient varies between 0.5 and 1.0 on the left side and between 1.0 and 1.5 on the right side. The model sensitivity parameter is the absolute residual mean (ARM).

The Hydraulic Conductivity sensitivity results show that ARM variation is less than 5 m for the most sensitive zone. Considering Recharge, Figure 12 shows that zone 1 (Ririba sub-basin) stands out in relation to the others. This is related to the higher hydraulic conductivity in this area. The remaining 7 zones, including all the rivers courses and the Teltele, Bulal and Megado sub-basins, display a sensitivity comparable to Hydraulic Conductivity (same variation magnitude of ARM). Although the Recharge and Hydraulic Conductivity knowledge available throughout the basin is still insufficient, a special effort should be made to improve the Recharge knowledge in the Ririba sub-basin.

3.5. Simulation of Groundwater Resource Exploitation

3.5.1. Simulation Using the Global Model

The calibrated model is then used to simulate the exploitation of the basaltic aquifer system. The network of exploitation wells is distributed on 7 wellfields (Figure 13). The cumulative discharge of these wells is on average 44,000 m³/d (15.8 Mm³/year).

The model water balance, when all wells are operating with the current pumping rates, is shown in

Table 6. Water balance of the Borena basin with all the wells being exploited.

	Inflow (m3/year)	Outflow (m3/year)
Inflow from the Eastern boundary	175 × 106	---
Diffuse Recharge	328 × 106	---
Indirect recharge through rivers beds	39 × 106	---
Outflow to the North	---	23.9 × 106
Outflow to the Kenyan border	---	502.6 × 106
Wells discharges	---	15.8 × 106
Total	542 × 106	542 × 106

. The whole yearly pumped volume of water represents about 4.3% of the total recharge of the basaltic aquifer. However, caution should be exercised before considering any increase in well pumping rates. This issue, i.e., the increase in pumping rates, is discussed below using the refined model of each wellfield.

3.5.2. Wellfield Simulations and Optimization

The exploitation of the wellfields should be carried out taking into account the sustainability of groundwater resources at the basin scale. There are two concepts, widely discussed in the literature, for defining the safely exploitable groundwater resources: Safe Yield and Sustainable Yield [66,67,68,69]. Safe yield, defined at the beginning of the 20th century [70], is now an obsolete concept. Roughly, Safe Yield can be equated to the natural recharge of groundwater systems. Obviously, the whole natural recharge of any groundwater system cannot be fully exploited, as the role of groundwater in the hydrologic cycle is essential. Any uncontrolled abstraction may not only affect the aquifer (excessive groundwater depletion) but also the groundwater-fed surface water (springs and base flow) and groundwater-dependent ecosystems (wetlands and vegetation). The concept of Sustainable Yield, on the other hand, reserves a fraction of natural recharge for the benefit of the environment (surface waters, ecosystems) [71,72,73]. There is, however, a lack of consensus about the percentage of natural recharge that should constitute sustainable yield of a basin. Average percentages may be around 30% to 40% [74,75], while in a less conservative approach, this percentage may reach 70%.

A global constraint in the optimization procedure of wellfield exploitation was determined in agreement with the end-users (MoWE), consisting of a reasonable compromise so as not to exceed a ratio of sustainable yield to recharge of 30% to 40% at the basin scale.

The wellfields under consideration, shown in Figure 13, are the following ones: Gelchet-Wobock, Liso-Sadeka, Mermero-Gocha, Megado, Hobok, Sarite and Elkune (or Teltele). The global model of the Borena basin is used to enhance the analysis of the potential of each wellfield. The objective is to optimize the total abstraction of each wellfield, while preserving the environment. The following steps are involved in this analysis: (i) Extract the wellfield from the Borena general model; (ii) Refine the wellfield model; (iii) Simulate and optimize pumping scenarios.

The wellfield models are extracted from the global model and refined at cell sizes of 50 m × 50 m. A gradual refinement is also applied around each well in order to locate each pumping well in a cell with sizes of less than 1 m × 1 m, so that the drawdowns simulated in the wells are more realistic.

The constraint used to maximize abstraction at each wellfield is related to acceptable drawdowns. The saturated thickness of the basaltic aquifer (Bulal basalts) is variable, but does not exceed 100 m (see Figure 4). There is a general consensus in the literature that excessive groundwater depletion should be avoided [32], which is in line with the Sustainable Yield concept. Setting constraints on groundwater levels to maximize extraction is a common practice in optimization procedures [76,77]. In the Borena basin, the groundwater level constraints are set according to the hydrogeological environment at each wellfield.

The wellfield optimization approach is illustrated using the example of the Hobok wellfield. This wellfield includes 3 wells (BTW11, NBTW6 and NBTW7). They are plotted in Figure 14, which also shows the piezometric map of this area when no well is operating. This map indicates that there is more concentrated groundwater flow in the SE part of the area, in relation to the general flow from Ririba rift valley to the Ethio-Kenyan border. The water balance of the wellfield is given in

Table 7. Water balance of the Hobock wellfield when no well is operating.

	In (m3/year)	Out (m3/year)
Wells	0	0
Recharge	41.8 × 106	0
Head Dependant Boundaries	169 × 106	210 × 106
Total	210 × 106	210 × 106

. Current discharge rates, simulated drawdowns and theoretical Specific Capacity are reported in

Table 8. Current discharges rates and drawdowns simulated by the model at Hobock wellfield. DD: drawdown; SC: Specific Capacity.

Well	X	Y	Q (m3/d)	Simulated DD (m)	Theoretical SC (m2/d)
NBTW6	296003	487267	2400	5.1	474
NBTW7	310241	474412	2400	3.7	647
BTW11	322212	482717	518	0.7	711

. The water balance of the wellfield with the current pumpings is given in

Table 9. Water balance of Hobock wellfield with the current pumping rates.

	In (m3/year)	Out (m3/year)	Ratio Abstraction/Recharge
Wells	---	1.94 × 106	4.6%
Recharge	4.18 × 107	---
Head Dependant Boundaries	1.69 × 108	2.09 × 108
Total	2.11 × 108	2.11 × 108

.

The results show that the potential of the Hobock wellfield is almost equivalent for all three wells (474 m²/d < SC < 711 m²/d). With the current exploitation, the ratio abstraction/recharge is quite low (4.6%), and drawdowns (DD) are small (1 m < DD < 6 m). Abstraction is maximized under the constraints that the ratio abstraction/recharge is maintained at a moderate value and the piezometric depression is limited. In agreement with the end-users (MoWR), a reasonable compromise was reached, limiting the drawdowns to 20 m at this wellfield. This constraint is acceptable given the average saturated thickness of the aquifer. The results of this scenario are given in

Table 10. Optimized discharge rates at the Hobock wellfield.

Well	x	y	Current Q (m3/d)	Optimized Q (m3/d)
NBTW6	296,003	487,267	2400	9073
NBTW7	310,241	474,412	2400	12,232
BTW11	322,212	482,717	518	9756
Total discharge			5318	31,061

and

Table 11. Water balance of the Hobock wellfield with optimized pumping rates.

	In (m3/year)	Out (m3/year)	Ratio Abstraction/Recharge
Wells	---	11.4 × 106	27.3%
Recharge	41.8 × 106	---
Head Dependant Boundaries	173.6 × 106	204 × 106
Total	215.4 × 106	2.15 × 106

.

Under this scenario, total discharge rates can be augmented up to 11 Mm³/year, while current pumpings amount to 1.9 Mm³/year. As drawdown is limited to 20 m, this scenario preserves the environment. Note also that the ratio abstraction/recharge remains moderate (27%).

The 6 other wellfields were optimized using the same procedure. The optimization work is given in the Supplementary Materials. An optimized scenario is proposed for each wellfield under constraints respecting the environment and in line with the Sustainable Yield concept. The results show that wellfield potential is spatially variable in each wellfield, being very different from one wellfield to another. Optimum abstraction at each wellfield is reported in

Table 12. Current and optimized abstraction (m³/year) at each wellfield.

Wellfield	Current Abstraction	Optimized Abstraction
Mermero-Gocha	4.56 × 106	4.78 × 106
Gelchet-Wobock	5.46 × 106	18.7 × 106
Hobock	1.94 × 106	11.4 × 106
Megado	2.01 × 106	4.0 × 106
Liso-Sadeka	3.98 × 106	69.5 × 106
Sarite	1.55 × 106	11.6 × 106
Teltele	1.30 × 106	0.98 × 106
Total	20.8 × 106	120.9 × 106
Ratio Abstraction/Recharge	4.6%	33%

. Compared to the current abstraction rate, Table 12 shows that at Liso-Sadeka, the exploitation can be greatly increased (×17 times). At other wellfields (Gelchet-Wobock, Hobock, Megado, Sarite), the increase in exploitation can be significant. At Mermero-Gocha, the current abstraction almost represents the predicted optimized abstraction. At Teltele wellfield, however, the current abstraction is significantly larger than the optimized one. Contrary to other wellfields, exploitation of the Teltele wellfield should be decreased in order to not cause excessive piezometric depression. The optimized total abstraction for all wellfields represents 33% of the basin recharge, which is fully in keeping with the sustainable exploitation of basin water resources. This issue, regarding the relation between available groundwater for abstraction and aquifer recharge, is widely debated in the literature [78,79]. The old concept of ‘safe yield’ has nowadays been abandoned, and has been replaced by approaches aiming at sustainable groundwater management, which should account for groundwater uses (drinking, agriculture, industry) and ecological water demands as well (baseflow to streams, etc.).

The optimized exploitation of the Borena basin’s groundwater resources, using the existing wellfields, amounts to 121 Mm³/y. The comparison with the present exploitation (21 Mm³/y) shows that the potential of the basin is significantly higher. The optimized abstraction volume, simulated by the model, must of course be considered as a framework. The model is a perfectible tool, and needs to be validated in the transient regime. To this end, additional data are recommended, specifically continuous monitoring of the groundwater, at least, at each wellfield. Natural phenomena, like the impact of climate change, should also be considered in the future exploitation of basin resources.

4. Conclusions

The overall objective of this modeling work was to propose to end users a sustainable framework for exploiting the Borena volcanic aquifer system and preserving this vital resource for the development of the region. A numerical flow model was successfully developed, using the Pilot Points calibration approach. The model successfully accounted for the strong heterogeneity of the volcanic system. An important outcome of the modeling task consisted of the estimation of the water balance of the basin. The total inflow to the basin amounts to 542 × 10⁶ m³/year, of which 367 × 10⁶ m³/year are provided by superficial recharge. The abstraction rates of seven wellfields exploiting the groundwater resources were optimized using the model and within the framework of the Sustainable Yield concept. Total optimized rates amount to 121 × 10⁶ m³/year, which represents 33% of the natural recharge and is fully in keeping with the Sustainable Yield concept.

The model laid the groundwork for the rational exploitation of this aquifer by optimizing the abstraction rates of 7 wellfields in order to preserve the environment. Recommendations were made to stakeholders to improve the database (e.g., groundwater level monitoring on all wellfields) and thereby improve the reliability of the model. It is expected that the model will efficiently assist the stakeholders with respect to the sustainable planning and management of groundwater resources in Borena. Ethiopia has major agricultural development projects for rural area planning. The need for water is immense. The work presented here is part of the national effort to assess and optimize water resources, specifically groundwater resources.