Numerical Groundwater Flow Modeling of Dijil River Catchment, Debre Markos Area, Ethiopia

Dijil River catchment is a sub-catchment of the Abay drainage basin and covers 138.28 km. This paper presents numerical groundwater flow modeling at steady-state conditions, in a single-layer aquifer system under different stress or scenarios. A numerical groundwater flow models represent the simplification of complex natural systems, different parameters were assembled into a conceptual model to represent the complex natural system in a simplified form. The conceptual model was input into the numeric model to examine the system response. Based on geologic and hydrogeological information, confined subsurface flow condition was considered and simulated using MODFLOW 2000. The model calibration accounts matching of 24 observation points with the simulated head with a permissible residual head of ±10m. The sensitivity of the major parameters of the model was identified during the calibration process. According to the simulated water budget in the model, the simulated inflow is found to be 1.2791870E+05 m/day which is nearly equal to the simulated outflow of 1.2791755E+05 m/day with the difference being only 1.1484375E+00 m/day. Water budget analysis reveals that outflow from river leakage accounts for 92.8 % of the total outflow and 14.1 % of the total inflow comes from the river leakage in the study area. Three scenarios of increased withdrawals and one scenario of altered recharge were used to study the system response. Accordingly, an increase in well withdrawal in scenario-I (existing wells pump simultaneously), scenario-II (existing drilled wells yield withdrawal increased by 30%), and scenario-III (additional eight wells having expected yield of 30 l/s drill and pump) resulted in an average decline of the steady-state water level by 1.06m, 1.68m, and 4.46m, respectively. They also caused the steady-state stream leakage to be reduced by about 2.93%, 4.58%, and 11.23%, and subsurface outflow by 9.41%, 14.67%, and 37.86%, respectively. A decrease in recharge by 25% and 50% results in a decrease of the head by 6.1m and 13.4m respectively, and a stream leakage decrease by 20.3%, and 40.3% respectively as compared to the simulated steady-state value. Therefore, adequate groundwater level monitoring wells should be placed in the catchment to control the total abstraction rates from the aquifer and fluctuations in groundwater levels.

