Watershed hydrological responses to land use land cover changes at Bilata watershed, Rift Valley Basin, southern Ethiopia Open Access

LULC has pronounced effects on ecological, environmental, and hydrologic systems and processes (Liyew et al. 2019). Urbanization, agriculture expansion, mining, and other land uses have exponentially changed the earth's environment. Most of these anthropogenic changes have a profound impact on the water cycle and, finally, on the availability and quality of water (WoldeYohannes et al. 2018). To reduce these impacts, environmental protection, sustainable development, and a study of the effects of land use change are needed. It is necessary to forecast the effect of land use changes, agricultural practices, afforestation and deforestation, and related activities on hydrology (Liu et al. 2022).

LULC have pronounced effects on ecological, environmental, and hydrologic systems and processes (Ndulue et al. 2015). Land cover change, combined with the expansion of agriculture, cattle raising, and urbanization, could have a profound impact on the hydrological process of a watershed (Terefe et al. 2021). Land use/land cover changes occur in the country, Ethiopia, as a whole and in the study area in particular due to the increasing population, which has almost doubled (Mathewos et al. 2019; Sulamo et al. 2021). The effect of land use land cover, especially in the mountainous part of the catchment has a significant effect on the rainfall–runoff response of the catchment (WoldeYohannes et al. 2018; Sulamo et al. 2021; Tamire et al. 2022).

Even though today, the spatial and temporal change of LULC and its impacts on hydrological processes are given attention in the country (Choto & Fetene 2019; Engida et al. 2021; Negese 2021; Tavares et al. 2021; Gebre et al. 2022) still it needs further study, particularly in the study area to protect reservoirs, lakes, and rivers. Land use and land cover change are important characteristics of the runoff process that affect infiltration, interception, erosion, and evapotranspiration. These changes have caused severe stress on forested water resources due to rapid development in the catchment, which is subjected to land use and cover changes causing the area to form impervious surfaces. Hydrological modeling and water resource management are highly related to the processes of the hydrologic cycle. This cycle can be affected by land use and land cover change. This LU/LC change is mainly caused by anthropogenic activities. Thus, it is essential to analyze the possible impacts of these changes on hydrological characteristics for environmental protection.

To evaluate LULC change impacts on the river flow characterization of the rainfall–runoff relationships of the watershed need both spatial and temporal data to observe the effect, the chosen model should incorporate both temporal and spatial parameters. Therefore, the physically based, semi-distributed, computationally efficient and public domain model, SWAT was used. Many scholars in the country have tried to quantify the impact of land use change on the hydrological process using the SWAT model (Belihu et al. 2020; Aragaw et al. 2021; Kuma et al. 2021; Ayele et al. 2023). Therefore, this research sought to assess the changes in land use and land cover (LULC) effects on the hydrological characteristics of the Bilata watershed. A deep understanding of the LULC change is important for the conservation of previous LULC changes, predicting future changes for adequate and sustainable streamflow and decreased sediment depositions, and developing sustainable land resource management techniques intended to protect essential landscape functioning. This knowledge gives crucial information for producing evidence that helps local communities, actors functioning within a given environment, spatial planners, and decision-makers to design the best policies and strategies. An understanding of the extent of the hydrological process in the watershed and knowing the changes in streamflow are necessary for an effective reservoir and basin management program. Also, understanding the hydrological process associated with LULC change is vital for researchers for further study and for non-governmental organizations who are interested in developing any measurable structures to protect the fast-growing degradation of the land and in improving human wellbeing.

Bilata watershed is characterized by diverse topographic conditions. The upper part of the watershed is characterized by mountainous and highly separated terrain with steep slopes and the downstream part is a gentle slope with elevation ranges from 3,318 m a.m.s.l. in the mountains area to 1,116 m a.m.s.l. in the lowland. The general topography of the catchments is undulating hills and flats. It gradually decreases in elevation to the south (Kuma et al. 2021).

According to the Ethiopian agro-climatic zone classification, the climate of the study area ranges from hot temperate (Kola) to cool temperate (Wurch) in the mountains and escarpment.

Rainfall distribution over the area is Bimodal, characterized by a short rainy season (Belg) that occurs between March and May and a long rainy season (kiremt) that occurs from June to July, with a dry season from December to February. The watershed area has a tropical climate with pronounced aridity in the southern part and a warm-temperate rainy climate in the central and northern highlands. It receives about 854–1,039 mm of rainfall annually (Sulamo et al. 2021). The other factor which characterizes the watershed is the temperature. The average value of the temperature in the watershed has been recorded as 20.6 °C in the summer season and 12.57 °C has been recorded in winter (dry season). The mean annual maximum and minimum temperatures are between 22.8 and 30.4 °C and between 10.7 and 16.8 °C (Kuma et al. 2022).

