Integrated Hydrological Analysis of Little Akaki Watershed Using SWAT-MODFLOW, Ethiopia

: In Ethiopia, groundwater is the main source of freshwater to support human consumption and socio-economic development. Little Akaki watershed is located in Upper Awash basin, known for its high annual rainfall and considered as the potential groundwater recharge zone. On the contrary, urbanization and industrial expansion are increasing at an alarming rate in the area. This became a concern threatening the groundwater resources’ sustainability. To address these challenges, integrated analysis of groundwater recharge and groundwater numerical simulations were made. For groundwater recharge estimation, SWAT model was used. The result indicated that recharge in the watershed mostly occurs from July to October with maximum values in August. On average, the estimated annual catchment recharge was 179 mm. For the numerical simulation and prediction of the groundwater ﬂow system, MODFLOW 2005 was used. The model simulations indicated that the groundwater head converges towards the main river and, ﬁnally, to the outlet of the watershed. The study indicated areas of interactions between the river and groundwater. The scenario examination result reveals increasing the present pumping rate by over ﬁfty percent (by 50%, 100%, and 200%) will surely cause visible groundwater head decline near the outlet of the watershed, and substantial river baseﬂow reduction. The recharge reduction scenario also indicates the huge risk of groundwater sustainability in the area.


Introduction
Groundwater is one of the main sources of fresh water to support human needs for domestic consumption and socio-economic development in Ethiopia. Sustainable use of the groundwater depends on understanding of the groundwater resources dynamics and its relationship with the other hydrological cycle components [1,2]. Most of the groundwater comes from rainfall that infiltrates down from the land surface. Infiltration of rainfall through the soil unsaturated zone is the source of water to the water table, groundwater recharge. The distribution of precipitation, which is the source of almost all freshwater in the hydrologic cycle is highly variable. Likewise, evaporation and transpiration rates differ significantly according to climatic conditions. This determines the recharge distribution, and the groundwater recharge-discharge dynamics, in addition to the hydrogeological feature of the aquifer [3,4].
Today, to sustainably manage the water resources, the consideration of groundwater and surface water as an interconnected and single resource has been increasing all over the world due to the increasing sustainability and environmental concerns [5]. The scientific concepts and aspects about the interaction of groundwater and surface water and the necessity of dealing with it in a unified way are widely addressed by USGS [6].
In such an interconnected water resources network, both surface water and groundwater resources need to be managed conjunctively to satisfy the current demand without compromising its sustainability for the future. To do that, a physically-based distributed hydrological models handling both the surface and groundwater hydrological processes can be used to estimate water availability under present and forthcoming conditions and to manage suitable integrated water management policies [7][8][9].
The central Ethiopia is the place where large number of populations lives (including the capital city of the country, Addis Ababa), and large number of industries are active. According to Kebede [9], the groundwater contribution to the surface water is high in the central Ethiopia, mostly in the Awash basin because of the high groundwater recharge contribution in the central highlands. The water use in the area are both from surface and groundwater sources. Groundwater contributes for more than 25% of the water supply for the capital city, Addis Ababa [10,11]. Due to the alarmingly increasing population of the capital city and the nearby cities, the Addis Ababa Water Supply and Sewerage Authority (AAWSSA) is expanding additional groundwater sources as sources of water supply for the city. In addition, significant amounts of groundwater are pumped by private industries, institutions and individuals without necessary scientific support.
To sustainably use and manage the groundwater resources in the area, spatial distribution, and amount of groundwater recharge, the mechanism of groundwater flow, occurrence, the interactions of groundwater with the river, and assessing the response of the system to different abstraction are the important issues to be addressed. These could be done with a reliable surface and groundwater hydrological model simulations, which also depend on thorough calibration and testing with the availability of quality hydro-climatic data [12].
As groundwater is vital to meet the growing water demand of population, dependable estimation of groundwater recharge is fundamental for the evaluation of the groundwater resources potential [3,13,14]. Since directly measuring the groundwater recharge like that of precipitation or streamflow is not possible, it needs to be estimated by indirect methods. In this study, the water balance method of groundwater recharge estimation is followed to estimate the catchment scale groundwater recharge of the Little Akaki watershed using the Soil and Water Assessment Tool (SWAT) model.