hypothetical flow situations and then gain a general understanding of what type of flow system persists (Anderson and Woessner, 1992).
The use of numerical groundwater flow models in the general groundwater resources management and in evaluating the spatial variations of the important aquifer parameters represents one of the most trustworthy methods of regional hydrogeological investigations (Anderson and Woessner, 2002). Geologic information like geologic maps, geological cross-sections, and well logs together with data on hydrogeologic properties are useful to describe hydrostratigraphic units for the development of a conceptual model (Anderson and Woessner, 2002).
Groundwater models have been applied to four general types of problems: groundwater flow, solute transport, heat flow and aquifer deformation (Wang and Anderson, 1982;Trescott et al., 1976;McDonald and Harbaugh, 1988;Pinder and Gray, 1977;Konikow and Bredehoefet, 1978;Prickett and Lonquist, 1971). Over the last decade, the concern of modeling groundwater flow and transport systems has grown promptly (Zhou and Li, 2011;Huizar-Alvarez et al., 2016;Havril et al., 2017;Mussa et al., 2020). The application of numerical simulations in 2-and 3dimensional models to distinguish the various types of groundwater flow systems is becoming increasingly important for assessing potential groundwater management strategies and simulating under the steady condition and any changes in groundwater budget components in transient conditions in complex basins (Cardenas and Jiang, 2010;Bresciani et al., 2016;Mussa et al., 2020).
Groundwater modeling is being used more frequently as a tool to help answer optimum water management questions because it can lead to a better understanding of how the real system behaves and it can be used to make predictions about the system's future behavior. Today, groundwater flow modeling is the most widely used to solve complicated hydrogeological issues.
The foremost importance of numerical modeling is its speed, accuracy of results and reliability of executed calculations (Batu, 2006;Kresic, 2007;Anderson et al., 2015). Groundwater evaluation and management in many countries can be easily handled using an interface called MODFLOW (Wang et al., 2008;Veiko et al., 2013). The easiness, simple structure and ability to handle with separate package for resolving several special hydrogeologic problems can be the comparative advantages of MODFLOW to be applied as a calculation program (Tenbus and Fleck, 1996).
The interaction between groundwater-surface water in the study area has been done by many researchers. Surface water from the Tana Lake percolates through its floor to aquifers of Tertiary volcanic and Mesozoic sedimentary rocks (WWDSE, 2007) and based on baseflow separation and numerical modeling techniques, groundwater contribution is less than 7% of the total inflow to the Lake Tana (Seifu, 2004;Getachew, 2008;Samson, 2010) and groundwater inflow to the lake is insignificant (SMEC, 2007). Abay and Tana basin have various groundwater types that were recorded as a result of hydrogeochemical and isotope studies (Seifu et al., 2005 and. Even though a number of studies have been undergone in the study area on local and regional levels, the groundwater flow system of the study area is not adequately characterized. Moreover the yield and static water level of the study area are decreasing from time to time. According to Debre Markos town water and sewerage office (2016 and 2017), the yield of and the water level decreases by 13l/s and 31 meters on average, respectively. Therefore, this study aims to carry out numerical groundwater flow modeling at steady-state condition, in a single layer (1-D) aquifer system under different stress or scenarios.

STUDY AREA
The study area, Dijil River catchment, is located in the northwestern part of Ethiopia within the Blue Nile/Abay basin (Fig 1). Geographically, it lays between 347122 mE to 362601 mE and 1139983 mN to 1161150 mN, having an area of 138.28 km 2 and a perimeter of 82 km. The study area is characterized by a cool climate with a mean annual temperature range of 10°C to 15°C and situated in the altitudinal range of 2300 masl to 3300 masl. The climate of the area can be categorized into two broad seasons: the dry season which covers the period from October to May and the wet season extending from June to September, with slight rainfall during autumn and spring (Tenalem and Tamru, 2001). The mean annual precipitation of the area calculated using Thiessen polygon is 1354.15 mm.

GEOLOGY AND HYDROGEOLOGY
The study area consists of six geologic formations: the Middle (Yujbe) basalts, the Ignimbrite and tuff deposits, the Debremarkos basalts, the Lumame basalts, the Robgebeya basalts and the quaternary eluvial deposits (Fig 2). These geologic formations are described from the oldest to the youngest. The Middle (Yejube) Basalt is dominantly composed of Olivine-plagioclase phyric basalt, plagioclase phyric basalt and at places, it also consists of pockets of pyroxene-plagioclaseolivine phyric basalt and interlayered with pyroclastic (Wolela, 2002;Getaneh, 1991;Ayalew et al., 2015). The Ignimbrite and Tuff are also found intercalated with the basaltic rock of the area. It is composed of tuff, agglomerate, and ignimbrite. The color of tuff shows great variation from light to pinkish-gray and it is fine-grained. The Debremarkos Basalt (Wolela, 2002;Getaneh, 1991;Ayalew et al., 2015) is aphanitic to porphyritic and at places it is vesicular olivine-plagioclase phyric basalt, plagioclase-olivine phyric basalt and olivine phyric basalt, and pyroxene phyric basalt, grading to each other; with minor trachytes and pyroclastic tuff interlayered with the basalt at the top and bottom of the formation. Its topmost part is the hill forming pyroclastic tuff. The Lumame Basalt is characterized by variable phenocryst proportion, olivine-plagioclase phyric basalt, plagioclase phyric basalt, pyroxene plagioclase phyric basalt, pyroxene-olivine phyric basalt and pyroxene phyric basalt (Wolela, 2002;Getaneh, 1991;Ayalew et al., 2015). The top part is overlain by the Choke shield volcano. The Robgebeya basalts are composed of varying phenocryst proportion; plagioclase phyric basalt, olivine phyric basalt, olivine-plagioclase phyric basalt, olivine-pyroxene-plagioclase phyric basalt and pyroxene-olivine phyric basalts. The Quaternary Eluvial sediments are reddish-brown to black clay soils and commonly cover the lowland plain part of the study area (Wolela, 2002;Getaneh, 1991;Ayalew et al., 2015). The sediments are derived from chemical weathering of underlying volcanic rocks or parent rock materials in the upstream of the catchment.
Most of the faults and fractures in the area are trending North-South, East-West, and Northeast-Southwest (Fig 2). The volcanic aquifers show both the primary and secondary porosity and permeability. The primary porosity and permeability of the volcanic aquifers include vesicles or gas cavities and porous flow textures. However, weathering and fracturing also play a great role in the groundwater potential of the volcanic aquifers. The volcanic aquifers are characterized by using twelve pumping test data from drilled wells (Fig 2). The yield of the wells mainly penetrating the volcanic aquifers ranges from 4 l/sec to 39 l/sec ( Table 2). The pumping test data analyses were made by using the aquitest-v4 software and the results show that the transmissivity ranges between 4.6 m 2 /day -768 m 2 /day (Table 2 and Fig 2). Based on the aquifer classification by Krasny (1993), the transmissivity values of most of the aquifers lay in the moderate to high transmissivity category (Tables 1 & 2  The partially penetrating boreholes in the sediments and the upper weathered and fractured aquifer and fully penetrating boreholes that penetrated into the volcanic aquifers do not show match difference in their depth to groundwater level indicating that the different litho-stratigraphic units can be considered as hydraulically interconnected. Therefore, a single layer, 2-D, steadystate condition, and confined aquifer layer was considered for the modeling practice in this study. The quantification of the groundwater recharge in the study area was made based on precipitation distribution, geologic character, soil type, land use/land cover and topography of the area by using water balance method (De Vries and Simmers, 2002;Lucas et al., 2012;Diniz Melo et al., 2015) and it is estimated to be around 267.02 mm/annum (7.3 x10 -4 m/day).

METHODOLGY
In this research, the appropriate computer code/ MODFLOW-2000 (McDoland andHarbaugh, 1988 as developed by USGS) was used to simulate the numerical groundwater flow system in the study area using Processing MODFLOW program. The MODFLOW is widely used, tested, and verified software. It is a water balance computation incorporated model that simulates all hydrologic features independently using its grouped packages. It is a deterministic model approach that assumes the stage or response of aquifer is predetermined by the help of physical laws governing the groundwater flow.
The Arc GIS 10.1 was also used to develop the location map and the geological and hydrogeological maps of the study area. Landsat images were used to delineate the study area boundary and conceptualize the boundary conditions. The surface water divide of the study area is assigned as a no-flow boundary, recharge, and well discharge are represented using specifiedflux boundaries, and rivers are represented as head-dependent flux boundaries. Recharge to the groundwater systems was estimated using the water balance method (De Vries and Simmers, 2002;Lucas et al., 2012;Diniz Melo et al., 2015). Assumptions are made to derive the water balance equation of the study area surface water divide coincides with groundwater divided hence the groundwater in the basin is in closed system, the change in the storage on annual basis is assumed to be zero surface and subsurface water exchange with neighboring catchment is assumed to be zero, assuming no artificial diversion from other catchments, since the computation is made on annual basis, a net change of the soil moisture and groundwater storage is assumed to be zero.
Aquifer properties were determined using Aquitest-4 software from the pumping test data.
The volume of water that leaves as surface runoff from Dijil river catchment is 396.11mm/year calculated by the Rational Method (Chow et al., 1988). A rational approach is to obtain the yield of a catchment by assuming a suitable runoff coefficient. (Eq. 1). The value of the runoff coefficient C varies depending upon the soil type, vegetation, and slope.
Where, A = area of catchment, P = precipitation, C = runoff coefficient.
Calibration of the model was done using the field measured heads, by adjusting parameters based on trial and error method. The Mean Error (ME), Mean Absolute Error (MAE) and Root Mean Squared Error (RMSE) Error were computed to express the average difference between simulated and measured heads are commonly used (Duffield et al., 1990;Ghassemi et al., 1989;Konikow and Bredehoefet, 1978;Anderson and Woessner, 1992).
ME is the mean difference between the measured heads (hm) and simulated heads (hs).
Where, n-is the number of calibration value at a given target (i).
MAE is the mean absolute value of the difference in measured and simulated heads.

Simulation Code
The numerical groundwater flow modeling of the Dijil river catchment, the general 3-D governing equation has been adjusted according to the prevailing field condition (Eq. 5). During the conceptual model development, it is identified that there is no significant variation in water level and recharge in the catchment. Therefore, the conceptualized model is a steady-state, two dimensional confined aquifer and a single layer with no possible flow in the Z direction; the equation can be simplified and rewritten into the following equation (Eq. 6). Where, W is a general sink or source basically positive to represent recharge and negative for withdrawal of groundwater from aquifer system.

Spatial Discretization and Grid Layout
A grid with a small number of nodes is preferred to minimize data handling, computer storage and computation time (McDoland and Harbaugh, 1988;Anderson and Woessner, 2002;Kresic, 2007;Anderson et al., 2015). Two sets of parallel lines form a grid and are orthogonal. These lines form blocks which are called cells and the node is the point at which the model measures hydraulics.
Hydraulic and hydrogeological properties are presumed to be consistent over a cell's length so that its node represents the cell (Batu, 2006;Anderson et al., 2015).
In this work, the model discrete with a uniform cell size of 100m by 100m with 211 row and 154 column arrays in a single layer. The model uses a total number of 32,494 cells in which only 13,828 are active and used to compute the hydraulic head and the remaining number of cells are inactive cells (Fig 4).

Boundary Conditions
Boundary conditions are mathematical statements specifying the dependent variable (head) or the derivative of the dependent variable (flux) at the boundaries of the problem domain (Anderson and Woessner, 1992).
The model under study considers the volcanic ridge surface water divide as the no-flow boundary and the same time the decreasing of the conductivity down the depth of the aquifer as bottom no-flow boundary. The model also considers specified head boundary which is simulated by setting the head at the relevant boundary nodes equal to known values. The southeastern and southwestern part of the study area is considered as a specific head boundary. Three types of boundary conditions were used to define the groundwater flow system in Dijil river catchment: no-flow boundaries, specified-flux boundaries, and head-dependent flux boundaries (McDoland and Harbaugh, 1988;Anderson and Woessner, 2002;Kresic, 2007;Anderson et al., 2015).
Geologic or hydrologic barriers to groundwater flow were simulated using no-flow boundaries. In the study area, the contact between the permeable groundwater flow system and nearly impermeable bedrock is an example of a no-flow boundary. Known or estimated hydrologic fluxes, such as recharge and well discharge, are represented using specified-flux boundaries (Fig 4).

Model Calibration
There are two ways of finding model parameters to achieve calibration; manual trial and error adjustment of parameters and automated parameter estimations (Kresic, 1997;Anderson, 1982, Anderson andWoesser, 1992). In this model, calibration was performed by the traditional trial and error processes in which model parameters and hydrologic stresses were adjusted manually within reasonable limits of existing data and field hydrogeological observation to achieve the best model fit. Additionally, hydrologic stresses literature review and point hydraulic conductivity and transmissivity data were used as initial guessing during calibration of hydraulic conductivity. The calibration criteria set by the modeler is to match simulated groundwater head and hydraulic gradient with an estimated one by taking into account the variation of horizontal hydraulic conductivity of aquifers in a small distance. Water level measurement data (primary and secondary data) of 24 water points have been used to calibrate the model.
The qualitative methods of calibration evaluation that are applied in this study are the comparison of the observed and simulated head contour, post errors in the contour and the scatterplot. A comparison between contour maps of measured and simulated heads (Fig 5) provides a visual, qualitative measure of the similarity between patterns, thereby giving some idea of the spatial distribution of error in the calibration.  (William, 2007). In the present study, the resulting model has a mean error of -2.1 m that is not too far to be a good calibrated model, an absolute residual mean of 6.43 m which is below the residual criteria set before the calibration process and has resulted 7.31 m root mean squared error of standard deviation which is a good range in calibrations. For an ideal calibration, the scatter plots where observed values are plotted against the value computed by the model, the points will fall on the straight line with a 45 degree slope that means the computed value equals the measured value. The correlation between observed and simulated values is significant and the determination coefficient is high: R 2 = 0.99 (Fig 6). Figure 6. Comparison of the model computed and observed hydraulic head data after model calibration.

Sensitivity Analysis
The purpose of sensitivity analysis is to quantify the uncertainty in the calibrated model caused by uncertainty in the estimated values of the aquifer parameters, stresses, and boundary conditions (Anderson and Weossner, 1992). It is the process of varying model input parameters over a reasonable range (range of uncertainty in values of model parameters) and observing the relative change in model response. The determination of a sensitivity analysis is to quantify the uncertainty in a calibrated model caused by uncertainty in the estimates of aquifer parameters, stresses and boundary conditions of the model area (Senthil and Elango, 2004).
In this study, following the calibration, the sensitivity analysis was carried out by changing (increasing and decreasing by 25%, 50%, and 75%) the calibrated model parameters such as hydraulic conductivity, recharge and river bed conductance (Fig 7).
Generally, the model is found to be more sensitive for decreasing of the three parameters (hydraulic conductivity, recharge and river bed conductance) than increment. And hydraulic conductivity is found as the most sensitive parameter (Fig 7).

Simulated Water Balance Result
The simulated water budget has been computed for the study area (Table 3). The simulated inflow of the model is 1.2791870E+05 m 3 /day which is nearly equal to simulated outflow 1.2791755E+05 m 3 /day with a difference of 1.1484375E+00 m 3 /day. Recharges from precipitation, river leakage, and groundwater enters the aquifer across the southeastern general head boundary are the major inflow components. River leakage simulated discharge holds 92.8% of the outflow. It also contributed as recharge into the aquifer that accounts for 14.1% of the inflow. Hence in this catchment, it can be generalized that the contribution of the groundwater to the rivers is higher than the contribution of the river to the groundwater (Fig 5).
Groundwater leaving the aquifer, across the southwestern general head boundary and abstraction by pumping for the domestic, industrial and agricultural purpose also contribute some proportions to the outflow from the Dijil river catchment. Groundwater leaving the aquifer across the general head boundary holds 4.3% of the total outflow and that enters to the aquifer accounts to 3.3% of the total inflow (Table 3). Detail water budget result is finding in the table 3.

Scenario Analysis
A calibrated flow model can be used as a tool to evaluate and compare the responses of an aquifer system to potential future stresses. In this study, the model was used to simulate the aquifer-system response to increased withdrawals and decreased recharge.
In all scenarios, other model parameters were kept to the steady-state values except the stress for which the projection was carried out and the resulting changes in water levels and fluxes are interpreted as responses of the system to the changes introduced on it.
In the study area, industrialization and urbanization is booming. Thus, it is logical to think that groundwater withdrawal will increase in the catchment and it is necessary to project its future effects on groundwater table and fluxes of the area so that appropriate water management practices that could mitigate the likely adverse effects of increased withdrawal can be proposed. To evaluate this condition, three scenarios were made: First, simultaneous pumping of the existing wells for 12 hours/day, second scenario is made by increasing the existing withdrawal (3672 m 3 /day) by 30% and the third scenario is made by drilling eight additional new wells with 30 l/s yield from each additional well (Table 4).
As it has been seen in the scenario analysis result, the increasing of well withdrawal has adverse effect on stream leakage and groundwater outflow but relatively less impact on groundwater head. System responses to increased groundwater withdrawal for the above scenarios are shown in table 4 and figure 7. The different scenario results shows that the change in drawdown, river leakage and sub-surface outflow in the third scenario is about four times than the first scenario. Therefore, the managed aquifer recharge (MAR) technologies can provide a variety of water resources management benefits by increasing the volume of stored water and decreasing the drawdown in the aquifer.  Figure 7. Trend of system response to increased groundwater withdrawal rates for river leakage and subsurface outflow.
The second scenario has been done by decreasing recharge to aquifers that may result from the expansion of agriculture, deforestation, urbanization, and decreasing precipitation. If the recharge is decreased by 25% and 50%, the simulation results showed on an average head decrease of 6.1m and 13.4m over the whole area. Moreover, a stream leakage decrease as compared to the simulated steady-state value, and the changes was about 20.3% and 40.3%. This indicates decreasing recharge has an adverse effect on the groundwater table and stream leakage.

CONCLUSION
The calibrated steady-state groundwater flow model of this study was able to reasonably simulate the hydraulic heads that fit the measured heads particularly for the lowland plain part of the catchment. Besides, the model simulated heads contour map shows that the general hydraulic gradient in the basin pursues the surface topography and the gradient is towards the north-south direction and it was in agreement with the groundwater flow system defined in the conceptual model. A minor local irregularity of flow pattern is attributed to the existence of geologic structures.
As the model was intended to study the response of the hydrologic system, three scenarios of increased withdrawals and one scenario of altered recharge were used to study the system response. The effects of the scenarios were evaluated with respect to changes induced on stream leakage, subsurface outflow and groundwater heads compared to the steady-state simulated values.
Accordingly, an increase in well withdrawal of scenario-I (existing wells pump simultaneously), scenario-II (existing drilled wells yield withdrawal increased by 30%) and Scenario-III (additional eight wells having expected yield of 30 l/s drill and pump) resulted in an average decline of the steady state water level by 1.06m, 1.68m and 4.46m, respectively and caused the steady state stream leakage to be reduced by about 2.93%, 4.58% and 11.23%, and subsurface outflow by 9.41%, 14.67% and 37.86%, respectively over the whole catchment. A second scenario was done by decreasing the recharge values shows that a decrease head and a stream leakage in the whole area of study.
An adequate groundwater level monitoring wells should be placed in the catchment to control the total abstraction rates from the aquifer and fluctuations in groundwater levels.

ACKNOWLEDGEMENTS
The first author wish to thank to Amhara Design and Supervision Works Enterprise, Bahirdar branch for financial support. We are also grateful to the two anonymous reviewers, whose suggestions significantly improved the quality of the paper.