Modeling the spatio-temporal flow dynamics of groundwater-surface water interactions of the Lake Tana Basin, Upper Blue Nile, Ethiopia

Water availability and water quality are major challenges of the 21st century due to a growing water demand. The challenge is twofold for developing countries like Ethiopia because of a limitation in technology and financial resources to invest in water resources development. Thus, groundwater is the only reliable water supply option for most rural areas in Africa to meet the dispersed demand (Hiscock 2011). However, extraction of water from the aquifer system is rarely regulated, resulting in over-pumping of the groundwater. The complete drying-up of Haramaya Lake in Eastern Ethiopia since 2005 is one example of the consequences of decreasing groundwater levels due to over-pumping for agriculture and household use (Abebe et al. 2014). Therefore, a comprehensive understanding of catchment hydrology is essential for sustainable water management in catchments with a strong groundwater (GW) and surface water (SW) interaction (Bailey et al. 2016; Dowlatabadi & Zomorodian 2016; Guevara Ochoa et al. 2020).

The GW-SW interactions have significant impacts on the quantity and quality of water in both, aquifers and surface water bodies (Sophocleous 2002; Lin et al. 2018). According to Winter et al. (1998), understanding the interaction between GW and SW is imperative to address conjunctive use of GW and SW, mitigate contamination of aquifer and surface water bodies, manage water rights and reservoirs, and integrate GW flow into watershed management and planning. As a result, the number of studies focusing on GW-SW interactions has increased substantially over the last two decades (e.g. Guevara Ochoa et al. (2020); Lin et al. (2018); Becher Quinodoz et al. (2017); Wondzell 2015; Elsawwaf et al. (2014); Ayenew et al. (2008); Winter et al. (1998)).

According to Winter et al. (1998) and Guevara Ochoa et al. (2020) climatic, geological, topographic, and biotic factors together with land use type and management are the most important variables that affect GW-SW interaction. Thus, process-based physical models that are capable of incorporating such variables are needed to understand GW-SW interaction.

Since the 1980s, the physically based, 3D, distributed numerical groundwater model MODFLOW (McDonald & Harbaugh 1988; Harbaugh 2005) is the most commonly applied model to simulate GW-SW interaction (Brunner et al. 2010). The Soil and Water Assessment Tool (SWAT) (Arnold et al. 1998; Arnold & Fohrer 2005; Neitsch et al. 2011) on the other hand is one of the most widely applied hydrological models globally. However, these two models have some limitations: (1) MODFLOW is limited with regard to the simulation of surface processes such as infiltration and surface runoff, nutrient cycling and transport, plant growth and impacts of management practices on agricultural systems (Bailey et al. 2016); and SWAT has commonly known weaknesses with regard to (1) non-spatial reference of hydrologic response units, (1) difficulty of calculating distributed groundwater parameters, and (2) its simplified groundwater subroutine (Arnold & Fohrer 2005; Nguyen & Dietrich 2018). This underlines that most of the physical standalone models are designed to simulate part of the hydrological processes (either the land surface or subsurface hydrological processes). For this reason, the use of coupled land surface and subsurface hydrologic models has increased (Guevara Ochoa et al. 2020).

A number of individual attempts to couple MODFLOW (McDonald & Harbaugh 1988; Harbaugh 2005; Niswonger et al. 2011) with SWAT (Arnold et al. 1998; Arnold & Fohrer 2005; Neitsch et al. 2011) have been published (Sophocleous & Perkins (2000); Kim et al. (2008); Guzman et al. (2015); Bailey et al. (2016); Ehtiat et al. (2018)). Bailey et al. (2017) and Park et al. (2019) have also created graphical user interfaces to facilitate the use of the SWAT-MODFLOW model (Bailey et al. 2016).

More recently, the SWAT-MODFLOW became an important modeling tool for many researchers worldwide to (1) investigate regional and catchment scale GW-SW interaction (Bailey et al. 2017; Semiromi & Koch 2019; Guevara Ochoa et al. 2020); (2) to study water management strategies, climate change impact, and abstraction scenarios (Aliyari et al. 2019; Chunn et al. 2019; Molina-Navarro et al. 2019); (3) to determine hydrologic system responses such as streamflow, groundwater level, and groundwater discharge to streams (Izady et al. 2015; Abbas et al. 2018; Wei & Bailey 2019) and (4) to simulate the spatio-temporal variability of water resources (Gao et al. 2019; Guevara Ochoa et al. 2020).

In Ethiopia GW is widely used by many sectors. Due to its proximity to the point of demand, GW provides 85% of the domestic water supply, 95% of the industrial use, and a small fraction of irrigation water (Khadim et al. 2020). Surface water is the major source of water supply for small scale irrigation (Awulachew et al. 2007). These national statistics are similar for different basins in the country. The Lake Tana Basin, where the largest fresh water lake of Ethiopia is located, has a similar structure of water use for different purposes as the national statistic. GW is the primary source of water supply for towns, industries and rural communities in the Lake Tana basin, and there is an attempt by the federal and regional governments to develop groundwater using deep boreholes and surface water using dams for medium- and large-scale irrigation in the basin (Mamo et al. 2016). However, the knowledge about groundwater recharge/discharge, GW-SW interaction, and conjunctive use of GW and SW is scanty due to the lack of relevant data.