For this study, GIS version 10.4 and SWAT 2012 version were used. The image classification was performed by the Earth Resources Data Analysis System (ERDAS) Imagine 2015 software package to generate the land use and land cover maps so as to detect the change in different categories of land use and land cover. Other software and materials like SWAT-CUP, Google Earth and GPS were used for model calibration and validation, homogeneity test, hydro-meteorological data analysis, ground verification, and collection ground control point's (GCP's), respectively.

Different data were collected from various sources and field observations. It includes primary data such as ground control points and secondary data including recorded meteorological data, hydrological data, digital elevation model (DEM), soil map, and land use land cover. After collecting the required data, the quality of the data was checked mainly through preliminary checking, plotting and removal of errors in order to ensure the quality of the data for further investigation.

The data which are considered as spatial data include DEM, land use land cover, and soil map. A 20 m × 20 m DEM were obtained from the USGS website, https://earthexplorer.usgs.gov.

Soil properties are one of the major input data required by the SWAT model of the watershed. The soil data covering the study area were obtained from the Digital Soil Map of the World (DSMW) on the FAO website, in Environmental System Research Institute (ESRI) shape file format.

Land use is the most important factors that affect runoff, evapotranspiration, and surface erosion in a watershed. These data are required by the SWAT model to describe the hydrological response units (HRUs) of the watershed together with slope and soil map. The land use land cover map of the study area for different years was obtained from the website.

Landsat imageries were used for change detection of land use and land cover in the Bilata watershed from 2000 to 2020. For this period, three Landsat images from 2000, 2010, and 2020 were downloaded from the United States Geological Survey (USGS) website (earthexplorer.usgs.gov) GeoTIFF file format.

Landsat 7 Enhanced thematic mappers (ETM), Landsat 8 operational land imager (OLI), and thermal infrared sensor (TIRS) were used for data acquisition. To avoid a seasonal variation in vegetation pattern and distribution throughout a year, the selection of dates of the acquired data was made in the same season of the acquired years. The images used in this study area were geo-referenced to a Universal Transverse Mercator (UTM) projection using datum WGS (World Geodetic System) 84 zone 37 N prior to image classification. So as to view and distinguish the surface features clearly, all the input satellite images were composed using the RGB color composition.

In order to use remote sensing data for land use land cover analysis, it needs to be preprocessed. This study utilized both radiometric and geometric processing techniques. Radiometric processing is necessary to minimize the effects that hinder the measurement of surface reflectance through remote sensing. Geometric correction, on the other hand, involves distorting an image to fit a specific map projection. However, this process results in a resampling of the image pixels which can lead to some degradation of the original data.

An image classification process is required for classifying all pixels in the form of digital images into several land cover classes, or ‘themes’. Image classification is typically carried out using ERDAS Imagine software. ERDAS Imagine is a software package designed for image processing that enables users to process geospatial and other types of imagery, as well as vector data. It can handle hyperspectral imagery and LiDAR data from various sensors, and includes a 3D viewing module (Virtual GIS) and a vector module for modeling. The native programming language used in ERDAS is called EML (ERDAS Macro Language). ERDAS is integrated with other applications such as Google Earth and remote sensing applications, and the imagery storage format can be read in many other applications. In this study, land use land cover was classified using a supervised classification system in ERDAS Imagine software. In supervised image classification a set of training points is developed as a pre-processing step for speeding up image retrieval in large databases and improving accuracy, or for performing automatic image annotation. The integration of Google Earth and ERDAS Imagine software was used, and the results were verified with survey data obtained from the field.

Ground truth or field survey is done in order to observe and collect information about the actual condition on the ground at a site to determine the relationship between remotely sensed data and the object to be observed. It is recommended to have a ground truth at the same time of data acquisition, or at least within the time that the environmental condition does not change.

Accuracy assessment is an important step in the image classification process. The objective of this process is to quantitatively determine how effectively pixels were grouped into the correct feature classes in the area under investigation. It is a process used to estimate the accuracy of image classification by comparing the classified map with a reference map. The most widely used classification accuracy is in the form of an error matrix which can be used to derive a series of descriptive and analytical statistics. The columns of the matrix depict the number of pixels per class for the reference data, and the rows show the number of pixels per class for the classified image. From this error matrix, the number of accuracy measures such as overall accuracy, mapping accuracy, producer's accuracy, user's accuracy, and Kappa coefficient of the agreement are estimated. The overall accuracy is used to indicate the accuracy of the whole classification (i.e., the number of correctly classified pixels divided by the total number of pixels in the error matrix).

