Application of coupled WetSpass-M and MODFLOW models to estimate spatial–temporal water balance components in the Chemoga watershed, Ethiopia

Groundwater is an essential global supply of fresh water for drinking, industry, agriculture, and the preservation of ecosystem services (Amiri et al. 2022). Water availability decreases as a result of increasing human intervention that modifies the hydrological process. Furthermore, groundwater levels are dropping as a result of pollution, climate change, and inadequate groundwater recharge, in addition to the world's rising population. Understanding the spatial extent and change of groundwater levels is vital to safeguarding available water resources, especially as a key source for drinking water (Salem et al. 2019; Amiri et al. 2022). Effective management of water resources and subsurface fluid and contaminant models depends on a deeper comprehension of the unique and temporal distributions of water balance components, especially groundwater recharge, especially as these resources become the primary source of drinking water (Salem et al. 2019; Amiri et al. 2022; Demissie et al. 2023).

The process of delivering water to the groundwater storage below the surface is known as groundwater recharge, and it is indicated by a change in the water table's level (Han et al. 2017). Understanding recharge processes and their quantification is vital for sustainable management and protection of groundwater resources (Xu & Beekman 2019). Since groundwater resources are sensitive functions of climatic factors, geological formation, topography, soil properties, and land use types (Gebremeskel & Kebede 2017), a thorough understanding of watershed physical and biological characteristics is important. Precisely assessing recharge is among the most difficult aspects of water balance since it varies greatly both in space and time and is influenced by numerous variables (such as terrain, climate, vegetation, soil, and geology) (Mekonen et al. 2023).

There is a strong groundwater recharge variability across the Ethiopian volcanic rock aquifers (Yenehun et al. 2022). The Chemoga watershed, which is one of the major basins in the Ethiopian volcanic plateau, has high hydrogeological and topographical variabilities. In addition, the region experiences brief, rainy summers and lengthy, dry winters. Therefore, it is anticipated that there will be significant regional and temporal variation in the groundwater recharge–discharge processes. The Chemoga watershed is an area of increasing competitive demand for water for agriculture and water supply. The groundwater recharge study of the area has not been based on physically varying methodologies for estimating the long-term average. Scientific research in the Chemoga watershed was not undertaken in accordance with the quantification and mapping of groundwater recharge areas in the sub-basin. The components of the water balance were not properly defined, and the hydraulic head distribution in relation to stress was not modeled. Lack of a good understanding of groundwater recharge was a serious concern for sound and suitable groundwater management in the watershed, given the high pace of population growth and increased reliance on groundwater. Groundwater recharge is typically quantified using a variety of methods, such as empirical methods, hydrological budget (HB), chloride ion mass balance method, regression curve displacement method (Rorabaugh method), statistical approaches, water table fluctuation (WTF), and numerical methods like water balance simulation.

However, a number of variables, like the features of the catchment, the accessibility, and the precision of field data, affect how reliable the approaches are (Mekonen et al. 2023). Although most approaches are unable to calculate the recharge precisely, they do not take temporal and spatial variability into consideration (Mekonen et al. 2023). The hydrological models such as SWAT (Yifru et al. 2020), SVAT (Doble & Crosbie 2017) and DREAM (Sheffer et al. 2010) becoming more popular and widely used for estimating recharge; however, they are data-intensive and require extensive parameterization and detailed hydroclimate and hydrogeological data (Gebremeskel & Kebede 2017; Mekonen et al. 2023). Groundwater models are often applied to simulate water flow in aquifer. One important feature in groundwater models is recharge. This paradigm has various drawbacks that should be taken into account despite its advantages. One of the primary drawbacks of these models is their unreliability due to the improper estimation of recharge and its spatiotemporal dispersion (Dowlatabadi et al. 2024).