The majority of hydrological modeling studies in the Lake Tana basin applied SWAT, and the studies focused on investigating the land surface hydrological processes. For example, Setegn et al. (2008) proved the suitability of the SWAT model to predict streamflow and the water balance of the Lake Tana Basin. The study showed that streamflow in the Lake Tana basin is dominated by groundwater discharge. In a more recent contribution, Tigabu et al. (2019) also confirmed that groundwater discharge dominates the total streamflow of the major tributaries of the Lake Tana basin (Gilgelabay, Gumara and Ribb catchments, 73, 62, and 60%, respectively). Setegn et al. (2011) and Dile et al. (2013) applied the SWAT model to investigate effects of climate change on the hydrology of the Lake Tana basin. According to Setegn et al. (2011), groundwater contribution to streamflow is expected to decrease for two future time periods (2046–2065 and 2080–2100), and Dile et al. (2013) reported that the mean monthly flow volume in the Gilgelabay catchment is expected to increase up to 135% during the 2080s. Another study by Dile et al. (2016) explored the effect of intensive water harvesting in the upstream part of the Lake Tana basin on downstream ecosystem services. This study revealed that the implementation of water harvesting ponds for irrigation increased the annual groundwater recharge by more than 2 mm. A study carried out by Woldesenbet et al. (2017) showed that the average annual groundwater discharge and percolation decreased gradually due to land use and land cover changes between 1986 and 2010. Recently, Teklay et al. (2019) investigated the impact of static and dynamic land use data on hydrological responses such as surface runoff, evapotranspiration, and peak flow in the Gummara catchment.

There are a couple of studies that focus on groundwater flow and GW-SW interaction in the Lake Tana Basin. Chebud & Melesse (2009) applied MODFLOW (McDonald & Harbaugh 1988) to investigate the temporal distribution of percolation in the Gumara catchment in particular, and in Lake Tana Basin in general. This study confirmed that groundwater contributes to the Lake Tana as a baseflow. Mamo et al. (2016) investigated the GW-SW interaction of the Lake Tana basin using baseflow separation, soil water balance, and chloride mass balance methods. The study concluded that a large part of the groundwater recharge (68.7%) discharged as base flow to the streams. Overall, the past SWAT and MODFLOW modelling studies in the Lake Tana basin confirmed that streamflow is groundwater dominated. Thus, understanding the interaction between the GW and SW is required to address the use of both resources.

Dile et al. (2018) reviewed the state of the art of hydrological research in the Upper Blue Nile basin and concluded that previous studies have been focused on single model applications to estimate a single output such as streamflow or sediment loss at the outlet of a catchment. Therefore, the application of coupled models to predict multiple outcomes across multiple spatial and temporal scales in the upper Blue Nile basin is recommended. So far, no studies have been published on the application of coupled hydrological models to investigate GW-SW interaction in the Lake Tana Basin, Ethiopia. Therefore, this study aims at investigating the flow dynamics of GW-SW interaction in three groundwater dominated catchments (Gilgelabay, Gumara, and Ribb) of the Lake Tana Basin. To this end, a model that is capable of simulating surface and groundwater hydrological processes as well as GW-SW interaction is required.

Thus, a coupled SWAT-MODFLOW model (Bailey et al. 2016, 2017) was applied to enhance the understanding of spatio-temporal interactions of GW and SW in the Lake Tana Basin, Upper Blue Nile, Ethiopia. The specific objectives of this study are:

• to investigate if the aquifer and stream network systems are hydraulically connected,

• to assess GW-SW exchange rates at river reach and sub-basin levels, and to identify the spatial extent of GW discharge and recharge areas,

• to determine if there are differences in GW-SW interactions between the three catchments,

• and to understand the spatial variability of the GW head.

The Lake Tana basin, which is the second largest sub-basin of the Blue Nile, is situated in the north-western part of Ethiopia (Figure 1). Its catchment area at the outlet of Lake Tana is 15,321 km². The lake has a surface area of 3,156 km² that covers about 20% of the total catchment area (Alemayehu et al. 2010; Tigabu et al. 2020). The catchment has a very diverse topography with an altitude ranging from 1,322 m to 4,111 m above sea level (Tigabu et al. 2020). Moreover, the lake is shallow with a maximum depth of 15 m and characterized by a steep slope at the borders and by a flat bottom (Kebede et al. 2006). It comprises about 50% of the country's fresh water (Costa et al. 2014). The Blue Nile River originates from Lake Tana. More than 40 rivers and streams flow into Lake Tana with a mean annual inflow of 158 m³/s (Alemayehu et al. 2010), but 86% of the water originates from three major rivers: Gilgelabay, Gumara, and Ribb (Alemayehu et al. 2010). The annual volume of surface outflow measured at the lake outlet is 4 billion m³ (Figure 1).