Study Area Description
Little Akaki catchment is one of the suppliers to the flow of the Akaki river. It is located in the west of Addis Ababa and within Oromia Region. It is located between UTM grids 453,334.5 m to 469,476 m East and 991,929 m to 1009,619 m north. The elevation of the catchment area varies from 2370 to 3368 m above sea level ( Figure 1). The physiographic components in this area are rolling plains, valleys, steep river banks, hills, and mountains.
The daily average temperatures in the area range from around 10 to 25 • C and the aerial average annual catchment rainfall is 1204 mm, as estimated at Addis Ababa Observatory. Like other parts of Ethiopian highlands, the climate of the Little Akaki catchment is described by two separate seasonal weather patterns. The main rainy season extends from June to September and contributing around 70% of the yearly total rainfall (Figure 2). A small rainy season from mid-February to mid-April contributes additional moisture to the region.
The geology of the Little Akaki watershed is a basic part of the growth and advancement of the Ethiopian Plateau and the rift system. The catchment system is protected by volcanic rocks covered by fluvial and residual soils, in which black cotton soils are the predominant, varying in thickness from a few centimeters to about 20 m. The leading lithologies include basalts, rhyolites, trachytes, scoria, trachy-basalts, ignimbrites, and tuff. The geological structures, fractures, and faults are mainly related to the rift tectonic extensions. These greatly weathered, broken lithologies favor the movement and storage of groundwater [11,15]. The aquifer properties and the hydraulic complexity of these volcanic rocks in the area are controlled by the litho-stratigraphy of the volcanic rocks and their complex spatial distribution, their different joint stratigraphic relationships, their substantial compositional, structural, and textural variability, and their diverse stages of tectonization and weathering [9,16].