User's accuracy is regarded as the probability that a pixel classified on the map actually represents that class on the ground or reference data, whereas producer's accuracy tells us how well a certain area can be classified. It is obtained by dividing the number of correctly classified pixels in the category by the total number of pixels of the category in the reference data. The producer's accuracy is also known as an omission error, which is the probability of reference pixels being classified correctly. It gives only the proportion of correctly classified pixels. The accuracy assessment of the classified map was done through the comparison of the classified image and the sampling points from the field survey and existing land cover maps. In this study, a total of 80 testing sample points for each land use class were selected randomly for the years 2000, 2010, and 2020.

The SWAT model requires weather data, hydrological and spatial data like watershed characteristics such as topography, land use/land cover, and soil.

SWAT input parameters are process-based and must be held within a realistic uncertainty range. The first step in the calibration and validation process in SWAT is the determination of the sensitivity parameters for a given watershed or sub-watershed. Sensitivity analysis is the way of determining the rate of change in model results concerning changes in model inputs (parameters) (Arnold et al. 2012).

The SWAT hydrological model was calibrated and validated using a static land use (SLU) (2000) setup. Calibration is an effort to better parameterize a model to a given set of local conditions, thereby reducing the prediction uncertainty (Abbaspour et al. 2017). The prediction of the uncertainty of SWAT model calibration and validation results will be analyzed by the SWAT calibration uncertainties program known as SWAT-CUP. It can integrate various calibration/uncertainty analysis procedures for SWAT in one user interface. It is a public domain program that links sequential uncertainty fitting (SUFI-2) (Abbaspour et al. 2007), particle swarm optimization (PSO) (Beven & Binley 1992), and parameter solution (ParaSol) (Beven & Binley 1992), and Marko Chain Monte Carlo (MCMC). For this study, SUFI-2 was used for the calibration of the model because its algorithms account for several sources of uncertainties.

Model validation is the process of describing that a given site-specific model is capable of making satisfactory simulations. It is the comparison of model results with an independent data set without further adjustment of the model parameters. Validation embraces running a model using parameters that were estimated during the calibration process and comparing the predictions to observed data not used in the calibration process.

Herein the performance of the model has been checked by statistical tests that can be used to judge the SWAT model. According to Coffey et al. (2004), the statistical tests that are used to judge the SWAT model include R², NSE, root mean square error (RMSE), nonparametric tests, t-test, objective functions, autocorrelation, and cross-correlation. For this study, the most widely used tests of Nash–Sutcliffe efficiency (NSE), per cent bias (PBIAS), and coefficient of determination (R²) are used which is recommended by Moriasi et al. (2007) (Table 1).

The accuracy assessment was conducted for all classified imageries using a standard method. The user's, producer's, overall accuracy, and the kappa coefficients were computed. In general, all the classifications met the recommended minimum 85% accuracy (Anderson et al. 1976). The result of the classification reveals that the overall accuracy was 90.97% and the Kappa coefficient of 88.94 for 2,000 land use/land cover classification, which indicates that the land cover classes were correctly classified.

Similarly, the overall accuracy and kappa coefficient for the years 2010 and 2020 are 91.8 and 86.8%, 90.52 and 86.21%, respectively (Table 2). A similar study by Mathewos et al. (2019) on one of the Bilata sub-watersheds on land use land cover dynamics shows that the overall accuracy and Kappa coefficient are greater than 85% which is almost identical to this study. The spectral property similarities among other land cover classes cause the class accuracy values the lowest which indicates misclassification.

Land use land cover maps were made from the data that have been generated from Landsat imageries of 2000, 2010 and 2020 years to describe land use land cover change patterns with time series and analyze its impacts on runoff and sediment yield of the Bilata watershed.

The satellite image of 2010 was classified similar with 2000 LULC class which includes settlement 0.631% (3,547.91 ha), water bodies 0.125% (701.13 ha), regularly flooded 0.1% (563.28 ha), crop land 58.73% (330,383.33 ha), grass land 11.40% (64,124.783 ha), forest 1.11% (6,268.63 ha), shrub land 6.5% (36,585.58 ha), and bare land 21.39% (120,314.97 ha) (Figure 4).