Similar to its hydrologic variability, the Lake Tana basin has a heterogeneous hydrogeology. According to Kebede et al. (2005), areas surrounding Lake Tana are covered by quaternary basalts and alluvial sediments. These geologic formations are characterized by high groundwater transmissivity that vary between 100 and 700 m²/day. The highland part of the basin is covered by basalts, where the groundwater transmissivity is lower compared to the quaternary basalt and sediments (Kebede et al. 2005). Being driven by the diversified parent geology, the soil structure, lateral and vertical extents of the soil profile, and hydraulic conductivity is spatially heterogeneous. The soils vary from hydrologic group B to group D, which represent infiltration rates from moderate to very low. Gilgelabay, Gumara, and Ribb are three perennial rivers that are characterised by a succession between bedrock types in their higher reaches, and alluvial types in their lower reaches and floodplains (Poppe et al. 2013). Depending on the lithologic units, different aquifer types such as fluvio-lacustrine sediment aquifer, intergranular aquifers, fractured and weathered upper basalts aquifers, highly fractured quaternary basalt aquifers, and vesicular-scoraceous aquifers can be found (Mamo et al. 2016).

The SWAT model is a continuous time, semi-distributed, process-based suitable model for small watersheds to river basin-scale developed to evaluate the effects of alternative management decisions on water resources and nonpoint-source pollution (Arnold et al. 1998, 2012). Since its development in the early 1990s, the model has constantly reviewed and updated (Neitsch et al. 2011). The model simulates the land phase of hydrological processes based on land use, soil, slope, and weather data. A catchment is partitioned into multiple sub-basins that are spatially connected, and further subdivided into hydrologic response units (HRUs), which are unique combinations of land use, soil and slope characteristics within a sub-basin. In SWAT simulation, HRUs represent a percentage of sub-basin that may not be contiguous or spatially identified (Arnold et al. 2012).

MODFLOW (McDonald & Harbaugh 1988) is a physically based, spatially distributed, three-dimensional and finite-difference groundwater flow model widely applied in groundwater flow modeling studies. It can simulate steady and non-steady flows in a saturated system, in which aquifer layers can be confined, unconfined, or a combination of confined and unconfined. Growing interest to use the model for different purposes has led to the development of multiple MODFLOW versions, like MODFLOW-2005 which was designed to simulate steady and non-steady flow in an irregularly shaped aquifer system (Harbaugh 2005), and MODFLOW-NWT which includes a Newton-Raphson formulation to improve solution of unconfined groundwater-flow (Niswonger et al. 2011).

SWAT-MODFLOW (Bailey et al. 2016) is a public domain model that links SWAT (Version 2012) with MODFLOW-NWT. In this integration, SWAT simulates land surfaces processes, crop growth, in-stream processes and soil zone processes, whereas MODFLOW-NWT simulates three-dimensional groundwater flow and all major groundwater sources and sinks (e.g. recharge, pumping, discharge to tile drains and interaction with stream networks). The source codes of the two models are combined into a single FORTRAN executable that runs independently (Bailey et al. 2016; Abbas et al. 2018).

While linking SWAT and MODFLOW-NWT, SWAT simulated HRU-recharge (i.e. water that exits the bottom of the soil profile) as recharge is passed to MODFLOW grid cells, and then groundwater-surface water fluxes calculated by MODFLOW are passed to the streams in SWAT for routing (Bailey et al. 2016; Bailey et al. 2017). Unlike other coupled models (Sophocleous & Perkins 2000; Kim et al. 2008; Chung et al. 2010; Guzman et al. 2015), SWAT-MODFLOW by Bailey et al. (2016, 2017) is capable of using SWAT and MODFLOW models of different spatial extents, has an efficient HRU-grid cell mapping scheme, and uses HRU disaggregation techniques to represent the recharge at each individual disaggregated hydrologic response unit (DHRU), and then map it to each MODFLOW grid cell according to the percent of area of the DHRU (Bailey et al. 2016, 2017; Guevara Ochoa et al. 2020). Moreover, the model has a graphical user interface that allows the SWAT-MODFLOW model to be set up and run (Bailey et al. 2017). The River package (RIV) of MODFLOW-NWT is used to simulate GW-SW interaction. Darcy's law is applied to calculate the volumetric flow rate of water through the cross-sectional area between the aquifer and the stream channel. A detailed model description including variable descriptions and mathematical equations is available in Bailey et al. (2016).

The SWAT model requires climate, soil, land use and topographic data. Three independent SWAT models were set up for the three catchments of Gilgelabay, Gumara, and Ribb based on available climate, land use, soil, and DEM data (Tigabu et al. 2019). Model setup included a 30 m DEM-SRTM (USGS 2016) for land surface topography, daily records of rainfall and temperature dataset (NMA 2016), basin land use and soil dataset (ADSWE 2016), and the streamflow data (MoWIE 2016). Dominant agricultural crops were considered and each crop was considered as a separate HRU. The simulations were carried out for the period of 1980–2014 with a 5-year warm period and the rest for calibration and validation. Measured streamflow data at the outlets of each catchment were used to test performances of each model. A multiple flow segment calibration approach using performance metrics and signature metrics was applied (Pfannerstill et al. 2014a, 2014b; Haas et al. 2016). Six thousand model runs were tested to end up with reasonable simulation outputs. Model performances were evaluated based on multiple objective functions such as Nash Sutcliffe efficiency (NSE), Kling–Gupta efficiency (KGE), percent bias (PBIAS), and standardized root mean square error (RSR) and were generally satisfactory for daily and at least good for monthly streamflow in the calibration and validation periods for the three catchments. Further information on calibration and validation periods, performance indices and employed data are found in (Tigabu et al. 2019).

