Application of GMS-MODFLOW to Investigate Groundwater Development Potential in River Meme Catchment, Kogi State, Nigeria

Groundwater modelling is an important tool that can be used to determine appropriate management strategies for groundwater conditions. It is useful in monitoring groundwater levels due to extensive exploitation as is presently done in the Otokiti area on the River Meme catchment, Lokoja, Nigeria. A conceptual model was developed using Groundwater Modelling System (GMS) software. The model was calibrated both in steady and transient states. The steady state calibration was done for March 2009, aquifer performance data while the transient model was calibrated for March 2016 aquifer data. Results from calibration showed values of hydraulic conductivity varying from 0.02 to 25.6 m/d while the recharge rates varied from 0.0001 to 0.0007 m/d. A predictive run was done from 2016 to 2026 where the model examined the response of the aquifer to abstractions under three different schemes. In Scheme A, the abstractions remained the same as that of the current year 2016. In Schemes B and C, abstractions were increased by 20% and 60% respectively over the 2016 rates. The results showed that there was very little decline in head for locations near the rivers at the eastern and western parts of the catchment. However, for locations in the central part of the study area, which were mostly residential, there was a gradual decrease in head of up to 6 m for scheme A, 6.5 m for scheme B and 6.8 m for scheme C. The implies that even if the rate of abstraction is increase by up to 60% i.e. at 122 m/d, the groundwater system would still be sustainable as a main source of supply for domestic consumption. The volumetric budget of the study area for the three abstraction schemes also showed that the groundwater abstraction was generally less than the annual groundwater recharge. Consequently, more water is available in the formation which can be abstracted for future development without any appreciable loss of head.


Introduction
Water is essential for the survival of mankind. Over the last years, the demand for water for agricultural, industrial and municipal usage has been on the increase. To meet these demands, surface and groundwater sources are required (Rastogi, 2007). In many regions, the groundwater resource is so huge that its occurrence and hydrological significance cannot be overlooked in the planning and management of water resources (Zakir et al., 2011). A major tool in the assessment of water resources management scenarios is the application of computer models, which due to their capability enable decision makers to compare alternative actions and take management decisions to achieve efficient goals without violating specified constraints (Olayinka, 2013). Prediction of subsurface flow, water table level, solute transport, and simulation of natural or human-induced stresses is necessary for groundwater management (Karamouz et al., 2003). Groundwater modelling is an important tool that can be used to determine appropriate management strategies for groundwater conditions particularly in the areas where the hydrological cycles is predicted to be accelerated because of climate change (Mall et al., 2006). An important step in the modelling process is the selection of modelling software in which all possible options are considered (Kumar, 2013). In many cases, researchers have adopted the use of MODFLOW, developed by the US Geological Survey (Harbaugh et al., 2000). Ahmed and Umar (2009) used Visual MODFLOW to simulate the behaviour of the flow system and evaluate the water balance of Yamuna-Krishni interstream. The principle of finite difference method and its applications in groundwater modelling was discussed by Igboekwe and Achi (2011). Their study used finite difference method to solve the equations that govern groundwater flow to obtain flow rates, flow direction and hydraulic heads through an aquifer. Rodr'lguez et al. (2006) developed a groundwater/surface water interaction model for the shallow alluvial aquifer of the Choele Choel Island in Patagonia, Argentina. The model utilized MODFLOW and its stream package, and was successfully calibrated for a historical irrigation season. Their modelling results indicated that drainage through streams is significantly higher than drainage through artificial drains. Inverse modelling of groundwater flow was the approach used by Abdelaziz and Baker (2012) to evaluate the potential of groundwater resources for the development of the Sinai Peninsula. GMS was used for the numerical modelling, and parameterisation, PEST was used in GMS. A steady state calibration of the model was done and the results used as initial values for a transient run. The model assumed the calibrated recharge was constant over the transient simulation from 1988 to 2006. At the end of their study, they came to the conclusion that the water balance is in a critical situation and needs an improved management to exploit the water resources.
Groundwater flow modelling at the source of River Ethiope, Delta State, Nigeria was undertaken by Okocha and Atakpo (2013). The watershed was modelled with a grid of 40 rows by 15 columns with 2 layers. Recharge from rainfall served as input while abstraction from boreholes served as output. A steady state groundwater simulation was carried out and calibrated using six target heads. The model converged after 150 iterations. The results from the modelling showed that abstraction was much more less than recharge. Therefore, further exploitation of the groundwater is possible. Igboekwe and Udoinyang (2011) carried out a study aimed at assessing the extent of the interaction between the Kwa Ibo River and the viable groundwater wells located within the premises of Michael Okpara University of Agriculture, Umudike. The watershed was modelled with a grid of 50 rows by 40 columns with two layers. Computed heads and porosity were used to compute groundwater velocity. The maximum groundwater velocity was 0.44m/d along the Kwa Ibo River. The predominant groundwater flow was found to be towards the south, but it was mostly towards Kwa Ibo River on both sides. Ujile (2013) used the principles of mass transfer to evaluate groundwater contaminants flow model in Yenagoa, Nigeria. The application of the model was carried out with the MATHCAD software systems. The model solution showed that the concentration levels of iron pollutants in the groundwater depends on the distance from source of pollution, time, process of transfer and flow regime. Overexploitation of the aquifers in Lagos Area, Nigeria led Abiola and Agbede (2012) to study and ascertain the extent of sea water intrusion in the area. The study used the two dimensional solute-transport in an incompressible porous medium and transport in the absence of sources and decay reactions. A computer program was developed to solve the 2-D dispersion problems. The model was calibrated in two steps, with the transient runs modelling chloride concentration for three aquifers. The results indicated that the area of chloride concentration increased from 13Km 2 in 1996 to 38Km 2 for the upper coastal plain sands, while the concentration increases from 68mg/L to 83mg/L . Edet et al. (2014), quantified the amount of exploitable aquifer in Akwa Ibom by a regional numerical groundwater flow modelling using MODFLOW. The model was calibrated under steady state conditions to determine the aquifer's hydraulic conductivity and recharge characteristics. Groundwater modelling system (GMS) software, a MODFLOW-based software was also used for this study.

