Groundwater Modelling of Motloutse Alluvial Aquifer, Eastern Botswana ()
1. Introduction
1.1. General
Botswana is generally an arid country with little surface water except in the far north. Much of the country depends on groundwater. Many regions of Botswana are characterized by low and erratic rainfall varying from an annual minimum of less than 300 mm in the southwest of the Kgalagadi District to over 650 mm in the northern part of the country near Kasane [1]-[3].
Most rural areas in Botswana are supplied with groundwater, which is used for multiple purposes. Sand rivers present the availability of this resource in the eastern parts of the country where this study was conducted. Rainfall over the Motloutse catchment area occurs in the summer, and the average monthly evaporation exceeds the mean monthly rainfall. Communities within the river catchment use the groundwater from the sand river as a water source for irrigation to sustain their agricultural needs.
Geologically, the area comprises granitic rocks overlain by alluvial sediments [4]-[8]. The alluviums consist predominantly of gravel and sand of variable thickness, with gravel and coarse sand in the middle of the river and fine-grained sediments occurring adjacent to riverbanks. According to Edwin [7] and Nord [9], the Motloutse sediment is made up of poorly graded gravelly sand. The alluviums in the flood plain are the aquifers that serve as a groundwater source for the community in the study area. The aquifer is phreatic and recharges from direct percolation and river flow.
Even though the phreatic aquifer serves as a water source for irrigation, the farmers in the study area rely on groundwater, and they do not know how much groundwater resources that are available. Moreover, the demand for more water is increasing yearly because of the increasing number of farmers involved in irrigation schemes. Furthermore, the hydraulic characteristics of this aquifer are not properly determined by considering the areal distribution and their thickness throughout the study area. To support the ongoing effort of the farmers in their irrigation activities and thereby improve their livelihoods, investigation of this phreatic aquifer potential is indispensable. This study was conducted to generate important information regarding the Motloutse phreatic aquifer that can be utilised to plan sustainable irrigation development in the area. The objective of this study was to provide a basic groundwater flow model to quantify and evaluate the resource potential on a one-kilometre river stretch scale.
1.2. Location
This study was conducted in Eastern Botswana along the Motloutse River, bounded between UTM coordinates of grid references 600,000 – 620,000 m E and 7,590,000 - 7,570,000 m S (Figure 1). It is downstream of the Letsibogo dam, near the Bobonong, Mmadinare and Tobane villages.
2. Methodology
The hydraulic properties of the alluvial aquifer system investigated previously in the Motloutse River catchment by the same authors were used as input parameters in this study. In their research, Edwin et al. [10] conducted a Ground Penetrating Radar Survey (GPR) to determine the thickness and porosity of the aquifer, used
Figure 1. Location map of the study area.
laboratory method to determine specific yield for the riverbed and riverbank sediments, utilized the Alyamani-Sen [11] empirical formula to determine the hydraulic conductivity of the riverbed sediments, and performed slug tests to determine the hydraulic conductivity of the riverbank sediments.
This study implemented a simple steady state, finite difference groundwater model to simulate and plot groundwater flow to improve the understanding of groundwater flow in the Motloutse sand river aquifer. The steady-state simulation assumes pre-development conditions (those prior to pumping) in the study area.
Groundwater modelling was undertaken using a three-dimensional groundwater flow code MODFLOW [12]. This modelling code solves the general governing equation (Equation 1) for heterogeneous and anisotropic unconfined aquifers [13] [14];
(1)
Where K is the Hydraulic conductivity; Ss is the Specific storage; h is the Piezometric head, and t is the Time.
Table 1. A summary of hydraulic properties of the Motloutse sediments [10].
No |
Parameters |
Determined values |
1 |
Thickness of the alluvium |
Varies from 9 – 12 m |
2 |
Porosity of sediments |
40% |
3 |
Riverbed sediments-specific yield |
13.68% |
4 |
Riverbank sediments-specific yield |
8.84% |
5 |
Riverbank sediments hydraulic conductivity |
26.43 m/day (3.059 × 10−4 m/sec) |
6 |
Riverbed sediments hydraulic conductivity |
160 m/day (1.856 × 10−3 m/sec) |
3. Results and Discussions
Groundwater Modelling
Conceptual Model
The model area (Figure 1) is 1000 m in length by 135 m, which is defined by the average river width at that river section under investigation. The vital hydrogeological processes in the study area include areal recharge from precipitation, stream inflow, and discharge, such as river outflow, evaporation, and evapotranspiration.
Based on the geology (Figure 2) and geophysical data [10], the model domain was conceptualised as a single-layer unconfined aquifer system (Figure 3). The productive aquifer, which is sand, does not have a uniform thickness throughout the model area. Its thickness ranges from 4.8 to 6.4m, with an average thickness of 6 m.
Figure 2. Cross section of Motloutse alluvial aquifer in the model area.
Figure 3. Conceptual model of Motloutse alluvial aquifer.
Model Discretization
A uniform grid is laid over the phreatic aquifer of the model area (Figure 4), and a single-layer model with an overall size of 50 columns x 54 rows was developed for the Motloutse River phreatic aquifer. To effectively define the resource the river boundaries were defined as accurately as possible, and the aquifer geometry is governed by top and bottom elevations incorporated from the results of the geophysical survey [10].
Figure 4. Model boundary conditions (X-axis and Y-axis are in UTM coordinate system representing easting and southing, respectively).
Boundary Conditions and Assumptions
Groundwater model boundaries effectively dictate the flow direction and influence the water balance of a numerical model; hence, the type of boundary used depends on the modelling scope and model purpose. As shown in Figure 4 above, the active areas of the model are bounded by no-flow cells. General Head Boundary (GHB) nodes were placed just inside a no-flow boundary representing lateral subsurface flow, with a no-flow boundary representing lateral subsurface flow, which is anticipated for such an ephemeral river. Consequently, and this will allow groundwater flow inside or outside the model domain as baseflow proportional to head differences.
Unnatural boundaries were allocated to the west and eastern boundaries for the purpose of modelling. The western boundary is located at the contact between gneiss formations and is conceptualised as a specified flux boundary. This upstream section was lined up with injection wells to simulate groundwater flow into the model, whereas the east is a general head boundary. General Head Boundary elevations were set to match available water level data. The water table is not fixed at this boundary and may change by some stresses in the basin hence General Hydraulic Boundary (GHB) was selected.
Recharge and Discharge
Rainfall - The Motloutse alluvial aquifer is an unconfined aquifer that contains fresh groundwater of recent age, recharged through precipitation almost annually by the episodic flow of the Motloutse River. The precipitation was derived from meteorological data for each month with rainfall. The higher conductivity and storativity of the Motloutse River aquifer unit allows it to be readily and rapidly saturated in the “wet” months. Recharge was spatially applied to the model domain or simulated using MODFLOW’s recharge package as 115 mm/yr. This estimated value is a quarter of the long-term mean annual precipitation of the area, representing the expected recharge on the model area for the 2014 wet season, during which the water levels were used for calibration.
River Inflow - Hydraulic gradients are low within the aquifer [9], resulting in low groundwater velocities and subsurface throughflow in the aquifer. This inflow was simulated by flux through ten injection wells. Darcy’s law was used to estimate the pumping rate of the injection wells. Riverbed sediments hydraulic conductivity was estimated as 160 m/day (Table 1). The estimated total pumping rate of each injection well was 21 cubic meters per day resulting in 210 m3/day for the ten injecting wells as per the number of cells in the western boundary.
River Outflow - Subsurface outflow at the downstream section has been represented by a general head boundary of which its hydraulic head will rise and fall according to the aquifer conditions in and out of the model domain.
Evaporation - After a river flow event or flood, the aquifer is discharged when the sand becomes unsaturated locally due to falling groundwater levels, mainly due to subsurface evaporation. Evapotranspiration presents another potential loss by the sparse riverine vegetation. However, there are no visible deep rooted plants growing within the channel, so an effective extinction depth was taken as 1 m [9] below which direct evaporation becomes insignificant.
MODFLOW’s inbuilt evapotranspiration (ET) package incorporates and simulates the effects of direct evaporation in removing water from the saturated groundwater regime [12]. This package requires 3 components: a maximum rate of evapotranspiration (ET) (L/T) is specified for each surface cell, ET surface elevation (L) and an extinction depth (L). When the water table depth is at or above the ET surface, evapotranspiration occurs at the maximum specified rate, decreasing linearly with a decline in the water table. When the water table lies below the extinction depth, zero evaporation occurs.
Evaporation rates are at a maximum when the water table is at the sand surface therefore, an estimated starting value of 635 mm/year from the MacDonald and Partners studies (Table 2) has been applied [15]. Measured pan evaporation is useful for estimating potential ET when the water table is near land surface, but actual maximum ET rates may be overestimated or underestimated. The ET surface was taken as the average sand elevation of 778 m a.m.s.l.
Table 2. Month and open water evaporation (mm).
Nov |
Dec |
Jan |
Feb |
Mar |
Apr |
May |
Jun |
Jul |
Aug |
Sep |
Oct |
209 |
211 |
215 |
185 |
183 |
145 |
118 |
93 |
105 |
139 |
180 |
217 |
Abstraction - Another obvious potential loss is abstraction by hand-dug wells within the sandy aquifer from which unknown quantities are harvested. Still, there was no evidence for this in the modelled area.
Aquifer Parameters
Due to the difference in the nature of the sediments in the riverbanks and in the riverbed, two hydraulic conductivity zones are recognized in the model area (Figure 5): K for zone 1 and K for zone 2. The hydraulic conductivity (K = 160 m/day) for the sediments in the riverbed (Table 1), which is zone 1, was used and assumed to be isotropic in all directions. The isotropic assumption was deemed reasonable for simplification as the aquifer mainly comprises sand. The hydraulic conductivity (K = 26 m/day) for the sediments in the riverbank (Table 1), which is zone 2, was used in this model. A porosity of 40% and an average specific yield of 13.68% were also used in this model (Table 1).
Storativity, or the storage coefficient is the volume of water released from storage per unit decline in hydraulic head in the aquifer per unit area of the aquifer [13] [14] [16]. For an unconfined aquifer, storativity is approximately equal to the specific yield, and ranges between 0.02 - 0.30. In this model, an initial value of 0.30 was used.
The groundwater table occurs at a depth of 0.3 to 0.8 m; ±0.2 m below ground level in the sand river during the wet season.
In summary, the aquifer was modelled as: unconfined with uniform aquifer properties; has an impermeable base; recharges during the rainy season and discharged during the dry season, i.e. recharge is zero during the dry season; and, it is sloping, and gravity is the driving force.
Figure 5. Spatial distribution of Hydraulic Conductivity on model domain.
Calibration
The initial estimates of recharge, hydraulic parameters, and boundary conditions used in the model were based on the above-mentioned values. Most of these values were adjusted accordingly within the range that would be expected in the Motloutse sand aquifer for similar conditions or materials. The model was iterated until the computed solution matched the observed values within an acceptable level of accuracy by trial and error.
The steady-state Motloutse River sand aquifer groundwater-flow model was calibrated using the 2014/2015 wet season climate conditions. This period represents slightly wet conditions as the area received above-average precipitation in 2014. However, most parameters employed were average and conservative. Therefore, all results and alternatives derived from the model simulations are an average representation of the system.
The hydraulic-head data used for calibration consisted primarily of water-level measurements in nine (9) hand-dug wells between November 2014 and January 2015. Five of these wells were dug in the riverbed, and the remaining four were in the riverbank. Throughout the calibration process, no adjustments were made that conflicted with the general understanding of the geology and hydrology. A plot of modelled against observed heads is presented in Figure 6 with the data points representing hydraulic heads. The 45˚ line is the reference line, with points lying exactly on this line, depicting a perfect fit. Data points above the line reflect that the model is over-predicting the hydraulic heads in the system, while those below the line reflect that the model is under-predicting the hydraulic heads in the system.
Figure 6. Comparison of Observed vs Simulated heads.
The calibrated model fits the data reasonably well, implying that the conceptual model of flow and the representation of hydrogeology in the groundwater flow simulation model are reasonable. However, the overall result of the model was comparable with the measured well data, and some observations, as in Figure 6 have poor convergence, being a bit further from the 45-degree line. This can be attributed to errors in water level data collection, or the model does not reflect the geology and hydrogeology in that location and may also be due to the small model scale. The output statistical parameters derived from the calibration are shown in Table 3.
Table 3. Errors of the calibrated Model.
Error (m) |
Mean Error |
Mean Absolute Error |
Root Mean Square Error |
Correlation Coefficient |
−0.008 |
0.089 |
0.126 |
0.997 |
The evaluation of the calibrated model result shows that the water balance discrepancy was almost zero (Table 4). The overall results of the groundwater model are comparable with the measured well data and in agreement with the conceptual model.
Figure 7 shows the output contour map of the hydraulic heads and the expected flow path direction along the sandy aquifer. The general hydraulic gradient in the Motloutse River follows the surface topography, and the gradient is west-east, which agrees with the area’s conceptual model.
Figure 7. Simulated head distribution and flow direction and paths.
Water balance
The groundwater budget can be quantified based on the calibrated model output, and the balancing of inputs and outputs is expected for any accurate groundwater flow model. A comparison of inflow and outflow components of the Motloutse alluvial aquifer groundwater flow system is shown in Table 4. Calibrated net aerial recharge and evaporation of the Motloutse alluvial aquifer is 172 mm/year and 120 mm/yr, respectively. Calibrated flux through the western boundary is 60 m3/day for each of the ten injection wells, and calibrated K values are 145 m/day and 11 m/day for the riverbed and riverbanks, respectively.
Table 4. Wet season water budget.
Flow term |
Inflow (m3/day) |
Outflow (m3/day) |
Recharge |
102.29 |
|
Flow through specified flux boundary |
540 |
|
Evaporation |
|
70.15 |
Head-dependent flow through the eastern boundary |
|
571.85 |
Total |
641.995 |
642.29 |
Per cent discrepancy |
0.05% |
Sensitivity Analysis
A sensitivity analysis enables an assessment of the overall model performance by quantifying the sensitivity of the model simulations in the calibrated model caused by uncertainty in the estimates of aquifer parameters, stress and boundary conditions [13]. The system’s response to variations in the calibrated recharge, hydraulic conductivity zones and evaporation was evaluated. The parameters varied by a factor of 0.25, and the resulting heads were used to compute the mean error (ME), mean average error (MAE), and root mean square error (RMSE) for each factor. These were compared against the calibrated and the departures plotted in Figure 8(a) to Figure 8(d). Figure 8(a) shows high hydraulic head errors when the hydraulic conductivity of zone 1 of the calibrated model is decreased, showing that it is sensitive and susceptible to a decrease in hydraulic conductivity than its increment. For K for zone 2, the same response is repeated with the mean error decreasing with K increasing (Figure 8(b)). Figure 8(c) shows that the model is sensitive to increases and decreases in recharge and generally generates a non-linear sensitivity response. However, this response is sharp near the calibrated value within a single factor of 0.25 and steadier as recharge is adjusted further from the calibrated value. The calibrated model is more sensitive to evaporation (Figure 8(d)) increment than its decrease. It is concluded that the model is more sensitive to decrease and increase in recharge than adjustments in both hydraulic conductivity and evaporation.
The model’s sensitivity to the applied boundary conditions was also tested. The general head boundary was replaced with a constant head boundary, and from the computed water budget, there was a negligible change of 0.5%, showing that the model is insensitive to the general head boundary or constant head boundary.
Figure 8. (a) Sensitivity plot of the calibrated model with respect to K zone 1; (b) Model sensitivity to K zone 2; (c) Plot of model sensitivity to recharge; and (d) Plot of model sensitivity to evaporation.
Sensitivity analysis on the modelled recharge indicates a similar response when increasing or decreasing the recharge value. This is expected since equally distributed recharge is applied over the entire model domain and the heads in the model respond linearly to increases or decreases in the recharge.
Resource Quantification
The estimated sustainable extractable volume for the Motloutse River 1 km stretch is 29,400 m3. Nord [9] has calculated the groundwater resources of Motloutse sand river aquifer on the basis of total volume of groundwater per kilometre length of river and found 29,000 m3, which compare well to the value obtained in this study.
Model Limitations
In groundwater modelling, some part of the data is used for model calibration, and the remaining part is used for model validation. It is may be impossible to validate a model because usually too short set of observed data is available, which is already required for calibration; therefore, model validation was not fully accomplished here.
4. Conclusions
The groundwater potential of the Motloutse water course has been determined successfully using the input parameters determined using different techniques. These values are considered realistic for a gravelly sand aquifer, with the calibrated groundwater model producing K values of 145 and 11 m/day, respectively. Porosity values and specific yield have been calculated as 40% and 13.68%, respectively. Calibrated net aerial recharge and evaporation of the Motloutse alluvial aquifer is 172 mm/year and 120 mm, respectively. The Motloutse alluvial aquifer yields a safe yield of 29,400 m3 for a kilometre of river stretch.
In the study area, areal recharge, evaporation, subsurface inflow and outflow are the main processes that define the water balance components of the sand riverbed aquifer. In addition, it has been shown that direct rainfall is likely to recharge the aquifer fully during a wet season.
The groundwater potential estimated in this study is the minimum water available per km stretch of the river at Tobane, and more water is expected from the riverbed a longer stretch of the river.
Acknowledgments
The University of Botswana sponsored the research work through the Office of Research and Development (ORD) and is greatly acknowledged. Anonymous reviewers are thanked for the constructive comments that shaped the manuscript