Three MODFLOW models were set up for the three catchments using the SWATMOD-Prep interface (Bailey et al. 2017). MODFLOW-NWT (Niswonger et al. 2011) solves the groundwater flow equation to simulate spatially-varying groundwater head subject to aquifer stresses (e.g. GW recharge, GW discharge, and interaction of GW-SW with the stream network) and system parameters (hydraulic conductivity, specific yield, specific storage). It has a Newton solver that can only be used with the Upstream Weighting (UPW) package to allow drying and re-wetting of grid cells in a transient simulation. Each model was discretized into finite difference MODFLOW grid cells with a lateral dimension of 210 m by 210 m. The boundaries of the catchments were considered as boundary conditions for active and inactive cells.

Only the first unconfined layer of the aquifer system was considered due to the limitation of SWATMOD-Prep in creating a MODFLOW model with a single layer (Bailey et al. 2017). A two-dimensional groundwater flow condition was considered assuming that heads do not vary in the vertical direction. While modeling an unconfined aquifer, the fundamental assumption is that groundwater flows horizontally and the groundwater discharge is proportional to the saturated aquifer thickness (Dupuit 1863; Forchheimer 1886). This implies that a horizontal water balance is applied to the entire saturated aquifer thickness.

The average aquifer thickness was estimated from the geologic log information from 36 drilled wells distributed throughout the Lake Tana basin (AWWCE 2016). Based on the information extracted from these wells, the depth to the impervious layer from the groundwater level values vary between 0.5 and 12.5 m. This information was not sufficient to produce a map of spatially distributed aquifer thickness. Thus, a single-layer aquifer with an average thickness of 6 m was used and all the wells were considered open only to the first layer. Spatially distributed hydraulic conductivity and specific yield values (Figure 2) were assigned to each soil unit based on previously published values (Morris & Johnson 1967), and river bed material K (default value of MODFLOW) were used as input data for the MODFLOW models.

Moreover, maps of initial GW head (Figure 2) were prepared by interpolating water level data from hand-dug wells and boreholes that were collected from Amhara Water Works Construction (AWWCE 2016) and Tana Basin Development Authority (TBDA 2016) using an inverse distance weighting approach. Chebud & Melesse (2009) had used the one time record groundwater head data as the initial condition for their MODFLOW modelling work in the Gumara catchment. Saha et al. (2017) have also applied the same approach in their modelling study of temporal dynamics of groundwater-surface water interactions for a watershed in Canada.

The three calibrated SWAT models (Tigabu et al. 2019) that were constructed for the catchments of Gilgelabay, Gumara, and Ribb River were integrated with the MODFLOW models. SWATMOD-Prep (Bailey et al. 2017) was used to set up coupled SWAT-MODFLOW. The interface is based on the coupled SWAT-MODFLOW model code of Bailey et al. (2016), and links SWAT to MODFLOW-NWT. Within this framework, MODFLOW grid cells are created, SWAT-HRUs are disaggregated into DHRU that are individual contiguous areas within a sub-basin to allow HRU calculations to be geolocated (Bailey et al. 2017). The DHRUs and the MODFLOW grid cells were intersected, and river segments were defined using polygons (called river cells) that were identified to pass data between SWAT and MODFLOW. Then, SWAT-MODFLOW input files were written that link HRUs to DHRUs, DHRUs to grid cells, and grid cells (i.e. river cells) to the SWAT sub-basins (Bailey et al. 2016). Recharge data from each HRU were taken directly from the calibrated SWAT model. Linking of the two models followed several step by step procedures (Figure 3). Finally, SWAT-MODFLOW simulations were carried out for the three catchments, and the spatial and temporal interactions of GW and SW were evaluated for the three catchments. The River package of MODFLOW was used to simulate the GW–SW interaction and the Darcy equation was used to calculate the volumetric flow of water through the cross-sectional flow area between the aquifer and the stream channel (Bailey et al. 2016). Due to the lack of documented information on water withdrawal, pumping rates were not considered. The pumping rates are for domestic wells only, and therefore are very small in magnitude. As surface water irrigation is dominant in the study area, it was assumed that pumping does not affect the groundwater and streamflow hydraulic heads significantly.

The model is capable of simulating the groundwater hydraulic head for each MODFLOW grid cell, deep percolation (called recharge in SWAT-MODFLOW) for each HRU, volumetric exchange rates (m³/day) between the stream network and the aquifer for each SWAT sub-basin and volumetric exchange rate (m³/day) between the stream network and the aquifer for each MODFLOW River cell. Our simulations for GW-SW interaction were carried out on a daily time step from 1985 to 2014 and the analyses focused on the monthly groundwater head, GW-SW exchange, and GW recharge. R codes were written to extract and analyse the large size simulated outputs. Spatial distributions were produced by mapping the model output to the HRU location. These maps presented the spatial extent of source and sink areas, the variation of the hydraulic head, and GW recharge to MODFLOW. Moreover, temporal variations of GW-SW interactions were shown with the help of time series plots for those sub-basins that include gauging stations.