Nowadays, numerous studies have emphasized the importance of utilizing a linked model for comprehensive research, particularly for surface and subsurface models, in order to achieve the benefits of diverse models integrated modeling. Because surface water and groundwater are related elements of hydrology, they should be examined together, for instance. In order to evaluate groundwater recharge in the Takelsa multilayer aquifer in northeastern Tunisia, Ghouili et al. (2017) developed a coupling model between WetSpass and MODFLOW. That indicates a WetSpass model with a modified groundwater model could better represent the water balance in the study area. Dowlatabadi et al. (2024) assess the spatiotemporal distribution of recharge and simulate groundwater flow in the Neyshabur watershed using WetSpass-M and MODFLOW models. These results are considered more appropriate than the simulated recharge by SWAT in the MODFLOW model in the Neyshabur watershed; additionally, (Dowlatabadi et al. 2024) shows the higher precision of the recharge calculation of the WetSpass-M model. Aslam et al. (2022) developed hydrological modeling of the aquifer's recharge and discharge potential by coupling WetSpass and MODFLOW for the Chaj Doab, Pakistan. These results indicate that the recharge caused a larger drop in the water table.

Recently, energy and water transfer among plants, soil, and the atmosphere under a quasi-steady-state (WetSpass) model (Salem et al. 2019; Amiri et al. 2022; Demissie et al. 2023) has been widely used for groundwater recharge assessment. Abdollahi et al. (2017) developed a WetSpass-M model by downscaling the seasonal resolution to a monthly scale. It has been demonstrated that the WetSpass model may be used to more accurately characterize recharge, including its variation across global geographic locations. For the Chemoga watershed to manage water resources sustainably and effectively, a deeper comprehension of the temporal and spatial fluctuations of water balance components, particularly actual evapotranspiration (ET), surface runoff, and recharge, is essential. The earliest version of the WetSpass model was only able to simulate the yearly and seasonal variations. For this study, the limitation was resolved in the WetSpass-M modified version for simulating monthly variation, which uses monthly climate data. Therefore, WetSpass-M is able to provide a more accurate evaluation of water balance components throughout time and space with the use of large monthly-scale data. Since climates tend to vary widely across topography, any hydrological model simulation at a monthly time step is more suited than one at the seasonal or annual scale (period and type) for the assessment of water resources.

The spatially distributed recharge output of the WetSpass-M model can improve the prediction of simulated groundwater level and the locations of discharge and recharge areas for a steady-state groundwater model (Aslam et al. 2022; Salem et al. 2023). Hence, WetSpass-M and the groundwater model run simulations consecutively while exchanging inputs for recharge and groundwater, respectively. MODFLOW was used to simulate the hydraulic head distribution using the groundwater recharge distributions acquired by WetSpass-M (Aslam et al. 2022; Salem et al. 2023). As a result, the groundwater level and discharge zones have a steady solution. Overall goal of this study was to assess the long-term spatial distribution of monthly, seasonal, and annual components of the water budget and spatial groundwater levels. This study represents the first work to assess the long-term monthly, yearly, and seasonal components of the water budget in the Chemoga watershed, and this information, together with aquifer geometry and other boundary conditions, were used to develop the groundwater model.

The Chemoga watershed, which is part of the Abay River Basin, covers an area of approximately 364 km² and is situated at the geographical coordinates of 10°18′N and 10°39′N latitude and 37°42′E and 37°53′E longitude (Figure 1). The Chemoga River originates from Choke Mountain and joins the Blue Nile in the East Gojam Zone, located in the Amhara Regional State of Ethiopia at an elevation of 2,401 m above sea level. The Chemoga watershed is one of the headstreams of the Blue Nile and contributes around 86% of its total annual discharge. In recent years, the amount of water carried downstream by the Blue Nile has declined significantly. The Chemoga watershed is a part of this degraded and degrading basin for the representative of conditions in large part of the temperature locally known as dega and wurch in the northern western highlands. The mean annual temperature is 14.5 °C, ranging between 13.2 °C in July and August and 17.3 °C in March. The average annual rainfall is 1,453 mm; more than 75% of the total precipitation occurs in the 4 months from June to September. The major types of soil groups in the study area are clay, and loam. Seven different types of land use have been identified for the Chemoga watershed. These include water, cultivated land, forest land, grassland, shrub land (bush land), urban settlement, and bare land. According to the Ethiopian geological survey, the lithologic units of the study area include Middle Yujbe basalts, Ignimbrite and tuff deposits, Debre Markos basalts, Lumame basalts, quaternary alluvial deposits, and others that fall in the choke shield volcano group like Rob Gebeya basalt, Arat Mekeraker basalt, and Kutye basalt (Teklebirhan et al. 2012).

