Introduction

Climate and LULC changes are increasingly recognized as significant threats to the sustainability of water resources, particularly within African watersheds1. These changes can lead to altered rainfall patterns, increased evapotranspiration, and changes in surface and groundwater availability, thereby affecting human populations and ecosystems. Watersheds in sub-Saharan Africa have also received limited research attention in terms of LULC changes, particularly in response to increasing anthropogenic pressure and disturbances2. Climate projections indicate that these anthropogenic drivers will significantly alter hydrological patterns by 2050, further exacerbating water resource challenges3. The combined effects of climate change and LULC change are expected to substantially modify seasonal streamflow patterns and overall water availability, presenting significant challenges in achieving Sustainable Development Goals (SDGs) for water resources in the foreseeable future4,5,6,7. Water yield, a crucial indicator of water resource availability within watersheds, is significantly influenced by climate and LULC changes8,9,10. It is essential to understand how these factors impact water yield, as their effects differ based on the unique geographic, climatic, and socioeconomic characteristics of each watershed11. However, there are still gaps in our understanding of the diverse impacts across various watersheds, especially concerning how watershed-specific factors influence the effects of these environmental changes12,13,14,15. In addition, understanding the long-term influence of LULC changes on water resources is essential for evaluating the current and future impacts on hydrological systems.

Moreover, numerous studies have examined the effects of climate and land use at different scales, and there is a significant gap in the literature regarding the simultaneous analysis of similarities and differences in their influence on water yields across multiple watersheds. Previous studies have noted that geographic and climatic characteristics play a vital role in shaping water yield responses to climate and land use changes16,17,18,19,20. Nevertheless, the current literature lacks a comprehensive analysis that concurrently examines the effects of climate change and LULC across multiple watersheds with varied characteristics21,22. This gap may be due to the limitations of models such as InVEST, which, although effective in calculating watershed water supply, often fails to sufficiently capture factors specific to each watershed. Furthermore, inconsistencies in data and standard concerns may hinder accurate water yield assessments.

Water yield is an important measure of ecosystem health and sustainable development23,24. To better understand watershed dynamics, it is essential to analyze the spatial and temporal variations in water yield. The Gilgel Gibe watershed is particularly significant, not only for its role in powering Ethiopia’s Gibe I, II, and III hydroelectric systems, but also for serving as a key water source for Jimma Town. This watershed faces significant environmental pressure owing to human activity, which affects its ability to provide essential ecosystem services. The main factors contributing to these pressures were changes in land use and climate. Prior studies using the InVEST model have assessed the watershed water production capacity and explored the effects of climate and land-use changes on water yield25,26,27. Despite extensive research, there is a significant gap in the literature regarding the simultaneous analysis of similarities and differences in the influence of climate and LUC on water yields across multiple watersheds. This inconsistency can be attributed to several factors. First, while the InVEST model is commonly utilized for estimating watershed water yield, it lacks a tailored approach for deriving watershed-specific parameters for various watersheds. Second, the collection and processing of extensive data from watersheds is essential because ensuring data quality directly affects the accuracy of water yield assessments. Various data types often require transformation into a uniform format that aligns with modeling requirements.

This study uniquely investigated the combined effects of climate and LULC change on water yield in the Gilgel Gibe watershed, a region of substantial ecological, agricultural, and ecosystem value. The combined effects of these factors have significant impacts on regional water resource availability and biodiversity conservation, as well as vital ramifications for ecological stability and sustainable management. How do changes in LULC, specifically deforestation and agricultural expansion, affect the water yield trends over time in the region? To what extent do climate variability factors, such as changes in rainfall and temperature patterns, amplify the impacts of LULC changes on water yield? What are the integrated impacts of altered water yield on regional water resources and biodiversity conservation, and what are the broader implications for ecological sustainability? This watershed plays a critical role in downstream irrigation, municipal water supply, and biodiversity conservation, making it an ideal case study to explore how these changes can jeopardize water resource sustainability. By utilizing integrated modeling frameworks and advanced approaches, including machine learning and remote sensing, this study aims to bridge critical knowledge gaps regarding watershed hydrology. The diverse land cover of the Gilgel Gibe watershed, encompassing forests, shrublands, grasslands, agricultural areas, and wetlands, is essential for examining the varying effects of LULC changes on hydrological processes.

Additionally, the diverse geographic and climatic conditions within the watershed, including fluctuating rainfall patterns and elevation gradients, further enhance its suitability as a model for studying the spatial variability in water yield and availability. The presence of hydropower projects in the area also ensures the availability of substantial hydrological data, including river flow, reservoir levels, and meteorological data, which are essential for this type of analysis. Recent studies by28,29 highlighted the significance of integrating these advanced methodologies to reveal the spatiotemporal patterns of LULC change impacts on the water availability. Moreover, integrating climate projections into these models enables the exploration of potential future scenarios, providing insights into the risks of water scarcity and the socioeconomic impacts of diminishing water availability30,31. Water yield declines in the Gilgel Gibe watershed might have a substantial impact on hydropower production, which is essential to Ethiopia’s energy supply. Reduced base flow and higher sedimentation are expected to reduce reservoir storage capacity and turbine efficiency, potentially resulting in significant revenue losses of millions of dollars each year. For example, in May 2019, the Ethiopian government and the Ministry of Water, Irrigation, and Electricity (MoWIE) estimated a 476-megawatt power shortage caused by climate change impacts on the Omo-Gibe dams in the Omo-Gibe River basin32. On the agricultural front, reduced water availability may limit irrigation capacity, resulting in significant yield decreases for staple crops, such as maize and teff, which are critical for local people’s survival. The Gibe River watershed includes three cascading dams, Gilgel Gibe I, Gilgel Gibe II, and Gilgel Gibe III, that support irrigation projects on 175,000 hectares of agricultural land33.

By integrating climate and land use scenarios, this study aims to provide practical insights for policymakers and stakeholders to facilitate more sustainable water management in response to changing environmental conditions. The primary objectives of this study were as follows: (1) to evaluate temporal variations in surface water yield under climate variability, (2) to analyze the spatial distribution of LULC change and its direct impact on water availability, and (3) to develop predictive models that integrate climate and land use projections to forecast future water yield fluctuations. Furthermore, the research findings provide essential hydrological insights, presenting transferable methods and lessons that may be applied to model-driven watershed studies, both locally and worldwide. These findings are significant for designing efficient water management strategies in the Gilgel Gibe watershed and understanding how environmental changes affect water yields in this critical region. Incorporating scenario-based analyses and methodologies from similar research could improve the study’s results, giving policymakers solid facts to address socioeconomic concerns and encourage sustainable water resource management practices.

Materials and methods

Study area

The Gilgel Gibe watershed is located 272 km southwest of Addis Ababa, in the Oromia region of Ethiopia (Fig. 1). The watershed lies between 7°19′N and 8°4′N latitude and 36°37′E and 37°53′E longitude. The Wabe River, which drains much of the West Gurage land, contributes to the Omo River, forming part of the Omo-Gibe River Basin. This basin encompasses both the Gibe and Omo Rivers, which drain the upper and lower parts of the watershed, respectively. These watersheds have a total catchment area of 5,105.22 km2, of which the Gibe watershed comprises 961.21 km2 (18.83%). This study assessed the potential impacts of climate and LULC changes on water yield in the Gilgel Gibe watershed, which was selected based on several criteria. First, the availability of reliable runoff data at the watershed outlets is a major consideration. Second, preference was given to watersheds in the upper reaches of river systems to minimize external influences from adjacent watersheds. Third, the watershed’s diverse topography, climatic conditions, and socio-economic activities, especially varying precipitation patterns, were factors that made it suitable for comprehensive analysis. These considerations ensure the ability of the study to accurately assess the watershed response to potential future climate and land-use changes.

Fig. 1
figure 1

Source: The maps were created using ArcGIS software (version 10.8).

Location map of the study area.

The Gilgel Gibe watershed exhibits a contrasting elevation range, varying from 1,100 m above the mean sea level (m.a.s.l.) in the lowlands to 3,341 m.a.s.l. in the highlands. Such variation could have significantly influenced precipitation patterns, with annual rainfall ranging from 990 mm in the southern lowlands to 2,500 mm in the highlands. The average annual rainfall across the basin was about 1,508 mm. The mean annual temperature within the basin varies with regional differences, ranging from less than 7.53 °C in the western highlands to approximately 34 °C in the southern lowlands34,35. Highland areas are distinguished by their steep slopes and rugged, dissected hills, whereas lowland regions are characterized by smoother, more gently undulating terrain. The primary soil classifications observed within the study area contained predominantly Nitosols (46.72%), followed by Vertisols (27.75%), and Luvisols (18.31%). Additionally, minor occurrences of Leptosols (3.32%), Cambisols (2.4%), and Xerosols (0.61%) were also recognized36.

Data source and preprocessing

Data preprocessing plays an essential role in LULC classification and has a substantial impact on the accuracy of classification results37. This involves several procedures performed before LULC classification to enhance the data quality and optimize the accuracy of the classification outputs38. Earth Observation Data (EOD), particularly satellite imagery, plays an important role in the wide use of Landsat data because of their availability through platforms such as the United States Geological Survey (USGS) and Google Earth Engine (GEE). In this study, Landsat-5, Landsat-7, and Landsat-8&9 surface reflectance Tier 1 data were selected for their high quality and accuracy. The use of multitemporal Landsat data presents several challenges, particularly owing to cloud cover and the occurrence of missing data. Cloud cover significantly affects the quality of satellite images. To ensure the reliability of the data, atmospheric corrections such as geometric, radiometric, and noise removal were applied using the Landsat-8–9 Surface Reflectance Code (LASRC). LaSRC algorithm enables greater flexibility in processing Landsat data compared with using pre-processed Level-2 products. This flexibility allowed us to address the unique environmental conditions in the Gilgel Gibe watershed, such as variations in atmospheric composition, topography, and seasonal dynamics, ensuring more accurate reflectance values.

Landsat Surface Reflectance (Level-2) is reliable for general applications and may not always be consistent for all periods or regions, particularly in areas with complex terrain and microclimates. Therefore, the use of LaSRC preprocessing techniques in this study is essential to effectively address atmospheric issues and ensure high-quality data tailored to the specific environmental characteristics of the study area. Moreover, a combination of software tools Arc Geographic Information System (ArcGIS Desktop 10.8), ERDAS Imagine Hexagon 2015, and Environment for Visualizing Images (ENVI) Exelis 5.0 was used for preprocessing tasks, including image registration, resampling, mosaicking, image classification, and extracting relevant features from the geospatial data39. Additionally, Invest 3.13.0 played an essential role in integrating the biophysical parameters and hydrological data, which provided a valuable context for our analysis. It provides spatially explicit assessments that quantify the impacts of land use and climate change on ecosystem services, including water yield and sediment retention, to inform sustainable resource management decisions. This integration allowed for more precise modeling of the interactions between land use, climate change, and water yield. Furthermore, advanced machine learning models were employed to model the combined effects of LULC changes and climate variability on water yield within the Gilgel Gibe watershed. These models, including Random Forest (RF), Support Vector Machine (SVM), Extreme Gradient Boosting (XGBoost), and Light Gradient Boosting Machine (LightGBM), provide a sophisticated framework for analyzing complex interactions and forecasting future hydrological trends based on historical data and projected scenarios. The land cover analysis was conducted for the years 1993, 2001, 2009, 2017, and 2023, providing a comprehensive overview of the changes in the study area over time (Table 1).

Table 1 Description of satellite image data.

A comprehensive analysis was conducted using satellite imagery from various Landsat missions to assess the integrated impact of climate change and LULC changes on the Gilgel Gibe water yield. Specifically, data from the Operational Land Imager (OLI) on Landsat 8–9, the Enhanced Thematic Mapper Plus (ETM +) on Landsat 7, and the Thematic Mapper (TM) on Landsat 5 were utilized. These datasets provide comprehensive insights into temporal and spatial changes in land cover over time. Each LULC class was systematically classified according to the description provided in Table 2. LULC classifications describe separate land surface types based on human activities and natural vegetation. Each class plays a crucial role in ecosystem function, and changes in these classes owing to land cover or climate change can significantly impact hydrological processes and ecosystem services. This study employed a supervised maximum likelihood classification method to analyze LULC in the Gilgel Gibe Watershed. This method is widely used and recognized in remote sensing because of its ability to classify satellite imagery based on predefined training samples. Supervised classification involves selecting representative sample areas within the imagery that are known to correspond to specific land cover types, such as forests, shrubland, grassland, wetland water bodies, croplands, and urban areas. These training samples are used to train a classifier algorithm, such as the maximum likelihood, which calculates the probability of each pixel in the image belonging to each predefined class.