Description of Study Area
The study area lies within latitudes 7º48'23''N and 7º48'9''N and longitudes 6º41'22''E and 6º40'23''E ( Figure 1). The area drains an area of about 1.7 Km 2 and ranges in altitude from 66 m to about 94 m above sea level. The topography of the area is typical of the Mid Niger basin which is characterized by undulating topography with limited flat terrain. This is typical of Mount Patti and Agbaja plateau which are the two main highlands in the area (Omali, 2014). The type of climate in Lokoja is that of the tropical hinterland of Nigeria. The area is characterized by two distinct seasons namely; the rainy season and the dry season. Average annual precipitation from 2001 -2010 for the area ranged between 1000 mm -1500 mm (Akpah, 2014). The study area is dominantly underlain by the Precambrian Basement Complex. However, part of the area is underlain by Cretaceous sediments which uncomfortably overlie the Basement Complex. The basement complex covers about 80-90% of the western and central areas of the Kogi state (Omanayin, 2009). River Meme is one of the tributaries of the Niger which flows through the study area. It is a seasonal river as it is totally dependent on rainfall. During the dry season around February, the river stops flowing with ponds of water at various locations. River Meme gains most of its water from several springs emerging from mount Patti and Agbaja Mountains. Other tributaries of River Meme are minor streams and small rivers around the Crusher village and Otokiti area (Akpah, 2014). Boreholes in Lokoja metropolis in recent years are seen to be a major supplement of the water needs of its residents (Omada et al., 2009). Their study also showed that groundwater is concentrated in the aquiferous conglomerate and sandstone horizons of the Lokoja and Patti formations in the sedimentary terrain, and in the weathered and or fractured basement rocks in the basement area. The hydrogeology is divided into two based on the geological terrain in the state; the basement complex and the sedimentary terrain (Omonayin, 2009).

Figure 1: Map showing study area extent (Google Maps, 2016)
The study area is within the Lokoja metropolis in Kogi state, Nigeria. The area has witnessed an upsurge in commercial and industrial activities and an increase in population since it became the capital of Kogi State in 1991. Hence, the demand for potable water for domestic, industrial and agricultural needs has grown rapidly (Omali, 2014). This has led to substantial increase in groundwater abstractions in the study area, hence it has become necessary to have a good understanding on the use of groundwater resource in the area. Thus the objective of the study is to investigate the effects of these abstractions on the groundwater system and the need for the development of effective management strategies and regulations.