Method of Recharge Estimation
SWAT model, a physically-based basin-scale spatially distributed continuous daily time step and computationally efficient hydrological model [15], is used for recharge estimation. It was developed to simulate the impact of land use (management practices) on quantity, as well as the quality of water, sediment, and agricultural chemical yields, in big complex watersheds [16]. ArcSWAT, the interface of the GIS divides the watershed into hydrological response units (HRUs) that contain similar characteristics of land use, soil, slope, and management. HRU is used to describe the spatial heterogeneity within a watershed [16,17]. The inputs data required for the GIS interface of ArcSWAT 2012 to simulate the stream flow includes Digital Elevation Model (DEM), soil and land use spatial data reclassified as per the SWAT format, and daily weather data from the representative stations.
A 30 m by 30 m resolution ASTER global DEM (ASTER GDEM) version 2.0 [17] was used. For the Land Use and Land Cover requirement, the open access global landcover map datasets (http://www.globallandcover.com/GLC30Download/download.aspx, accessed 15 January 2020), high-resolution (30 m × 30 m), released by the government of China [18] was used.
The soil data from the Harmonized World Soil Database (HWSD) v 1.2 [19] was extracted and used for this study. The spatial resolution of the data is about 1 km (30 arc seconds × 30 arc seconds).
The input meteorological data used to run the SWAT model were the conventional weather data from 1988 to 2004 from five stations in and around the Little Akaki sub-basin (Addis Ababa Bole, Addis Ababa Observatory, Holota, Ginchi, and Teji) obtained from the National Meteorological Agency of Ethiopia. The hydrological flow data, daily streamflow data for the Little Akaki gauging station was obtained from Ethiopia Ministry of Water, Irrigation, and Electricity department of Hydrology.
The observed daily streamflow data from 1992-2003 was used for the model calibration and validation. The calibration and validation analysis, the sensitivity, and uncertainty were used in the SUFI-2 algorithm in the SWAT-CUP [20]. The SWAT model simulation output was imported to SWAT-CUP for analysis. The SUFI-2 algorithm computes the 95% prediction uncertainty (95 PPU) range through an iterative procedure and attempts to capture the observed streamflow data within the band.

Groundwater Modelling Method and Procedures
MODFLOW-2005 was used to develop a 3-dimensional model grid used to characterize a two-dimensional groundwater flow through a single aquifer layer. For the use of finite difference calculations, a grid is overlaid over the 30 m × 30 m digital elevation model (DEM) of the study area, and aquifer hydraulic parameters necessary to solve the flow equation are been in an average of over the area of cell and given at a node at the center of the block. In the block-centered preparation, the nodes for which water levels are simulated placed at the midpoint of the grid cells. These cells are the smallest volumetric units over which the hydraulic properties are assumed constant.
The extent of the model area comprises the boundaries of the regional flow system of the Little Akaki River above the Little Akaki river gauging station. The model covers an area of about 134 km squares.
The model conceptualization process is an essential process to control the modeling technique to be used. The next important step, after building the necessary record of the database required for the modeling, is the development of a conceptual model. This is the most fundamental part of the modeling procedure [21,22]. The important data necessities in the course of conceptualization are the hydraulic characteristics of the aquifer, available surface water bodies, recharge value and zones, and area of outflow zones. The conceptual model for the study area was developed using GMS software (Aquaveo 2016). The coverages for the model boundary, source and sinks, recharge rate, and hydraulic conductivity zones were created.

Model Boundary, Local Source/Sink Coverage
The initial step in creating a conceptual MODFLOW model is to define the boundary of the model. The areal boundary was developed from a shapefile made by using the ArcGIS and then imported into GMS as GIS layer. The boundary of all the other layers was made by duplicating the boundary of the model and specifying the properties of each in the Coverage Setups.
After defining the boundary, the next step was to define the sources and sinks coverage. Under the source/sink coverage, the boundary conditions of the study area, such as wells, specified head boundary general head boundary, and river, were defined.

Initial Hydraulic Heads
Initial hydraulic head conditions refer to the groundwater level distribution in the system at the commencement of the simulation [23]. Distributed initial static hydraulic heads data collected from 11 observation wells of the current study are created by interpolating from the head values. The interpolated hydraulic heads distribution was then imported to MODFLOW at the start of the model run, and every cell in the model was designated its head value (Harbaugh et al., 2000). The land surface elevations at particular model cell are derived from 30 m × 30 m DEM of the area, and then depth to static water level values from the wells database is subtracted from each DEMs obtained. The results of this calculation are used as initial prescribed hydraulic heads in each cell for the initial specification of head values.

Hydraulic Conductivity
Groundwater flow rate in the saturated zone depends on the hydraulic conductivities of the aquifer. For the current study, the horizontal hydraulic conductivity values were obtained from the calibration process in PEST [24], considering that there is no available data. Zonation for hydraulic conductivity area was based on the available geologic and geomorphologic information. Then, the hydraulic conductivity map obtained from the calibration is overlaid on the model grid, and then respective average values are assigned to each model cells.

Boundary Conditions
Because of the boundaries mainly determine the flow pattern, an appropriate selection of boundary conditions is a crucial action in model design [23]. In groundwater flow, the system physical boundaries, including impermeable geologic formations and surface water bodies or hydraulic boundaries that include groundwater, divides, and flow lines, can be used as boundary conditions (Figure 3). The employed boundary condition in this study includes specified flow boundaries, head-dependent flow boundaries, river, and drain packages. A general head boundary package was used to simulate groundwater outflow through the outlet of the basin. Specific head boundaries were defined for Lake Gefersa (lake level of 2585 m is used). The catchment boundaries in the western, northern, and eastern parts of the catchment are taken as a no-flow boundary. General Head boundaries usually used at the edge of the model to enable groundwater to flow into or out of the model under a regional gradient. In the General-Head Boundary package, the flux is always proportional to the difference in head. The tributary rivers are considered with a drain package.
The recharge zone is obtained from the SWAT sub-basin classification, and the subbasin average values estimated was used in the study area. Using the estimated recharge as an input, a conceptual groundwater model for the area was developed. The recharge zone in the MODFLOW model is the same as the sub-basin of the SWAT model.

Hydraulic Conductivity Zone
The aquifer hydraulic conductivity is variable with the aquifer types. Multiple polygon zones were used to define the hydraulic conductivity. Since there is no reliable recorded data of hydraulic conductivity data in the area, it was estimated during calibration with the PEST automatic parameter estimation tool for the zones created based on geo-morphology of the area.

MODFLOW Grid
With the coverage's setup, the MODFLOW grid was created. The number of cells in the x and y dimensions were generated automatically because the number of rows and columns and the location of the cell boundaries were controlled by the refine point data at the wells. The grid dimensions of the model automatically assigned due to the refinement assigned for the vicinity of the pumping wells. The model divided into 228 rows and 232 columns with the size of each cell being varied from 40 m near the wells to 140 m. The grid cells outside and inside the model domain (Little Akaki catchment) were designated as "inactive" and as "active", respectively. The total number of active cells of the model were 29,840.

MODFLOW Model Calibration
A steady-state calibration was carried out for the model based on the groundwater level data available. The MODFLOW model was calibrated using an automatic calibration tool, PEST for hydraulic conductivity using the observed water table and estimated baseflow at the out late of the watershed. It is assumed that the stream baseflow (separated from the runoff component of the river streamflow) is sourced from the groundwater flow to the river. Based on this assumption, the groundwater flow model, MODFLOW is calibrated with baseflow obtained from streamflow, in addition to calibration with few available groundwater observation data, to improve the model performance.
The primary effective calibration practice for the fine-tuning of the hydraulic conductivity field in the model was to initially define fewer conductivity zones and then progressively raise the number of zones based on the available geomorphology of the area. The hydraulic conductivity values were continually adjusted during calibration in each hydraulic conductivity zone.

Recharge Estimation
The SWAT model performance for monthly streamflow simulation was evaluated based on basic statistics, such as Nash-Sutcliff (NSE), percent bias (PBIAS), Correlation coefficient (R 2 ), and root mean square error divided by the standard deviation (RSR). The p and the r-factors measures the fraction of observed data contained in the 95 PPU and the ratio of the average width of the 95 PPU band, respectively (Table 1). Using the best parameter values from the calibration process, monthly water balance components were simulated using the SWAT model for the year 1988 to 2004 (Figure 4). Average monthly variations of basin values for the various water balance elements estimated during the simulation period display average gain and losses on the monthly base from the watershed with the change in soil water storage ( Figure 5). Actual evapotranspiration (AET) from these components contributed a more considerable amount of water loss from the watershed.

Basin Average Recharge Rates
The average monthly and annual long-term groundwater recharge in the Little Akaki watershed was simulated by using the SWAT hydrological model. The average monthly recharge distribution indicates clearly that the recharge amount and timings are correlated with the precipitation amount and months. In the watershed, most of the groundwater recharge occurs during the wet/rainy season, from June and September of the summer season ( Figure 6). The average annual recharge estimated for the Little Akaki watershed is 179 mm (Figure 7), about 15% of the annual precipitation of the area (1204 mm). Based on that, the annual average volumetric groundwater recharge in Little Akaki watershed was estimated at about 24.3 Million cubic meters. The average annual simulated groundwater recharge is in line with the average groundwater recharge reported by Reference [9], 150 to 250 mm/year, less than the 265 mm recharge estimated by Reference [10] from selected points in the Akaki Catchment, and by far greater than the regional groundwater recharge estimation for Ethiopia, 39 mm/year [25].

Sub-Basin Average Recharge Rates
The spatial distribution indicates that the annual average sub-basin recharge for the 21 sub-basins are in the range of 96 to 240 mm ( Table 2). Regarding the spatial variability of the recharge estimated, as clearly seen in Figure 8, the estimated recharge amount is significantly low in urban areas, and high recharge values were obtained in areas dominated by agriculture and forest. In sub-basins 13, 10, 9, 6, and 15 (Table 2), the land uses are dominated by urban land use type, and the estimated recharge in those areas is lower than in other land use classes.

(a) MODFLOW Model Calibration
The MODFLOW model was calibrated using automatic calibration tool, PEST for hydraulic conductivity in the range of 0.001 to 25 m/day using observed water table and estimated baseflow at the out late. The best fit simulation results were obtained while the model region was divided into different hydraulic conductivity zones through trial and error considering the hydrogeology and geomorphology of the area. The hydraulic conductivity zones and results of hydraulic conductivity values from the steady state calibrations are shown in Figure 9. The hydraulic conductivity result from calibration shows higher values below the Gefersa Lake which could be due to the presence of geological faults [11]. The plot of simulated steady-state hydraulic head values versus the measured water levels at the observation wells shows fit by comparison to a 1:1 line in Figure 10.

(b) MODFLOW Model Performance Evaluation
The simulation results of the calibration would be assessed both qualitatively and quantitatively [23]. The calibrated outcomes were appraised based on the calibration objective and calculation of the mass balance of the system. A scatter plot of measured heads against simulated heads showed the calibration performance ( Figure 10). The scatter plots are observably examined to see how much the points in a plot differed from the straight line. Additionally, the calibrated model results were evaluated statistically by evaluating with the three common methods of error quantifications known as Mean error, Mean absolute error, and Root Mean Squared Error [26].
where n denotes to the number of observations considered, h s to the simulated head [L], and ho to the observed head [L]. The mean error refers to the variation between the observed and simulated heads. The mean absolute error indicates the mean differences of the absolute values of observed and simulated heads.
The calibration of the model was evaluated by applying all these objective functions (Table 3) highlighting the RMSE as it is the most reliable evaluation method as compared to the other two methods [26]. Therefore, the model performance evaluation to simulate the groundwater level shows statistically acceptable and unbiased. The water table varies from 2370 to 2710 m above sea level. Figure 11 shows the predicted groundwater contours using the calibrated model. Flow is generally to the southeast toward the Akaki well field. The simulated heads consistently indicated that the groundwater flows directions are toward the stream/river. It shows that flow converges from all sides towards the rivers. The simulated and the measured groundwater levels indicate that the water table from the surface tends to be deeper in elevated areas. The model estimated groundwater head distribution slightly higher than the regional groundwater head made for the Akaki basin [11]. The gains to and losses from streams indicate the interaction of the stream with an aquifer. The interactions were simulated using river and drain processes in the MODFLOW. Constant river stage was considered because the model is steady-state. The interaction between rivers and groundwater indicates that flow occurs from the rivers into aquifers above the Gefersa Lake in river 1 (R1) and above an elevation of 2375 masl in the second river (R2) (Figure 11), and from the aquifer to the rivers in the downstream courses of the rivers.

(d) Sensitivity Analysis
In groundwater modeling, there are unavoidable limitations that result in uncertainties in aquifer parameter estimation. Sensitivity analysis helps to understand the level of uncertainty in the calibrated model caused by those controls in the calculations of aquifer parameters. Aquifer hydraulic conductivity is one of the input parameters the groundwater models are sensitive to. Little variations in most sensitive parameters could result in a large change in model simulation outputs.
Parameter sensitivity analysis ( Figure 12) shows that the hydraulic conductivity value is highly sensitive in hydraulic conductivity zone HK_5, which is located below the Gefersa reservoir. The hydraulic conductivity value near the center of the watershed, below the Gefersa Lake, could be due to the presence of significantly different geological features, like faults, in those areas [11]. In this study, hydraulic conductivity is the only parameter obtained through the calibration process. The MODFLOW numerical model response to the deviation from the calibrated hydraulic conductivity during the simulation, while keeping other parameters constant was used as a tool to evaluate the sensitivity of the model to that particular parameter.
To check the sensitivity of the model to the parameters, continuous simulations were made by changing the calibrated values of hydraulic conductivity both in increases and decreases by 10%, 25%, and 50%. After the continuous run of six times, the corresponding root means squared error (RMSE) for head changes (%) in observation wells from the calibrated value were plotted in Figure 13. It was noted that the model was less sensitive to a small change in parameter value and more sensitive to the big changes of the parameters from the calibrated values.

(e) Groundwater Balance and Scenario Analysis
One of the fundamental ways to quantitatively estimate the flow of groundwater within an aquifer system is through a calculation of the water budget. The basic water budget equation is that the difference between inflows and outflows equals the change in storage of the groundwater system. The water balance is established based on the modeling water budget tool in MODFLOW. The water balance of the baseline scenario (Table 4) shows that inflow and outflow of the total Source/Sink are equal. The baseline groundwater pumping recorded considered was the static pumping rate obtained from 17 pumping wells (Figures 14 and 15). Because of the increasing groundwater withdrawals as a result of the increasing population and development activities in the area, five pumping scenarios were adopted to test the response of the aquifer system under variable groundwater abstraction rates. In addition to the pumping scenarios, increasing and reducing groundwater recharge scenarios are considered because of the recharge in the area is significantly dependent on land use type (reported in the recharge estimation section).  The considered pumping scenarios are the current pumping scenario (Scenario 1) and the increasing of the current abstraction by 10% (Scenario 2), 25% (Scenarrio 3), 50% (Scenario 4), 100% (doubling the current pumping, scenario 5), and 200% (increasing fourfold, Scenario 6) from the existing pumping wells (Table 5). Figures 14 and 15 summarizes the results of the groundwater response to the pumping scenarios. Simulation results of stream base flow, water table elevations (drawdown), and subsurface outflows during the scenarios were analyzed with model estimated results. The pumping scenarios result clearly indicates that increasing the pumping rate considerably altered groundwater dynamics and fluxes. Aquifer outflow from the head-dependent boundary and groundwater flux from the wells is the biggest change observed in response to the increased abstractions. Equal percentage increasing to pumping wells has a different effect on every well because of the abstraction rates are different and well proximity to constant head boundary and rivers. In some wells having a high rate of abstraction at the baseline scenario, imposing additional increasing scenarios significantly change the flow system, and a noticeable decline of groundwater level and drawdown were observed.
Considering the possibility of the recharge rate can be increased (if managed well) or reduced due to the changing climate and human effect on the recharge area, both increasing and decreasing recharge scenarios were considered. The considered recharge scenarios are increasing the current recharge rate by 10%, 25%, and 50%, and if the current recharge rate reduced by 10%, 25%, and 50% from the recharge rate estimated ( Table 5). The result indicates that increasing and reducing the recharge rate considerably changed groundwater dynamics and fluxes in a reverse way.
The result indicated that increasing pumping significantly affect the baseflow out of the watershed, which is equal to aquifer outflow from the head-dependent boundary. Similarly, low recharge also significantly affects the stream baseflow at Little Akaki gauging site.

Conclusions
In the Littlie Akaki watershed, areal recharge was estimated using SWAT model. The SWAT model simulates average monthly and annual groundwater recharge at sub-basin scale for 21 sub-basins in the Little Akaki watershed. SWAT model is a good tool to estimate the spatial and timely variability of the recharge. Recharge in the watershed mostly occurs in the months of July to October with maximum values in August. On average, the estimated annual catchment recharge was 179 mm.
In general, the groundwater recharge estimated was lower for sub-basins dominated by urban land covers. As the cities of Addis Ababa and other cities of the Oromia region in the area are located in a high rainfall area (potential recharge area), the expansion of the cities could significantly affect the recharge of the groundwater. This could be alarming for the sustainability of the increasing groundwater exploitation in the area as a result of a fresh water shortage.
The groundwater recharge estimated was used as an input for the MODFLOW model simulation. The model simulations indicate that groundwater converges towards the outlet, near the Akaki well field. The groundwater contours indicated that the water table of the groundwater varies with the topography. The hydraulic conductivity value near the center of the watershed, below the Gefersa Lake, shows high level of sensitivity to the water table and the base flow simulation during calibration, indicating there could be significantly different geological feature in those areas. The study clearly indicated areas of interactions between the river and groundwater.
The scenario examination result reveals a substantial risk of over-abstraction of the groundwater as a result of the increasing freshwater demand in the area. Increasing the present pumping rate by over fifty percent (50%, 100%, and 200%) will surely cause noticeable groundwater head decline near the outlet of the watershed, and substantial river baseflow reduction which is estimated from the aquifer outflow from the head-dependent boundary. The recharge reduction scenario also indicates the huge risk of groundwater sustainability in the area.

Conflicts of Interest:
The authors declare no conflict of interest.