Effects of land use/land cover change on hydrological responses of a watershed in the Central Rift Valley of Ethiopia

Land use/land cover (LULC) changes, driven by agricultural expansion, urbanization, and deforestation, significantly impact hydrologic processes and water resources. These changes have mostly resulted from diverse anthropogenic practices, which are performed to satisfy the immediate needs of human beings, arising from the rapid population dynamics. Understanding the hydrological impacts of LULC change is crucial for sustainable planning and management of water resources (Alemu et al. 2022), as well as for addressing human-induced natural resource degradation.

Several studies reported that changes in LULC alter the basic water balance components, like groundwater recharge, interception, infiltration, mean flow, and evaporation of any watershed (e.g., Puno et al. 2019; Aghsaei et al. 2020; Bekele et al. 2021; Chaemiso et al. 2021; Husen et al. 2023). This alteration is mainly through interception of the hydrological cycle such as through evapotranspiration (ET), surface runoff (SURQ), and soil infiltration, thereby affecting the availability of the entire water resources (Wakjira et al. 2016). Welde & Gebremariam (2017) indicate that an increase in bare land and agricultural land caused an increase in the annual and seasonal streamflow of the watershed. Engida et al. (2021) demonstrate that the streamflow trend increased while groundwater revealed a decreasing trend due to LULC changes in the Upper Baro Basin during the past decades. Moreover, the impacts of LU/LC change on water resources of the Tana Subbasin, Ethiopia, were assessed recently by Bewuketu et al. (2023) and they showed that the annual runoff and water yield (WYLD) increased by 16.16 and 12.6%, respectively, while base flow (or groundwater flow, GWQ) and ET decreased at a rate of 9.93 and 13.49%, respectively. However, the hydrologic effect of LULC changes varies from watershed to watershed, and is also a time-variant phenomenon. This necessitates exploring the significance of its effects at regional or local scales.

Estimating water balance components using models is challenging because of the nonlinear land surface processes and complex hydrological relationships. Selecting the right model for a catchment and research purpose is crucial for reliable results. Catchment modeling is used to describe the relationship between the hydrological processes and natural and environmental factors, such as land use, climate, or other management practices. Depending on input data availability, the purpose of the study, and desired accuracy, different methods may be used to assess the effects of LULC changes on watershed hydrology and water balance components (Li et al. 2019; Tadese et al. 2020; Husen et al. 2023). Both fully distributed and semidistributed models are included in the physically based distributed model types. However, fully distributed models are difficult to employ at a large catchment because of their extensive parameterization, which makes them comparatively high data-intensive models, while the semidistributed models are practical alternative tools in most previous studies (e.g., Anand et al. 2018; Jaiswal et al. 2020; Bekele et al. 2021). Based on the free availability of the model, the possibility of spatial and temporal distribution of input data, as well as the efficiency of the model in different topographic and climatic regions, the soil and water assessment tool (SWAT) model (Arnold et al. 1998) was chosen for this study to simulate the water balance and hydrologic flux of a tropical watershed. Furthermore, several studies used the SWAT model to assess the hydrological impacts of LULC changes on water resources and hydrological fluxes in different watersheds in Ethiopia (e.g., Gessesse et al. 2015; Shawul et al. 2019; Bekele et al. 2021). For example, Wakjira et al. (2016) used the SWAT model to assess the effects of LULC changes on the hydrological process of the Gilgel Gibe catchment. Welde & Gebremariam (2017) on their part evaluated the impacts of LULC changes on the hydrological response of the Tekeze Dam watershed, in northern Ethiopia. A study by Engida et al. (2021) utilized the SWAT hydrological model to estimate the effects of LULC changes on hydrological processes, and they found that SWAT is effective in simulating water balance components under LULC changes in the Upper Baro Basin during the past decades. Moreover, the recent study by Bewuketu et al. (2023) showed that the annual runoff, WYLD, GWQ, and ET were estimated by the SWAT model. The outcomes of these studies are generally valuable and encouraging, as they indicate the usefulness of the SWAT model in simulating watershed responses to LULC changes under time-invariant conditions. However, the hydrologic effect of LULC changes varies from watershed to watershed, and it is also a time-variant phenomenon. This necessitates exploring the significance of its effects at regional or local scales.

Catchment models have a large number of parameters, which tend to change under soil type, topography, and land use types. However, in most previous studies, the SWAT model was employed to be calibrated using a given land use map situation and then used the same calibrated parameters for new land use scenarios (Welde & Gebremariam 2017; Shawul et al. 2019; Engida et al. 2021, among others). Unless the watershed experienced significant human interventions during the simulation period, this approach is not feasible, because as the LULC conditions of a certain catchment change, the calibration parameters and model results also change. Such an approach examines the effect of the changes in the hydrologic response unit (HRU) composition, rather than the effect of the LULC changes on the hydrological components (Guzha et al. 2018). Moreover, a study by Deng et al. (2016) described that for highly human-modified watersheds, assuming the parameters under one LULC scenario as a suitable approach for another LULC condition is implausible. This may lead to underestimating parameter values, and to model outputs that cannot rationally describe the existing condition (Li et al. 2019). To this end, a more suitable approach that should be employed for catchments that have been significantly disturbed by human activities is to allow model parameters to vary in time, rather than assuming a constant optimal value under different LULC scenarios. Therefore, intensive parameter sensitivity analyses need to be carried out to obtain improved model parameterization as well as reasonable model results due to LULC changes (Eckhardt et al. 2003; Breuer et al. 2006; Guzha et al. 2018).

