Challenges of Groundwater Flow Model Calibration Using MODFLOW in Ethiopia: With Particular Emphasis to the Upper Awash River Basin ()
1. Introduction
Water-centered development is explicitly seen as the entry point for growth and improved livelihoods in Ethiopia. Especially in the last few decades, together with surface water, the large potential role and contribution of groundwater are recognized [1] . Hence, focus on sustainable groundwater use is larger than assumed previously and the knowledge is evolving. Along with the increasing groundwater development, there is growing awareness that management is needed to ensure the sustainability of the resource. In spite of the awareness, institutions in managing groundwater are still weak in the country. Thus, recording and monitoring of basic information on groundwater levels or total number of functioning wells for instance, has not been done to the required level. Furthermore, the complex hydrogeological setup of Ethiopia demands detailed hydrogeological investigation along with proper monitoring plans. As a result, designing tools, for example groundwater models, which would be used as basis for groundwater management decision support system, has become very difficult task. The problem is clearly in vicious circle: there is a need for more intensive use of groundwater, resulted from ambitious water centered development plans. This obviously requires a commensurate development of groundwater management institutions. However, groundwater management institutions are weak leading to lack of basic information important to design well fields and manage the resources properly.
Even with such insufficiency of data, however, there has been an increase in the usage of groundwater flow models, mainly MODFLOW numerical code to address a wide range of water-related problems in Ethiopia [2] . Due to lack of monitoring well networks, or representative monitoring of the groundwater systems, all the groundwater models are, however, simulated and calibrated under steady state condition, and there are only a couple of attempts of transient groundwater flow simulations [3] [4] . It is impossible to calibrate a transient model without time series data, and hence choice of steady state groundwater model development over transient simulation seems unquestionable. Otherwise, transient models are more important tool than models simulated under steady state conditions for groundwater management that leads sustainable utilization of the resource.
The steady state calibration of the groundwater flow modeling practice in the country have been done with the common approach that, because there are no monitoring wells, the equilibrium state groundwater models can be calibrated with static water level (SWL) records of water supply boreholes measured at the time of drilling, before pumping test. There are some assumptions involved in this calibration practice, which are commonly not thoroughly justified. These assumptions are elaborated in this work.
Attempt was made to introduce the major challenges of numerical groundwater flow model calibration practices in Ethiopia. Particular emphasis is given to looking at ways of calibrating groundwater flow models under steady state without the need for having monitoring wells in selected basins. The assumptions made in the modeling procedures to use SWL collected over many years from water supply boreholes, to calibrate steady state models is elaborated. Detailed discussion on the calibration procedure and the implications of the assumptions have been made by considering SWL records of bore holes drilled in the Upper Awash river basin as a working example. This paper also indicates how the model calibration approach should be geared in terms of logical assumptions and reasonable representation of groundwater systems to better contribute to sustainable use of groundwater resources in Ethiopia.
Though may not be comprehensive, the authors think that this paper provides general assessment of researches on calibration of models done with MODFLOW which is Modular Three Dimensional Finite Difference Groundwater Flow Model [5] . A detailed account of limitations and guidelines of groundwater model calibration has been presented by [6] . This guide line has been taken as a major reference in evaluating model calibration procedures in Ethiopia.
2. Overview of Numerical Groundwater Flow Modeling
In the most general terms, a model is a simplified representation of the real world. Groundwater flow models are simplified representation of groundwater system with mathematical equations solved by a computer program [6] [7] . And yet, the simplifications should be based on logical assumptions and without compromising reasonable representation of complex hydrgeological dynamics.
For instance, all natural systems fluctuate in response to climatic variations that can be diurnal, seasonal, annual, decadal or longer. As a result, true steady- state conditions do not exist in natural systems, especially when the time taken in modeling is longer. Hence, assumption is made so that a system can be represented by a state of dynamic equilibrium or an approximate equilibrium condition. This can be a reasonable assumption, only if the system approximate equilibrium during the period of interest and objectives of the investigation do not require information on the time it takes for a system to respond to new stresses. If the system is not at a period of equilibrium or if there is a need to know how the system respond to new stresses during the periods of interest, then a transient simulation is required. The selection of steady state simulation is, therefore, never a random selection [6] . The appropriateness should be justified with adequate information as to whether the system is reasonably close to steady state condition or not, within the selected simulation period.
Equation (1) describes the three dimensional groundwater flow process in porous medium [8] . The equation can be solved for equilibrium (steady state) or changing (transient) conditions. For steady state simulation, the simulated head is assumed to be constant within the model time domain, and so the storage term is zero. For transient simulation, there is a change in simulated head resulted from time dependent stresses on the system for each stress period.
(1)
where Kx, Ky, and Kz are values of hydraulic conductivity in the x, y and z directions along cartesian coordinate axes, which are assumed to align with principal directions of hydraulic conductivity (LT-1), h is hydraulic head (L), W is a volumetric flux per unit volume and represents sinks and/or sources (T-1), Ss is the specific storage of the porous material (L-1), and t is time (T).
Because any hypothetical system can be simulated with a combination of hydraulic parameters, artificial boundary conditions and simplified conceptualizations, simulation is not the final phase of the modeling procedure. It is only through calibration that a real system can be logically depicted. Therefore, models simulated either for steady state or transient conditions should be calibrated with observed heads and/or flows, and it is one of the most important phase of numerical groundwater flow modeling procedure.
Generally, without calibrating the model with real time head and/or flow observations, it is impossible to investigate real groundwater dynamics, or evaluate recharge, discharge and existing aquifer storage processes, or use the model as management decision support tool. Model calibration in its most limited meaning is, the modification of model input parameters for the purpose of making the model more closely match observed heads and/or flows [9] [10] [11] . Adjustment of parameters can be done manually or automatically by using nonlinear regression statistical techniques [8] .
But, this is by no means entails only matching the observation and simulated head is enough [11] [12] . Multiple calibration of the same system can yield different combinations of boundary conditions and aquifer properties due to non- unique characteristics of the calibrated models [7] [13] . Non-uniqueness of calibrated parameters is partly caused by complexity of real system coupled with lack of sufficient data [14] . Although a good match in itself does not prove the validity or adequacy of the model [15] , a poor match provides evidence of errors in the model.
In urban and peri-urban environment where there is continuous change in land use and associated recharge processes and with widespread groundwater utilization practices, steady state models may not represent the dynamics of the groundwater. Therefore, in such an environment transient models represent the real world much better than the steady state models. The focus of this work is to address this problem by considering the Upper Awash river basin as working example. The basin is characterized by many cities including Addis Ababa where groundwater is being used for municipal, industrial and locally agricultural purposes. Hence groundwater dynamics is high that demands development of a more robust transient groundwater model or careful calibration of steady state conditions.
3. Numerical Groundwater Flow Modeling Practice in Ethiopia
3.1. Brief Review
Numerical groundwater flow models have been used in developed countries since 1970’s. However, its importance and application was not well understood in Ethiopia until the 1990’s [2] . Afterwards, there has been an increase in the usage of groundwater models, especially MODFLOW to address a wide range of water-related problems in Ethiopia. The great majority of works have been made for academic exercises rather than using the model as decision support tool.
A more or less similar researches, in terms of aquifer characterization, setting appropriate pumping rate, describing surface and groundwater interactions and demonstrating potential consequences of a groundwater system developed in unregulated manner and under natural stress conditions include: [3] for the Adelle-Haromaya dry lakes catchment; [20] for Mihone sub basin in Raya valley; [21] for the Dire Dawa area; [22] for Lake Hawassa basin [23] and [24] for Upper and middle Awash river basin respectively; [25] and [26] of the Lake Tana basin, [27] and [28] of the Kobo Valley; [29] for Meki River catchment; [30] for Raya Valley; [31] for the Upper Awash and south eastern flank of Blue Nile basin. These modeling exercises helped in understanding the groundwater flow system of the different catchments and enlightened the importance of groundwater models in understanding groundwater flow under different hydrogeological environment and in studying groundwater and surface water interactions.
3.2. Commonly Used Calibration Approach
Proper calibration of numerical groundwater flow models can be made, both under steady state and transient condition if appropriate well distributed data set is available. For steady state calibration, since there is no change in head with time, a single simulated heads or flows is calibrated against one time groundwater level or flow measurements. But, for the transient model calibration, time series data is needed. Because the head simulation will vary as the stresses change, head or flow observations for each stress period to calibrate models should be available.
Often, groundwater level measurements are made using deep meters in a onetime field campaign. Recently there is attempt to monitor few wells using data loggers in Upper Awash and Raya valley (northern Ethiopia). However, the data generated has not been used for developing robust transient models that can be readily used for groundwater management. Reference [3] and [4] who attempted transient simulation to see seasonal variations of groundwater dynamics (Table 1). Development of detailed transient predictive models, which can readily be used as a decision support tool for development and implementation of sustainable water resources utilization practices is more or less nonexistent. Table 1
Table 1. Summary groundwater modeling done and water supply wells used.
WSW = water supply wells, MW = monitoring wells, SS = steady state, T = transient.
summarizes the number of observations, type of wells used to collect groundwater level, model simulation type and data collection period for the different groundwater flow models calibration practices in Ethiopia.
The steady state calibration of the groundwater modeling practice in the country have been done with the commonly used approach that the equilibrium state groundwater models can be calibrated with SWL records from water supply wells collected at the time of drilling. The exception is steady state numerical groundwater flow model for the case of Akaki well field [2] , in which case four monitoring wells were used in addition.
The number of groundwater head observations from water supply wells used to calibrate respective models in all the materials reviewed in this paper range from (32 to 530) (Table 1). Though not mentioned in most works, looking at the number of measurements, it can be anticipated that the SWL records used for the steady state calibration have been collected over many years. However, for how many years the SWL data have been collected, seems to be not a problem. At least, not mentioned as a major problem to calibrate a model for equilibrium condition. Except [16] , well head measurements taken within the model temporal domain from few open wells is used to calibrate the model for steady state condition; [4] [20] [29] , years of data records were specified; what seems taken enough is “calibration” of the model with SWL data of water supply wells, regardless of when the head observation data have been collected. Under extreme data scarce conditions, high discharge spring emanation points have been used as static groundwater levels to fill data gaps [16] . Such an approach requires detailed understanding of the nature and type of springs that is, whether such springs are connected with regional groundwater levels or not.
3.3. Considerations in Model Calibration
The “conventional” equilibrium condition model calibration practice in Ethiopia has been done based on the following assumptions: 1) water supply wells can be used as monitoring wells; 2) all SWL data have been collected within temporal model domain; 3) the groundwater level remain constant within the model simulation period. Though not well explained with reasonable justification, [17] also assumed that error introduced into the result due to heads measured at different times is lower than errors introduced due to uncertainties in recharge, hydraulic conductivity or other model input parameter estimations.
For the first assumption, since the data were collected at the time of drilling, before pumping started, the systems were not disturbed. Hence, to use the water level records from water supply boreholes for calibration seems to be reasonable. But, the second and third assumptions imply that, if a surface map is constructed using the groundwater level readings regardless of years of groundwater level records, it would give water level surface. Nonetheless, water table or potentiometric surface should be made using groundwater level measurements taken within short period of time [32] . Otherwise, the groundwater level can change with time and would be illogical to just assume that the piezometric surface remains constant for years.
4. Illustrative Example of Challenges of Model Calibration: The Case of Upper Awash River Basin
The rationale of selecting Upper Awash river basin as a good working example for groundwater flow model calibration lies in the availability of relatively better well distributed groundwater level records (although not continuous) as compared to many groundwater basins of Ethiopia. Furthermore, most of the above mentioned modeling practices have been conducted in the same basin and its major tributary catchments.
4.1. Upper Awash River Basin
The Upper Awash river basin is located in central Ethiopia at the western margin of the Main Ethiopian Rift (MER). A large portion of the basin is within the central Main Ethiopian Rift floor. In this study the most urbanized upper part of the basin is considered. The capital Addis Ababa is located at the northern end of the study area. The area considered in the analysis is confined within the limits of 8˚23'25.15" to 9˚18'42.15"N latitude and 37˚59'9" to 39˚04'12.69"E longitude (Figure 1). The Upper Awash basin is adjacent to the Blue Nile river basin to the north. The total area considered is 6735 km2.
Any deterministic groundwater flow modeling requires good conceptual hydrogeological model. Reasonably good works have been done in terms of
Figure 1. Location map of the Upper Awash river basin.
conceptualization of the lateral and vertical extent of the aquifers and the mechanism of groundwater flow [33] [34] [35] [36] . There are also numerical groundwater flow models done either for the entire Upper Awash or limited sub catchments of the Upper Awash river basin using MODFLOW [2] [18] [23] [36] .
Until very recently, groundwater level measurements have been done in a field campaign using deeper. But, very recently few groundwater level data loggers were installed in deep wells around Addis Ababa. However, no data from these monitoring stations have been availed. But, there are about 600production wells with SWL data collected at the time of drilling by different authors (Figure 2). It is very important to mention that, out of these production wells, only 273 SWL data are with known drilling time record (Figure 2). Based on the record, the SWL data from the 273 wells were collected from 1937 to 2007. However, since most of the wells were drilled since 1975, 220 wells are considered for demonstration in this study.
4.2. Previous Groundwater Flow Model Calibration Practices
Figure 3 show the piezometeric map constructed based water level records from more than 500boreholes. Taking only the 220 production wells with known time of drilling into account, the groundwater level data were collected over 32 years period. If such “water level” data is used to calibrate a model under steady state condition, it basically means, no matter what happened over 32 years in the natural and anthropogenic stresses, the groundwater level remains constant and the groundwater system approximate equilibrium state. In reality, with time the water level changes.
Aside from anthropogenic influences, climate influences the piezometric surface. Therefore, even static water level may change depending up on seasonal
Figure 2. Distribution of production wells in Upper Awash river basin.
Figure 3. Piezometeric surface constructed using SWL records collected over 32 years.
and internal variability. This can be easily understood from the average rainfall of the Addis Ababa station over 45 years (Figure 4). Due to changes in rainfall over the basin, the static water level changes from year to year. Therefore, even static water level of a given well should be averaged out over many years. In fact this is practically challenging unless records are taken systematically over a number of years. Therefore, errors could be introduced when one time static water level measurements are taken as representative piezometric surface of a certain locality.
Rainfall variability affects directly groundwater recharge and hence piezometric surfaces. To illustrate this three years climatic data have been considered to establish the distributed groundwater recharge in the same basin (Figure 5(a)- Figure 5(c)). In establishing the recharge, the most important time variant parameter is recharge followed by evapotranspiration.
The recharge distributions are established by interpolation of HYDRUS 1D, point recharge simulations [37] [38] , using ordinary kriging tool using ArcGIS software. The spatial variability over the three years (2004, 2005, and 2006) indicates that there has been significant spatial variation of recharge in the basin. This brings about significant variations in piezometric surfaces. In addition to the spatial distribution of recharge, it was tried to see the point recharge variability (Figure 6) for the years 2004, 2005, and 2006 at 18 simulation sites in Upper Awash.
Figure 4. Average annual variability at Addis Ababa meteorological station.
Figure 5. Spatial distributions of directrecharge ((a) HYDRUS simulation points; (b) 2004; (c) 2005 and (d) 2006).
Over the years, many boreholes have been drilled in the basin. These boreholes were used for municipal, industrial and to a limited extent irrigation purposes. This change certainly brings about significant variability in static water levels. This means at a particular point a borehole drilled 30 years ago will have different static water level when it is drilled currently. The dramatic increase of production boreholes with time is shown in Figure 7. Especially in the last three decades the increase in productive boreholes is very high.
Besides, the variation of the 220 water levels with time has been plotted with topography in order to see how significant time was in casing groundwater dynamics (Figure 8). Topographic variation is plotted together with groundwater level measurements taken for years to take into account the possible variation of groundwater level due to variation in space. It can be clearly observed that, even with random selection of wells at any location in Upper Awash basin, the variation pattern of groundwater level due time elapse is not in line with variation of water level due to topographic variation. It is evident that the variation in groundwater level is not just because of location, but also time, especially in the last couple of decades when significant number of boreholes introduced in the system.
Figure 6. Three years direct recharge variability at 18 selected sites in Upper Awash basin.
Figure 7. Increase in water supply boreholeswith time in Upper Awash basin.
Figure 8. Groundwater level variation with time vs topographic variation.
The above mentioned factors certainly demonstrate challenges of groundwater model calibrations under variable topographic conditions, growing water uses and one time groundwater level observations. This demands careful observations in developing deterministic groundwater flow models that can be used as decision support tools.
4.3. The Way Model Calibration Should Be Made
In the case of the Upper Awash river basin, there are well distributed groundwater level records from water supply boreholes. The problem lies in how and when water level observations have been collected over many years. Though there are many wells in the basin, it is only 220 wells with known date of drilling. Based on the boreholes with known date of drilling, the 220 wells have been drilled in the past 32 years. The groundwater level records are onetime SWL measurements collected during drilling of the boreholes before pumping tests were conducted. Monitoring wells to continuously measure basic information is almost nonexistent. Hence, using either all the wells or the 220 wells water level readings, which are records collected over more than 30 years, to construct “average water table” or “potentiometric surface” is contradicting with the meaning of water level surface [32] . Therefore, all the data cannot be used to calibrate the model and end up with wrong interpretation of the dynamic system.
Even constrained under such situation, there are some alternatives much better than the way numerical groundwater flow models calibration have been commonly done in Ethiopia. One of the reasonable options could be to calibrate the model for steady state simulation. From the 220 wells with groundwater level data with known date of drilling, it is possible to extract fair distribution of water level data collected from 2004 to 2006 (Figure 9). Instead of using all the available data collected for more than 32 years in the area of interest for calibration
Figure 9. Potentiometric surface map constructed using three years groundwater level records (2004-2006).
(Figure 3), just because their spatial distribution is good, it is better to rely only on the average water table of head observations of the three years record to calibrate the model.
The three years water level readings are relatively short period records, within which average water table can be assumed [32] , which is incomparably better than the assumption “head remains constant over 32 years”.
The other plausible calibration alternative is, instead of taking average water table of the three years records, to use only water level data of boreholes constructed within the model temporal domain to calibrate steady state models. For instance, in the earlier case, steady state model can be simulated using recharge distribution and pumping for 2004 and calibrate the model with 2004 water level readings. Similarly, simulation for steady state condition for 2005 or 2006 can be done and calibrate the models with their respective head observations. The limitation with this approach is, it may not be possible to map the water table due to lack of fair distributions of boreholes. Therefore, visual matching of groundwater level surfaces is difficult. But, calibration is still possible with comparison of head simulations and observations from the available water level records.
If the need is to see how the system responds to different stress conditions, pseudo transient simulation is also possible, provided reasonable distribution of head data for consecutive years is available. For the case of Upper Awash river basin, the three years (2004, 2005 and 2006) water level can be used (Figure 9) with recharge distributions and the pumping rates as stresses for the same period and simulate transient model with three stress periods. Calibration of such pseudo transient model can be done using water level observations of the same years. Even though the groundwater level are not time series observations from the same boreholes, the different year records of the different boreholes can still be used to calibrate the pseudo transient model.
In fact in its real sense, this may not represent the conventional transient model and the approach is methodologically different from the standard transient model calibration procedure. However, essentially it is more or less the same except this method only considers time not space. In the conventional procedure, time series observations should be recorded from monitoring wells, but in this particular approach, water level readings from boreholes drilled in different times are utilized to calibrate transient model simulations. The preconditions to this approach include: 1) Stresses should be measured for those years, observations to calibrate the model are available. 2) The groundwater level observations records used for calibration should be of consecutive years. In that way the variation of the groundwater system in time can be fairly controlled.
5. Conclusion and Recommendation
5.1. Conclusions
Modeling is a step-by-step procedure through which logical assumptions and simplifications are made, which would possibly lead to realistic representation of a system. There are assumptions in conceptual model development phase, not to overcomplicate the environment which is difficult to comprehend. There is regionalization during parameterization, in order to avoid excessive heterogeneity. Artificial boundary conditions are created with degree of uncertainty, to give the model a certain level of degree of freedom. These approaches are simply inherent to any modeling practice which can never be avoided. However, additional simplifications on top of all these uncertainties, through improper model calibration procedure lead to unrealistic simulation of the reality. This issue is well illustrated by taking the time variant fluxes and water levels in the Upper Awash basin.
A groundwater table mapped using groundwater level records collected over significant time periods, cannot realistically produce a water table or potentiometric surface. Unless it is well justified or tested, to just assume the head remains constant for many years is too much of an assumption. Hence, it leads to oversimplification of the groundwater model, which in turn results in ground water flow models, which are more of hypothetical than depicting the reality.
The assumptions of steady state condition can be reasonable, at least if the observation head values were collected within short period of time, within which it is logical to assume that there is no significant change in all the model stresses and therefore, the model assumes equilibrium state. If the need is to see how a groundwater system responds to time dependent stresses, SWL records of sequential years can be used to calibrate transient conditions.
Multiple calibration of the same system can yield different combinations of boundary conditions and aquifer properties due to non-unique characteristics of the calibrated models. The assumptions made at the beginning regarding appropriateness of the models, in terms of conceptualization, parameterization and definitions of boundary conditions should always be seen seriously.
5.2. Recommendations
Generally, both scarcity and the nature of the available data put pressure in the model simulation and calibration process. In most of the cases, there is shortage of data, and it is obviously a bottleneck for model simulation and calibration of any system. But still, it shouldn’t lead us to illogical simplifications and unreasonable representations.
The over simplifications in the model calibration procedures in Ethiopia have been justified by data scarcity or unavailability of appropriate data. The selection of steady state groundwater model simulation is also a forced choice due to lack of monitoring wells. In order to meet the ambitious water centered development plan and use badly needed groundwater resources sustainably, groundwater management institutions should be strengthened so as to establish monitoring wells fitted with modern data loggers.
Mechanisms should be put in place to build monitoring well networks, at least for well fields which are being intensively exploited. It is only with time series data from monitoring wells, researchers will be able to simulate models under transient conditions, which are by far more powerful groundwater management tools than steady state models and better contribute to sustainable use of groundwater resources of Ethiopia. Otherwise, regardless of the extent of efforts exerted, model will never play the expected key role, as it should, in solving practical groundwater management problems, but rather lay foundation for non-ending future studies.
Competing Interests
The authors declare that they have no competing interests
Authors’ Contributions
The first author, Behailu Berehanu is currently conducting his PhD research on integrated water supply systems modeling using WEAP and MODFLOW models, a case of Upper Awash river basin. This work is a major component of the research in fulfilment of one of the major objectives. Hence, he has made substantial contribution in every aspect of the work. He is conducting the research under very close supervision and technical support by Tenalem Ayenew, who is professor of Hydrogeology at Addis Ababa University. Professor Tenalem involved from conception and design to analysis and interpretation of data. Tilahun Azagegn, assistant professor working for the same university involved in drafting the manuscript or revising it critically for important intellectual content.
Acknowledgements
The authors thank Addis Ababa Science and Technology University, the Swedish International Development Cooperation Agency (Sida), Addis Ababa Water and Sewerage Authority (AAWSA), Water Works Design and Supervision Enterprise (WWDSE)and Addis Ababa University, School of Earth Sciences for their support and provisions.