Conventionally, performances of integrated surface and sub-surface hydrological models have mostly been tested using measured streamflow and groundwater head values (Semiromi & Koch 2019). The current study was carried out in extremely data scarce area where there is no time series records of groundwater head data. However, the SWAT models were calibrated by considering parameters that affect both surface and subsurface processes (eg. deep aquifer percolation fraction, groundwater delay, alpha baseflow factor, saturated hydraulic conductivity, and groundwater revap coefficient); further information on list of model parameters are found in Tigabu et al. (2019).

Therefore, the performance of the coupled SWAT-MODFLOW models for the three catchments were evaluated against measured streamflow data on a monthly basis for calibration and validation (Table 1). Using streamflow data to evaluate the performance of the three models is reasonable as the streamflow is groundwater dominated in the study area (Setegn et al. 2008; Tigabu et al. 2019). Various performance indices like the Nash-Sutcliffe efficiency (NSE), (KGE) Kling-Gupta Efficiency, Coefficient of determination (R2), percent bias (PBIAS), and the ratio of the root mean square error to the standard deviation of measured data (RSR) are used to evaluate the goodness-of-fit of the models. In addition, flow duration curves (FDCs) are used to assess the model performance in different segments (high flow, middle flow, and low flow portion) of streamflow values.

Moreover, the mean simulated groundwater head values of the 36 wells were validated using observed static groundwater heads recorded during the course of drilling (AWWCE 2016). Use of GW head (static groundwater head) records of water supply boreholes measured at the time of drilling, before pumping, is a widely applied approach in Ethiopia to calibrate groundwater flow models due to the lack of monitoring wells (Berehanu et al. 2017). Here, the static groundwater head recorded during the time of drilling is assumed to be representative as it is recorded during an undisturbed state.

The models sufficiently reproduced the measured monthly time-series of streamflow. Figure (4) shows a good agreement between the measured and simulated monthly streamflow for the entire simulation periods in the three catchments. Only some peak flow events show smaller deviations. Deviations could result from the coarse spatial representation of rainfall in the input data. In addition, the under estimation of peak flows may be attributed to the SCS-curve number method used in the SWAT model, which does not consider the duration and intensity of rainfall (Tigabu et al. 2019). For all the three catchments, the statistical values of calibration and validation periods for streamflow indicated a good (Moriasi et al. 2007) temporal performance (Table 1). Simulated results indicate that the streamflow systems of the catchments are dominated by groundwater discharge (Table 2). The FDCs of the simulated streamflow show good agreement with the measured ones (Figure 5) especially for middle and low flow segments, providing additional confidence in the representation of groundwater flows.

Comparing the static groundwater head values of the 36 wells against the mean long-term simulated groundwater head values (Figure 6), very good results are achieved for Gilgelabay and Ribb catchments, whereas for the case of the Gumara catchment the statistical values are lower (Table 1). It has to be noted that the statistical values computed for the groundwater head are influenced by the fact that the mean simulated value for the entire modelling period was compared to an observed value measured during the time of drilling. Thus, this evaluation indicates a reasonable spatial model performance in the long term.

In this study, we identified source and sink areas of aquifer and stream systems. The GW-SW interactions within the three catchments vary spatially and temporally throughout the simulation period of January 1985 to December 2014. This is demonstrated by the volumetric exchange rate (m³/day) between the stream network and the aquifer for each MODFLOW river cell, and between the stream network and the aquifer for each SWAT sub-basin. As expected from the classical river hydraulics, the stream systems in the highland areas are mainly a source for the aquifer indicating that water is flowing from the river reach to the aquifer, and that the aquifer is a sink for the stream system. The spatio-temporal variations in the GW-SW exchange rates depend on the head difference between the aquifer and the stream and on the aquifer properties such as hydraulic conductivity, specific storage and initial groundwater head (Kim et al. 2008).

In the Gilgelabay catchment, the long-term mean annual volumetric GW exchange rate for each of the SWAT sub-basins is positive indicating that the aquifer system is a source area for the stream systems. The net annual mean contribution from the aquifer to the stream network in the sub-basin at the catchment outlet is 170 m³/day. In contrast to Gilgelabay, the mean annual volumetric GW-SW exchange rates between the stream network system and the aquifer for each SWAT sub-basin of Gumara and Ribb are not unidirectional (Figure 7). The aquifer system in the Gumara catchment is mainly controlled by the surrounding SW system. However, the flow system at the downstream reaches (sub-basin at catchment outlet) which is found in the floodplain area (34 km²) is from the aquifer to the river network system. This indicates that the floodplain is serving as a storage area for both surface and groundwater. The current result confirms the findings of Dessie et al. (2014) who reported that the floodplain of Gumara catchment was recharged by the groundwater during the rainy season. Chebud & Melesse (2009) also reported that 0.09 billion cubic meters of water flow out from the Gumara floodplain aquifer to Lake Tana during the wet season. Besides the downstream reaches (floodplain area), there are a few sub-basins from the middle and upstream part of Gumara catchment in which the flow system is from the aquifer to the stream network (Figure 7).