The Modjo watershed is one of the major tributaries of the inland Awash basin of Ethiopia. The watershed is one of the notable examples where there is insufficient information regarding the LULC changes and their effects on the water resources and the hydrological flux. Little information has been found in the existing literature that assesses hydro-climatic trends (Zena et al. 2022) and the effects of changes in LULC on SURQ generation and soil erosion (Gessesse et al. 2015). These studies have greatly contributed to the current understanding regarding the hydro-climatic trends and soil erosion over the watershed. Some of the key challenges of the watershed are inadequate resource development to support the rapidly increased population; loss of soil and land because of deforestation and agricultural land expansion; improper and insufficient soil conservation practices for watershed management; and river bank and bed degradation. As a result, water resource degradation, sediment deposition in channels, production of toxic substances in lakes, and discharging of point sources from industries and factors are the most commonly experienced challenges in the catchment. Additionally, soil erosion and sedimentation resulting from anthropogenic activities mainly in land use practices and management are reported as major environmental concerns within the catchment (Gessesse et al. 2015). Against these backgrounds, therefore, this study was designed to assess the impacts of LULC changes on the hydrological components of the Modjo watershed, using multivariate statistics and hydrological modeling approaches. The specific objectives of the study are as follows: (i) to evaluate the LULC changes in the Modjo watershed during the past decades; (ii) to estimate water balance components of the watershed using three different LULC conditions (i.e., 2000, 2010, and 2020), and (iii) to evaluate the impacts of multidate LULC changes on the water resource components, and suggest future research direction. This research explores the relationship between LULC changes and water balance components of a highly human-modified watershed in Ethiopia's Central Rift Valley, arguing that LULC changes significantly influence the hydrological components in the region. The findings of the study will provide conceptual information on the changing properties of hydrological regimes for stakeholders while planning and making decisions regarding water resource management under environmental changes.

The land cover of the watershed is typically agricultural (as it covers more than 85% of the total area), with small coverages of shrubs, woodlands, and grasslands in different parts of the watershed. The agricultural land is covered by mixed farming (Zelalem & Kumar 2018). Intensively cultivated dry land in the watershed is currently used for the production of crops like teff, wheat, lentils, haricot beans, and chickpeas under rain-fed agriculture. But some root crops (carrot, potato, onion, chili pepper, etc.) and vegetables (cabbage, tomato, chili pepper, etc.) are produced using small-scale irrigation, along the Modjo River. Vertisols and Cambisols are the dominant soil types in the study area. In general, the land surface is broken with numerous underground cracks, causing a severe underground loss of runoff. The significant changes in the surface characteristics as well as changes in water resources have made this watershed an ideal place to study how anthropogenic activities, particularly the LULC changes, affected the hydrological process and the availability of water resources in the recent past.