Table 2 Categories of LULC changes in the study area.

Data types

This study used hydroclimatic data from multiple sources. Specifically, data on minimum and maximum temperatures as well as rainfall from 1993 to 2023 were acquired from the Ethiopian Meteorological Institute (EMI) and National Aeronautics and Space Administration (NASA) power access (https://power.larc.nasa.gov/data-access-viewer/). This hydroclimate dataset was presented in a gridded format with a resolution of 4 km × 4 km. It combines observed station data from the national network managed by EMI with satellite-based rainfall and temperature estimates from reliable sources such as the National Oceanic and Atmospheric Administration (NOAA), the Coupled Model Intercomparison Project Phase 6 (CMIP6), and the US NASA40,41,42. This gridded dataset serves as high-quality data, filling the spatial and temporal gaps in Ethiopia’s national observation network43.

Additionally, this acquired precipitation, LULC, potential evapotranspiration (PET), soil moisture, and watershed boundary data from the Google Earth Engine and the Ethiopian Space Science and Geospatial Institute (SSGI). Furthermore, a biophysical table containing attributes for each LULC type, along with the seasonality factor (Z), was obtained. The biophysical table aided the model in differentiating between wet and dry season water yield, which is essential for managing water resources effectively. Daily rainfall and temperature data from four meteorological stations within and surrounding the Gilgel Gibe watershed were obtained from EMI for 30 years, from 1993 to 2023. To address spatial climate variability in the Gilgel Gibe watershed, this study integrated station data with gridded datasets (CHIRPS and NASA) and applied spatial interpolation techniques to estimate climatic variables across the watershed, thereby minimizing the biases associated with station density. To generate a continuous spatial layer of annual average precipitation, this study utilized ordinary Kriging interpolation in ArcGIS. Kriging was selected to analyze rainfall changes due to its advanced spatial interpolation capabilities and ability to account for precipitation data variability and spatial autocorrelation. This method provides more accurate rainfall estimations, particularly in areas with complex topographies and diverse microclimates across the Gilgel Gibe watershed. In this study, the effectiveness of the Kriging, multivariate linear interpolation (MvLI), and elevation-based interpolation (SbLI) techniques in generating grid-based rainfall datasets for the Gilgel Gibe watershed was evaluated. Based on this assumption, ordinary Kriging was deemed to be the most suitable technique for this study. Figure 2 was created using CHIRPS satellite imagery and DEM data obtained from USGS Earth Explorer (https://earthexplorer.usgs.gov) and NOAA, processed and analyzed using ArcGIS 10.8

Fig. 2
figure 2

Source: The maps were created using ArcGIS software (version 10.8).

Spatially interpolated rainfall distribution in the gilgel gibe watershed.

Hydrogeological data

This study extensively examined the water yield of the Gilgel Gibe Watershed by integrating hydrogeological, remote sensing, and climate data. These data are essential for understanding geological properties, groundwater flow patterns, and aquifer characteristics that significantly influence surface water availability and hydrological processes44. It also accounts for the buffering effect of groundwater on surface water fluctuations caused by climate change and LULC dynamics45. This integrated method enables a holistic evaluation of future water availability by considering surface and subsurface water resources46. In addition, it enhances the predictive accuracy of machine learning and remote sensing models by providing a contextual understanding of the interactions between climate variables, LULC changes, and surface water yield. In the water yield models, essential variables such as precipitation, vegetation, runoff coefficients, and evapotranspiration were calibrated to maximize the prediction accuracy. Sensitivity analysis was used to determine the most influential parameters, while observed hydrological data over time enhanced the model performance across different climatic conditions and flow regimes. Previous studies have emphasized the importance of incorporating hydrogeological information into watershed-scale hydrological modeling, emphasizing its role in improving model reliability and robustness47. By combining hydrogeological data with advanced analytical techniques, this study intends to improve our understanding of the complex relationship between environmental factors shaping water resource dynamics in the Gilgel Gibe watershed (Fig. 3).

Fig. 3
figure 3

Source: The maps were created using ArcGIS software (version 10.8).

Hydrogeological maps of the gilgel gibe catchment.

Climate trend analysis

Previous research has studied temperature and precipitation patterns in Ethiopia at several geographical and temporal scales. This study improves watershed-level analysis by focusing specifically on precipitation, temperature, and potential evapotranspiration (PET) to examine their effects on surface water recharge. Climate models and trend analysis tests, along with WorldClim data, were utilized to identify trends in time-series climatic variables. The Mann–Kendall test and Sen’s slope estimator were used to detect and quantify trends in annual rainfall, maximum and minimum temperatures, average yearly temperature, and potential evapotranspiration (PET) at a 95% confidence level48,49. These methods, combined with data from the WorldClim database, ensured a robust and comprehensive analysis of the climate trends in the Gilgel Gibe watershed. Based on daily climate data from 1993 to 2023, consistent annual trends were observed for solar radiation, maximum and minimum temperatures, and wind speed, indicating stability over this period. In particular, precipitation showed a significant anomaly in 2000, deviating from typical patterns (Fig. 4a). The analysis of yearly variations (Fig. 4b) highlighted distinct trends across all variables, reflecting the dynamic nature of these climatic factors influenced by natural variability and anthropogenic impacts. Such variability emphasizes the importance of continuous monitoring and analysis for a better understanding of long-term climatic patterns and their impact on Gilgel Gibe’s water resources and to inform water yield management and planning activities. Analysis of climatic data trends for the Gilgel Gibe Basin revealed prominent variability in climate variables, reflecting significant climatic changes that affect water yield. Changing precipitation and temperature patterns directly influence the basin’s water yield, with changes in precipitation leading to periods of excessive runoff or reduced water availability, thereby altering the water balance.

Fig. 4
figure 4

a Daily time series of key climate variables in the gilgel gibe watershed. b Plots of daily climate data points over years.

Methods

This study used advanced remote sensing techniques and water yield modeling to assess the integrated impacts of climate and LULC changes on surface water yield in the Gilgel Gibe Watershed. These methodologies are systematically categorized according to their structural composition and predictive approaches, creating a comprehensive framework for assessing environmental changes50,51. By integrating various datasets and sophisticated modeling techniques, this study aimed to elucidate the complex interactions among climate variables, LULC changes, and their impacts on regional water resources (Fig. 5). The study used RF, a robust ensemble learning algorithm that enhances prediction accuracy by aggregating outputs from multiple decision trees trained on random data subsets52,53,54,55. Additionally, an SVM is used to model the complex relationships between climate factors, land-use patterns, and surface water yield, leveraging historical data to forecast future trends while accommodating additional spatial and temporal variables56,57. Furthermore, this study encompasses XGBoost and LightGBM, both of which are recognized for their efficiency and accuracy in processing large datasets. XGBoost iteratively constructs decision trees to refine predictions, whereas LightGBM applies novel strategies to enhance the training efficiency58,59,60. This provides a realistic and robust assessment of the models’performance on time-dependent datasets. The models were executed on an MITRO V15 device equipped with a 13 th Gen Intel(R) Core (TM) i7-13700H processor operating at 2400 MHz, with 10 cores and 16 logical processors. The implementation utilized Scikit-Learn and the Keras API within TensorFlow (version 2.10.0), offering a powerful and optimized framework for model training and testing. Through the integration of these advanced algorithms, this study aims to develop robust predictive models that illustrate the complex interaction between climate variables, land-use dynamics, and surface water yield, thereby providing essential insights for sustainable water management (Fig. 5).

Fig. 5
figure 5

Methodological flowchart of the study.

Hyperparameter optimization for ML models

In this study, hyperparameter optimization was vital for enhancing the accuracy and generalization of machine-learning models for predicting water yield under diverse climate and land-use scenarios. The hyperparameters of the four machine learning models, Random Forest Regressor, XGB Regressor, LGBM Regressor, and SVR were optimized using BayesianSearchCV with threefold cross-validation. For the Random Forest Regressor, the number of trees (n_estimator) was tuned between 50 and 500, while the maximum tree depth (max_depth) was adjusted between 5 and 20. These settings were evaluated over 50 iterations using threefold cross-validation. The XGB Regressor was optimized for n_estimator, max_depth, and the learning rate (ranging from 0.01 to 0.1), also employing 50 iterations and threefold cross-validation. Similarly, the LGBM Regressor’s n_estimator, max_depth, and learning rate (ranging from 0.01 to 0.1) were tuned using the same procedure. For the SVR, the regularization parameter (C) was adjusted between 0.1 and 1000, the epsilon value ranged from 0.01 to 0.1, and the kernel was selected from options’linear’and’rbf.'This systematic optimization ensured that each model was finely tuned for high predictive accuracy while mitigating overfitting. For time-series data, walk-forward cross-validation was implemented to evaluate models sequentially over yearly increments. This approach incrementally expands the training dataset with each iteration, testing the model on the next time period. Such a method simulates real-world prediction scenarios and ensures that predictions are based only on past data, avoiding data leakage. This provides a realistic and robust assessment of the models’performance on time-dependent datasets. Models were trained and tested using a rolling-forward approach. The initial training period was 1993–2016, with testing performed in 2017. This procedure was repeated iteratively, with training periods of 1993–2017 (tested on 2018), 1993–2018 (tested on 2019), 1993–2019 (tested on 2020), and 1993–2020 (tested on 2021). Independent cross-validation was then conducted using climate data (1993–2017 training, 2018–2022 validation) and hydrological data (1995–2016 training, 2017–2021 validation). To avoid overfitting and assess the generalization capabilities of each model, an 85/15 data split was employed, where 85% of the data were used for training and 15% for testing. Cross-validation using a k-fold approach was applied, which aided in preventing overfitting and ensured the robustness of each model. Parameter tuning was conducted using BayesianSearchCV to optimize hyperparameters specific to each model: depth, minimum samples split, and estimators for Random Forest; kernel type and regularization for SVM; learning rate, maximum depth, and estimators for XGBoost and LightGBM. Validation metrics included mean absolute error (MAE), root mean square error (RMSE), Nash–Sutcliffe efficiency (NSE), Willmott’s Index of Agreement (d), and R-squared (R2) values to comprehensively evaluate the model performance.

InVEST water yield model

The InVEST model was used to determine the water yield in the Gilgel Gibe watershed. This model operates on a spatially explicit framework utilizing the cell size of the Digital Elevation Model (DEM) to assess water yield dynamics. The InVEST Water Yield model was selected based on its capacity to estimate water yields across 10 representative watersheds within the Gilgel Gibe watershed, as documented in previous studies61,62. The InVEST model was used for water yield simulations in this study because of its distinctive advantages in simulating ecosystem services in complex hydrological settings63,64,65. The InVEST model is widely used for assessing ecosystem services; however, it has significant limitations, particularly the assumption of uniform soil properties over large geographic areas. This generalization can lead to inaccuracies when modeling hydrological processes in diverse landscapes, such as the Gilgel Gibe watershed, where soil characteristics, such as texture, permeability, organic matter content, and water-holding capacity, vary significantly. To address these limitations, this study integrates the InVEST model with high-resolution, localized soil datasets collected by the Ethiopian Agricultural Transformation Institute (ATI).

Additionally, complementary tools, including hybrid machine learning models, are employed to enhance the predictive accuracy in regions with insufficient data. These models combine physical principles with data-driven techniques. The integration of remote sensing, GIS, and advanced climate models further improves the capacity of the model to evaluate the impacts of LULC changes on runoff, evapotranspiration, and groundwater recharge. These methods aim to enhance the accuracy and reliability of hydrological modeling, providing a comprehensive understanding of the interactions between water, land, and climate. This model, which is widely used for ecosystem service valuation and land-use planning, requires calibration using various input parameters, some of which are either unknown or subject to significant variability. To manage uncertainty during calibration, sensitivity analysis was employed to evaluate the influence of key input variables, such as LULC, slope, and precipitation, on the model results.

Furthermore, uncertainty analysis and error propagation were used to quantify the overall uncertainty in model predictions by accounting for the variability in each input and evaluating the cumulative effect on the model’s results. To further validate the model, it was calibrated against independent data sources by comparing model predictions with independent datasets derived from field observations and remote sensing products. InVEST is suitable for integrating land use and climatic factors into hydrological outputs, which aligns with perfectly the study’s objective of understanding how LULC changes and climate variability influence water yield. The model used remote sensing and GIS data to generate quantitative water yield outputs at the watershed scale, which were then compared with the observed runoff data for calibration and validation. Based on the Budyko curve and the average annual precipitation, the annual water yield (WY) for each location (pixel) was calculated using Eq. 1, as outlined by66,67.

$${\text{Y}}\left( {\text{x}} \right) = \left. {\left( {1 - \frac{{{\text{rAET}}\left( {\text{x}} \right)}}{{{\text{P}}\left( {\text{x}} \right)}}} \right.} \right){{{\text{*}}}}\,\,{\text{P}}\left( {\text{x}} \right)$$
(1)

where AET(x) and P(x) represent annual actual evapotranspiration (mm) and precipitation (mm) at pixel x, respectively. This study employs the Budyko equation to estimate actual evapotranspiration (AET) at each pixel or location within the landscape Eq. 268,69,70.

$$\text{AET}(\text{x}) /\text{ P}(\text{x}) = 1 + (\text{PET}(\text{x}) /\text{ P}(\text{x})) - \{ 1 + [(\text{PET}(\text{x}) /\text{ P}(\text{x})) *\upomega ] \}{^\wedge}(1/\upomega )$$
(2)

where AET(x) represents the annual actual evapotranspiration (mm) and P(x) is the precipitation (mm) at pixel x. Additionally, a non-physical parameter (\(\upomega\)) Is incorporated to account for the partitioning of precipitation between runoff and evapotranspiration. Watersheds with elevated ω values demonstrated an increased ability to transform precipitation into evapotranspiration, whereas those with lower values exhibited a reduced capacity for this process.

$$\text{PET}\left(\text{x}\right)=\text{ Kc}(\text{lx}) *\text{ ET}0(\text{x})$$
(3)

The PET for each pixel was calculated using Eq. 3. This equation incorporates a vegetation coefficient (Kc(lx)) that varies based on LULC type and ranges from 0 to 1.5. In addition, the reference evapotranspiration (ET0(x)) was included, which was computed using a modified version of71,72 Azzam et al. (2023) and Hargreaves and Samani (1985).

$$\text{ET}0=0.0013 \times 0.408 \times \text{ RA }\times (\text{Tav }+ 17) \times (\text{TD }- 0.0123\text{P})^\wedge0.76$$
(4)

In Eq. 4, RA represents solar radiation, Tav is the average of the daily maximum and minimum temperatures, TD is the temperature difference between the daily maximum and minimum, and P is the monthly precipitation (Li et al., 2022).

$$\upomega \left(\text{x}\right)=\text{ZAWC}(\text{x}) /\text{ P}(\text{x}) + 1.25$$
(5)

The \(\upomega\) parameter is an experimental value requiring further investigation for site-specific calibration. The InVEST model uses a computational formula for ω (Eq. 5), based on global data73,74,75. The Z parameter describes the spatial distribution of precipitation events. Areas with more precipitation events (for the same total annual precipitation) typically have higher Z values76. Budyko’s theory states that a higher Z value reduces the sensitivity of the model to this parameter77. Table 3 presents an in-depth evaluation of the water yield in the Gilgel Gibe Watershed, emphasizing the effects of LULC variations and climate variability.

Table 3 Input data for the InVEST water yield model.

The sensitivity ranges for various input data reflect the inherent uncertainties and variability. LULC classification values (± 20–30%) account for misclassification and changes over time, whereas precipitation data (± 10%) capture uncertainties from satellite and ground-based measurements. Temperature data (± 1 °C to ± 2 °C) allow for minor fluctuations, and watershed and sub-watershed boundaries depend on the vector data and DEM resolution. The evapotranspiration range (± 10%) reflects climatic variability and soil moisture measurements (± 10%) address spatial and temporal data uncertainties. The ± 10% Z factor was selected to address potential variability in slope steepness and account for uncertainties in terrain data, which limitations may influence in the resolution or quality of topographic datasets. By applying a ± 10% range, the calibration process ensures that the model accurately depicts the impact of terrain on hydrological outcomes without overfitting localized data that may not capture broader trends. The Kc factor (± 10%) accounts for crop-specific variations, and the discharge data may have uncertainties owing to measurement errors or gaps. These ranges ensure that the model calibration accommodates data variability.

Spatially explicit modeling of water yield

This study employed a spatially explicit model to assess water yield dynamics within the study area. The model operates on a gridded map where each pixel represents a specific location. This method allows the integration of spatial heterogeneity in key water yield drivers such as soil type, precipitation, and vegetation type78,79. The model calculates the average annual water yield (Yxj) for each grid cell (usually 300 m x 300 m) using factors such as the average annual precipitation (Px), annual reference evapotranspiration, soil depth, plant available water content (PAWC), root depth, and land-use characteristics. This equation-based calculation enabled the estimation of the water yield across the entire study area, accounting for spatial variations in these driving factors. The modeling process was complicated by uncertainties in the high-resolution input data, including LULC classification variability, which was addressed through quality control and cross-validation. The uncertainty in the precipitation data was decreased by combining satellite observations with ground-based measurements and spatial interpolation, whereas temperature variations were handled by integrating historical and satellite-derived data. The elevation data resolution was improved for a more precise watershed boundary delineation, and evapotranspiration and soil moisture uncertainties were minimized using multi-source data and interpolation. Seasonal adjustments and literature-based estimates were used to account for the Kc factor’s variability. Finally, discharge data gaps and seasonal variations were addressed by cross-referencing with the model predictions, ensuring that the model effectively handled the uncertainties in each variable to provide reliable results.

Climate and LULC change scenarios

Scenario analysis is a powerful tool used to forecast future outcomes based on projected conditions, assuming that current observed phenomena or trends continue. Understanding the association between climate change and LULC change is essential for predicting and managing water resources. The water yield of the Gilgel Gibe Watershed was estimated using the InVEST water yield module for 1993, 2001, 2009, 2017, and 2023. Three distinct scenarios were developed using scenario analysis methodologies to examine the effects of climate change and LULC on water yield. The scenarios considered a range of variables, including probable temperature changes, precipitation patterns, and LULC change trajectories. For instance, scenarios might include variations in urban expansion, deforestation rates, and agricultural practices, along with different climate projections, such as increased rainfall variability or higher average temperatures. Combining LULC changes with climate change projections provides a more precise and holistic understanding of the water resource dynamics, as noted by80. Similarly81, emphasized the necessity of scenario analysis to anticipate and mitigate potential adverse effects on water yield and ecosystem services. Finally, the scenario analysis of climate and LULC change in the Gilgel Gibe watershed underlines the potential impacts on water yield and the significance of integrated management methods that encompass land-use policies and climate adaptation measures.

Evapotranspiration data

This study used Google Earth Engine satellite data, including the Global Land Data Assimilation System (GLDAS), Moderate Resolution Imaging Spectroradiometer (MODIS), and Penman–Monteith, to assess the impact of evapotranspiration on the water yield of the Gilgel Gibe Watershed. The FAO-56 Penman–Monteith (PM) equation is widely recognized as the benchmark method for estimating reference evapotranspiration (ET₀), ensuring consistency and comparability across diverse geographic regions. MODIS and GLDAS datasets were selected for their extensive spatial coverage, high temporal resolution, and suitability for the hydrological processes of the study area. MODIS provides vital data on vegetation indices, land cover, and temperature, essential for understanding evapotranspiration (ET) dynamics (Heck et al., 2019). The daily evapotranspiration data series were obtained by temporally interpolating values adjacent to MODIS 16-day Terra climate products, yielding data that corresponded to satellite overpass times (Singh et al., 2024). Various datasets were used for the water yield analysis, including evapotranspiration (ET), Latent Heat Flux (LE), Potential ET (PET), Potential LE (PLE), Latent Heat layers (LE and PLE), and annual quality control (ET_QC). The resolution of the MODIS products, ranging from 1 km to 500 m, was selected as the final resolution for ET calculations82. The GLDAS air temperature products were downscaled to a 1 km resolution using a series of steps: (a) For each GLDAS grid (0.25° × 0.25°), the air temperature was considered as the temperature at the grid’s average elevation. (b) The sea level temperature data were then resized to a resolution of 1 km using bilinear interpolation. (c) The air temperature at the original elevation was calculated from the resized data, effectively reversing step (a).

In addition, GLDAS offers key meteorological variables such as precipitation, temperature, and solar radiation, which are ideal for modeling water balance and ET83. GLDAS data are particularly valuable in areas with limited ground-based measurements. The daily ET estimates are derived through a process that includes (1) data assimilation, (2) model calibration, (3) potential evapotranspiration (PET) calculation using atmospheric variables, and (4) conversion of PET to actual ET by considering soil moisture and vegetation. These steps ensured that the GLDAS provided accurate and reliable ET estimates for hydrological assessments (Table 4). The Penman–Monteith equation is a widely recognized method for estimating evapotranspiration (ET), which refers to the water transfer from the land surface to the atmosphere84. This process occurs through evaporation from the soil and water bodies, as well as transpiration from vegetation. However, its accuracy is highly dependent on the availability of high-quality meteorological data, including solar radiation, wind speed, temperature, and humidity. In data-scarce regions, the lack of reliable meteorological inputs can introduce significant errors, leading to misclassifications in ET estimations. To enhance the accuracy of this study, we incorporated high-resolution spatial datasets along with ground-based meteorological observations. This integration provides a more robust assessment of ET₀ variability and mitigates the uncertainties associated with missing or low-quality data. However, the PM equation is computationally demanding and requires specialized meteorological stations and remote sensing datasets for precise parameter estimation. Given these challenges, the application of the Penman–Monteith equation necessitates careful data handling and methodological refinement to ensure reliability. An accurate estimation of ET₀ depends on precise incident solar radiation (SRD) measurements, either directly obtained from meteorological stations or derived through gap-filling techniques85. The selection of the PM method in this study was validated by comparing its ET₀ estimates with MODIS-derived ET data and ground-based observations, thereby confirming its suitability for assessing regional evapotranspiration patterns. The FAO Penman–Monteith method for calculating reference evapotranspiration (ETo) is expressed as follows (Eq. 6):

$$\text{ET}=\frac{0.408\Delta \left(\text{Rn}-\text{G}\right)+\frac{\upgamma 900}{\text{T}+273}\text{u}2(\text{es}-\text{ea}) }{\Delta +\upgamma (1+0.34\text{u}2)}$$
(6)

where ET is evapotranspiration (mm/day), Δ is the slope of the saturation vapor pressure curve (kPa/°C), Rn is the net radiation (MJ/m2/day), G is the soil heat flux density (MJ/m2/day), γ is the psychrometric constant (kPa/°C), T is the mean air temperature (°C), u₂ is the wind speed at 2 m height (m/s), es is the saturation vapor pressure (kPa), and ea is the actual vapor pressure (kPa).

Table 4 Satellite data used in ET calculations.

Model scenario settings

This study developed detailed and immersive simulated scenarios, advancing beyond the mere comparisons found in the literature. Previous studies on ecosystem service simulations typically created two to four development scenarios, such as economic priority, ecological protection, and sustainable development, by limiting land-use expansion86,87,88. In contrast, this study adopted a novel approach that integrates historical climate and land-use data to create diverse scenarios that assess the relative impacts of climate change and LULC change on water yield services. Table 5 outlines the three main types of scenarios explored in this study. Initially, the actual water yield conditions in the Gilgel Gibe River Basin were computed using the InVEST model’s water yield module for 1993, 2009, 2019, and 2023. In the climate change scenarios, LULC change was held constant to isolate the effects of climate variation on ecosystem water services. This study segmented the period from 1993 to 2023 into five intervals (1993–2009, 2009–2023, 2009–2019, 2019–2023, and 1993–2023) corresponding to scenarios 1 to 3. For instance, scenario 1 utilized climate data from 1993 (multi-year average from 1996 to 2006), whereas the LULC data remained unchanged from 1993. This allowed for the assessment of climate change impacts on water yield compared to the actual conditions in 1993. Similarly, the LULC change scenarios kept climate elements constant and divided the study period into scenarios 2–3 to assess the impact of LULC changes on regional water yield over three time periods. Comparing these ten sub-scenarios with actual conditions revealed the specific contributions of climate change and LULC change to water yield across different phases.

Table 5 Scenarios for climate change, LULC change, and surface flow.

Accuracy assessment

Assessing accuracy is crucial for determining how precisely maps can be used for accurate and effective data usage. Inadequate precision levels often indicate poor satellite data classification for the LULC categories89,90. The classification accuracy was assessed using precision metrics derived from the error matrices. An error matrix was created using the validation points and classified data. The Kappa coefficient (KC) was used to evaluate classification accuracy. Per-class evaluation improves the Kappa coefficient by resolving the accuracy limitation caused by the amount of correctly identified cases in each class, often known as user accuracy91. In LULC classification, where some classes may be under or overrepresented, the Kappa coefficient is essential for managing class imbalance and introducing bias into accuracy assessments92. It corrects for such imbalances by considering classification consistency across all classes rather than just overall accuracy, which is especially essential for rare or fragmented land classes that may distort conventional measures93. A confusion matrix was used to examine metrics such as overall accuracy (OA), producer accuracy (PA), and consumer accuracy (CA)94. Using a random sampling process, 700, 936, 679, 800, and 750 Ground Control Points (GCP) were collected for 1993, 2001, 2009, 2017, and 2023, respectively. The Kappa coefficient was calculated using the following formula:

$${\text{K}}^{ \wedge } = \frac{{{\text{n}}\mathop \sum \nolimits_{{{\text{i}} = 1}}^{{\text{r}}} {\text{xii}} - \mathop \sum \nolimits_{{{\text{i}} = 1}}^{{\text{r}}} \left( {{\text{xi}} + { }.{\text{ x}} + {\text{i}}} \right)}}{{{\text{N}}^{2} - \mathop \sum \nolimits_{{{\text{i}} = 1}}^{{\text{r}}} ({\text{xi}} + { }.{\text{ x}} + {\text{i}})}}$$
(7)

In the error matrix, represents the kappa coefficient, r denotes the number of rows, Xii indicates the number of observations on the main diagonal, Xi +  is the total number of observations in row i, X + i is the total number of observations in column i, and N is the total number of observations within the matrix.

Change detection

Change detection in LULC is crucial for understanding the ecosystem of a watershed and its impact on water yield. Change detection matrices, which provide a structured approach for quantifying and interpreting LULC changes over time, were used in this process95. The Gilgel Gibe watershed in Ethiopia has experienced significant LULC changes over the years due to factors such as agricultural expansion, urbanization, unsustainable water use, and climate change96,97. This study used a change detection matrix based on remote sensing data and geospatial analysis to analyze LULC changes and their impact on water yield. The change detection matrix identifies transitions between different land cover types, allowing for a detailed understanding of the landscape dynamics98,99. The following equations were used to quantify changes in LULC:

$$\text{Loss}=\text{Row total}-\text{diagonals of each class unchanged}$$
(8)
$$\text{Net change}=\text{Gain}-\text{Loss}$$
(9)
$$\text{Gain}=\text{Column total}-\text{diagonals of each class}(\text{unchanged})$$
(10)
$$\text{Net persistence}(\text{NP})=\frac{\text{Net change}}{\text{diagonals of each class}(\text{unchanged})}$$
(11)

Model calibration and validation

The model was calibrated using data from January 1, 1993, to December 31, 2023. The prediction accuracy of the model was validated by segmenting the study period. Calibration data were used from 1993 to 2018, whereas validation data were set aside from 2019 to 2023. This division enabled us to assess the model’s capacity to generalize novel, previously unseen data. This involved minimizing the differences between the simulated and observed streamflow data at the three stream gauges in the model area. Calibration was accomplished by fine-tuning the essential parameters to eliminate variances in various streamflow characteristics, such as volume error, mean flow, minimum flow, maximum flow, runoff, and seasonal volume (summer and winter). The performance of the InVEST model is determined by parameter Z, which is an empirical constant representing regional precipitation patterns and hydrogeological properties100,101. The best Z value is often established by comparing the simulated water yield of the model to the observed runoff data102. This study evaluates the accuracy of the InVEST model by analyzing runoff data from three hydrological stations in the Gilgel Gibe Basin: Gilgel Gibe, Bulbul, and Bidru Awana, from 1995 to 2023. The Ethiopian Ministry of Water and Energy published data obtained from the monitoring systems at each station. Five statistical metrics were used to evaluate model performance: Nash–Sutcliffe Efficiency (NSE), Willmott’s Index of Agreement (d), Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and Coefficient of Determination (R2)103,104,105,106. The NSE index, which ranges from negative infinity to one, assesses the prediction accuracy; a value of one indicates a perfect match between the observed and simulated data. The NSE can be defined mathematically as

$$\text{NSE}=1-\frac{{\sum }_{ i=1}^{n}\left(yi-\widehat{yi}\right){ }^{2 }}{{\sum }_{ i=1}^{n}\left(yi-\widehat{yi}\right){ }^{2}}$$
(12)

Willmott’s Index of Agreement (d) measures the alignment between model predictions and observed data, ranging from 0 to 1, where 1 signifies perfect agreement. Higher numbers imply better model performance. The calculation is performed using the following formula, where n represents the number of observations, yi denotes the actual value of the i-th observation, and y ̂i indicates the predicted value for the i-th observation:

$$\text{d}=1-\frac{{\sum }_{ i=1}^{n}\left(yi-\widehat{yi}\right){ }^{2 }}{{\sum }_{ i=1}^{n}\left(yi-\widehat{yi}\right){+\left(yi-\widehat{yi}\right)}^{2}}$$
(13)

where yi represents the observed values, \(\widehat{yi}\) denotes the simulated values, and \(\widehat{yi}\) bar \(\left|\right.\widehat{yi\left|\right.}\) Is the mean of the observed values. The Root Mean Square Error (RMSE) measures the average magnitude of errors between expected and observed values. This method calculates the standard deviation of the residuals, which reflects the degree of divergence from the true values. The formula for RMSE is given by:

$$\text{RMSE}=\sqrt{\frac{1}{\text{n}}\sum_{ i=1}^{ n} \left(yi-\widehat{yi}\right){ }^{2}}$$
(14)

where: yi​ represents the observed values, \(\widehat{y}\) i​ Denotes the predicted values and n is the number of observations. Mean Absolute Error (MAE) measures the average size of mistakes in a dataset by computing the mean of the absolute differences between observed and anticipated values. This metric measures the average magnitude of the residuals and provides insights into the correctness of the model. The MAE was calculated using the following formula:

$$\text{MAE}=\frac{1}{n}{\sum }_{ i=1 }^{n} \left| \right.yi- yi \left|\right.$$
(15)

The Coefficient of Determination (R2) indicates the variance in the dependent variable that can be explained by a linear regression model. It is scale-independent, meaning that it is not affected by the size of the dependent variable. R2 ranges from 0 to 1, where a value of 1 indicates a perfect fit.

$${\text{R}}^{2}=1-\frac{{\sum }_{ }^{ } \left(\text{yi}-\widehat{\text{yi}}\right){ }_{ }^{2 }}{{\sum }_{ }^{ } \left(\text{yi}-\widehat{\text{yi}}\right){ }_{ }^{2 } }$$
(16)

Linear regression analysis was conducted to evaluate the accuracy of the simulated water yields from the InVEST and ML models. The results shown in Fig. 6 reveal a strong correlation between the simulated and observed datasets. This strong correlation validates the efficacy of the climate and hydrological models in accurately simulating the water yield and is closely aligned with the observed runoff patterns.

Fig.6
figure 6

Correlation matrix for various features in the datasets.

Results

LULC change analysis

The results from our remote sensing classifier provide insightful predictions for LULC classes across five distinct periods: 1993, 2001, 2009, 2017, and 2023 (Table 6). These results underscored the dynamic nature of the various LULC categories over time. For instance, shrubland, forest, grassland, and wetland covers showed notable fluctuations. Shrubland decreased from 1108.37 km2 (21.54%) in 1993 to 295.22 km2 (5.74%) in 2023, reflecting a significant coverage area loss of 849.15 km2 between these two periods. The extent of water bodies also changed by increasing from 12.51 km2 (0.24%) in 1993 to 41.57 km2 (0.81%) in 2023, largely due to the construction of the Gilgel Gibe hydroelectric dam starting in 2003. This study showed that the Gilgel Gibe Dam has a significant impact on LULC, resulting in both direct and indirect environmental changes. Since its construction in 2003, the dam has altered local hydrology, comparing data from 1993, leading to changes in vegetation, soil properties, and land use in the surrounding areas. These changes are primarily driven by the creation of the reservoir, altered water availability for agriculture, and dynamic socioeconomic activities. The dam has caused land inundation, converting terrestrial areas into water bodies, which has destroyed natural habitats, such as forests, grasslands, and wetlands, leading to habitat loss and fragmentation.

Table 6 LULC changes from 1993 to 2023.

In addition, it has contributed to the expansion of irrigated agriculture. The forested areas declined from 626.73 km2 (12.18%) in 1993 to 534.18 km2 (10.38%) in 2023. The main reason for forest conversion to agricultural land in the Gilgel Gibe watershed is the rapid population growth and increasing demand for arable land. Local households depend on agriculture as their primary source of livelihood, clear forests to meet their subsistence, and commercial farming needs, further exacerbated by limited access to alternative income sources. Government initiatives such as villagization and large-scale farming projects also contribute to agricultural land expansion. Additionally, the fertile soils in the watershed make them highly attractive for agriculture, leading to further encroachment on forests. To address the loss of ecosystems, stakeholders should develop strategies that balance agricultural development with forest conservation, and promote sustainable land use and ecosystem health. Similarly, grasslands and wetlands showed a dramatic reduction, with grasslands shrinking from 941.67 km2 (18.30%) in 1993 to 184.63 km2 (3.59%) in 2023 and wetlands decreased from 454.98 km2 (8.84%) to 302.35 km2 (5.88%) over the same period. Grasslands and wetlands are vital ecosystems that support biodiversity, regulate the climate, and provide essential ecological services. These changes can be attributed to complex interactions between natural and anthropogenic factors. Climate cycles and long-term trends may have contributed to the reductions in grassland, wetland, forest, and shrubland areas. Shrinking grasslands fragment habitats, reduce food sources, and breeding grounds for species such as herbivores, predators, and pollinators. Similarly, wetland loss disrupts critical breeding and migration routes for birds, fish, and amphibians, leading to a population decline in the study area. Built-up areas have expanded consistently, driven by increasing demand for housing and infrastructure development. Urban areas expanded consistently from 31.02 km2 (0.60%) in 1993 to 150.78 km2 (2.93%) in 2023. This growth reflects the progression of urban development shaped by resource availability and prevailing environmental conditions. Conversely, from 2017 to 2023, the urban landscape of the study area has changed significantly owing to government policies targeting illegal housing. A demolition policy implemented in 2022 and continuing to 2023 reduced built-up areas by removing illegal houses and settlements. This effort addressed the unregulated urban expansion and promoted urban greenery. Cultivated land also showed a substantial increase, from 1969.58 km2 (38.28%) in 1993 to 3636.17 km2 (70.68%) in 2023. This unprecedented expansion was driven by population growth, improved agricultural practices, and favorable irrigation conditions. These results highlight the critical need to implement sustainable land management strategies that harmonize developmental goals with environmental conservation. Targeted policies and conservation efforts are essential to mitigate adverse effects and promote sustainable development in the Gilgel Gibe watershed (Fig. 7).

Fig. 7
figure 7

LULC map of 1993 to 2023 Source: The maps were created using ArcGIS software (version 10.8).

Generally, in the study area, land use changes, particularly agricultural expansion, significantly affect water yield. Clearing vegetation for crops reduces soil permeability and increases runoff, thereby disrupting water infiltration, groundwater recharge, and streamflow. For example, large-scale farming near the Gilgel Gibe Dam has led to higher runoff during rainfall. The conversion of wetlands and grasslands also exacerbates this by removing natural water regulation, leading to lower water yields in dry periods and increased flooding during rainy periods. To mitigate these effects, sustainable land management practices like agroforestry, water-efficient irrigation, and soil conservation should be prioritized.

Accuracy assessment

The accuracy assessment of LULC classification for 1993, 2001, 2009, 2017, and 2023 shows a significant and consistent improvement in predictive performance over time. The overall accuracy of the LULC classification progressively improved, reaching 94.94% in 1993, 95.95% in 2001, 95.44% in 2009, 97.22% in 2017, and 96.71% in 2023 (Table 7). Correspondingly, The Kappa coefficient, indicating the agreement between predicted and actual categories, ranged from 0.94 to 0.96 across the study period. This metric is highly useful in remote-sensing-based LULC classification, particularly when ground truth data are acquired on handheld GPS devices. These findings demonstrate the excellent precision of our LULC classification system, which exceeds the widely accepted accuracy threshold in land cover investigations. Particularly, in 2017, the classification performed admirably, with an overall accuracy of 97.22% and a Kappa value of 0.96, indicating a high agreement between predicted and actual land cover types. In 2023, the kappa coefficient showed a slight decline to 0.9615 compared with previous years. This decrease can be attributed to rapid and fragmented LULC changes, including agricultural expansion and climate variability, within the study area. Agricultural expansion, characterized by the conversion of natural ecosystems into agricultural land, results in increasingly dynamic and unpredictable LULC patterns. This transition presents challenges in accurately classifying land cover, particularly when agricultural plots are small and dispersed or when there are shifts in crop types that may not be easily detected by classification models. The results of our classification method can be attributed to the use of advanced classification algorithms, high-resolution satellite imagery, and extensive field data. Finally, this study indicated that high-resolution satellite imagery combined with comprehensive ground-truth data can improve the precision of LULC mapping.

Table 7 Presents the accuracy assessment results for the classified images from 1993 to 2023.

Change detection matrix

A change detection matrix was developed to analyze LULC changes in the Gilgel Gibe watershed between 1993 and 2023. The analysis revealed significant transitions, with over two-thirds (65.9% or 3390.14 km2) of the watershed experiencing a change in LULC classes during this period (Table 8). The most significant change was from forest to agricultural land, accounting for 51.33% of the transformation, mainly in areas far from urban centers. Forests also transitioned into shrubland (114.55 km2) and grassland (14.32 km2), threatening biodiversity and ecosystem services Interestingly, about 5.98 km2 and 321.55 km2 of Cultivated land were converted to water bodies and forests, respectively, largely due to the construction of Gilgel Gibe I hydroelectric power dam and the reservoir created and government initiatives promoting private and community forests, particularly for coffee growers.The conversion of forests into agricultural lands significantly disrupts natural water cycles, worsening water scarcity. These changes affect the water yield by reducing groundwater recharge, increasing surface runoff, and contributing to soil erosion, ultimately leading to lower water quality and availability. For example, the absence of tree cover increases evaporation, soil moisture diminishes rapidly, and surface runoff intensifies, further reducing the ability of land to replenish groundwater recharge. Comparatively, more agricultural land transitioned to built-up areas (14.42 km2), with minimal reverse transitions. These findings align with previous studies that highlighted the expansion of urban areas into natural landscapes, necessitating sustainable urban management and ecosystem-based adaptation strategies.

Table 8 Change detection matrix of 1993 to 2023.

The results of the image difference and LULC change analysis for the Gilgel Gibe watershed from 1993 to 2023, as shown in Fig. 8, indicate substantial changes across various classes. The largest change was observed adversely, affecting 3484.62 km2, which constitutes 67.73% of the total area. By contrast, there was an increase of 615.49 km2 (11.96%). Additionally, there were minor decreases and increases, impacting 42.58 km2 (0.83%) and 736.40 km2 (14.31%), respectively. Notably, 266.04 km2 (5.17%) of the study area remained unchanged. Overall, changes that accounted for a total area of 5145.14 km2 indicated significant land cover dynamics in the watershed during the study period. This analysis highlights significant transitions in land cover types in the Gilgel Gibe watershed, emphasizing the need for sustainable land management practices. The observed dynamics, indicating that water yield changes are more attributable to loss mechanisms (68%), indicate significant implications of increasing cropland at the expense of decreasing green cover, such as shrubs and forest areas, on basin water yield. As cropland expands, it typically replaces areas of natural vegetation, which has substantial effects on the hydrological cycle. Reduced green cover diminishes the land’s ability to absorb and retain rainfall, leading to increased surface runoff and decreased groundwater recharge. This change in land use can exacerbate the impacts of observed changes in rainfall, evapotranspiration, and temperature. For instance, increased cropland often correlates with higher levels of evapotranspiration due to intensified agricultural activities, potentially leading to reduced water availability in the basin. Concurrently, decreased green cover lowers the transpiration rates and disrupts local climate regulation, which can further alter temperature patterns and precipitation dynamics. This reduction in natural vegetation also diminishes the basin’s capacity to stabilize water yields, making it more susceptible to extreme weather events and variability in rainfall. Overall, the LULC dynamics marked by increased cropland and decreased green cover can thus have profound implications for water yield, exacerbating the effects of changing climatic variables and reducing the basin’s resilience to hydrological fluctuations.

Fig. 8
figure 8

Change matrix and image difference map of 1993 to 2023 Source: The maps were created using ArcGIS software (version 10.8).

Climate change scenario analysis

Figure 9 shows the annual variability in rainfall, revealing notable fluctuations in monthly precipitation levels across different years. In 1993, the mean monthly maximum rainfall reached 87.54 mm, the highest recorded during the study period, while the minimum was 20.26 mm. The elevated rainfall levels observed in this study are likely driven by strong monsoonal activity and positive climate anomalies, such as the Indian Ocean Dipole (IOD) or El Niño-Southern Oscillation (ENSO). Similar research highlights the influence of atmospheric circulation patterns, such as the Intertropical Convergence Zone (ITCZ), in generating high rainfall variability107.The year 2009 also saw substantial rainfall, with mean monthly values ranging from a peak of 68.54 mm to a low of 26.97 mm. These values suggest enhanced precipitation, possibly linked to global climatic phenomena, such as La Niña events, which are known to bring wetter-than-average conditions to parts of East Africa. The year 2017 was particularly significant because of the severe precipitation conditions. During this year, the mean monthly rainfall ranged from 12.64 mm at its highest to 7.24 mm at its lowest, coinciding with a major drought affecting much of Ethiopia. Severe drought coinciding with low rainfall, attributed to Indian Ocean warming and altered atmospheric circulation, affects Ethiopia, causing widespread agricultural losses, water scarcity, and food insecurity108. In 2001 and 2023, there were notable variations in the mean monthly rainfall. In 2001, rainfall ranged from a high of 26.54 mm to a low of 9.28 mm. In 2023, the maximum was slightly lower at 16.84 mm, and the minimum was 7.95 mm. These variations illustrate the influence of climatic conditions on rainfall patterns over time. Rainfall variability during these years likely reflects regional land use changes and moderate climatic disturbances. Declining trends have been attributed to urbanization and global warming-driven changes in precipitation patterns109.

Fig. 9
figure 9

Spatial distribution of precipitation and temperature variations and their influence on surface water yield in the gilgel gibe catchment.

Figure 10 shows the seasonal temperature and precipitation patterns. Temperatures are highest from March to June (median ~ 21 °C in May) and lowest in December and January (~ 18 °C), with constant variations. Precipitation is low from December to February and increases significantly from June to September, peaking in July and August with medians over 2,000 mm and a maximum exceeding 3,500 mm. These patterns emphasize various climatic cycles that are important for agricultural and water resource management.

Fig.10
figure 10

Correlation matrix of climate variables from 1993 to 2023.

The correlation analyses presented in Table 9 further highlight the performance of this model. For the period from 1993 to 2001, the correlation value was 0.9957, indicating a high agreement between the model simulations and observed data, reflecting the model’s reliability during this period. The correlation improved to 0.9986 from 2001 to 2009, suggesting enhanced model accuracy and performance, which is likely due to improved calibration or data quality. From 2009 to 2017, the correlation remained high at 0.9983, demonstrating consistent model performance. The most recent period, from 2017 to 2023, showed the highest correlation value of 0.99997, indicating an almost perfect alignment between the simulated and observed data, possibly due to advancements in modeling techniques or data quality. Overall, the correlation value for the entire study period from 1993 to 2023 was 0.9959, highlighting the model’s strong performance and ability to accurately simulate water yield dynamics under diverse climatic and land-use conditions.

Table 9 Correlation analysis of long-term climate data.

Comparative evaluation of IDW and kriging interpolation models

Spatial interpolation was used in this study to estimate climate variables (precipitation and temperature) across ungauged locations in the Gilgel Gibe watershed for water yield modeling. This study compared two widely used methods: Inverse Distance Weighting (IDW), which is deterministic, and distance-based weighting. Ordinary Kriging: Geostatistical method incorporating spatial autocorrelation (Table 10).

Table 10 Comparison of annual mean precipitation interpolation (IDW and Kriging) metrics.

The comparative analysis of annual mean precipitation interpolation methods demonstrated that Ordinary Kriging with an exponential variogram model (including the nugget effect) outperformed IDW (power parameter p = 2) across all validation metrics. Kriging achieved a lower RMSE (89.7 mm/yr vs. 112.3 mm/yr) and MAE (68.4 mm/yr vs. 87.6 mm/yr), representing a 20.1% and 21.9% improvement in accuracy, respectively. Additionally, Kriging exhibited a stronger correlation with the observed data (R2 = 0.81 vs. 0.72), confirming its superior ability to capture spatial precipitation patterns, particularly in regions with complex terrain or sparse station coverage. These findings align with those of previous studies in similar climatic regions, reinforcing Kriging as the preferred method for hydrological applications requiring high-precision precipitation estimates.

Evaluation of ML algorithms for climate data analysis

The climate data for the Gilgel Gibe watershed indicate a significant temperature range, with maximum values reaching 34 °C and minimum values decreasing to 14 °C. Precipitation data revealed that the region experiences monthly rainfall ranging from 0 to 50 mm. Solar irradiance levels in the area vary between 6 and 7.75 kilojoules per square meter (kJ/m2), while wind speeds fluctuate between 0 and 8 m per second. The relative humidity in the watershed ranges from 40 to 90%, reflecting the region’s role as a condensation point for air masses that condense into liquid water at constant atmospheric pressure. This climatic variability throughout the year has a notable impact on water yield in the watershed (Fig. 11a, b, and c).

Fig. 11
figure 11figure 11figure 11

(a, b, and c) Comparison of observed and predicted values for climatic variables using different machine learning algorithms.

The climate model performance was assessed using data from two precipitation stations (Gilgel Gibe and Jimma) and two temperature stations (Omo Nada and Deneba). The performance of the SVM algorithm was notably high across the analyzed climate stations. For the Gilgel Gibe station, the SVM model achieved the following performance metrics: MAE = 0.36, RMSE = 0.45, R2 = 0.81, NSE = 0.81, and Willmott’s Index of Agreement (d) = 0.94. These results indicate that the SVM model effectively minimized the prediction errors and was closely aligned with the observed data, demonstrating a robust capacity for capturing complex climatic patterns. Similarly, at Jimma station, the SVM model performed with MAE = 0.32, RMSE = 0.44, R2 = 0.74, NSE = 0.74, and d = 0.91. Table 11 and Fig. 12 provide a comprehensive overview of the model performance, showing that all models exhibited strong fits to the observed data, as evidenced by metrics approaching 1. This alignment indicates that the models are highly accurate with minimal inconsistencies between the predicted and actual values. Finally, the analysis underscores the robustness of our climate models, revealing significant correlations between climatic variation and water resources in the Gilgel Gibe watershed. The results confirm the accuracy and reliability of the model in assessing climate change impacts on water yield, offering a solid foundation for effective water resource management and climate adaptation strategies.

Table 11 Comparison of machine learning model performance for monthly climate data in the gilgel gibe watershed.
Fig.12
figure 12

Comparison of performance metrics for RF, SVM, XGB, and LGBM on climate data.

Generally, the Gilgel Gibe watershed is a crucial hydrological resource in Ethiopia and faces water yield challenges driven by climate change, land use dynamics (deforestation, agricultural expansion, and urbanization), and socio-economic pressures. These activities disrupt the hydrological cycles, leading to increased surface runoff, reduced infiltration, and diminished groundwater recharge. Over the past 30 years, the study area has experienced a loss of 276.51 km2 (5.36%) of forest cover. The removal of forest vegetation leads to increased surface runoff and sediment yield, ultimately reducing water availability. Although deforestation may initially lead to a temporary increase in water yield, it disrupts the hydrological cycle by accelerating soil degradation and reducing groundwater recharge capacity. Findings from the Gilgel Gibe sub-basin indicate that deforestation has resulted in a 15% increase in annual runoff while causing a 20% reduction in dry-season baseflow due to diminished groundwater recharge110. This study confirmed that deforestation is a primary factor driving LULC fluctuations within the Gilgel Gibe Catchment. Forests are being cleared for agriculture, construction, firewood, charcoal production, and timber extraction; however, replanting and conservation efforts remain inadequate. Consequently, the region faces multiple environmental challenges, including increased soil erosion, water scarcity, heightened flood risks, rising temperatures, and more frequent droughts. Ethiopia’s climate, characterized by alternating wet and dry seasons influenced by the Inter-Tropical Convergence Zone (ITCZ), exacerbates these challenges, resulting in erosion, sedimentation, and altered water yields. Similar trends have been observed in other studies, such as111, who explored topography-based methods for downscaling sub-grid precipitation in Earth system models, as well as in watersheds such as the Awash River Basin112, which revealed reduced base flow and increased runoff due to land cover change. To address these issues, actionable recommendations include strengthening climate-resilient infrastructure (check dams, terraces, reservoirs), as demonstrated in the Upper Tana Basin113; implementing community-based water use planning, such as water-efficient irrigation training; incorporating ecosystem-based approaches such as reforestation and controlled grazing114 and integrating local actions with national (Ethiopia’s Growth and Transformation Plan (GTP) and global (Sustainable Development Goals (SDG) policies, particularly within frameworks such as the Nile Basin Initiative. These strategies, which focus on infrastructure, community engagement, and ecosystem preservation, can inform broader sustainable water management across Ethiopia and similar regions.

Machine learning-based streamflow modeling

In this study, four machine learning models, RF, SVM, XGB, and LGBM, were used to simulate streamflow. The daily streamflow predictions of these models are illustrated in Fig. 13 (a, b, and c), offering a comprehensive view of the simulated streamflow patterns. The dataset was partitioned into two phases: the training phase, which encompassed January 1995 to December 2019 (85% of the data from January 1995 to) and the testing phase, which spanned January 2020 to December 2022 (15% of the data). The effectiveness of RF, SVM, XGB, and LGBM in streamflow simulations was evaluated using several statistical indicators. The RF model showed the highest accuracy in predicting observed discharge values compared to the SVM and XGB models, as shown in the regression graphs and the GitHub repository (https://github.com/amanuel2025/Climate-and-Hydrology-Data-Analysis-using-ML). This superior performance can be attributed to RF’s capability of RF to capture complex nonlinear relationships in the data, making it particularly adept at handling complex streamflow dynamics. Both the RF and XGB models excel in processing sequential data, whereas SVM and LGBM are better suited for dealing with fuzzy and uncertain rule-based systems.

Fig.13
figure 13figure 13figure 13

(a, b& c) of Observed and predicted streamflow for RF, SVM, XGB, and LGBM models during training and testing phases.

Fig. 13 (a, b, and c) present the analysis of the RF, SVM, XGB, and LGBM models, examining the relationship between the five input parameters (daily mean, daily flow, discharge, and runoff) and the output parameter (mean daily discharge). The model performance was evaluated using five statistical indicators: RMSE, MAE, R2, NSE, and Willmott’s Index of Agreement (d). During training, the RF model achieved MAE, RMSE, R2, NSE, and d values of 0.38, 0.61, 0.76, 0.76, and 0.93, respectively, indicating an average deviation between the predicted and observed discharges. During testing, these values improved slightly, with MAE, RMSE, R2, NSE, and WIA values of 0.25, 0.34, 0.82, 0.82, and 0.94, respectively, showing a minor increase in prediction error compared to the training phase. The XGB model showed a moderate correlation between predicted and observed discharge values, with R2 values of 0.72 during training and 0.80 during testing. This indicated that the model explained 72% and 80% of the discharge variance in the training and testing phases, respectively. These results highlight the effectiveness of the model in predicting streamflow dynamics using daily flow, maximum and minimum discharge, and runoff data from Gilgel Gibe water yield (Table 12 and Fig. 14).

Table 12 Evaluation of machine learning models for daily streamflow simulation in the gilgel gibe watershed.
Fig.14
figure 14

Statistical evaluation of RF, SVM, XGB, and LGBM models on hydrological data.

Machine learning cross validation results for climate and hydrological data

The assessment of climate and land use changes on water yield within the Gilgel Gibe watershed necessitates the integration of diverse spatial and non-spatial datasets to identify key influencing factors. Although the selection of cross-validation methods may appear secondary when a model exhibits satisfactory performance on a given dataset, our findings underscore the profound impact that cross-validation strategies have on the model’s capability to accurately evaluate climate and land use change effects on water yield. To enhance the model’s generalizability and mitigate overfitting, hyperparameter optimization was conducted using BayesianSearchCV with the hyperparameter’c'configured to 3. Furthermore, a walk-forward cross-validation approach was employed to incrementally validate the model’s performance. Initially, the dataset was partitioned into 90% for training and 10% for testing. The validation process followed a sequential progression, where the model was trained on the 1993–2016 data and tested on 2017, then retrained on 1994–2017 and tested on 2018, continuing this pattern incrementally up to 2021 for climate. This rigorous approach facilitated a more reliable evaluation of the model’s predictive accuracy over varying temporal spans. The evaluation of model performance was conducted using a comprehensive suite of statistical metrics, including MAE, RMSE, R2, NSE, and d, as detailed in Table 13 &Table 14. The cross-validation results of SVM model demonstrated strong performance in predicting streamflow components. Mean streamflow predictions achieved consistently high accuracy during training (R2 > 0.99; 2018: 0.995, 2021: 0.994). Maximum streamflow also exhibited excellent performance (R2 > 0.97 for most years, 2018: 0.980, 2017: 0.974). Stream runoff predictions yielded R2 values between 0.82 (2021) and 0.99 (2018). Minimum streamflow prediction accuracy varied, with high R2 values in 2018 (0.986) and 2019 (0.995) but a lower value in 2020 (0.733). Overall streamflow prediction accuracy also varied interannually, with higher R2 values in 2018 (0.956) and 2019 (0.895) and a lower value in 2021 (0.708). These results suggest that while the SVM model performs exceptionally for mean and maximum streamflow, minimum streamflow and streamflow for some periods require further refinement. Overall, the SVM model shows strong potential for accurate water yield predictions, particularly for critical metrics such as mean and maximum streamflow, providing reliable insights for hydrological management in the Gilgel Gibe watershed. These metrics provide a robust framework for quantifying predictive accuracy and model reliability in capturing the underlying patterns of climate and hydrological data.

Table 13 Best-performing ML models across different training periods for hydrological data prediction/validation.
Table 14 Best-performing ML models across different training periods for climate data prediction/validation.

The cross-validation results of the SVM model for climate data in the Gilgel Gibe watershed (Table 13) revealed variable performance across different climate factors. For Maximum Wind Speed, the model achieved R2 values ranging from 0.60 to 0.95, with the highest performance observed in 2018 (R2 = 0.95). The model consistently demonstrated strong performance for Minimum Temperature, with R2 values spanning from 0.91 in 2021 to 0.97 in 2019 and 2022. Maximum Temperature displayed particularly impressive results, with R2 values nearing 0.99 in both 2020 and 2022. Solar Radiation maintained high accuracy across all years, with R2 values exceeding 0.97, peaking at 0.99 in 2020. Precipitation also showed strong predictive performance, particularly in 2020 (R2 = 0.97). Lastly, Relative Humidity exhibited a range of R2 values from 0.89 in 2018 to 0.96 in 2022. Overall, the SVM model demonstrated robust performance in predicting climate factors, with consistently high R2 values for most variables.

Evaluation and validation of the invest model

The InVEST model estimates water yield using Budyko’s hydrological framework, and its accuracy is highly dependent on proper calibration. The Z calibration constant is a key parameter in the model that controls the partitioning of precipitation between evapotranspiration (ET) and water yield. A higher Z value (> 10%) leads to an overestimation of ET, reducing water yield predictions, whereas a lower Z value (< 10%) increases the estimated water yield but may underestimate ET losses. In this study, an initial Z value of 10 was used, which was aligned with the global average for mixed land cover (Table 15).

Table 15 Calibrated Z values metrics.

The calibration of the Z factor in the InVEST model significantly influences the accuracy of water yield estimation. The results show that adjusting Z improves model performance by reducing runoff overestimation in tropical regions and refining the balance between evapotranspiration (ET) and water yield. For the Gilgel Gibe Watershed, a Z value of 11.2 resulted in NSE improvement from 0.59 to 0.70 and R2 = 0.65, with a + 5.8% bias, suggesting a better representation of mixed agricultural and forested land cover. When the Z value was set to 12, NSE improved to 0.68, R2 increased to 0.71, and bias reached + 8.2%, indicating a more stable estimate but slightly higher error in dry-season flow estimation. However, the initial Z value of 10, which aligns with the global average for mixed land cover, showed weaker performance (NSE = 0.55, R2 = 0.58, and bias = + 15%), demonstrating overestimated runoff in tropical environments. This highlights the necessity of local calibration to refine hydrological predictions. The machine learning-based optimization of Z was used for these purposes to reduce bias to + 3.1%, demonstrating its potential in fine-tuning water yield models for regional hydrology.

Evapotranspiration impact

Evapotranspiration (ET) is key to understanding global hydrological budgets and water yield modeling. Its precise estimation is essential for comprehending the water balance and formulating effective water resource management plans for the Gilgel Gibe Watershed. However, regional-scale estimations using meteorological data from the GLDAS, MODIS 1 km and 500 m, and the Penman–Monteith method, an equation used to compute water evaporation from vegetated surfaces, are less commonly employed (Fig. 15a.,15b and Fig. 16).

Fig. 15
figure 15

a.The daily potential evapotranspiration rates from 2000 to 2023 illustrate temporal variations and trends. b. Annual potential evapotranspiration rates from 2000 to 2023 illustrate temporal variations and trends.

Fig. 16
figure 16

Evapotranspiration maps from 2000 to 2023.

Analysis of potential evapotranspiration (PET) trends from 2000 to 2023 revealed varied spatial distributions without a distinct overall pattern. On a daily scale, PET trends (Fig. 15a) ranged from 1 to 4 mm, except for MODIS 500 m, which showed a slightly decreasing trend close to 0.9 mm. Annually (Fig. 15b), all evapotranspiration (ET) graphs consistently show increasing PET trends, ranging from 100 mm to 350–400 mm. Figures 15a and b illustrate the magnitudes of these detected PET trends across various meteorological stations on daily and annual scales, respectively. During summer, a seasonal increase in PET was evident, ranging from 1 to 4 mm. The Penman–Monteith-Leuning method for ET components highlights distinct behaviors in different evapotranspiration processes. Canopy Evaporation (Ec) showed values greater than 1.5 mm, reaching up to 3 mm. This represents the evaporation of water from wet canopy surfaces after rainfall, where water on the leaves and stems evaporates back into the atmosphere.

Interception Evaporation (Ei) ranges from 0 mm to 0.5 mm, indicating the water intercepted by the canopy and evaporated without reaching the ground, including the water caught on leaves and stems during rainfall. Soil Evaporation (Es) values ranged from 1 to 2 mm throughout the study period, highlighting the evaporation of water directly from the soil surface. This process is influenced by the soil moisture content, soil type, temperature, and other environmental factors. These findings are consistent with those of previous studies on evapotranspiration, such as115, which indicated that increasing temperatures and changing precipitation patterns strongly affected PET trends. The consistent increase in PET, particularly in summer, highlights the impact of climate change on water availability. The varying behavior of the canopy, interception, and evaporation of ET components offer valuable insights into the ecosystem’s water balance, stressing the need for adaptive water management strategies in response to these changes.

The spatial trend of evapotranspiration, as depicted in Fig. 16, showed that the highest ET was recorded in the years 2000 and 2017, with values of 357.34 mm and 355.24 mm, respectively. Ethiopia faces several significant natural hazards, including droughts, with different spatial and temporal patterns. In 1999/2000, the country experienced a short-season rainfall failure that affected approximately 10 million people. In 2017, spring rain failed, while the main rainy season brought heavy rainfall and flooding across the country. By June, the number of districts needing emergency food aid rose from 192 to 228, with approximately 461 districts facing severe difficulties. The Gilgel Gibe watershed was one of the sites affected by this drought in 2000 and 2017. Owing to inadequate spring rains, the number of individuals requiring food assistance has increased to 7.8 million. Additionally, linked to the El Niño weather phenomenon, the 2015–2017 drought was severe and impacted large areas of Ethiopia. This has resulted in widespread crop failure, livestock deaths, and acute food insecurity, affecting millions of people across the country116,117.

In addition, our study indicated that high ET values also reflect high water demand for crops, which is essential for maximizing agricultural yields. Ecosystems are intrinsically linked to water cycles. Rising ET values signify healthy and productive ecosystems with active vegetation growth. However, these ET values are influenced by and contribute to climate change. Increased temperature and altered precipitation patterns can elevate ET rates, creating response mechanisms that exacerbate warming. The Penman–Monteith equation proved to be the most effective model for estimating ET in the region. This model incorporates factors such as temperature, humidity, wind speed, and solar radiation. The results showed consistently high ET values in the Gilgel Gibe watershed, which was evident through both periodical and spatial variations in its distribution. Periodical variation of the changes in ET values observed over different time periods, reflecting seasonal or annual fluctuations driven by factors such as changes in vegetation, rainfall patterns, and temperature. Spatial variation denotes the differences in ET values across different geographic locations within the watershed, influenced by variations in land cover, topography, and soil characteristics.

In addition, this study investigated how climate change and land use affect the water yield. Correlation analysis showed strong associations between climate factors and water yield, while regression quantified land-use impacts. Machine learning models, including RF, SVM, EGBM, and LGBM, were trained on historical data and cross-validated to analyze future water yield under various scenarios. The combined effect of these periodical and spatial variations provides a comprehensive understanding of the dynamics of evapotranspiration in the Gilgel Gibe watershed. In summary, high evapotranspiration values have profound implications across various domains, from water resource management and agriculture to ecosystem health and climate regulation. Understanding and managing these high ET rates through advanced measurement, modeling techniques, and sustainable practices are crucial for addressing the challenges posed by changing climates and growing water demands.

Assessment of quantitative sensitivity analysis results

The sensitivity study evaluates the relative impact of land-use change, climate variability, and observed discharge trends on water yield dynamics in the Gilgel Gibe watershed (1993–2023).The quantitative sensitivity analysis determined that land-use change had a greater impact on surface runoff than climate variability, with deforestation causing considerable changes in water production. Precipitation variability was identified as the primary driver driving hydrological responses, followed by variations in temperature (Table 16).

Table 16 Quantitative assessment of LULC change, climate change, and observed discharge trends.

The analysis revealed that land-use change was the dominant driver of water yield alterations in the Gilgel Gibe watershed (1993–2023), with deforestation increasing quick flow by 22% and reducing baseflow by 15% (p < 0.01), exerting 2.7 times greater influence than climate factors. Climate change contributed to an 8% yield reduction primarily through ET increases (12%) and precipitation declines (5%), with particularly strong dry season impacts (r2 = −0.81). The observed discharge data from the three MoWE stations confirmed these trends, showing a net 14% yield increase with strong model validation (NSE = 0.76–0.88) while highlighting hydrological shifts: 18% higher peak flows and 25% lower baseflows, underscoring the watershed’s declining water regulation capacity. These quantified results demonstrate the urgent need for land use management policies to mitigate hydrological degradation.

Integrated impact of climate change and land use dynamics on surface flow and water yield in the gilgel gibe watershed

The Gilgel Gibe Watershed in Ethiopia has undergone notable climate, land use, and hydrological changes that have impacted surface flow and water yield. Over the past three decades, the region has experienced a significant increase in the mean annual temperature, exceeding the global average rate of warming. Climate data for the watershed revealed temperature variations ranging from a maximum of 34 °C to a minimum of 14 °C, with monthly precipitation fluctuating between 0 and 500 mm. Additionally, solar irradiance levels range from 6 to 7.75 kJ/m2, and wind speeds vary between 0 and 8 m/s. These climatic factors significantly affect the hydrological cycle, including evapotranspiration rates and precipitation patterns. Increased temperatures have been linked to increased evapotranspiration rates, which, along with changes in precipitation, affect water availability. Studies by118 emphasized the role of temperature variability in shaping hydrological processes, noting that higher temperatures accelerate snowmelt in high-altitude areas, altering streamflow timing and magnitude. These findings highlight the complex interactions between climatic variables and water yield dynamics in the Gilgel Gibe Watershed. Land use change further complicates this scenario. The conversion of natural vegetation into agricultural land and urban areas has significantly modified the hydrological processes of the watershed. Deforestation and urbanization have intensified surface runoff and decreased surface water flow, leading to a reduction in water yield. For instance, the expansion of irrigated agriculture, particularly in upper catchment areas, has increased water abstraction, which in turn has decreased surface runoff and water availability downstream. Data from the Gilgel Gibe Watershed show a decrease in surface runoff from 15.78% in 2021 to 15.28% in 2022 and a reduction in water yield from 1.22% in 1993 to 0.83% in 2022, with the lateral flow decreasing from 8.55% in 1993 to 7.0% in 2022. These trends underscore the impact of climate and LUC on hydrological dynamics.

The synergistic effects of land use and climate change have increased fluctuations in water yield. Changes in land cover not only alter surface runoff but also influence local microclimatic conditions, affecting precipitation patterns and evapotranspiration rates. This interaction leads to complex, nonlinear responses in water yield, with potential threshold effects under certain land-use and climate change scenarios. Changes in water yield significantly affect ecosystem services, particularly freshwater provision, agricultural productivity, and aquatic habitats. Decreased water yield has led to reduced water availability for irrigation, threatening food security and affecting aquatic habitats due to lower river discharge. These changes have further affected biodiversity by altering the ecological balance, particularly in wetlands and riparian zones, which depend on stable water levels for their functioning. The ecological health of the Gilgel Gibe watershed is increasingly threatened by climate change and land-use changes, posing serious long-term risks to both human and environmental well-being.

Effective water management strategies are crucial to mitigate these impacts. Ethiopia’s National Water Policy, which guides the country’s management of water resources, has several gaps that limit its effectiveness in addressing emerging challenges related to water yield and ecosystem health, especially in the context of climate change and land use changes. One significant issue is inadequate integration of climate change adaptation strategies. Climate change leads to altered precipitation patterns, rising temperatures, and more frequent extreme weather events such as droughts and floods, yet the policy does not sufficiently address these impacts on water resources or include measures to enhance climate resilience in water management. Additionally, the policy lacks a comprehensive approach for managing ecosystem services that are essential for regulating water quality, preventing soil erosion, and maintaining water flow. While it acknowledges the importance of water conservation, there is little emphasis on how land-use changes, such as deforestation, urbanization, and agricultural expansion, affect water yield. The policy does not fully embrace ecosystem-based approaches that protect and restore vital ecosystems such as wetlands, forests, and grasslands, all of which are crucial for improving water security.

Moreover, Integrated water resource management (IWRM) approaches that incorporate sustainable land management practices, climate-resilient infrastructure, and community engagement are crucial for mitigating the impacts of land use and climate variability on water yield. Ethiopia’s National Water Policy and Water Sector Development Program (WSDP) aims to enhance water resource management and infrastructure development. When compared to global standards such as the United Nations Sustainable Development Goals (SDGs), Ethiopia’s National Water Policy reveals several discrepancies. For instance, SDG 6 (Clean Water and Sanitation) emphasizes the need for integrated water resources management, equitable distribution, and the protection of water ecosystems. Although the policy highlights water access and infrastructure, it does not adequately promote sustainable water resource management across sectors, particularly regarding climate resilience and ecosystem conservation. Furthermore, SDG 13 (Climate Action) stresses the importance of addressing climate change impacts; however, the policy does not sufficiently consider the long-term implications of climate change on water systems. Incorporating climate adaptation strategies could help Ethiopia align its water management framework with global climate goals and improve its resilience to climate impacts such as extreme weather events and shifting rainfall patterns.

Additionally, SDG 15 (Life on Land) advocates for sustainable land-use practices that protect and restore ecosystems vital for maintaining water yield. The policy’s lack of focus on ecosystem-based approaches and land-use changes contradicts SDG 15, as ecosystems play a critical role in regulating water resources119. evaluated the WSDP’s impact on water availability in the Gilgel Gibe Watershed, finding improvements in water supply infrastructure but highlighting ongoing challenges in equitable distribution and conflict management. International agreements, such as the United Nations Sustainable Development Goals (SDGs), also influence water resource management in watersheds120. examined Ethiopia’s alignment with SDGs, emphasizing the need for policy coherence and effective implementation strategies to achieve water-related targets. Understanding the combined effects of climate change, land use change, and hydrological discharge on surface flow and water yield is essential for creating effective adaptation strategies. Adopting integrated water management practices aligned with Ethiopia’s national policies is key to ensuring sustainable water use and safeguarding livelihoods in the Gilgel Gibe Watershed.

Discussion

The combined effect of LULC changes and climate change on water yield in the Gilgel Gibe Watershed has been extensively studied, revealing complex interactions that have a considerable impact on hydrological processes. These collective effects have been shown to intensify water yield fluctuations in watersheds121. utilized hydrological modeling techniques to assess integrated impacts, revealing that increased deforestation and expansion of agricultural areas, coupled with rising temperatures and altered precipitation patterns, have led to increased hydrological variability. This intensification effect was particularly evident during extreme weather events, where increased runoff and decreased infiltration rates resulted in increased flood hazards and a reduced baseflow contribution to streamflow. Furthermore122, emphasized the importance of considering land–atmosphere interactions to understand the synergetic effects on water yield. Changes in land cover composition affect surface runoff and soil moisture dynamics as well as local microclimatic conditions, influencing precipitation patterns and evapotranspiration rates. The response mechanisms that produce nonlinear water yield patterns are frequently caused by the complex interactions between LULC changes and climatic variability.

These interactions can result in threshold effects, in which small changes in one element have a disproportionately large impact on water yield. For example, deforestation and agricultural intensification cause a gradual increase in surface runoff, compromising soil moisture retention; however, beyond a certain threshold, they can result in severe erosion, sedimentation, and loss of infiltration capacity, thereby reducing groundwater recharge. These nonlinearities can be enhanced in areas with existing climate pressures, such as prolonged dry seasons or extreme rainfall events, where reaction loops exacerbate these effects123. demonstrated that land-use changes negatively impact ecosystem services, especially the trade-off between provisioning and regulating services. While expanding provisioning services may seem beneficial in the short term, continuing current land use practices threaten the long-term supply of essential regulating services. This aligns with our findings for the Gilgel Gibe Watershed, where climate and land use changes have significantly degraded critical ecosystems, including flora and fauna. Addressing these land use dynamics is crucial for ensuring the sustainability of ecosystem services in watersheds.

Moreover, studies by124 in the Ganges–Brahmaputra-Meghna (GBM) Basin in Bangladesh and Nepal further elucidated the intricate relationship between LULC change and climate variability in water yield. In the GBM Basin, deforestation, agricultural expansion, and urbanization, combined with climate fluctuations, have led to increased surface runoff, reduced groundwater recharge, and altered water yield, affecting downstream ecosystems and livelihoods. Similar patterns were observed in the Gilgel Gibe Watershed, where human-induced changes in land use and climate variability exacerbated runoff, reduced groundwater recharge, and influenced water availability. These findings emphasize the need for integrated watershed management practices, including afforestation, reforestation, and water-efficient irrigation, to mitigate land degradation and enhance water yield.

Moreover, modeling and monitoring hydrological changes, as suggested in the GBM Basin study, can aid in predicting water yield fluctuations in the Gilgel Gibe Watershed, helping develop adaptive strategies for sustainable water management. This study underscores the importance of addressing both environmental and socioeconomic factors to ensure the long-term sustainability of water resources in the Gilgel Gibe region, benefiting both the ecosystem and local communities125. highlighted the direct consequences of reduced water yield on agricultural productivity, which is a primary source of livelihood for many in watersheds. Decreased water availability for irrigation owing to altered hydrological regimes negatively affects crop yields, leading to food insecurity and economic losses for farmers. Moreover, changes in water quantity and quality affect the ecological integrity of the watershed, influencing biodiversity, habitat suitability, and ecosystem functioning126.

The ecological consequences of water yield in the Gilgel Gibe Watershed are multifaceted and have significant implications for the health and functioning of the local ecosystem. Several studies have investigated these consequences, shedding light on the intricate relationships between water availability, LULC dynamics, and ecological processes. One of the primary ecological consequences of altered water yield in the Gilgel Gibe Watershed is its impact on vegetation dynamics. Research conducted by127 found that changes in water availability driven by both LULC changes and climate variability directly influence vegetation distribution and composition. Reduced water yield resulting from increased land degradation and deforestation can lead to the loss of vegetation cover and biodiversity, particularly in riparian zones and ecologically sensitive areas.Riparian zones, which rely on consistent water availability, have suffered severe biodiversity loss, with native plant species being replaced by invasive or drought-tolerant species, resulting in ecosystem homogenization. In addition, climatic variability has exacerbated these difficulties by introducing erratic rainfall patterns, prolonged droughts, and increasingly intense weather events, all of which have decreased water availability and accelerated soil erosion.

Adaptive strategies for water resource management in the Gilgel Gibe Watershed are essential to counter the effects of LULC change and climate variability on water yield. Key approaches include adopting sustainable land management practices to reduce soil erosion, improve infiltration, increase the water retention capacity, and ensure long-term water availability120. highlighted the effectiveness of terracing and agroforestry in reducing soil erosion and increasing groundwater recharge in similar watersheds. By adopting these practices, farmers can mitigate the loss of fertile topsoil and optimize water-use efficiency, thereby sustaining agricultural productivity and maintaining water yield. National water resource management policies play a crucial role in shaping the water yield dynamics of the Gilgel Gibe Watershed. Several studies have investigated the implications of these policies on regional water resources. For instance128, assessed the effectiveness of Ethiopia’s Water Sector Strategy (WSS) in managing water resources in the Gilgel Gibe Watershed. WSS emphasizes integrated water resource management (IWRM), highlighting stakeholder involvement, water distribution, and ecosystem conservation. However, this study identified challenges in implementing the WSS locally, particularly in addressing competing water demands from agriculture, hydropower, and domestic users.

Furthermore, the Ethiopian government implemented the Water Sector Development Program (WSDP) to enhance water resource management and infrastructure development129. evaluated the impact of the WSDP on water availability and access in the Gilgel Gibe Watershed. The findings revealed improvements in the water supply infrastructure, such as reservoirs and irrigation systems, which have increased water availability for agricultural activities. One notable example is the Tana-Beles Integrated Watershed Management Project, supported by the World Bank, which has successfully introduced community-based sustainable land management practices. These efforts have resulted in a significant 20% increase in groundwater recharge within pilot areas, demonstrating the project’s success in enhancing water resource sustainability130. This finding aligns with our study of the Gilgel Gibe watershed and other hydrological zones where Integrated Water Resource Management (IWRM) strategies have been employed to balance competing water uses and improve water yield modeling. Additionally, the Ministry of Water and Energy has actively engaged in watershed rehabilitation and afforestation programs aimed at improving water yield and alleviating the impacts of land-use changes. The Tana River Basin in Kenya serves as a relevant case study, illustrating that conservation agriculture and reforestation policies have increased surface water retention and enhanced dry-season flows131. This emphasizes the importance of integrating land management policies with hydrological models, which is consistent with our findings on sustainable water management. We recommend applying similar approaches to other watersheds across East Africa in order to bolster sustainable water and land management. Furthermore, a study conducted in the three-river headstream region of China utilized the InVEST model and machine learning to evaluate water conservation efforts132. The results indicated that land use change improved water yield; however, these improvements were less significant than the impacts of climate. The study also provides policy recommendations, including the implementation of stricter watershed protection laws and incentives for land conservation. Despite these advancements, several challenges hinder policy implementation, especially in the ethnically diverse regions of Ethiopia and similar areas. Even when effective policies are established, their execution is often hindered by economic, social, and institutional factors. A primary limitation is the lack of financial and technical resources, along with insufficient funding and expertise, which impede large-scale watershed rehabilitation projects. To address these challenges, it is advisable to leverage international climate funds, such as the Green Climate Fund, and foster public–private partnerships to finance sustainable watershed management.

Another issue is institutional fragmentation, as overlapping responsibilities among federal, regional, and local agencies can delay their implementation. To resolve this problem, it is crucial to strengthen cross-sectoral coordination and ensure that hydrological data are shared among the relevant stakeholders. Additionally, community resistance to land use policies poses significant challenges. Policies that restrict deforestation or enforce land conservation may face opposition from the local communities that rely on these resources. Therefore, implementing community-based watershed management programs, in which local participation has successfully led to small-scale reforestation efforts, is vital for effective execution133.

Moreover, the valuable insights provided by datasets such as Landsat, MODIS, and climate data from NASA POWER, challenges persist in ensuring an equitable water distribution and managing conflicts among different water users. One of the key limitations of this study lies in the accuracy and resolution of the remote-sensing data. Although these datasets offer important information, their spatial and temporal resolutions may have limitations in the in-depth capture of localized variations in water yield, particularly in heterogeneous landscapes. For example, Landsat’s 30-m spatial resolution has some limitations in detecting rapid land-use changes such as urbanization, and climate data with a 4 km resolution may miss microclimates within complex watersheds. Addressing these limitations requires the development of standardized frameworks for data integration and the incorporation of socioeconomic factors into hydrological models. This would not only improve the reliability of predictions but also enhance their applicability to sustainable water resource management. Finally, this study highlights the importance of adopting integrated watershed management practices to ensure the sustainable use of water resources and enhance water yield in the Gilgel Gibe Watershed. By combining ecological, social, and economic strategies, such as soil conservation, reforestation, community engagement, and strong governance, these approaches can effectively tackle the challenges of LULC changes and climate variability, enhancing the resilience and sustainability of the watershed’s water yield.

Impacts of land use and climate change on ecology

The Gilgel Gibe I watershed, an important hydrological resource in Ethiopia, faces substantial ecological challenges as a result of land use change and climate variability. Deforestation, agricultural expansion, and urbanization have significantly affected the ecosystems of the region, resulting in habitat fragmentation and the extinction of native forest species. For example, deforestation, which is frequently driven by the demand for cropland and fuelwood, disrupts habitat connectivity and endangers animals, resulting in population decrease and, in some cases, extinction. Increasing temperatures exacerbate these stresses by increasing evapotranspiration, worsening water scarcity during dry spells, and putting additional strain on native plant and animal species. Similar studies in the Awash Basin have found a 25% decrease in forest cover over two decades, resulting in the extinction of some plant species essential to local biodiversity134.

The aquatic habitats and riparian biodiversity are also threatened. Reduced plant cover increases surface runoff, accelerates soil erosion, and reduces groundwater recharge, which results in sedimentation in water bodies. Sedimentation not only degrades water quality but also jeopardizes hydropower production, irrigation, and household water supply. For instance, sedimentation in the Blue Nile Basin has reduced the reservoir capacity by 15% over ten years135. Meanwhile, altering rainfall patterns and escalating droughts and floods caused by the Inter-Tropical Convergence Zone (ITCZ) destabilize wetland ecosystems and water availability. Wetland-dependent plant species in the Omo-Gibe Basin lost 40% of their coverage during drought years, demonstrating the severe ecological impact of climate variabilityo136. Addressing these difficulties requires the implementation of comprehensive mitigation methods. Reforestation, sustainable land management, and community-based conservation can all contribute to ecosystem restoration and protection of native species. Furthermore, adopting water-efficient irrigation techniques, encouraging rainwater gathering, and minimizing urban water usage are critical for saving water resources. Adaptation methods, such as introducing drought-tolerant crops and improving flood control systems, can help reduce the vulnerability of both ecosystems and human societies. Understanding the complex relationships among land use, climate change, native species, and water resources might help in designing more effective and sustainable management techniques for the Gilgel Gibe watershed.

Future direction

The combined influence of LUCC and climate change on water yield modeling is increasingly being recognized as a critical area of research, owing to its profound impact on hydrological processes within watersheds. Although LUCC plays a significant role in altering water yield and ecosystem services, these changes are not solely driven by land-use modifications. Variations in land use can affect water yield and runoff; however, these impacts are often intertwined with changes in surface runoff patterns, making it challenging to attribute hydrological changes exclusively to LUCC and climate change. This study employed the InVEST model’s annual water yield module, which is based on the water balance principle, owing to its simplicity and effectiveness. However, the accuracy of water yield estimates depends on several factors. For instance, employing a longer time frame (e.g., 30 years) to align climate data with land-use maps can enhance the model precision. Input variables, such as plant-available water content and potential evapotranssaspiration, were calculated using empirical formulas, and the choice of these formulas can significantly influence the results. The seasonal constant “Z,” calibrated to local water resource conditions, further optimizes the model but may require region-specific adjustments. The"Z"calibration constant plays a key role in enhancing the accuracy of hydrological models by adjusting parameters such as evapotranspiration and soil moisture based on observed data. In this study, the"Z"constant was used to reconcile inconsistencies between the model predictions and real-world observations, ensuring that the model better reflects local conditions. This calibration improves the reliability of the predictions for water yield, runoff, and other hydrological variables. Additionally, spatial interpolation methods used for inputs, such as precipitation and evapotranspiration, can introduce errors, emphasizing the importance of selecting optimal empirical formulas and interpolation methods tailored to local conditions. Although parameter calibration reduces some errors, discrepancies between the observed and simulated runoff remain. In this study, Kriging interpolation was employed due to its ability to capture spatial heterogeneity, which is crucial in hydrological and environmental analyses with complex spatial relationships. By incorporating variograms to model spatial autocorrelation, Kriging provides high accuracy, particularly in areas with dense and well-distributed data. This method not only offers precise predictions but also quantifies the prediction uncertainty. However, it is computationally intensive and requires expertise in spatial statistics to select suitable models and parameters for obtaining optimal results. Future research could integrate alternative methodologies, such as complementary or total differential methods, to address these gaps. These approaches may provide deeper insights into the complex interactions between LUCC, climate variability, and hydrological responses, thereby advancing the precision of water yield modeling and improving watershed management strategies.

Conclusions

An assessment of the combined impacts of LULC changes and climate variability in the Gilgel Gibe watershed from 1993 to 2023 revealed significant alterations in water yield. This study highlights key trends, including an increase in mean temperature to 34 °C and a reduction in monthly precipitation from 500 to 100 mm. These changes are closely associated with land-use changes, particularly the expansion of cultivated land and reduction of vegetation cover. These findings show that declining precipitation and rising temperatures increase evaporation rates, resulting in a significant decrease in water yield. Importantly, LULC changes have shown a more pronounced impact on surface runoff than climate variability, influencing both groundwater recharge and the volume of surface runoff, as well as the amount of water lost through evapotranspiration. In general, this research emphasizes the critical need to integrate both climate and land use dynamics into water resource planning. Understanding the role of human-induced LULC changes is crucial for accurately assessing hydrological responses and for formulating effective watershed management strategies. The calibrated model developed in this study provides a robust framework for future assessments, enabling the evaluation of surface water flow and storage under various scenarios and guiding decision-making to ensure the long-term sustainability of water yield in the Gilgel Gibe watershed. Additionally, our previous research on climate change impacts in the Gilgel Gibe watershed, which incorporated the Shared Socioeconomic Pathways (SSP) and Representative Concentration Pathways (RCP), revealed that future precipitation is projected to decline under RCP 8.5 (−32.2%) and SSP 4.5 (−88.8%) for the period 2024–2049, with further reductions expected between 2050 and 2099. Conversely, temperature projections indicate an increase under RCP 4.5 (TMAX + 0.67 °C, TMIN + 1.11 °C) and RCP 8.5 (TMAX + 0.25 °C), with higher emissions under SSP 4.5 and SSP 8.5 further intensifying warming137. While RCP-based scenario modeling is not the primary focus of the present study, our previous findings support the likelihood of future precipitation declines and rising temperatures, both of which have significant implications for water yield trends in the region.