From 2020 Landsat image, LULC of Bilata watershed was classified as settlement 0.71% (4,004.49 ha), water bodies 0.11% (616.43 ha), regularly flooded 0.11% (637.01 ha), crop land 69.89% (393,110.9 ha), grass land 4.86% (27,311.83 ha), forest 1.08% (6,084.53 ha) shrub land 6.7% (37,702.65 ha), and bare land 16.54% (93,015.62 ha). In this time period, the grassland was decreased by approximately 50% from 2010 (Table 3), mainly it may be converted to arable land and settlements. A study on land use land cover change and their effects in the Abaya-Chamo landscape including the Bilata watershed by WoldeYohannes et al. (2018) reveals that between 1985 and 2010, natural grassland, shrubland, and heterogeneous agricultural areas were mainly converted to arable land.

The individual class areas and change statistics for the three periods are summarized in Table 3. A similar study in the previous period shows that the agricultural land settlement increased and the forest land was decreased. Recently, a study by Sulamo et al. (2021) in the Bilata watershed reports that the settlement and cropland was increased by 14.38 and 11.73% from 1989 to 2015, respectively. On the contrary, the forest land decreased by 4.3% from 1989 to 2015.

The results reported here are in line with the recent work of Kenea et al. (2021) in Fincha watershed to evaluate the response of hydrological components under LULC change. The result reports that there is a significant reduction of forest land (−32.3%) and shrubland (−9.6%) between 1994 and 2018, mostly due to the creation of new agricultural zones (+27.2%) and settlement (+3.7%) to answer to an increasing population that needs new means of subsistence.

The sensitivity of the model parameter for the given watershed was estimated based on the t-stat value and p-value. We can say the model parameter is sensitive when the t-stat value of the parameter is the higher and minimum value of the p-value.

The result of the global sensitivity analysis reveals that, from 20 model parameters, 16 parameters (Table 4) were found to be sensitive under the category of high to low sensitivity based on the value of p-value and t-stat. The higher the absolute value of the t-stat and minimum p-value, the more the sensitivity of the model parameter. From the 16 parameters, curve number (R_CN2), available water capacity of the soil layer (R_SOL_AWC), depth from soil surface to bottom of layer are the top five sensitive parameters (R__SOL_Z), baseflow alpha factor (days) (V__ALPHA_BF.gw), and threshold depth of water in the shallow aquifer (mm) (V__REVAPMN.gw) are the top five sensitive parameters. The result of the sensitivity analysis matches with the study by Beza et al. (2023) in the Awash basin at Jewuha watershed to model and assess the surface water potential using a combined SWAT model and spatial proximity regionalization technique. The result of the study reveals that curve number (R_CN2), saturated hydraulic conductivity (R_SOL_K), groundwater delay (days) (V_GW_DELAY), Manning's ‘n’ values for overland flow (R_OV_N), and available water capacity of the soil layer (R_SOL_AWC) are the most five sensitive parameters.

Another research conducted by Setegn et al. (2008) studied hydrological modeling in the Lake Tana basin watershed, Ethiopia tested the performance and feasibility of the SWAT model for the prediction of streamflow. During sensitivity analysis ESCO, CN2, ALPHA_BF [days], REVAPMN.gw [mm H₂O], [days], SOL_AWC[mm of H₂O], GW_REVAP, CH_K [mm/h], and GWQMN.gw [mm H₂O] are the most sensitive parameters which are similar to our study due to the physical similarity of the watershed except ESCO and ALPHA_BF due to the SWAT model performing better to simulate the groundwater in Lake Tana basin watershed than in Bilata watershed and also the study considers the impact of sub-basin discretization which resulted in a better representation of the hydrological processes and produced streamflow yield which had a better model efficiency in comparison to those who are not considered in basin discretization.

The SWAT (Soil and Water Assessment Tool) model is a widely used hydrological model that simulates the impacts of and management practices on water resources at the watershed scale. Calibration of the SWAT model is a crucial step in ensuring the model accurately represents the hydrological processes and provides reliable predictions. It involves adjusting the model parameters to minimize the difference between observed and simulated values of various hydrological variables, such as streamflow, nutrient concentrations, sediment yield, and groundwater levels (Abbaspour et al. 2017). The goal is to find the optimal set of parameter values that best represent the observed data (Abbaspour 2005; Beza et al. 2023).

The result of the model performance evaluation for the Bilata watershed at Bilata Tena gauge station showed good results of p-factor 0.9 and 0.6, and r-factor 1.31 and 1.5 for streamflow and sediment calibration, respectively (Table 5). Similarly, in the case of validation, the p-factor was 0.85 and 0.55 and r-factor was 0.74 and 1.5, respectively. The p-factor of 0.9, which means the model captured 90% of observed streamflow in the 95PPU band. A smaller p-factor and a larger r-factor could be acceptable range (Beza et al. 2023; Dinku & Kebede 2023).