Annual mean exchange fluxes show temporal variability for all the three catchments. In the Gilgelabay catchment, however the net annual mean flow is from the groundwater to the channel, and the year to year variability during 1985–2014 reaches about 30%. The temporal pattern of GW-SW exchange resembles some features of the observed streamflow measured at the outlet point (coefficient of correlation r = 0.20, Figure 7). The mean annual GW-SW interaction rate of the Gumara catchment at its outlet showed that the aquifer system was gaining a significant amount of water from the stream before the year 1994, and the flow magnitude gradually decreased afterwards (Figure 8). This decrease of the flow magnitude of water from the river network to the aquifer over time could be related to the increase of runoff, which has been reported by Tigabu et al. (2020), and less seepage of water from the stream. Likewise, the flow dynamics of the Ribb catchment are mainly controlled by the surrounding SW system for the period before 1990, whereas after 1990 the flow direction was changed from the aquifer to the stream network system (Figure 8). This might be due to a change in the number of rainy days and rainfall duration and intensity (Teshome 2016; Tigabu et al. 2020). As the total rainfall amount did not change significantly, the rainfall became more intense and its duration became shorter which very likely results in more surface runoff generation during the wet season, lower stream water level, and less seepage to the aquifer system. This is one reason for the fact that the stream systems have become gaining and perennial.

Moreover, the simulation indicates that the flow direction of GW-SW interaction is more variable at the daily scale compared to the annual scale (Figure 9). The prevailing flow direction is bi-directional (i.e. from the stream network to the aquifer system and vice-versa). The strong dynamics of the daily GW-SW interaction could be related to the large variations in the water table during a hydrologic year, whereas the changes in the flow domain may be attributed to the variation in topography and precipitation. The large variation in the water table during the year indicates that a stream segment could receive water from groundwater for a portion of the year and lose water to groundwater at other times (Figures 10). A comprehensive discussion was given by Sophocleous (2002) on how GW-SW exchange rate and flow direction vary with precipitation and hydraulic head.

Similar to the daily exchange fluxes, flow patterns are highly dynamic on a monthly basis for all the catchments. These are reflected by the considerable seasonal differences on the magnitudes and direction of exchange fluxes. Exchanges fluxes are higher in magnitude during the wet months (July, August, and September) and gradually diminish for the rest of the months, which signifies the seasonal variations of aquifer heads. Kim et al. (2008) found similar results to the Musimcheon Basin in Korea when they applied the integrated SWAT-MODFLOW model to simulate spatio-temporal distribution of GW recharge rates and GW evapotranspiration. Both an exchange of water from the stream network to the aquifer and vice-versa can be observed. The former especially occurs during rainy months indicating that the water levels of the stream systems respond fast (Figure 10). A similar study conducted in Pateira de Fermentelos, Portugal proved that the surface water level changes fast in response to aquifer discharge during rainy months (Sena & Condesso de Melo 2012). This might also result from a large variation in the total daily rainfall (Hassan et al. 2014).

Moreover, the magnitude of the exchange rate spatio-temporal variability between the dry and wet seasons depends on the geographic location. In Gilgelabay catchment, the seasonal flow pattern is from the aquifer to the stream system for both dry and wet seasons. During the wet season the mean contribution of the aquifer to the stream network is about 390 m³/day, whereas the dry season contribution is 60 m³/day. The lowland sub-basins showed strong dynamics of GW-SW exchange especially in the rainy season, whereas in the highland sub-basins the prevailing flow direction is mainly from the channel to the groundwater. The annual mean exchange rate of both dry and wet seasons is positive for each SWAT sub-basin and negative for each MODFLOW river cell indicating that water is leaving the aquifer to the stream channel. As stated previously, bi-directional flow patterns are observed in the Gumara and Ribb catchments. Flow patterns vary between dry and wet seasons. There is mainly stream network system loss water during the wet season and water gain from the aquifer system during the dry season (Figure 11).

The different flow patterns between the dry and wet seasons are related to the aquifer response to rainfall events. The only input variable to the hydrologic system of the three catchments is rainfall, and the hydrologic year starts in June. July and August are months when peak rainfall occurs and streams are at their full stage. During this time, the aquifer systems are mainly sinks to the stream system especially for the catchments of Gumara and Ribb. A very interesting finding is the existence of a time lag between peak rainfall event and the groundwater head. As stated above, the heavy rainfall events occur in June, July, and August, while the groundwater heads start to increase in July (one month later) and attain a peak one or two months after the peak rainfall (Figure 12). This indicates that aquifer reactions to rainfall changes are very slow compared to surface water systems. For this reason, the stream systems are gaining during the dry months and losing during the rainy months.

The spatio-temporal variations of the GW-SW interactions among the three catchments could be linked to the heterogeneity of the catchments with respect to topography, soil, geology, as well as to the seasonal variations and the temporal pattern of the precipitation in each catchment. The magnitude of GW-SW interaction is highest during the rainy (summer) season compared to other seasons. The high degree of variability on the daily and seasonal GW-SW exchange rates might be attributed to a large variation in daily rainfall (coefficient of variation ranging from 100 to 200%). Moreover, the GW–SW exchange rate between the stream network and the aquifer for each SWAT sub-basin varies spatially within each catchment and among the three catchments (Figure 7). On the one hand, mean values of GW-SW exchange rates are primarily positive for Gilgelabay catchment with daily discharge values from the aquifer system ranging from 1 to 2,000,000 (m³/day) per sub-basin, whereas seepage values from the river network system to the aquifer vary between 2,050 and 124,000 m³/day per sub-basin, indicating that discharge from the aquifer system is higher than the seepage from the stream network system. On the other hand, the mean values of GW-SW exchange rates for Gumara and Ribb catchments are mainly negative, i.e. seepage from the stream network is greater than the discharge from the aquifer system. One reason for this could be the higher influence of the maximum elevation especially for the Gumara catchment which has a maximum elevation of 3,650 m above sea level. This indicates that there are differences in the GW-SW interaction among the three catchments and within each catchment.