A digital elevation model (DEM) of the study catchment, with 12.5 m by 12.5 m spatial resolution, was derived from the Alaska satellite facility website (available at https://vertex.daac.asf.alaska.edu/). The obtained DEM was pre-processed using the fill and sink procedures in GIS software using a raster editor. We used the DEM for watershed delineation and to determine topographical parameters, such as slope gradient, aspect, and stream network. To improve the hydrographic segmentation and subbasin boundary delineation of the watershed, the stream network data obtained from the Ministry of Water Resources and Electricity (https://www.mowe.gov.et) were superimposed onto the DEM.

The soil database is needed by SWAT to describe the characteristics of the Modjo catchment. In the present study, the Harmonized World Soil Database (HWSD_v1.2) of the Food and Agricultural Organization (FAO) digital soil database (FAO et al. 2011) was obtained from the Ministry of Water Resources & Energy of Ethiopia (https://www.mowe.gov.et). The GIS-based soil data analysis reveals that there are 11 soil units and 9 dominant soil types (Table 1) in the Modjo catchment, including Pellic Vertisols, Vertic Cambisols, Luvic Phaeozems, Chromic Luvisols, Eutric Fluvisols, Dystric Nitisols, and Lithic Leptosols, based on the FAO classification system. Based on the data, Vertisols (42.02%) and Cambisols (33.56%) soil types cover a large area of the catchment, followed by Phaeozems (12.62%) and Luvisols (7.91%) as listed in Table 1. To append the necessary data into the SWAT database, first, the various parameters of the model (e.g. soil erodibility (K) factor, the capacity of the available soil water,) were estimated from the soil–plant–air–water (SPAW) model (Saxton & Rawls 2006). Additionally, the percentages of clay, sand, and silt were estimated by averaging the associated SoilGrids250m digital data (available at https://soilgrids.org).

To examine the land use change, three Landsat images (i.e., 2000, 2010, and 2020) were downloaded from the United States Geological Survey (USGS) website (available at https://earthexplorer.usgs.gov). The Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI) for 2000, 2010, and 2020 imageries, respectively, were used to assess the impacts on the hydrology of the study catchment.

Daily precipitation, air temperatures (maximum and minimum), relative humidity, wind velocity, and solar irradiance are the meteorological parameters needed for daily water balance calculations and WYLD components of the SWAT model. Generally, these meteorological parameters are provided by the National Meteorological Service Agency of Ethiopia (https://www.ethiomet.gov.et). Typically, the SWAT model assigns precipitation levels to subbasins based on precipitation at the nearest weather station to the center of each simulated subbasin (Arnold et al. 1998). However, due to the limited number of stations within the study catchment and surrounding areas, the meteorological parameters from six stations (Aleltu, Bishoftu, Chefe dons, Edjere, Koka dam, and Modjo) were used in this study.

All these data sets were checked for missing records, inconsistency, and unrealistic values (outliers) before being used for further analysis. As a result, the missing values were less than 3% in all the stations; hence, linear regression was used for the rainfall data while a simple average method was used to fill the air temperature data. Inhomogeneity meteorological data were checked at the station level using the standard normal homogeneity test (Zena et al. 2022). The result indicates all data are consistent at a 5% significant level, except rainfall at the Modjo station (result not given here). The rainfall data of this station are adjusted using the graphical double mass curve (DMC) analysis method before using it in the hydrological modeling. Moreover, the solar radiation, relative humidity, and wind speed input data were available only for the Bishoftu and Modjo stations. However, the data from these stations have missing values and there are inconsistencies in data recording. Therefore, solar radiation, relative humidity, and wind speed data were generated based on the weather generator installed in the SWAT model.

The daily river discharge measured at the outlet of the Modjo watershed was obtained from the Ministry of Water and Electricity (MoWE) of Ethiopia (MoWE 2020). To screen the outliers (which are greater than a threshold value that can affect the detection of the inhomogeneity of the hydro-climatic data sets), we applied the Tuckey fence method. Furthermore, streamflow missed records were filled by using the Markov chain Monte Carlo (MCMC) method, which uses multiple imputation algorithms by using the XLSTAT statistical software package. The XLSTAT statistical software package is one of the powerful Excel data analysis add-ons that are used to analyze, customize, and share results within Microsoft Excel. Table 2 summarizes the overall materials and data used in the study, including the DEM, soil data, the Landsat images for 2000, 12010, and 2020, and meteorological data for the period from 1981 to 2020, collected from different sources.

In the present work, pre-processing procedures including color composition, geometric rectification, atmospheric correction, radiometric calibration, topographic correction, and extraction of the study area were applied before image classification. The geometric correction was performed based on control points, and topographic maps obtained from MoWE (https://www.mowe.gov.et). The satellite images were already geo-referenced and geocoded to the World Geodetic System (WGS_84_UTM_Zone_37N). The removal of sensor noise, and missing lines due to solar position, radiometric, and geometric corrections was applied to the multidate images (Hilker et al. 2012). Finally, the study area was extracted from the obtained multispectral and multidate images, following the boundaries of the Modjo watershed.

A supervised image classification method coupled with the ‘maximum likelihood classification algorithm’ was applied to classify the three Landsat images used in the present study (Ginjo et al. 2022; Mahamba et al. 2022). Identification and categorization of LULC classes were defined by using Google Earth images and GPS-based ground control points during field visits for a recent image of the year 2020. As the reference data were not available for the 2000 and 2010 images, additional data were obtained from local community interviews regarding the cropland and vegetation cover changes and conversions during the past decades. In addition, the normalized difference vegetation index (NDVI) analysis technique was considered in the study to observe the historical changes and support the correctness of the classified images (Singh 1989). In general, considering both local characteristics and the global land use classification system, the obtained Landsat images were classified into seven suitable LULC categories as presented in Table 3.

In the classification process, mixed classes, such as forest areas and croplands, which are placed at higher altitudes, grassland fields, and shrub lands as well as a forest with shrub land, were separated with the aid of DEM data. However, due to the coarse resolution of the available satellite images, the present study could not distinguish between irrigation and rain-fed crop types precisely. Therefore, the agricultural land use class was considered a rain-fed cropland. Additionally, the irrigational land of the study watershed is scattered and relatively small in area coverage to be classified as a separate LULC type. In the present study, the 2000 LULC map was used as the baseline for the current study. The results of the model simulation from the other two LULC maps can be used to evaluate hydrology processes by comparing the reference scenarios.

Here, a_ij is the number of observations in row i and column j; a_i+ are the marginal totals of row i; a_+i are the marginal totals of column i; r is the number of rows in the matrix; and N is the total number of observations.

In this study, the SWAT model (Arnold et al. 1998) was used for long-term hydrological modeling of the study watershed under multiple land use scenarios. SWAT is a physically based, semidistributed, and continuous time-step model developed to estimate the effects of land management practices on water, sediment transport, and agricultural chemical yields in large complex watersheds with varying soils, land use, and management conditions over a long period (Neitsch et al. 2011). The model has been widely used in different parts of the world and proved to be an effective tool to examine the hydrological responses to land use and climate changes (e.g., Guo et al. 2016; Polanco et al. 2017; Anand et al. 2018; Gashaw et al. 2018; Guzha et al. 2018; Li et al. 2019; Bewuketu et al. 2023). Moreover, the model a public domain with SURQ, groundwater flow (or base flow, GWQ), percolation, lateral flow (LATQ) to streams, ET, transmission losses, reach routing, nutrient and pesticide loading, and water transfer components at each HRU (Arnold et al. 1998; Neitsch et al. 2011).

Surface runoff (SURQ) occurs when the rainfall rate exceeds the infiltration rate over catchments. In the SWAT model, the Soil Conservation Service curve number (SCS-CN) method (USA_SCS 1972) and the Green-Ampt Mein Larson Infiltration (GAML) method (Mein & Larson 1973) are used to estimate SURQ for each HRU. The SCS-CN is a widely used empirical model. Inputs for this method include soil, land use, elevation, as well as daily rainfall data sets. Whereas the GAML is a physically based model that uses spatial coverages in the same way as the SCS-CN method. Inputs for the GAML method include more detailed soil information, as well as subdaily rainfall data. Inputs for the GAML include soil and elevation data, as well as rainfall records (Arnold et al. 1998; Gassman et al. 2007; Neitsch et al. 2011).

Here, S is the potential maximum retention after runoff begins for a given day (mm), S_max is the maximum possible retention value that can be achieved on any given day during drought (mm), SW is the effective soil moisture of the entire profile excluding water held at wilting point (mm), W₁ & W₂ are form or shape coefficients.

To simulate the potential evapotranspiration (PET), three methods including the Penman–Monteith method (Monteith 1965), the Priestley–Tailor method (Priestley & Taylor 1972), and Hargreaves method (Hargreaves & Samni 1985) are available in the SWAT model. Based on data availability, this study relies on the Hargreaves–Samani method to estimate the daily PET (i.e., the collection of all processes by which the water is removed by evaporation from water bodies and soil as well as transpiration from plants, crops, etc.). This method uses the minimum and maximum temperatures to calculate the PET value from the watershed. Following the calculation of the daily PET parameter, SWAT quantifies actual ET (which is composed of surface evaporation and transpiration through plant cells). In the model, evaporation is estimated first from the canopy, second from the soil surface, and finally sublimation from snow (if any), at the HRU level (Belihu et al. 2020). The study used the variable storage method for flow routing in the channel.

In which ϕ_sh,i, ϕ_sh,i-1, W_rch,sh, Q_gw, W_revap, and W_pump,sh represent the amount of water stored in the shallow aquifer on a given day (mm), the amount of water stored in the shallow aquifer on day i − 1 , the amount of recharge entering the shallow aquifer on a given day (mm), the GWQ into the main channel on a given day (mm), the amount of water moving into the soil zone in response to water deficiencies on a given day (mm), and the amount of water removed from the shallow aquifer by pumping on a given day (mm), respectively.

In which W_{perc, sl} is the amount of water percolating out of the lowest layer of the soil on a given day (mm) and W_p represents the amount of water that flows past the lower boundary of the soil attributable to a bypass flow on a given day (mm).

Here, Q_LF represents the LATQ rate discharged to the main channel on a given day (mm H₂O), 0.024 is, is the amount of water about to flow out of the saturated zone (mm), K_sat is the saturated hydraulic conductivity of the soil (mm/h), S is the slope of flow (m/m), is the total porosity of the soil layer, and is the slope length (m).

The SWAT model has a large number of parameters, many of which cannot be measured directly. Therefore, sensitivity analysis was accomplished to identify parameters with a significant impact on model output variance. In the study, three model simulations were run to identify sensitive parameters using streamflow data recorded at the Modjo River gauging station. In this study, 32 SWAT parameters were selected for sensitivity analysis, by considering 11 GWQ parameters, 5 soil parameters, 7 SURQ parameters, 3 channel parameters, 3 ET parameters, and 3 parameters related to geomorphology. Sensitivity analysis was performed by considering the lower and upper bound values derived from the SWAT model database, and other related literature in the area (Zelalem & Kumar 2018). The sensitivity analysis was performed under three land use local area plans (LAPS) (i.e., 2000, 2010, and 2020 LULC maps), concerning similar climatic, soil, and topographic conditions to obtain three parameter sets associated with each LULC scenario. The t-test (which provides a measure of sensitivity) and p-value (which determines the significance of sensitive parameters) were used to identify the sensitivity parameters (Abbaspour et al. 2015). In general, a p-value close to zero has more significance, and larger absolute t-values are more sensitive in the process (Dahiru 2008; Sime et al. 2020).

In the present study, three SWAT model setups were conducted separately for 2000, 2010, and 2020 LULC conditions to simulate the various hydrological components of the study watershed. For these conditions, three model calibration and validation processes were undertaken in the study. The first model calibration and validation were undertaken for the period from 1981 to 2000, using a 2000 LULC map, with 3 years (1981–1983) of model warm-up. The second model calibration and verification were for the period from 2001 to 2020 using the 2010 and 2020 LULC maps. As the LULC changes were different during the second period, we further divided it into two equal parts of 2001–2010, and 2011–2020, with two LULC maps, of 2010 and 2020, respectively. In the process, the first halves of the two subperiods (i.e., 2001–2005 and 2011–2015, respectively) were used for calibration, while the second halves (2006–2010 and 2016–2020, respectively) were considered for model validation.

The sequential uncertainty fitting version-2 (SUFI-2) algorithm that analyzes all sources of uncertainties (such as model structure, measured data, driving variables, and input parameters) was applied to evaluate SWAT model estimation results (Abbaspour et al. 2015; Sime et al. 2020). The SWAT calibration and uncertainty programs (SWAT-CUP) are used for sensitivity analysis, calibration, validation, and uncertainty analysis of the SWAT model simulation results (Abbaspour et al. 2015). The two commonly used factors (i.e., p-factor and r-factor) were used to evaluate model uncertainty. The degree to which all uncertainties are accounted for is calculated using the p-factor, which is part of the measured data that are encircled by the 95PPU (95% prediction uncertainties) band, using Latin hypercube sampling. Conversely, the ratio of the average width of the 95PPU band and the standard deviation of the measured variable is determined by the r-factor.

The main purpose of using the partial least squares regression (PLSR) statistical method was to establish an association between independent and dependent variables during the analysis period (1981–2020), and to give a preliminary result about the land use types responsible for the changes in water resource components. Independent variables constituted the changes in the LULC classes, while the dependent variables were the changes in the hydrological components.

Before using the PLSR model, normality checks and collinear tests were performed on the independent variables (i.e., LULC classes) with the help of the Shapiro–Wilk (S–W) normality test (Shapiro & Wilk 1965). It is found that the Shapiro–Wilk (S–W) = 0.671, at α = 0.05 level of significance, with p-value = 0.003. As the computed p-value is lower than the significance level considered in the study (i.e., α = 0.05), the independent variables (i.e., LULC types) were not from a normal distribution population. As a result, the affected independent variables were log-transformed to achieve a normally distributed population (Yan et al. 2013). The multicollinearity test on the independent variable was conducted by considering the tolerance and variance inflation factor (VIF), and the result (not shown here) revealed collinear characteristics of the LULC data (because the VIF values are greater than 10, and tolerance is near zero) indicating the presence of collinearity in the data set. Hence, for the data that show collinearity, the PLSR model is useful in exploring the impacts of individual LULC types on the hydrological components (Abdi et al. 2014; Godoy et al. 2014).

Image classification accuracy for the year 2020 is presented in Table 4. The results indicated a valid supervised image classification, with a minimum of 90.1% of the user's accuracy (UA), which is calculated by dividing the row total by the total number of valid classifications for the shrub class, and 91.5% of the producer's accuracy (PA), which indicates how often the classified map accurately describes actual ground features for forest land class classification (Table 4). Conversely, the water body showed the highest (i.e., 100%) UA and PA values during image classification in the Modjo watershed, followed by cropland both for UA value (98.8%), and PA value (98.3%). In general, the overall classification accuracy (OA) was 96.2%, and the Kappa coefficient was 0.95, which denotes an almost perfect agreement for the classified map of the year 2020 (Yang & Lo 2002).

The cropland, built-up, and barren land covers maintained increasing trends throughout the periods, while the forest, shrubs, and grassland decreased between the different periods (Table 5). In the year 2000, about 83% of the watershed area was under agricultural land, followed by shrub land (7.1%), forest (6.85%), built-up area (1.1%), and the rest covered by water, grass, and barren land use types. After a decade (i.e. in the year 2010) cropland has increased to 86.7% while forest and shrub land reduced to 5.76 and 3.62%, respectively, at the cost of an increase in cropland and built-up area (Table 5). The increase in cropland area comes primarily at the expense of shrubs and grassland covers, which would lead to the degeneration of the ecological environment. The built-up area expansion was dynamic, with percentage proportions of 76.2 and 62.5% from 2000 to 2010 and 2010 to 2020, respectively (Table 5). During the period 2000–2020, the areas under cropland and built-up area have grown by 6.77 and 186.4%, respectively. By contrast, forest plantations, shrubs, and grassland covers decreased by 41, 67.1, and 36.9% in that order (Table 5). The areas under water bodies, conversely, have evolved with less pronounced variations. The cumulative decrease in vegetation cover from 2000 to 2020 was attributable to the indiscriminate expansion of farmlands and settlement areas (Table 5), presumably associated with the highest population growth (both naturally and by migration) in the area. Moreover, land-related policies that give better incentives during the construction of road infrastructure, industries, estates, factories, etc. also contributed their shares for the built-up area expansion in the Modjo watershed.

The Modjo watershed is one of the areas in central Ethiopia dominantly engaged in agricultural activities, urbanization, and high population pressure. The boundary of the catchment is also one of the areas with relatively high industries and urbanization in the country. Therefore, continuous pressure on the available land resources and natural vegetation cover resulting from the need for expanding settlement areas, cultivated lands, and developmental activities were expected in the area. Previous studies evidenced that improper agricultural-based economy and population increase are the main causes of LULC changes in developing countries, such as Ethiopia (Wakjira et al. 2016). The result obtained regarding the increase of built-up area coincides with the findings reported by Getachew & Melesse (2012) on the Angereb watershed.

SWAT has a large number of streamflow parameters, but all these parameters are not expected to be sensitive at a certain catchment. Therefore, sensitivity analysis is required to identify the influential parameters. By assuming that different LULC changes have an effect on parameter sensitivity and model outputs, flow sensitivity analysis was carried out separately for the three reference times that considered classified LULC maps of 2000, 2010, and 2020 at the subbasin number 27. Based on the principle of global sensitivity analysis (GSA), the parameters representing the SURQ, groundwater, and channel properties are more sensitive than others in all three model experiments (Table 6). This shows an accurate estimation of these parameters is important for streamflow simulation in the study watershed. Surface runoff is extremely sensitive to parameter SCS runoff CN for soil moisture condition II (CN2), and decreasing the CN2 values results in decreased runoff and increased soil infiltration, GWQ, and recharge. In the three model setups, the CH_K2 (effective hydraulic conductivity in the main channel), and CH_N2 (Manning's n value for the main channel) were found to have high sensitivity order as depicted in Table 6.

As can be seen in the given table, for the same DEM and soil types the sensitive parameters and their values have not been the same for the three different model setups (Table 6). This is mainly attributable to the differences in the climatic and LULC conditions we considered during the simulation. As a whole, the differences among these model parameters and model setups consistently reflect the differences in hydrological processes under different LULC information on the Modjo watershed. Hence, the consideration of different LULC data in the present study during parameter sensitivity analysis is justified to account for all uncertainties, parameter values, and reasonable model outputs.

The top sensitive parameters from each of the model simulation results (Table 7) were used to calibrate the streamflow separately. To examine the impact of parameter variations on the model dynamics, the SWAT model was calibrated and validated under different LULC conditions (i.e., 2000, 2010, and 2020). In this case, both the LULC map and climate data were from the same period; hence, the three scenarios can be considered baseline scenarios (Kumar et al. 2022). The curve number (CN2) and CH_K2 (a parameter that controls the transmission loss from the stream bed) were decreased between scenarios (Table 7). Whereas the GWQMN (the parameter that is needed for return flow), TRNSRCH (which indicates the fraction of transmission losses from the main channel that enters the deep aquifer), and CH_K2 (Neitsch et al. 2011) increased during the LU2020 scenario (Table 7).

It is important to note that the sensitive parameters have clear physical meanings in the study watershed. For example, the parameter value for CN2 is different for different LULC types. As indicated by Guse et al. (2013), the decrease in the CN2 value between periods is an indication of lower runoff generation from the watershed. In other words, more water remains in the soil for subsurface water contribution in the watershed. The CH_N2 and CH_K2 parameters can be related to the surface condition of the main channel, and will decrease/increase when impervious surface areas increase/decrease in the watershed, respectively (Du et al. 2013). As a whole, the variations among the three model setups consistently reflect differences in the hydrological processes in the investigated catchment.

The calibrated model was evaluated against a different set of measured streamflow data for the period 1984–1993, 2001–2005, and 2011–2016 with LULC maps of 2000, 2010, and 2020, respectively. Table 8 summarizes the statistical performance results of each calibration and validation period. The values for Nash–Sutcliffe Coefficient (NSE) are >0.50, the coefficient of determination (R²) is >0.60, and the PBIAS values are <± 15% during the calibration phase (Table 8). This indicates a satisfactory performance index (Moriasi et al. 2015). Moreover, the RSR value is between 0.6 and 0.7 both during the calibration and validation phases, suggesting satisfactory model simulation results (Moriasi et al. 2015). The lower the RSR, the better the model simulation performance (Swain & Patra 2019). Similarly, for the validation phase, except for the PBIAS value, all the statistical performance evaluation criteria of the SWAT model reflect satisfactory performance (Moriasi et al. 2015). Based on Table 8, we can say that the obtained model performances were ‘satisfactory’ with an R² of 0.69, NSE of 0.67, as well as an RSR of 0.57 (Moriasi et al. 2015), while they were ‘not satisfactory’ with regard to the PBIAS value of 15.7% during validation of the LULC2020 model simulation (Moriasi et al. 2015).

The impacts of LULC change on water balance components of the Modjo catchment were simulated for three separate model setups of LULC of 2000, 2010, and 2020, while other inputs were held constant throughout the simulations. The details of the water balance components of different LULC changes are summarized in Table 9. The water balance components considered in the study include SURQ, groundwater GWQ, flow LATQ, WYLD, SW, and ET. Accordingly, Table 9 illustrates that SURQ and total WYLD amounts decreased with LULC from 2000 to 2020, and increased from 2010 to 2020, with an overall decline from 2000 to 2020. In contrast, GWQ, ET, LATQ, and SW showed increasing trends with land use from 2000 to 2020, and decreased from 2010 to 2020, with an overall increase from 2000 to 2020, except the ET value.

Overall, the change analysis in annual water balance components showed 6.1, 2.9, and 9.3% declines in annual SURQ, WYLD, and actual ET, respectively (Table 9). This could be attributable to the expansion of cropland as well as the reduction of shrubs and forest lands. Although barren land, built-up, and settlement areas increased in the period (see Table 5), these changes would not have significant impacts on the SURQ and total WYLD of the Modjo River watershed, because their spatial coverages were insignificant, as the impermeability of these LULC classes has not affected the water balance components in recent decades. Therefore, the decline in annual SURQ and total WYLD components was mainly due to the expansion of cultivated areas and the removal of natural vegetation. Moreover, the decreased average annual SURQ associated with expanding cultivation is probably due to the increase in the infiltration of soil. However, the SURQ and total WYLD increased in 2020 compared with 2010 (Table 9), which is attributable to the increase in the area of both settlement and barren land during this period (Liu et al. 2023). The annual SURQ and WYLD declines depicted by the model results are indeed consistent with the findings presented by Li et al. (2019) and Liu et al. (2023). A comparison of the variations in the annual ET shows that the catchment experienced an increased ET value by 70.8 mm/year (15.6%) from 2000 to 2010, and then decreased by 112.9 mm/year (21.5%) from 2010 to 2020 (Table 9). An overall decrease in actual ET of 42.1 mm (9.3%) was observed at the end of 2020 over the Modjo watershed. The spatiotemporal changes in ET are directly related to the changes in forest land area and water bodies. Due to decreasing transpiration rates related to forest cover reduction, and the shrinkage of water bodies, ET values decreased over the study area. In general, significant expansions of cropland and reduction of forest, shrub land, and water bodies were observed in the watershed, which could cause ET to decline during the past decades.

By contrast, the annual LATW, SW, and GWQ values were increased, respectively, by about 150.4, 55.9, and 19.1% from 1981 to 2020, due to changes in the LULC of the period (Table 9). The inclining trend in the annual average of GWQ, LATQ, and SW parameters may be attributed to a decrease in SURQ, an increase in percolation rate, and high soil infiltration resulting from agricultural activities, as well as to soil conservation practices in the watershed. Some regional reports indicate there is a high water resource abstraction in the area for small irrigation activities, domestic use, and industrial consumption. However, this situation may not significantly affect, particularly the groundwater resource, as the watershed is supported by conservation activities and availabilities of natural lakes such as Hora, Bishoftu, and Arsede that support the groundwater recharge. The result depicted in this study is also partially consistent with the findings of Liu et al. (2023) that report an increase in LATQ, due to the cropland expansion in a watershed.

Using LULC maps of 2000, 2010, and 2020 scenarios, the seasonal water balance components of the Modjo watershed were also estimated and compared. Accordingly, similar changes in the wet seasons (summer: June–September and spring: March–May) water balance components have been observed as that of the annual time throughout the analysis period (1981–2020). Accordingly, a decrease of 18.2, 36.1, and 10.2% in SURQ, ET, and total WYLD, and an increase in GWQ, LATQ, and SW of 90.9, 235, and 14.3%, respectively, were found during the summer season. During the small rainy (spring: March to May) season, decreases in SURQ, actual ET, and total WYLD values, and increases in GWQ, LATQ, and SW, have been observed in the study watershed (Table 10). The decrease of total WYLD over these periods is, therefore, due to the decrease of SURQ in the wet months (i.e., Jun–Sept and Mar–May). Conversely, the reduction of vegetation covers and the increase of cropland and built-up area have reduced the dry season SURQ and ET values by 0.22 and 42.4%, respectively.

The effects of LULC changes on the increase of total WYLD value are clearly seen in the dry season (October–February), which is mainly attributed to the increase of subsurface flow components in these months (Table 10). The seasonal ET values show that a decrease in ET value was observed during both the summer and spring seasons during the analysis period. The ET values during the summer season decreased by 17.4 and 22.7% in the periods from 2000 to 2010 and from 2010 to 2020, respectively (Table 10). An overall decrease in the summer season ET of 36.1% was observed at the end of 2020 over the watershed. During the spring season, the catchment followed a similar trend as that of the major rainy (Ethiopian summer) season (Table 10).

To sum up, the water balance components of the Modjo watershed showed both positive and negative responses to the long-term LULC changes during the past decades. In the watershed, most of the cropland area is plowed, hence, as a result, SURQ and WYLD decrease as plowing increases the soil infiltration rate and groundwater recharge down to the soil layers. The interactive effect of LULC change with direct human interventions (e.g. surface water abstraction and groundwater exploration) could also have unavoidable impacts on the hydrological response of the studied watershed. In the study, expansions of agricultural, built-up, and barren land areas, and reductions of shrub and forest land areas, had led to a decrease in the SURQ, WYLD, and actual ET components of the period. By contrast, the GWQ, LATQ, and the SW contents of the watershed increased under the LULC changes observed between 1981 and 2020. Although the barren and built-up areas increased during the analysis period, the percentage change in these LULC types was comparatively small. However, the reduction of forest and shrub land covers could not imply increases in SURQ in our study watershed, as there might be other direct human impacts during the study period. As a result, the combined influence of LULC changes was attributed to SURQ and reduced WYLD. In addition, compared with 2000, ET decreased in 2020 as a result of reduced plant covers.

The findings of this study are consistent with those from some of the previous studies on the Ethiopian highlands. For example, a study by Rientjes et al. (2011) found significant decreasing trends in annual discharge with forest cover loss in the Gilgel Abbay watershed. The decrease in the dry season (October–February) flow in the Blue Nile basin was associated with the conversion of vegetation cover into agriculture and grasslands over large areas of the basin (Gebremicael et al. 2013). Getachew & Melesse (2012) reported that the 46% decrease in dry average monthly flow was because of the conversion of forest to agricultural land area in the Angereb watershed. Likewise, a study by Woldeamlak & Sterk (2005) showed that minimum daily flow shows a declining trend in the Chemoga watershed, which was partially explained by changes in LULC that involved agricultural land expansion, overgrazing, and an increase in the area of Eucalyptus plantation as well as watershed degradation that involves natural resource destruction. A study elsewhere by Liu et al. (2023) found that the changes in ET parameters were mostly caused by forest expansion.

A bivariate correlation analysis was conducted between the major LULC types and the water balance components at the annual time scale (Table 11). The correlation result indicated that almost all the LULC classes have a strong association with the major hydrological components considered in the present study. For example, the percentage change of forest land cover had an insignificant positive association with the decrease in SURQ (r = 0.689, with p = 0.130), and a significant direct correlation with total WYLD (r = 0.831, with p = 0.041) and ET (r = 0.949, with p-value = 0.004) at the 5% level of significance (Table 11). In the same pattern, the shrub has a significant inverse correlation with the annual flow (r = −0.875, with p-value = 0.022) and WYLD (r = −0.945, with p-value = 0.004), and a direct relationship with actual ET (r = 0.945, with p-value = 0.015). Correspondingly, grassland has a significant direct correlation with annual SURQ (r = 0.888, with p-value = 0.018) and WYLD (r = 0.920, with p-value = 0.009). It is also observed that cropland has a significant inverse correlation with the annual SURQ (r = −0.816, with p-value = 0.046), WYLD (r = −0.947, with p-value = 0.004), and ET (r = −0.970, with p-value = 0.001) (Table 11). However, the built-up area has no significant association (p > 0.05) with any of the parameters, while barren land has a significant correlation only with ET (Table 11).

The significant correlation of almost all LULC classes with most of the hydrological components suggests that almost all LULC classes are very important for the changes in the hydrology of the Modjo watershed. Changes in total WYLD were significantly attributed to the interplay between agricultural land expansion and the decline of forests, grassland, and shrubs. Similarly, changes in SURQ can effectively be attributed to the interplay between the decline of shrubs and grassland as well as the expansion of agricultural land areas as revealed by significant correlation coefficients (p < 0.05) and large Pearson's correlation coefficients (r > 0.80). In this regard, some previous studies found that the expansion of agriculture has a greater influence on SURQ (Nie et al. 2011; Yan et al. 2013). The negative correlational effects of built-up, barren land areas on ET (Table 11) suggest that an increase in impervious surfaces decreases ET in the study area, which is consistent with the findings reported by Zheng et al. (2020). Although the GWQ, LATQ, and SW contents experienced an overall increase across the watershed (see Table 10), all are insignificantly correlated (p > 0.05) with the forest, shrub, and grassland areas. This indicates that all these water balance components were not explained by the proportion of changes in LULC types over the Modjo watershed.

The summary of the PLSR model developed separately between the selected water balance components, such as SURQ, total water yield (WYLD), GWQ, LATQ, SW, and ET, and the LULC types (Table 12). The obtained values for R² and are all above 0.50, suggesting a good prediction of the PLSR models. For SURQ, the first three components explain the model (99.9% of the variation). The addition of other factors does not considerably improve the prediction capacity of the model (Table 12). Table 13 shows the variable importance for projection (VIP) values and PLSR weights (w*) for the major components of the water balance. The first SURQ PLSR model was dominated by cropland, forest, shrub, grassland, and built-up areas on the positive side; the second model, by barren land and grassland on the positive side, while the third SURQ model was dominated by built-up areas and cropland on the positive, and negative sides in that order (Table 13a). To differentiate between the relative importance of predictors, grassland, and shrubs (VIP > 1.0) is more important, followed by cropland in this model (Table 13a). Regression coefficients also indicate that shrubs and grassland positively affect annual flow, while the effect of cropland is negative.

For the GWQ model, the minimum cross-validated root mean squared error (RMSECV) was obtained using five components, but the first three explained 99.9% of the variance (Table 12). As can be seen in Table 13b, the changes in barren land appeared to dominate on the negative side of the first and second components, cropland and forest changes dominated on the positive side of the first component, whereas grassland dominated on the negative side of the second component of the GWQ PLSR model. In addition, higher VIP values were observed with the changes in the barren land, farmland, and forest (VIP > 1), followed by shrubs, grassland, and urban lands. Barren land had the highest negative regression coefficients, while farmland and forest had the next highest positive regression coefficients (Table 13b).

Similarly, the change in the WYLD was the same as the change in the outcomes of the SURQ regression model result. For the WYLD parameter, the top three components accounted for approximately 100% of the variance (Table 12). The addition of other factors does not significantly improve the prediction capacity of the model and does not substantially enhance the percent variation explained by predictors. The first WYLD PLSR model was dominated by shrubs, agriculture, and grassland on the positive side while barren land and built-up areas are on the negative side (Table 13c). Thus, changes in the WYLD were dependent on cropland and shrub areas (on the positive side), and barren and built-up areas on the negative side. Overall, the cropland, shrubs, and grassland areas have relatively greater importance (VIP > 1.0) in the WYLD PLSR model (Table 13c). A study by Gashaw et al. (2018) also showed the greater importance of shrubs, grassland, and cropland areas in the Andassa watershed.

In the LATQ and SW PLSR models, the first three components explained 99.9% of the variances of the two models (Table 12). Adding the fourth factor does not significantly improve the predictive capacities of the model and does not increase the percent change explained by the predictors. The changes in barren land appeared to dominate on the first, second, third, and fourth components of the LATQ PLSR model, built-up area dominated the second, third, and fourth components, whereas the cropland appeared to dominate in the first and third components (Table 13d). For the LATQ model, a higher VIP value was found with changes in farmland, forest, and barren lands (VIP >1). For the SW model, changes in barren and cropland seemed to dominate in the first, second, and third components, while the latter was dominant in the fourth component. In addition, changes in the forestland appeared to be predominant in the first SW component (Table 13e). Higher VIP value was observed with the changes in the farmland and barren land areas (VIP > 1) on the positive and negative sides, respectively (Table 13e). For the LATQ and SW, the direction of similar influence of the predictors is also reported in the regression coefficients (Table 12).

Finally, as has been presented in Table 12, the minimum RMSECV was obtained with three components, accounting for 99.9% of the total variance for the ET model. Among the LULC types, cropland appeared to be dominant on the first, second, and third components, grasslands on the first and second components, while forest and shrub land were dominant on the first component (Table 13f). A higher VIP value (VIP > 1.0) was found for the variation in cropland, forest, and shrubs, which indicated that these LULC types had a more significant role in the ET PLSR model.

Increased human intervention has resulted in rapid transitions in LULC, and adversely affected the hydrological cycle and water resources in the Modjo catchment. Therefore, it is very important to analyze the changes in LULC, and assess their impacts on the water balance components for the issue of sustainability. The objective of this study was to analyze the LULC changes and evaluate their impacts on the water resource components during the period 1981–2020. To achieve this objective, the SWAT model was used to generate the water balance components using the multidate Landsat images of 2000, 2010, and 2020. Then, the effects of each LULC class on the hydrologic components were examined using the PLSR model. Overall, the findings of the study are summarized as follows:

• The detailed analysis of remote sensing-based land cover mapping of the watershed for the years 2000, 2010, and 2020 shows that the shrub, forest, water body, and grassland areas were diminished by 67.12, 40.98, 19.71, and 36.88% of the total area, respectively, while built-up, cropland, and barren land areas increased by about 186.36, 6.77, and 30.55% in that order. Taking the internal conversion of various LULC classes into account, an overall trend from 2000 to 2020 was a change from shrubs, grassland, and forest land areas to crops and built-up areas.

• In the calibration and validation processes, good agreements have been identified in all scenarios. Therefore, the SWAT model was found to be useful in identifying the effect of LULC changes on the hydrology of the investigated watershed. The modeling result showed that the water balance components experienced positive and negative responses to LULC changes during the period 1981–2020. An overall decrease of 16.1, 2.9, and 9.3% in SURQ, WYLD, and ET, and an increase of 19.1, 150.4, and 55.9% in the GWQ, LATQ, and, SW content, respectively, were found resulting from LULC changes between 1981 and 2020.

• The major LULC types that affected the SURQ and WYLD components were cropland expansion and the shrinkage of shrub and grassland areas (with VIP values greater than 1, and larger regression coefficients). Conversely, expansions of cropland areas (with VIP > 1.0, w* > 0.30, and larger regression coefficients) and reduction of forest cover (with VIP > 1.0, w* > 0.30, and larger regression coefficients) were the major contributors to actual ET during the analysis period. In addition, the GWQ of the catchment was affected more by the barren land, cropland, grassland, shrubs, and forest land covers (with a VIP value > 1.0), and less by the built-up area expansion and water bodies (with VIP < 1.0).

The study provides insight into the hydrological response to current LULC changes in a highly human-dominated watershed. When analyzing the mean values of SURQ, ET, and WYLD, a decreasing trend was found with decreasing changes in the percentage of shrubs, grassland, and forest land cover, while increasing percentage proportion in cropland, built-up, and barren land areas. Thus, we assert that the impact of LULC changes on hydrological variables cannot be detected only by considering the LULC changes because of the more complex interacting nature of the watershed's specific characteristics such as climate change. Considering the growing water demand and agricultural activities, therefore, the study suggests further research on how the combined climate and LULC changes affect the hydrological responses of the watershed. The results of this study can be beneficial to the authorities, decision-makers, and planners for sustainable water resource management in the face of environmental changes

The authors are thankful to the Ministry of Water and Energy (MoWE) and the National Meteorology Agency (NMA) of Ethiopia for providing the raw data used in the study. The authors would like to acknowledge Jimma and Assosa Universities for their support with educational materials. We also thank the anonymous reviewers for their invaluable comments and constructive suggestions that improved the quality of this manuscript.

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

The authors declare there is no conflict.