Groundwater Flow Equations
Darcy's law states that "the saturated flow of water through a column of soil is directly proportional to the head difference and inversely proportional to the length of the column" (  It is an established hydraulic principle that groundwater moves from areas of higher potential i.e. recharge areas (higher elevation or higher pressure/hydraulic head) to areas of lower pressure or elevation. This implies that direction of flow of groundwater ideally follows the topography of the land surface. Cracks, inter-connected pore spaces make a rock material permeable. Some permeable materials may allow fluid to move several metres in a day; while some may move a few centimetres in a century. In the real subsurface, groundwater flows in complex 3D patterns. Darcy's law in three dimensions is analogous to that of one dimension. Combination of Darcy's flow equation with continuity equation, which represents the conservation of fluid mass, yields the following partial differential equation describing the three-dimensional movement of groundwater through porous media: Where, Kx, Ky, Kz = hydraulic conductivity along the x, y ,z axes which are assumed to be parallel to the major axes of hydraulic conductivity; According to Rausch (2010) the data required to formulate and solve such models include: Geometry (the shape of model area, elevation of aquifer bottom and top, and thickness of aquifer), aquifer parameters (hydraulic conductivity/transmissivity and specific storage/storage coefficient), inflows/outflows (well discharge/recharge, groundwater recharge by precipitation, boundary flows, and ex-/infiltration from surface water bodies), prescribed heads and observed discharges (e.g. spring discharges etc.).

Data Collection
The following data were collected to facilitate the modelling process.

Model conceptualization
The conceptual model for the study area was developed using the GMS software. GMS was chosen due to its versatility and user-friendly MODFLOW interface. The coverages for the model boundary, source and sinks, recharge rate and hydraulic conductivity were created in the model conceptualization as follows: (a) Boundary The first step in developing a conceptual model is to define the model boundary. The boundary was developed from ArcGIS using the .kml file created from Google Earth and then imported into GMS.
(b) Local source/sink coverage The next step was to define the sources and sinks coverage ( Figure 2). In this coverage, the boundary conditions of the study area were defined as well as the location of the wells.

(a) Boundary conditions
Specific head boundaries were defined along the northern and eastern sides of the model while the southern part of the model was designated as no-flow boundary. These conditions were chosen based on the fact that the northern (river Meme) and eastern sides are bounded by rivers and the southern part by hills.

(b) Creating the wells
A text file containing the Well information was imported into GMS. A well refinement was assigned to each well. This was necessary because a well represents a point of convergence in the groundwater flow and causes steep gradients in the head near the well. Therefore, in other to accurately model the flow near the wells, the grids are refined in the vicinity of the wells.

Figure 2: Wells and Boundaries (c) Recharge zone
The next step in building the conceptual model was to create a coverage for the recharge zone. The recharge was assumed to be uniform in the study area. In GMS, the recharge rate has been set to be calculated from: ℎ = ( ) 0.05 365 (3) The annual rainfall for 2009 was 1763mm (Akpah, 2014). Therefore the recharge rate was calculated to be 0.000242 m/day. This value was used for the steady state run and the initial recharge rate during the steady state calibration.

(d) Hydraulic conductivity zone
The next step was to set up the hydraulic conductivity (K) coverage. For the steady state run, a single polygon zone was defined for K. Initial value was set at 0.25 m/day (Ayenigba, 2016). The specific yield was set within this coverage. The value was set at 0.08 (Akpah, 2014).

(e) MODFLOW grid
With the coverages now 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. In total, the grid was divided into 30 rows and 65 columns.

(f) Layer elevation
The final step in completing the conceptual model is to create the layer elevations. A single-layer model for the study area was chosen due to the fact that the geometry of the fractured zones was unknown. It was clear which aquifer zones were penetrated by the abstracting well. The top elevation was determined based on the DEM data obtained. This data was interpolated into the MODFLOW layer of the model. The top elevation varied from 60 -90m. The bottom elevation was interpolated into the MODFLOW layer based on the depth of the boreholes. Starting head values were also interpolated into the MODFLOW layers. N In many cases, calibration can be achieved much more rapidly with an inverse model. An inverse model is a utility that automates the parameter estimation process. The inverse model systematically adjusts a user-defined set of input parameters until the difference between the computed and observed values is minimized. GMS contains an interface to an inverse model called PEST or automated parameter estimation (Aquaveo, 2014a).