Overall, the GW-SW interaction in the Gilgelabay catchment differs from the catchments of Gumara and Ribb catchments due to the differences in the hydraulic conductivity and specific yield values assigned based on the soil distribution as well as to the difference in the rainfall, river gradient, and hydraulic head. For instance, in Gilgelabay catchment the hydraulic conductivity values vary between 0.1 and 0.35 m/day and the specific yield varies from 10% to 18%, whereas in Gumara and Ribb the hydraulic conductivity values are between 0.1 and 1.25 m/day, and the specific yield varies from 12% to 28%. These relatively smaller hydraulic conductivity and specific yield values in Gilgelabay catchment caused a smaller GW-SW interaction rate compared to Gumara and Ribb catchments.

In general, from the above analyses one can certainly conclude that the GW-SW systems in the study area are hydraulically connected, and the GW-SW interactions differ between the three catchments.

The coupled SWAT-MODFLOW model simulates the amount of recharge in mm to MODFLOW for each HRU (called deep percolation in SWAT-MODFLOW) and the volumetric flow rate of recharge (m³/day) to each MODFLOW grid cell. The average monthly deep percolation during the period 1985–2014 varies from 0 to 90 mm, 0 to 75 mm, and from 0 to 60 mm in the Gilgleabay, Gumara, and Ribb catchment, respectively. The volumetric flow rate of recharge (m³/day) provided to each MODFLOW grid cell shows the spatial variation of recharge to the aquifer (Figure 13). Higher recharge rates correspond to the highland areas and low recharge rates are associated with the lowland areas. This is due to the fact that most of the time highland areas have groundwater heads below stream stage and lowland areas have groundwater heads above stream stage. Bailey et al. (2016) also found higher recharge rates in the highland areas of the Sprague River Watershed in southern Oregon, USA, by applying the same model.

Simulated cell-wise annual average recharge rates provided to each MODFLOW grid cell range from 10 to 140 m³/day in the Gilgelabay catchment, from 6 to 105 m³/day in the Gumara catchment, and from 4 to 90 m³/day in the Ribb catchment (Figure 13). Additionally, there are considerable differences between the wet and dry season recharge rates. The vast majority (75–91%) of recharge takes place during the wet (summer) season as the main source of recharge to the groundwater system is rainfall. Whereas only 9–25% occur during the dry season. This indicates that the regional groundwater flow system is seasonally limited. The current result is in agreement with the result of Woldeamlak et al. (2007) for the Grote-Nete catchment, Belgium that got 86% of its recharge during the wetter (winter) season (but different climatic zone). The recharge areas are mainly associated with the highland areas and with areas dominated by soils with a coarse texture like Eutric Leptosols and Eutric Regosols. GW discharge areas mainly correspond to the stream network and floodplains, which is in agreement with Sear et al. (1999). Aquifer systems of such a type (GW discharge areas) are mainly recharged locally and depend on the amount of precipitation, evaporation and temperature. However, at times of high river flows (wet season), water moves from the stream system into the groundwater system. This is clearly revealed by the GW-SW interactions between the SWAT sub-basins and stream network systems during the wet and dry seasons (Figure 11).