In this study, datasets such as physical-geographical, climatological, and hydrogeological data were collected and used by the model. The model made use of climatological and physical-geographical variables, such as temperature, wind speed, rainfall, and potential evapotranspiration (PET), as well as the digital elevation model (DEM), slope, soil texture, and LU/LC.

The hydrogeological data such as GWD was also used. The DEM for the watershed was acquired from ALOS-PALSAR. It was downloaded from the Alaska Satellite Facility (ASF) (https://www.asf.alaska.edu) at a spatial resolution of 12.5 m. The satellite images were downloaded from the United States Geological Survey (https://earthexplorer.usgs.gov), and used to classify the LULC thematic map (LULC). Time series climatological data with durations (1990–2021) were obtained from the National Meteorological Agency (NMA). From the daily precipitation the average monthly, annually and seasonal precipitations were derived. The stream flow data (1990–2012) was collected from the Abbay basin authority and used as to validate the WetSpass-M model. Groundwater level measurements were taken from accessible wells of the study area and used as to validate the MODFLOW model. Different software's and tools such as ERDAS Imagine 2014, ArcGIS 10.4, wetSpass-M model, MODFLOW, zonal statistics as table tools in ArcGIS, dep meter and GPS were used in this study. The input data of WetSpass-M model are listed in Table 1.

Two types of input data are needed for the WetSpass-M model: parameter tables and GIS-grid maps (Amiri et al. 2022; Yenehun et al. 2022). Weather information such as rainfall, PET, average temperature, and wind speed are all included in GIS grid maps, as well as GWD, slope, soil type, and topography. Land use, and soil characteristic was specified as parameters tables in the model (Batelaan & De Smedt 2007; Gidafie et al. 2016; Gebremeskel & Kebede 2017; Dereje & Nedaw 2019; Amiri et al. 2022; Demissie et al. 2023; Dowlatabadi et al. 2024). These tables were connected to the maps as attribute tables. The land use attribute table include parameters related to land use type; such as rooting depth, leaf area index (LAI), vegetation height, canopy cover among many others. The soil parameter table contains soil parameters for each textural soil class such as field capacity, wilting point, permeability. The runoff attribute is considered to be universal, because a certain combination of slope, land use, and soil type will produce a certain fraction of runoff independent of location (Amiri et al. 2022). The data have to be transformed using a spatial analysis technique appropriate for the type of data because of the many data sources and formats that have been found. The interpolation approach is employed because it can be quickly and easily developed for a specific purpose using limited and sparse measurement data.

For this study, ordinary Kriging interpolation methods were used in ArcGIS in a spatial analysis toolbox. Kriging is frequently utilized in many disciplines, especially in a heterogeneous setting, and is most suitable when a spatially associated distance or directional bias in the data is known.

Before using the data for further analysis, it is important to make sure that the data are homogenous, consistent and correct. Inconsistence may have resulted from changes in observation producers, changes in exposure to the gauge, changes in land use that make it unreasonable to maintain the gauge at the old location, and the occurrence of frequent instrument damage. Therefore, in this study, the accuracy and reliability of the input data were checked before further analysis. The inconsistency of the data were checked using double mass curve analysis and Alexanderson's SNHT test was also used to test the homogeneity for monthly rainfall.

Initially, the model input grid maps were prepared at the spatial resolution of the coarser biophysical factor, 12.5 × 12.5 m. The grid maps at 12.5 m resolution were quite sparse, making it difficult to see the distinct spatial variability of the hydro-meteorological observations in the watershed (Demissie et al. 2023). Therefore, the model input grid maps were further resampled to a spatial resolution of 30 m, and this spatial scale was found to produce coherent outputs during model simulation. All hydro-meteorological and biophysical maps were prepared, then resampled at a resolution of 30 × 30 m, and then converted into ASCII files so that they could be used in the model.

Land use/land cover (LULC) type has a major effect on groundwater recharge or infiltration (Amiri et al. 2022; Siddik et al. 2022). Another useful use of the LULC is estimating the values of vegetative parameters such as the LAI and evaporative zone depth. Surface evaporation and transpiration are both regulated by the LAI parameter (Amiri et al. 2022). One of the most important regulating elements in basin hydrology is LULC (Amiri et al. 2022). The LULC data for the watershed was downloaded from United States Geological Survey Global Visualization Viewer website (https://earthexplorer.usgs.gov). Satellite images for the study area were acquired at Path 166 and Row 53 on February 8, 2020.

Obtaining an accurate LU/LC thematic map through image classification is crucial in many remote sensing applications, such as detecting climate changes, monitoring and managing the environment, and tracking hazards and urbanization expansion (Abbass et al. 2023). The selection of a reliable classifier approach is essential to achieving an accurate estimation of the LU/LC thematic map (Abbass et al. 2023). In remote sensing, one of the most known types of classification that can be used to achieve classified outputs is called the Supervised Classification technique (Abbass et al. 2023). Using the standard ERDAS IMAGINE supervised image classification method, five different types of land use have been identified for each watershed.

To quantify the accuracy of land cover classification, the commonly applied confusion matrix and kappa coefficient (Defourny et al. 2019) were used. A confusion matrix was prepared, considering 200 randomly selected reference ground truth points for evaluation. Using Google Earth as a reference, 200 randomly selected points were selected and compared with each corresponding classification ground control points measured by global positioning system (GPS) for validation. There are seven LULC types such as agricultural (cultivated land), forest land, grass land, shrubs (bush land), and settlement (urban), water, and barren land have been identified for the watershed These are agricultural (cultivated land), forest land, grass land, shrubs (bush land), settlement (urban), water, and barren land. Each class accounts for 60% moderately cultivated, 15% grassland, 10% afro- alpine, 5% settlement, and 2% dominantly cultivated, respectively (Figure 3(c)).

Soil is a main physical characteristics that control runoff and recharge. Soil permeability is a function of soil infiltration capability and determines the soil's storage capacity as well as the degree of hostility with which water flows into deep layers. The soil's texture has a major impact on the infiltration capacity of the soil. The highest rates of infiltration are found in sandy soils, whereas heavier clay and loamy soils demonstrate lower rates of infiltration and greater surface runoff (Amiri et al. 2022). The soil texture map used as an input to the WetSpass model was collected from the Agriculture Bureau in Bahir Dar as a shapefile that considers the physical properties of soils, including texture and accessible water content, for every type of soil, including bulk density, hydraulic conductivity, and organic carbon content. The major types of soil groups in the study area were loam and clay. They cover 81% of loam and 19% of clay (Figure 3(d)).

This data includes four parameters (rain, temperature, wind and evaporation) as inputs to the Wetspass-M 1.3 model. Twelve raster maps were created for each climate parameter on a monthly basis; the total number of these maps was 48. All these rasters were converted to American Standard Code for Information Interchange (ASCII) files, in order to be input correctly to the model. These rasters have a cell size (30 m × 30 m) in order to run the model conveniently.

The identification of hydrographic stratigraphic units are important in defining the number of layers controlling groundwater flow within the aquifer system (Flores et al. 2023). Hydrographic stratigraphic units are frequently used to categorize geologic strata or formations for hydrogeological studies. In this study area, the geologic formations in the area can be broadly classified into two main aquifers: the unconsolidated sediment aquifer and the volcanic aquifers (Ketemaw et al. 2021). The unconsolidated sediments can be found in the lowland plain portion of the study area as alluvial fans along riverbanks and at the base of mountains. Alluvial deposits such as sand, gravel, pebbles, and boulders make up this aquifer. The results of geophysical surveying and the geological logs from boreholes indicate that the average thickness of this aquifer is around 30 m. As understood from lithological logging data, the major aquifer of groundwater is moderately fractured and weathered basalt and river gravel. Lithological logs show variable aquifer lithology with different degrees of fracturing and weathering. However, due to their permeability, the different hydro-stratigraphic units are hydraulically interconnected. Therefore, a single layer, steady-state condition, and unconfined aquifer layer were considered for the modeling practice in this study.

ModelMuse was used to develop the groundwater flow model. ModelMuse is a graphical user interface (GUI) designed to create model input files for MODFLOW (Akter & Ahmed 2021). In order to develop a conceptual model that is suitable for the modeling aim, a set of reasonable assumptions that minimize the actual situation are used. The following presumptions applied to the modeled area: (i) It was believed that the system operated in a steady-state throughout the year; and (ii) that the extent of the problematic geological formations was horizontal. ModelMuse was used to create a one-layered MODFLOW model using the basin's grid, which consists of 440 rows and 203 columns. The lower layer's bottom was aligned with the bedrock elevation, while the upper layer's top surface was set to match the elevation of the groundwater surface.

To construct a model, you need the following two input packages: (i) model properties and (ii) model boundary conditions. Among the inputs for the model properties are the aquifer parameters and initial heads. Since the groundwater flow model was single-layered, only horizontal hydraulic conductivities were important. In order to generate initial heads for the whole model, starting heads were measured directly from existing boreholes and interpolated within the model. Using the inverse distance weighting (IDW) method, the observation heads were interpolated. Boundary conditions have a great influence on the computation of flow velocities and heads within the modeled area. The numbers of boundary conditions are available in MODFLOW, such as general head, recharge, constant head, river, and no flux boundary conditions. In this study, the recharge and evapotranspiration boundary conditions are considered as model inputs. The recharge package is designed to simulate distributing recharge to the groundwater system (Dong et al. 2012). The recharge value represents the amount of water that goes into the groundwater system rather than the amount of precipitation. The evapotranspiration (ET) package is used to simulate the effects of plant transpiration and direct evaporation by removing water from cells during a simulation. The spatially distributed recharge and evapotranspiration are estimated using the WetSpass-M model, which can be imported in the form of a shape file into ModelMuse for the simulation of the flow of groundwater. The upper boundary of the system was simulated as specified flux as recharge was applied to the water table. The lower boundary was assumed to be a no flow boundary. Since the catchment's remaining boundaries were not specified as distinct boundary conditions, MODFLOW automatically regards them as inactive. The value ‘0’ was assigned to cells external of such boundaries to make them inactive, and no flow occurs across such boundaries.

Hydraulic conductivity is the most essential aquifer parameter that determines the flow system of a model (Vogelgesang et al. 2020). It is a measure of the ability of fluids to move through aquifer media. It is dependent on the properties of both the porous medium and the fluid. In this model, the hydraulic conductivity and specific yield were obtained from the pump test analysis, the study of Ketemaw et al. (2021), the geologic series, and a literature review. The hydraulic conductivity of an aquifer has a directional value, and in this model, the model area is conceptualized as an isotropic, single-layer unconfined aquifer. It has no vertical flow and has the same value in the x and y directions. The initial hydraulic conductivity value used as an input was gradually fine-tuned and updated during the calibration process. The hydraulic conductivity values in the study area range from 0.2 to 25.6 m/day with an average value of 6 m/day. The specific yield values range from 0.02 to 0.25.

The WetSpass model is essentially limited in how it handles groundwater flow due to its lumped structure. The primary inputs to the groundwater model, scattered groundwater recharges, are difficult for MODFLOW to identify (Aslam et al. 2022). MODFLOW and WetSpass ensure data interchange until the rates of recharging and hydraulic heads stabilize (Aslam et al. 2022). With a range of input data, the WetSpass-M model was used to execute the first simulation. Using the estimated groundwater recharge, which is supplied into WetSpass-M, MODFLOW was utilized to model groundwater head. The groundwater level was determined and transferred to WetSpass-M, and the hydraulic conductivity of the aquifer, hydraulic head of aquifers, aquifer thickness, and WetSpass-based groundwater recharge were provided as input data in MODFLOW. The research region's spatiotemporal features were accurately portrayed. Using a coupled model of WetSpass-M and MODFLOW the present and future circumstances of the hydrology (elements of the water balance, system water budget, and groundwater level) were investigated. Integrating WetSpass-M with MODFLOW is necessary to comprehend the groundwater flow characteristics that will be utilized as a decision support system (DSS) for water resource management (Aslam et al. 2022; Bezabih & Alemayehu 2022).

The model is calibrated by comparing the observed and simulated groundwater recharge (base flow) and direct runoff. The model was calibrated manually using the trial and error method by changing the parameters such as interception parameter (a), evapotranspiration coefficient (α), surface runoff (LP), runoff delay factor (x), and base flow (β) which were optimized according to the goodness-of-fit between the simulated runoff and the runoff from observed discharge at station. The calibrated model shows sensitivity to a parameter in which small changes in its values affect the model significantly. The sensitivity of a given parameter was determined by fixing all calibration parameters at their calibrated values except for the selected parameter and run the model by incrementally increasing and decreasing its value by some percent from its calibrated value. Based on the results, the model displays the highest and lowest sensitivity to the model parameters. The model is sensitive to both parameters. However, the model is relatively highly sensitive to decrease of the evapotranspiration coefficient (α), and surface runoff (LP) and less sensitive to increasing of calibrated evapotranspiration coefficient (α), and surface runoff (LP).

The main outputs of the WetSpass-M model are raster maps of monthly groundwater recharge, surface runoff, and actual evapotranspiration for the period 1990–2021. Each pixel on these maps corresponds to the water budget component's magnitude (in millimeters). By adding the evaporations from open water, impermeable surface area, bare soil, interception of vegetated area, and transpiration of the vegetative cover, a WetSpass-M model determines the total real evapotranspiration per pixel (Bezabih & Alemayehu 2022). The spatial mean monthly, seasonal, and annual evapotranspiration simulated by the WetSpass-M model is presented in Tables 2 and 3 below.

In this watershed, the annual evapotranspiration was very high in afro alpine and grassland (in the northern and southern parts of the watershed).

According to the annually simulated surface runoff of the catchment (Figure 7(c)), the central, western and eastern parts of the watershed have the highest surface runoff due to the presence of rural settlements and agricultural land use types where the probability for the formation of surface runoff is high. On the other hand, the southern and northern parts have less surface runoff. This is caused afro alpine and grass coverage of the area, which hinder surface runoff formation.

The average yearly recharge of groundwater is calculated using monthly simulated data. The maximum, minimum, and mean values of annual groundwater recharge for the whole period are 581, 70.8, and 169 mm, respectively. The average recharge accounts for 11.63% of the total average annual rainfall (1,453 mm). About 70.4% of the annual groundwater recharge of the watershed occurs during the wet season (summer), and the remaining 29.6% occurs in the dry season (winter). The groundwater recharge potential was highly variable from season to season; this high variation occurs due to variation in amounts of precipitation.

The groundwater recharge was fairly distributed in the wet season throughout the watershed; this is due to the existence of precipitation in all parts of the watershed. But the degree of distribution fluctuates from location to location depending on soil type, rainfall intensity, infiltration characteristics of the soil, land use, and slope. The south and northern parts of the Chemoga watershed that receive high amounts of precipitation have higher annual and seasonal groundwater recharge. The central parts of the Chemoga watershed. On the other hand, the eastern and central parts accounted for a lower rate of annual and seasonal groundwater recharge, attributed to the presence of settlements and agricultural land with less permeable clay soil.

Groundwater recharging is one of the hydrological cycle's process elements, and this study conceptualized it as the water balance residual of the hydrological cycle. It was defined, in other words, as the water that remains in the groundwater system after surface runoff, precipitation, and actual evapotranspiration have all been taken care of (Demissie et al. 2023). Because the environmental system is complex, it is typically impossible to measure these different components of the hydrological cycle directly. Consequently, in order to enhance our comprehension of the elements of water balance and underlying physical processes, modeling of the watershed is frequently necessary (Demissie et al. 2023). The relationship between surface water and groundwater resources is significantly influenced by groundwater recharge. To assess groundwater recharge and look into how surface water and groundwater resources interact, a variety of models are available (Demissie et al. 2023). Nevertheless, implementing the majority of these models in the research area is highly costly and data-intensive. For places with shallow aquifers and limited data, the WetSpass model is determined to be the most effective and appropriate modeling technique. The study's findings may offer priceless insights into the creation of an integrated master plan for the equitable and healthy use of water resources as well as sustainable management.

The estimation of groundwater recharge with the adoption of a spatially distributed water balance model has provided a new hydrological insight for the Chemoga watershed. Because of their interdependence, it is typically preferable to comprehend and talk about the numerous elements of the water balance of the area. Hence, besides groundwater recharge, which is the focus of this research, the other two important components (AET and surface runoff) of the water balance were discussed as well. The mean evapotranspiration is 405 mm/year, which accounts for 27.87% of the annual precipitation loss in the watershed (1,453 mm). The more water is able to infiltrate into the soil and percolate to recharge the groundwater system in the wet months, particularly in July, coincides with the month of highest rainfall in the basin.

The annual surface runoff in the Chemoga watershed ranges from 400 to 1,050 mm/year. The mean value represents 879 mm/year, which accounts for 60.5% of the annual precipitation loss in the watershed (1,453 mm). As expected, there is more surface runoff in the wet months compared to the dry months. The percentage of runoff from rainfall, however, varies from month to month depending on the rainfall amount, rainfall intensity and the antecedent soil moisture. The high surface runoff in July can be attributed to the antecedent soil moisture of the preceding months. This may not have anything to do with the soil becoming totally saturated from the preceding months' rains, which reduces or eliminates infiltration. Due to the expectedly dry soil matrix throughout the dry months, especially January (2.9 mm), the little rainfall that did fall was rapidly absorbed, producing little to no surface runoff.

The maximum, minimum, and mean values of annual groundwater recharge for the whole period are 581, 70.8, and 169 mm, respectively. The average recharge accounts for 11.63% of the total average annual rainfall (1,453 mm). The Wetspass-M model predicted higher groundwater recharge in the wet months of June to September, while less recharge was in the dry months October to May. From October to May, there was low rainfall and the watershed was dry; therefore, the available rainfall first satisfies the moisture deficit in the catchment before contributing to recharge. A higher proportion of recharge (over 70.4% of the annual recharge) occurred in June, July, August, and September during periods of high rainfall. This result is comparable to many studies that estimated water balance components using the WetSpass model in Ethiopia with similar climatic conditions and reported that groundwater recharge accounted for 14.1% of rainfall in the Gerado basin (Gidafie et al. 2016), 22% of rainfall in the Lake Tana basin (Yenehun et al. 2022), 4.2% of rainfall in the Werii watershed (Gebremeskel & Kebede 2017) in the Tekeze River basin, and 9.4% of rainfall in the Upper Bilate catchment (Dereje & Nedaw 2019) in the Rift Valley Lake basin. The novelty of this study from the previous studies is that for this study, the limitation was resolved in the WetSpass-M modified version for simulating monthly variation, which uses monthly climate data and is integrated with the groundwater model to analyze the effect of recharge on hydraulic head distributions. Furthermore, the results from this study can also be used to provide inputs to a three-dimensional groundwater model to analyze the available resources against the growing water demands for future groundwater availability and management.

The primary factor influencing the spatiotemporal variance in water table elevation is recharge, which is determined by transferring the recharge values from WetSpass to MODFLOW. Thus, the results from spatiotemporal drawdowns should be correlated to the results from spatiotemporal recharge. As seen in Figure 9(b), the decline of groundwater head occurs in dry/winter period, due to the relatively low recharge in dry/winter period. The southern and northern parts of the Chemoga watershed are characterized by high groundwater recharge; this results in the groundwater head rise. On the other hand, the central parts accounted for a lower rate of annual and seasonal groundwater recharge which indicate the groundwater head decline. Recharge values are comparatively high in high-elevation terrain. Similarly, several scholars (Putthividhya & Laonamsai 2017) in the Upper Chao Phraya River Basin show that the highest water table elevation occurs in the mountainous area of the Yom River and a low groundwater level of 181 m (MSL) occurs along the Yom River and tributaries (Chung et al. 2010) in Mihocheon watershed, South Korea.

The groundwater recharge in the Chemoga watershed was evaluated by applying the Coupled WetSpass-M and MODFLOW models, which are crucial for the management and planning of water resources for sustainable development.

• i. The WetSpass-M model estimates the annual actual evapotranspiration of the basin for the period from 1990 to 2021 at 366 and 440 mm as minimum and maximum values, respectively. This represents 27.87% of the annual average precipitation (1,453 mm).

• ii. The minimum and maximum values of annual runoff in the Chemoga watershed are 400 and 1,050 mm, with a mean value of 879 mm, which accounts for 60.5% of total rainfall (1,453 mm).

• iii. Annually, simulated groundwater recharge ranges from 70.8 to 581 mm, with a mean of 169 mm, which represents 11.63% of the annual precipitation (1,453 mm).

• iv. The groundwater level in the Chemoga watershed was studied under various stress conditions (dry season, wet season, and annually).

• v. The groundwater head due to the wet/summer stress period (recharge) varied from 4 to 41 m. While in the dry/winter stress period (recharge) groundwater head varied from 3.5 to 39.8 m, and also the groundwater head due to the annual stress period (recharge) varied from 3.7 to 40 m, with a correlation coefficient of 0.99, which is an acceptable fit between the simulated and observed heads in steady-state for all stress periods (summer, winter, and annual recharge).

Future groundwater resource development plans in the valley must balance groundwater recharge with anticipated abstraction rates for domestic and agricultural water supply in order to maintain the groundwater resource's long-term viability. The study reveals that the result obtained from the coupled model would be useful to the planners and decision makers in order to ensure sustainable water resources development in the watershed. The model that was created and the results that were obtained in this work are general and could be used globally to inform sustainable water management decisions. Although the application is demonstrated for an aquifer system in Ethiopia, the proposed model and the methodology described in this work could be successfully implemented or repeated in other areas.

The research encountered several limitations, including an inadequate distribution of stations, observation wells, and the collection and restoration of data. These limitations highlight the need for improved infrastructure and data collection methods to conduct more thorough and precise research in the feature. In light of these limitations, the following recommendations are proposed:

There is 11.63% recharge which is the lowest value in the water cycle. A major portion goes in evaporation and runoff. Therefore, to harvest this excess water, it could be advantageous to practice flood control dams (artificial recharge). In low potential areas, the groundwater potential map can be used to find ideal locations for artificial recharge. Artificial recharge facilities including infiltration basins, recharge pits, and agricultural ponds are required to raise the local water table. This is helpful in one way to reduce soil erosion and in the other way to enhance more recharge to groundwater.

• Design a monitoring groundwater level network for the identification and establishment of new observation wells or the elimination of inappropriate, but available, observation wells;

• Upgrade the WetSpass-M model on a daily basis for the estimation of daily recharge and to improve groundwater management decisions;

• WetSpass-M and MODFLOW models are internally coupled and presented as a single Python program.

I am also grateful to Ethiopian Water Resources Ministry, National Meteorological Service, Geological survey, Mapping Agency and Water Works Design Enterprise (WWDE) for providing me meteorology data, and relevant documents, which helped me to carry out my research work.

T.A. contributed to conceptualization, methodology, software analyses, writing – original draft preparation, review conceptualization, supervision, writing – review and editing.

This research received no external funding.

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

The authors declare there is no conflict.