Model calibration
The first step in setting up the inverse model is to "parameterize" the input. This involves identifying the parts of the model input, the inverse model utility is to optimize. There are two approaches for parameterization, zonal and pilot points. The zonal method involves identifying polygonal zones of hydraulic conductivity and recharge, marking the zones as parameters, and assigning a starting value for each zone. PEST will then adjust the K or recharge values assigned to the zones as it attempts to minimize the residual error between computed vs. observed heads and flows (Aquaveo, 2014a). Pilot points can be thought of as a 2D scatter point set.
Instead of creating a zone and having the inverse model estimate one value for the entire zone, the value of the parameter within the zone is interpolated from the pilot points. Then the inverse model estimates the values at the pilot points. Using pilot points, values will vary from cell to cell. When the inverse model runs, the values at the pilot points are adjusted and re-interpolated to the grid cells until the objective function is minimized (Aquaveo, 2014b).
After running MODFLOW, a calibration target appears at the observation wells in the grid. The components of a calibration target are illustrated in Figure 3. The centre of the target corresponds to the observed value. The top of the target corresponds to the observed value plus the interval and the bottom corresponds to the observed value minus the interval. The box represents the error. In MODFLOW, if the box lies entirely within the target, the box is drawn in green. If the box is outside the target but the error is less than 200%, the box is drawn in yellow. If the error is greater than 200%, the box is drawn in red (Aquaveo, 2014c). Predictive model The transient model, having been successfully calibrated, was used to perform predictive simulations so as to access the sustainability of the groundwater system. A predictive model can be used to forecast up to twice the calibration period (Mondal, et al., 2011). It was decided to simulate ten years into the future due to the uncertainty in future hydraulic conditions and stresses to the groundwater system. Hence, predictive simulations were performed from January 2016 to January 2026.

(g) Recharge
The generation of future rainfall data would require extensive statistical data analysis which is not within the scope of this study, therefore it was decided to take future rainfall values by simply replicating the past and present rainfall pattern. This approach was also used by Mondal et al. (2011). The rainfall values replicated were from 2005 to 2015. The recharge was calculated as in Equation 3.

(h) Abstractions
Since the objective of this study is to investigate the effect of future groundwater abstractions, predictive simulations were performed under three different abstraction schemes. a. Scheme A: Abstraction remained the same as current year 2016. b. Scheme B Abstraction was increased by 20% over 2016 level c. Scheme C: Abstraction was increased by 60% over 2016 level

Model calibration
The results from the steady and transient state calibrations are shown in Figures 4 and 5. The result of the sensitivity analysis showed that the model was sensitive to hydraulic conductivity, recharge and specific yield. The value of the specific storage had no effect on the model. Sensitivity analysis showed that a slight change in any of these parameters, especially the values of hydraulic conductivity and recharge, would affect the results of the model calibration. Thus, the model was more sensitive to hydraulic conductivity and recharge.

Groundwater level
The groundwater head distribution for the three abstraction schemes are shown in   Figure 13(a) shows the water table for 2016 while Figure 13(b) shows the water table for 2026. From this, it was observed that a difference in head of about 10m occurred during that duration. Figure 14(a) shows the water table for 2016 while Figure 14(b) shows the water table for 2026.
From this, it was observed that a difference in head of about 10.5m occurred during that duration.
The time series plot for scheme C is shown in Figure  15.

Volumetric budget
The volumetric budget for the various schemes are shown in Tables 1 to 3.   The outflow of groundwater into the rivers was 311.66 m 3 /day which is comparatively higher than the inflow (48.37 m 3 /day) of groundwater from the rivers. For scheme B, groundwater abstraction was 91.8 m 3 /day which is about 31.6% of annual groundwater recharge (290.66 m 3 /day). The outflow of groundwater into the rivers was 300.34 m 3 /day which is comparatively higher than the inflow (51.77 m 3 /day) of groundwater from the rivers. For scheme C, groundwater abstraction was 122.40 m 3 /day which is about 42.11% of annual groundwater recharge (290.66 m 3 /day). This implies that the recharge through rainfall was much more than the groundwater being abstracted. The outflow of groundwater into the rivers was 277.90 m 3 /day which is comparatively higher than the inflow (58.69 m 3 /day) of groundwater from the rivers. Generally, for all the schemes, the recharge through rainfall was more than the groundwater being abstracted. The higher values of outflow from the rivers in comparison to the inflow would suggest that the rivers are gaining more groundwater than they are adding to the catchment.

6. Conclusion
Groundwater flow model for Otokiti area of the River Meme catchment, Lokoja has been developed. This was achieved with the aid of the Groundwater modelling system (GMS) software and other tools such as GIS and Google earth. The model was calibrated in both steady and transient states with the available data. Results from the calibration in both states showed that the model calibration was successful as indicated by the calibration targets and the plot of computed vs observed heads.In order to simulate the model, input and output stresses were applied. The heads distribution would suggest that the groundwater flow in the area is towards the northeast direction from south to north. This is known due to the fact that higher heads occurred around the southwest while the lower heads were around the north-eastern parts of the model.
A predictive run was done from January 2016 to January 2026 where the model was examined for its response to abstractions under different conditions. The results from these scenarios showed that if the rate of abstraction increases by 60% over the 2016 rates, the groundwater system would still be sustainable as a source of water supply to the community. This is further buttressed by the time series plots of flow heads at various locations within the catchment which showed little or no further decline beyond 2026. Hence, there is room to provide additional wells in the area for water supply needs of the people.