In hydraulically connected GW-SW systems, the groundwater head is an important variable that controls the resulting exchange between the stream network and the aquifer system. The main driving force of GW flow is the groundwater head potential that is influenced by the topography in the case of unconfined aquifers (Fan et al. 2007). As stated in Fan et al. (2007), the rate and direction of groundwater flow at a given location is driven by the hydraulic gradient, which in turn is determined by the boundaries of the system and the location and strength of recharge and discharge (topographic control). Simulated mean daily groundwater head values for the given modeling period are shown in Figure 14. Both the highest (3,640 m) and the lowest (1,783 m) groundwater head values were estimated for Gumara catchment (Figure 14). The temporal groundwater head differences for MODFLOW grid cells at the outlets for dry and wet seasons vary between 1 and 2 m. The groundwater head values show a similar pattern as the topography for all the three catchments indicating that the shape of the water tablet was highly influenced by topography. This is in agreement with research findings of Sena & Condesso de Melo (2012), who produced a water table map from water table and surface water level data measured in the field and reached the conclusion that the shape of the water table maps was greatly influenced by the topography and the streams. The depth from the ground surface to the water table varies between 1 and 30 m in the case of Gilgelabay and from 1 to 15 m for Gumara and Ribb. It is shallow in the lowlands (1–2 m and deeper in the highlands and at the water divides ranging from 1 to 30 m.

As stated in Fan et al. (2007): The rate and direction of groundwater flow at a given location is driven by hydraulic gradient, which in turn is determined by the boundaries of the system and the location and strength of recharge and discharge (topographic control).

As the Lake Tana basin is considered one of the development corridor sites by the Government of Ethiopia, understanding GW-SW water interaction and implementing combined use of GW and SW is the top priority for water managers to contribute to sustainable development in the region. The present study focuses on the investigation of GW-SW interaction of three catchments in the Lake Tana basin. Our findings provide insights for water managers and policy makers for future planning. We found that GW and SW are hydraulically connected and the stream network systems are dominated by GW flow in Gilgelabay catchment with annual net aquifer discharge of 170 m³/ day. In Gumara and Ribb catchments, the annual flow patterns are bi-directional, i.e. both aquifer discharge and seepage from the stream network system are observed. Thus, sustainable water resources management will be achieved if the future water resource development plan puts a particular focus on the combined use of groundwater and surface water. The groundwater head and GW-SW interaction showed considerable spatial and seasonal variations in all three catchments. In Gumara and Ribb catchments, the wet season flow pattern is predominantly from the river to the aquifer system, whereas the dry season flow is to the river system. This implies that GW and SW withdrawal e.g. for drinking water use affect this interaction differently in the wet and dry seasons. Spatial patterns of the groundwater head should be considered to define GW recharge and discharge zones prior to water withdrawal from the aquifer system. Moreover, understanding the GW-SW interaction is important to characterize the general ecosystem status as the streamflow in all the three catchments is connected to the GW. Hence GW is critical with regard to water availability as well as with regard to the control of environmental flow.

The aim of this study was to investigate GW-SW flow dynamics in the Lake Tana basin by considering three catchments as case studies (Gilgelabay, Gumara and Ribb) using the coupled SWAT-MODFLOW model. Results from this study indicated that SWAT-MODFLOW can simulate the surface and subsurface hydrologic processes in the study areas and can provide insights into GW-SW interactions. The capability of the model to produce reliable results was reflected by the good match of the simulated and observed streamflow data as well as by the comparison between observed hydraulic head values measured at water wells during the course of drilling and mean simulated hydraulic head values. This coupled hydrologic modeling approach of the dynamic GW-SW interaction provides a basis for further study of groundwater and surface water interaction and combined use of groundwater and surface water in the Nile and other catchments in the sub-Saharan region. The key findings specific to the study catchments are:

• 1.

  The aquifer and the stream network systems are hydraulically connected in all catchments in both wet and dry seasons. This was proven by the simulated volumetric GW-SW exchange rates between each SWAT sub-basin (aquifer) and stream network, and the stream network and the aquifer for each MODFLOW river cell. The vast majority of water flow is from the groundwater to the channel, particularly in Gilgelabay catchment.

• 2.

  In all catchments, the groundwater and surface water flow system is more variable on the daily time scale than on the monthly time scale. This was shown by the simulated results of GW-SW exchange rates, water recharge to the groundwater head and discharge from the aquifer system.

• 3.

  Due to the variability of rainfall, the groundwater head values vary on a monthly basis. The groundwater head reaches its peak one to three months later than the rainfall. In the studied catchments, the peak rainfall occurs in the month of July or August, leading to an increase in the simulated groundwater head in the month of July that reaches its peak in October to December. This result illustrates the slow groundwater response to rainfall changes.

• 4.

  Despite the fact that daily GW-SW exchange rates are highly variable in the Gilgleabay catchment, the annual volumetric exchange is from the aquifer system to the stream network system indicating that stream systems are gaining.

• 5.

  In both Gumara and Ribb catchments, the annual volumetric GW-SW exchange rates are more dynamic in time compared to Gilgelabay catchment. This was suggested by the GW-SW exchange rates between each SWAT sub-basin and stream network system (Figure 6). During the wet seasons, the flow is predominantly from the stream system to the aquifer, while it is from the aquifer to the stream network system during the dry season (Figure 12).

• 6.

  The annual GW-SW exchange rates at the outlet sub-basins show a change in the flow direction of water for Gumara and Ribb after 1995 and 1991, respectively (Figure 8). Before the mentioned years, the flow system was from the channel to the groundwater and after these years the flow direction changed from the groundwater to the channel. This might be related to the decline with the number of rainy days (especially the dry season) since 1990 while the total annual rainfall was not significantly changed.

• 7.

  The modelling results have been checked for their reliability using measured and simulated streamflow time series and groundwater head values measured during the drilling of wells. However, the model could be improved if regularly measured time series of water table elevation and stream stage values were available.

• 8.

  This study primarily focuses on the setup of a coupled SW and GW model to understand the hydrodynamic condition of the GW-SW system in the study area. The model is valuable because it can be used to assess climate change impacts on groundwater resources in the future.

We would like to acknowledge funding for a doctoral study grant from the Federal State of Schleswig-Holstein, Germany, through the Landesgraduiertenstipendium of Kiel University. We are also thankful to the Ministry of Water, Irrigation, and Electricity, the National Meteorological Service Agency of the Government of Ethiopia, and the Amhara Water Works Construction Enterprise for their support by providing hydrological, climate, and groundwater well data. The authors are very grateful to the three anonymous reviewers for their thorough reviews and constructive comments which greatly improved the manuscript.

All relevant data are included in the paper or its Supplementary Information.