The effect of land use land cover on the hydrology of the watershed was analyzed by keeping the input parameters constant and changing LULC data. The change in the annual average surface runoff, evaporation and transpiration, lateral flow, return flow, percolation to the shallow aquifer, and recharge to the deep aquifer was presented in Table 6. As shown in Table 6, the reduction in return flow and percolation to a shallow aquifer between 2000 and 2010 could be due to factors such as changes in land use, increased urbanization, or alterations in precipitation patterns leading to less water reaching the aquifer. On the other hand, the increase between 2010 and 2020 might be attributed to improved water management practices, increased recharge efforts, or changes in climate patterns resulting in more water infiltrating into the aquifer during that period.

The streamflow was increased in the wet season, the expansion of agricultural area throughout the study period had a direct influence on streamflow and a decrease of forest/bush land to some extent increased a generating runoff during rainfall. As shown in Figure 9, the mean monthly streamflow increased by 4 and 42 m³/s in the wet season and the mean monthly streamflow decreased by 1 and 20 m³/s in the dry season in the period between 2000 and 2010, 2010, and 2020, respectively.

LULC change has substantial impacts on the hydrological process within a watershed. This study evaluates the effect of LULC change on the hydrological process of Bilata watershed. Alterations in land cover, significantly influence the hydrological processes. Landsat 7 and 8 in the years 2000, 2010, and 2020 are used for land use land cover classification assessment using a supervised classification system. The overall accuracy of the classified image was found to be above the recommended minimum of 85% for all years. The kappa coefficient was also high, indicating a good agreement between the classified image and reference data. The result of the classified land use/land cover map showed that the LULC of the study area includes eight major classes such as settlement, shrubland, forest, bare land, water bodies, regularly flooded/vegetative aquatic areas, grassland, and agriculture. It is visible that the Bilata watershed has been subjected to significant land use/cover change from 2000 to 2020. It reveals that the maximum change was observed in shrubland, grassland, and cropland class. The shrub and grassland decreased by 57,774.93 and 50,027.79 ha from the year 2000 to 2010. However, crop land increases by 46,891.04 ha from 2000 up to 2010. This is due to the increase in population growth that causes the increase in demand for cultivation land for different agricultural production. This might be because of the deforestation activities that have taken place for the purpose of agriculture and urban expansion. The remaining class type bare land, forest, water bodies and regularly flooded decreased by 34,731.45, 4,443.71, 23.19, and 363.99 ha, respectively. On the contrary, the settlement increased by 416.7 ha between 2000 and 2010. The watershed was modeling using the SWAT model. Analysis of sensitive model parameters reveals that, from 20 model parameters, 16 parameters were found to be sensitive under the category of high to low sensitivity. The calibration and validation of the model was held between 1990 to 2005 and 2006 to 2015 for streamflow. The performance of the model shows very good results which have R² 0.9, NSE 0.7, and PBIAS −10.8 for calibration and R² 0.8, NSE 0.6, and PBIAS −7.4 during validation. Alterations in land cover can impact surface runoff patterns and infiltration capacity of the soil, leading to changes in streamflow regimes. Due to intensive agricultural expansion and settlement, the surface runoff increases from 172 to 259.55 and 265.05 mm in 2000, 2010 and 2020, respectively. Similarly, evaporation and transpiration increased from 174.4 mm in 2000 to 185.1 mm in 2010, and 199.3 mm in 2020. Lateral flow has reduced from 38.41 mm in 2000 to 22.4 mm in 2010, and to 17.43 mm in 2020. Also, the average streamflow increased from 132, 134, and 150 m³/s between 2000, 2010, and 2020. The rapid expansion of agricultural land, deforestation and high population growth in the area resulted in a high rate of hydrological change in the watershed therefore a landscape intervention is needed to prevent the rate of hydrological change that risks the capacity and quality of Lake Abaya. Therefore, implementing conservation practices, such as terracing, contour ploughing, or reforestation, can help mitigate the impacts of LULC change on streamflow.

This work was not supported financially by any institute or organization.

A.M., M.B., E.T., and M.C. conceptualized the whole article; A.M. and M.B. developed the methodology; M.B. arranged the software; M.B. validated the data; A.M. and M.B. rendered support in formal analysis; A.M. M.B., E.T., and M.C. investigated the data; A.M. M.B., E.T., and M.C. arranged the resources; A.M. and M.B. rendered support in data curation; M.B. rendered support in writing the original draft; A.M., E.T., and M.C. wrote the review and edited the article; and A.M., M.B., E.T., and M.C. visualized the data. All authors read and approved the final manuscript.

The submitted work is original and has not been published elsewhere in any form or language (partially or in full).

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

The authors declare there is no conflict.