# Land Use/Land Cover Change Impact on Hydrological Process in the Upper Baro Basin, Ethiopia

Authors: Tewodros Getu Engida; Tewodros Assefa Nigussie; Abreham Berta Aneseyee; John Barnabas
Journal: Applied and Environmental Soil Science
 (2021)

---

## Abstract

Understanding the hydrological process associated with Land Use/Land Cover (LU/LC) change is vital for decision-makers in improving human wellbeing. LU/LC change significantly affects the hydrology of the landscape, caused by anthropogenic activities. The scope of this study is to investigate the impact of LU/LC change on the hydrological process of Upper Baro Basin for the years 1987, 2002, and 2017. The Soil Water Assessment Tool (SWAT) model was used for the simulation of the streamflow. The required data for the SWAT model are soils obtained from the Food and Agriculture Organization; Digital Elevation Model (DEM) and LU/LC were obtained from the United States Geological Survey (USGS). The meteorological data such as Rainfall, Temperature, Sunshine, Humidity, and Wind Speeds were obtained from the Ethiopian National Meteorological Agency. Data on discharge were obtained from Ministry of Water, Irrigation and Electricity. Ecosystems are deemed vital. Landsat images were used to classify the LU/LC pattern using ERDAS Imagine 2014 software and the LU/LC were classified using the Maximum Likelihood Algorithm of Supervised Classification. The Sequential Uncertainty Fitting (SUFI-2) global sensitivity method within SWAT Calibration and Uncertainty Procedures (SWAT-CUP) was used to identify the most sensitive streamflow parameters. The calibration was carried out using observed streamflow data from 01 January 1990 to 31 December 2002 and a validation period from 01 January 2003 to 31 December 2009. LU/LC analysis shows that there was a drastic decrease of grassland by 15.64% and shrubland by 9.56% while an increase of agricultural land and settlement by 18.01% and 13.01%, respectively, for 30 years. The evaluation of the SWAT model presented that the annual surface runoff increased by 43.53 mm, groundwater flow declined by 27.58 mm, and lateral flow declined by 5.63 mm. The model results showed that the streamflow characteristics changed due to the LU/LC change during the study periods 1987–2017 such as change of flood frequency, increased peak flows, base flow, soil erosion, and annual mean discharge. Curve number, an available water capacity of the soil layer, and soil evaporation composition factor were the most sensitive parameters identified for the streamflow. Both the calibration and validation results disclosed a good agreement between measured and simulated streamflow. The performance of the model statistical test shows the coefficient of determination (R2) and Nash–Sutcliffe (NS) efficiency values 0.87 and 0.81 for calibration periods of 1990–2002 and 0.84 and 0.76 for the validation period of 2003 to 2009, respectively. Overall, LU/LC significantly affected the hydrological condition of the watershed. Therefore, different conservation strategies to maintain the stability and resilience of the ecosystem are vital.


---

## Body

## 1. Introduction

The response of hydrologic circulation is closely related to land use planning and management [1, 2]. Land Use/Land Cover (LU/LC) change is one of the factors that directly affect the watershed hydrological cycle [3–5]. It has been widely accepted that human activities cause LU/LC change and have imposed a great impact on the hydrological processes and water resources of the watershed [6] For example, a study [1] showed that the water provision and the hydrological process decreased as consequence LU/LC change that was aggravated by significant increasing population pressure and development in Pennar Basin in India. This can be either by increasing the water yield or by eliminating the flow in some circumstances, thereby increasing the sediment load and decreasing the groundwater [7, 8].LU/LC changes and population growth are the most common problems in developing countries such as Ethiopia since their economic development mainly depends on agriculture [9]. In the past 40 years, the increases in human activities resulted in the expansion of agricultural land, extraction of timber, and urbanization, thereby resulting in deforestation in southwest Ethiopia [10, 11], while these changes in land cover affect the local watershed hydrological cycle and flood vulnerability of various subwatersheds [12]. The modifications/conversions of natural vegetation and physical soil conditions are usually the principal cause of changes in rainfall-runoff characteristics of the local catchment, which consequently change the river flow regimes [13, 14]. Several studies show that the changes in vegetation cover, i.e., deforestation, lead to an increase in water yield and sedimentation [15, 16].Watershed studies can provide direct evidence of LU/LC change impacts on runoff [17]. For these reasons, hydrological models are becoming important for studying LU/LC change effect on the hydrological cycle in a landscape [18, 19]. Several studies were carried out on the hydrology of the watershed by utilizing LU/LC data in different regions of Ethiopia using the Soil and Water Assessment Tool (SWAT). They indicated an increase in the wet seasonal flow, surface runoff, and water yield [20–23]. However, studies on the impact of the LU/LC dynamic on the hydrology regime of the Upper Baro Basin are limited.The Upper Baro Basin is hurt from human-induced degradation mainly due to large-scale agricultural expansion, which leads to LU/LC dynamism. The impacts of such changes on local hydrology are poorly understood. Hence, systematic analyses and understanding the impacts of LU/LC changes on the watershed hydrology are important for planning water resources management and maintaining the sustainable flow of water to the Baro-Akobo River. This study also helps decision-makers know where to invest to conserve the hydrological services.Therefore, the objectives of the study are as follows: (1) to assess the LU/LC change over the past 30 years on the hydrological response of Upper Baro Basin and (2) to check the applicability of SWAT model for the study area.


## 2. Materials and Methods

### 2.1. The Study Area

The Upper Baro Basin is found in southwestern Ethiopia located at 766 km away from Addis Ababa, the capital city of Ethiopia, between latitudes 7.44° and 9.41°N and longitudes 34.52° and 36.31°E, covering an area of 23462 km2. The altitude ranged from 390 m to 3266 m above sea level (Figure 1).Figure 1
Map of the study area.The Upper Baro basin has a long rainy season in the months of May-September. The area indicates a monomodal, bimodal, and triple-modal rainfall pattern. The annual average rainfall is varying from about 1163.03 mm to over 2258.36 mm with a monthly maximum rainfall record of 302.79 mm in August, while relative mean monthly rainfall for dry seasons occurs between December and March (Figure2).Figure 2
Average monthly areal rainfall distribution throughout 1987–2017.The annual rainfalls show pronounced annual and seasonal fluctuations (Figure3) and the maximum temperature of Upper Baro Basin ranges between 27.7°C (in the wet season) and 37.2°C (in the dry season) and occurs in June to August, while the minimum falls between 10.8°C and 13.7°C in the year.Figure 3
Annual areal rainfall distribution (mm) for the period of 1987–2017.The topography of the Upper Baro Basin is characterized by variation from flat to mountainous terrain and falls in the altitude range of 390 m to over 3000 m a.m.s.l. with 42% in between 1000 m and 2000 m a.m.s.l. The eastern two-thirds of the basin area lies between 1000 m and 2400 m a.m.s.l. and a gently sloping plain lies in the west between 380 m and 500 m a.m.s.l. Dystric nitisols and dystric gleysols are the dominant soil type in the study area.


### 2.2. SWAT Model Data Types and Analysis

#### 2.2.1. Meteorological Data

The SWAT model requires daily meteorological (Weather) data such as daily precipitation, minimum and maximum air temperature, relative humidity, wind speed, dew point, and daily sunshine. The standard deviation of each parameter was calculated. Thiessen method is used to adjust for nonuniform meteorological station distribution by proportion to the area which is closer. Eleven meteorological stations were selected (Table1) based on the quantity, quality, period, consistency, homogeneity, and their uniform distribution within and around the Upper Baro Basin. These data were obtained from the National Meteorological Agency of Ethiopia (NMAE) over the period (1987–2017). The missing values were estimated from other stations around the missed record station by using both arithmetic mean method and normal ratio method and also consistency test which shows double mass.Table 1
Summary of the meteorological stations.
StationCoordinateAltitude (m)Station classPeriod of dataLongitudeLatitudeGambella34.58338.25430C21983–2017Mettu35.566678.28331711C21981–2017Dembi Dolo34.88.51671850C21982–2017Ayira35.559.11555C21987–2017Gatira36.27.9832358C21987–2017Masha35.46677.752282C11981–2017Gore35.53338.13332033C11981–2017Seko Humbi34.98338.71661860C21987–2017Gimbi35.78339.16661970C21983–2017Bure35.18.23331750C21981–2017Alge35.66678.53331880C21987–2017The SWAT weather generator model (WGEN) was used to match the access table and lookup table. The daily precipitation and temperature of all gauging stations were prepared in comma-delimited (.csv) format. Solar radiation, relative humidity, and wind speed data were available only for principal Gore and Masha stations (synoptic stations).


#### 2.2.2. Soil Data

Seventy soil types have been identified at Upper Baro Basin (Table2). The spatial soil data were obtained from the Ethiopian Ministry of Water, Irrigation, and Electricity (MWIE), which is prepared according to FAO soil classification. The most dominant soil in the basin is Dystric nitisols and Dystric gleysols.Table 2
Soil type’s area coverage and SWAT code.
ValueSoil typesArea (km2)Area covered (%)SWAT code1Dystric nitisols14151.7660.32DYNITISOLS2Dystric gleysols5423.5223.12DYGLEYSOLS3Orthic acrisols326.871.39ACRISOLS4Chromic luvisols120.450.51CHLUVISOLS5Orthic solonchaks778.063.32SOLONCHAKS6Calcic xerosols363.201.55CLXEROSOLS7Leptosols255.951.09LEPTOSOLS8Iatric cambisols1041.604.44EUCAMBISOLS9Eutric nitisols108.430.46EUNITISOLS10Eutric fluvisols276.991.18EUFLUVISOLS11Dystric cambisols58.400.25DYCAMBISOLS12Gypsic yermosols43.390.18GYYERMOSOLS13Dystric fluvisols231.440.99DYFLUVISOLS14Cambisols2.540.01CAMBISOLS15Calcic fluvisols92.400.39CLFLUVISOLS16Calcic cambisols167.830.72CLCAMBISOLS17Calcaric flubisols19.170.08CLFLUBISOLSSWAT model requires physical and chemical soil properties in the simulation of a hydrological component such as available water content, soil texture, hydraulic conductivity, bulk density, and organic carbon content for the different layers of each soil type where shown. The major soils of the Upper Baro Basin are as shown below table and a lookup table was prepared for each type of soil.


#### 2.2.3. DEM

DEM with a resolution of 30 m was downloaded from United States Geological Surveyhttps://earthexplorer.usgs.gov/. The DEM is one of the essential inputs required by SWAT to delineate Upper Baro Basin (Figure 4), define the stream network, and determine sub-basin parameters such as slope area and slope length.Figure 4
DEM map (meters a.m.s.l.) of Upper Baro Basin.


#### 2.2.4. Slope Classes in Subwatershed

To develop the hydrological response unit (HRU) in the SWAT model, the slope is essential. It generates from the resolution of 30 m ∗ 30 m DEM for the study area. Slope classification may be single class or multiclass. For this study, the slope option (an option for considering different slope classes for HRU definition) was selected. Therefore, the slope class in this study was classified into four classes. According to [22], slope classification was used to account for lower range in hydrological modelling. Depending on the slope of 0–5%, 5–10%, 10–15%, 15–20%, and >20% were selected for HRU creation of the study Upper Baro Basin. Finally, HRU definition analysis in SWAT helps to load LULC and soil type projects.


#### 2.2.5. Land Use/Land Cover (LU/LC) Analysis

Landsat Satellite images of the years 1987 (Landsat 5), 2002 (Landsat 7) and 2017 (Landsat 8) were obtained from United States Geological Survey (USGS) (https://www.glovis.USGS.gov) (Table 3). The multidate multisensor satellite imagery was collected for successive cropping for each study year. Dataset selections were in dry season fixed with the lowest percent or zero monthly cloud cover. Therefore, all images are cloud-free (0%) in the study area.Table 3
Landsat data information for the study area.
Reference dataPath/rowAcquisition dataImageSensorResolution1987170/05422/01/1987Landsat 5TM30 m ∗ 30 m170/055171/0542002170/05423/01/2002Landsat 7ETM+30 m ∗ 30 m170/055171/0542017170/05427/01/2017Landsat 8OLF15 m ∗ 15 m170/055171/054The images were already terrain corrected projected to the Universal Transverse Mercator (UTM) of WGS84 zone 37 N. They were radiometrically calibrated to reflectance value and the Quick Atmospheric Correction algorithm was applied in ERDAS software. Radiometric and geometric corrections were applied to the removal of sensor noise, haze, correction data loss, and missing lines due to a solar position, satellite calibration [24, 25]. A supervised image classification method was employed using Maximum Likelihood Algorithm classification. For supervised classification, the Region of Interest (ROI) was made a signature for different LU/LC categories. The reference data collected using a GPS receiver were first converted to a vector file, then the region of interest (ROI) was made in ERDAS software. Using the ROI, the spectral signature of each LU/LC type has been extracted. Further investigations were carried out by taking these sample points in the land use types using a handheld GPS receiver. Based on this analysis and the existing land use distribution information, the area was categorized into eight broad LU/LC classes: cultivated land, woodland, forest land, grazing land, bare land, shrubland, water bodies, and built-up (Table 3).Overall accuracy and Kappa coefficient, the classify maps will check classification accuracy. According to [26, 27], a good overall classification map accuracy is accepted if and only it possesses an accuracy of at least 85%. The LU/LC was an input of the SWAT model to describe the HRU of the subwatershed. The SWAT Model has predetermined four-letter codes for each LU/LC category (Table 4). These codes were used to link the LU/LC map of the study area with the SWAT land database. To make it compatible, LU/LC with the SWAT model lookup table was prepared for each type of LU/LC.Table 4
LU/LC class SWAT code.
ValueLU/LC classSWAT databaseSWAT code1Forest landForest evergreenFRSE2WoodlandForest-mixedFRST3GrasslandRange-grassesRNGE4ShrublandRange-brushRNGB5SettlementResidential-medium densityURMD6Agricultural landAgricultural land-genericAGRL7Bare landBarrenBARR8WaterbodyWaterWATR


### 2.3. SWAT Model Description

SWAT is a physically based semidistributed continuous time-scale hydrological model, which works on a daily time step. This model can simulate runoff, sediment, nutrients, pesticide, and bacterial transport from agricultural watersheds [28]. It simulates the hydrological cycle parameters based on the water balance represented in equation (1) within the watershed [29]. More details on the information about the SWAT model can be found in the SWAT user manual [29].The hydrologic cycle simulated by the SWAT is based on the water balance equation [29] as follows:(1)SWt=SWo+∑i=1tRi−Qi−ETi−Pi−QRi,where SWt and SWo are the final and initial soil water content, respectively, (mm), Qi = daily surface runoff (mm), Ri = daily rainfall (mm), ETi = daily evapotranspiration (mm), Ri = daily lateral flow (mm), and Pi = daily percolation (mm).


### 2.4. Peak Runoff Rate

The peak runoff rate is a result of the erosive power of a storm and is used to predict sediment loss. SWAT calculates the peak runoff rate with a modified rational method for each HRU using the following equation:(2)Qpeak=αtc∗Qsurf∗Atconc,where Q peak is peak runoff rate (m3/s), αtc the fraction of daily rainfall that occurs during the time of concentration, Qsurf is the surface runoff (mm), A is the sub-basin area (km2), and tconc is the time of concentration (hr).


### 2.5. Surface Runoff Generation

The CN method was initially developed for small agricultural watersheds and the CN varies nonlinearly with the moisture content of the soil. It drops to zero as the soil approaches the wilting point and increases to be near 100 as the soil approaches saturation, with higher CNs associated with higher runoff potential watershed. This method is widely used [28].According to USDA, 1985 [30], this can be mathematically expressed using the following equation:(3)Qs=Ra−Ia2Ra−Ia+S,ForRa>Ia,Q=0forRa≤Ia,where Qs is the accumulated runoff (mm), Ra is the rainfall depth for the day (mm), and(4)S=25.4∗100CN−10,where CN is the curve number for the day and its value is the function of land use practice, soil permeability and soil hydrologic group. Conceptual framework showing the components and relationships have been used as a framework for the analysis of hydrology in research (Figure 5).Figure 5
General conceptual methodology used for modelling hydrology in the Upper Baro Basin.


### 2.6. SWAT Model Simulation, Sensitivity Analysis, Calibration, and Validation

The simulation result cannot be directly used for further analysis. Instead, the ability of the model to sufficiently predict the constituent streamflow should be evaluated through sensitivity analysis, model calibration, and model validation [31].


#### 2.6.1. Sensitivity Analysis

Fifteen hydrological parameters related to streamflow were selected for sensitivity analysis in the study area. The sensitivity of each parameter was identified using theT test and P values. P values were used to determine the significance of the sensitivity. A value close to zero has more significance. T test provides a measure of sensitivity and the most sensitive parameters usually appear with P values less than the alpha level of 0.05 (larger absolute values are more sensitive) [32].


#### 2.6.2. Model Calibration and Validation

The SUFI-2 algorithm was applied and frequently used in this study of calibration for the SWAT model at a large scale due to its easy application and the reduced number of model runs needed to achieve a decent prediction [33, 34]. In this study, automatic calibration was made every month from January 1, 1990, to December 31, 2002, until the average simulated value came closer to the measured value. Automatic calibration makes use of a numerical algorithm to increase the performance of the model and to optimize the numerical objective functions. After the simulation result for the calibration, the period had fulfilled the above statistical criteria, validation was performed for an independent period of records from January 1, 2003, to December 31, 2009, and the results were compared against an independent set of measured Upper Baro Basin.


#### 2.6.3. Model Performance Evaluation

Model performance was carried out to verify the robustness of the model to simulate hydrological processes. A model performance framework proposed by [35] was used in this study.Nash and Sutcliffe simulation efficiency (NSE) indicates the degree of fitness of observed and simulated data, and it is calculated using the following equation:(5)NSE=1−∑i=0NQobs−QSim2∑i=0NQobs−Q¯obs2,where N = number of compared values, Qobs is the observed data, Q¯obs is the observed mean, and Qsim is the simulated data.The NSE indicates how well the plot of observed versus simulated value fits the 1 : 1 line. The closer the model efficiency is to 1, the more accurate the model, and if it is found between 0 and 1, it indicates deviations between measured and predicted values. If NSE is negative, predictions are poor, and the average value of output is a better estimate than the model prediction [36].Coefficient of determination (R2) measures the ability of a model to predict or explain an outcome in the linear regression setting (Equation (6)). R2 ranges from 0 to 1, with higher values indicating less error variance, and typically values greater than 0.5 are considered acceptable [35, 37, 38].(6)R2=∑Xi−XavYi−Yav∑Xi−Xav2∑Yi−Yav2,where Xi is measured value, Xav is average measured value, Yi is simulated value, and Yav is the average simulated value.Percent bias (PBIAS) measures the average difference between the simulated and measured values for a given quantity over a specified period (usually the entire calibration or validation period) and it is calculated using the following equation [39, 40]:(7)PBIAS=∑Yi−∑XiXi,where Xi is measured value and Yi is simulated value.Model evaluation classifies as follows:(i)
0.75 < NSE < 1.00, PBIAS < ±10—better predictor(ii)
0.65 < NSE ≤ 0.75, ±10 < PBIAS < ±15—good(iii)
0.50 < NSE ≤ 0.65, ±15 < PBIAS < ±25—satisfactory(iv)
(NSE ≤ 0.50, PBIAS >±25)—unsatisfactory


## 2.1. The Study Area

The Upper Baro Basin is found in southwestern Ethiopia located at 766 km away from Addis Ababa, the capital city of Ethiopia, between latitudes 7.44° and 9.41°N and longitudes 34.52° and 36.31°E, covering an area of 23462 km2. The altitude ranged from 390 m to 3266 m above sea level (Figure 1).Figure 1
Map of the study area.The Upper Baro basin has a long rainy season in the months of May-September. The area indicates a monomodal, bimodal, and triple-modal rainfall pattern. The annual average rainfall is varying from about 1163.03 mm to over 2258.36 mm with a monthly maximum rainfall record of 302.79 mm in August, while relative mean monthly rainfall for dry seasons occurs between December and March (Figure2).Figure 2
Average monthly areal rainfall distribution throughout 1987–2017.The annual rainfalls show pronounced annual and seasonal fluctuations (Figure3) and the maximum temperature of Upper Baro Basin ranges between 27.7°C (in the wet season) and 37.2°C (in the dry season) and occurs in June to August, while the minimum falls between 10.8°C and 13.7°C in the year.Figure 3
Annual areal rainfall distribution (mm) for the period of 1987–2017.The topography of the Upper Baro Basin is characterized by variation from flat to mountainous terrain and falls in the altitude range of 390 m to over 3000 m a.m.s.l. with 42% in between 1000 m and 2000 m a.m.s.l. The eastern two-thirds of the basin area lies between 1000 m and 2400 m a.m.s.l. and a gently sloping plain lies in the west between 380 m and 500 m a.m.s.l. Dystric nitisols and dystric gleysols are the dominant soil type in the study area.


## 2.2. SWAT Model Data Types and Analysis

### 2.2.1. Meteorological Data

The SWAT model requires daily meteorological (Weather) data such as daily precipitation, minimum and maximum air temperature, relative humidity, wind speed, dew point, and daily sunshine. The standard deviation of each parameter was calculated. Thiessen method is used to adjust for nonuniform meteorological station distribution by proportion to the area which is closer. Eleven meteorological stations were selected (Table1) based on the quantity, quality, period, consistency, homogeneity, and their uniform distribution within and around the Upper Baro Basin. These data were obtained from the National Meteorological Agency of Ethiopia (NMAE) over the period (1987–2017). The missing values were estimated from other stations around the missed record station by using both arithmetic mean method and normal ratio method and also consistency test which shows double mass.Table 1
Summary of the meteorological stations.
StationCoordinateAltitude (m)Station classPeriod of dataLongitudeLatitudeGambella34.58338.25430C21983–2017Mettu35.566678.28331711C21981–2017Dembi Dolo34.88.51671850C21982–2017Ayira35.559.11555C21987–2017Gatira36.27.9832358C21987–2017Masha35.46677.752282C11981–2017Gore35.53338.13332033C11981–2017Seko Humbi34.98338.71661860C21987–2017Gimbi35.78339.16661970C21983–2017Bure35.18.23331750C21981–2017Alge35.66678.53331880C21987–2017The SWAT weather generator model (WGEN) was used to match the access table and lookup table. The daily precipitation and temperature of all gauging stations were prepared in comma-delimited (.csv) format. Solar radiation, relative humidity, and wind speed data were available only for principal Gore and Masha stations (synoptic stations).


### 2.2.2. Soil Data

Seventy soil types have been identified at Upper Baro Basin (Table2). The spatial soil data were obtained from the Ethiopian Ministry of Water, Irrigation, and Electricity (MWIE), which is prepared according to FAO soil classification. The most dominant soil in the basin is Dystric nitisols and Dystric gleysols.Table 2
Soil type’s area coverage and SWAT code.
ValueSoil typesArea (km2)Area covered (%)SWAT code1Dystric nitisols14151.7660.32DYNITISOLS2Dystric gleysols5423.5223.12DYGLEYSOLS3Orthic acrisols326.871.39ACRISOLS4Chromic luvisols120.450.51CHLUVISOLS5Orthic solonchaks778.063.32SOLONCHAKS6Calcic xerosols363.201.55CLXEROSOLS7Leptosols255.951.09LEPTOSOLS8Iatric cambisols1041.604.44EUCAMBISOLS9Eutric nitisols108.430.46EUNITISOLS10Eutric fluvisols276.991.18EUFLUVISOLS11Dystric cambisols58.400.25DYCAMBISOLS12Gypsic yermosols43.390.18GYYERMOSOLS13Dystric fluvisols231.440.99DYFLUVISOLS14Cambisols2.540.01CAMBISOLS15Calcic fluvisols92.400.39CLFLUVISOLS16Calcic cambisols167.830.72CLCAMBISOLS17Calcaric flubisols19.170.08CLFLUBISOLSSWAT model requires physical and chemical soil properties in the simulation of a hydrological component such as available water content, soil texture, hydraulic conductivity, bulk density, and organic carbon content for the different layers of each soil type where shown. The major soils of the Upper Baro Basin are as shown below table and a lookup table was prepared for each type of soil.


### 2.2.3. DEM

DEM with a resolution of 30 m was downloaded from United States Geological Surveyhttps://earthexplorer.usgs.gov/. The DEM is one of the essential inputs required by SWAT to delineate Upper Baro Basin (Figure 4), define the stream network, and determine sub-basin parameters such as slope area and slope length.Figure 4
DEM map (meters a.m.s.l.) of Upper Baro Basin.


### 2.2.4. Slope Classes in Subwatershed

To develop the hydrological response unit (HRU) in the SWAT model, the slope is essential. It generates from the resolution of 30 m ∗ 30 m DEM for the study area. Slope classification may be single class or multiclass. For this study, the slope option (an option for considering different slope classes for HRU definition) was selected. Therefore, the slope class in this study was classified into four classes. According to [22], slope classification was used to account for lower range in hydrological modelling. Depending on the slope of 0–5%, 5–10%, 10–15%, 15–20%, and >20% were selected for HRU creation of the study Upper Baro Basin. Finally, HRU definition analysis in SWAT helps to load LULC and soil type projects.


### 2.2.5. Land Use/Land Cover (LU/LC) Analysis

Landsat Satellite images of the years 1987 (Landsat 5), 2002 (Landsat 7) and 2017 (Landsat 8) were obtained from United States Geological Survey (USGS) (https://www.glovis.USGS.gov) (Table 3). The multidate multisensor satellite imagery was collected for successive cropping for each study year. Dataset selections were in dry season fixed with the lowest percent or zero monthly cloud cover. Therefore, all images are cloud-free (0%) in the study area.Table 3
Landsat data information for the study area.
Reference dataPath/rowAcquisition dataImageSensorResolution1987170/05422/01/1987Landsat 5TM30 m ∗ 30 m170/055171/0542002170/05423/01/2002Landsat 7ETM+30 m ∗ 30 m170/055171/0542017170/05427/01/2017Landsat 8OLF15 m ∗ 15 m170/055171/054The images were already terrain corrected projected to the Universal Transverse Mercator (UTM) of WGS84 zone 37 N. They were radiometrically calibrated to reflectance value and the Quick Atmospheric Correction algorithm was applied in ERDAS software. Radiometric and geometric corrections were applied to the removal of sensor noise, haze, correction data loss, and missing lines due to a solar position, satellite calibration [24, 25]. A supervised image classification method was employed using Maximum Likelihood Algorithm classification. For supervised classification, the Region of Interest (ROI) was made a signature for different LU/LC categories. The reference data collected using a GPS receiver were first converted to a vector file, then the region of interest (ROI) was made in ERDAS software. Using the ROI, the spectral signature of each LU/LC type has been extracted. Further investigations were carried out by taking these sample points in the land use types using a handheld GPS receiver. Based on this analysis and the existing land use distribution information, the area was categorized into eight broad LU/LC classes: cultivated land, woodland, forest land, grazing land, bare land, shrubland, water bodies, and built-up (Table 3).Overall accuracy and Kappa coefficient, the classify maps will check classification accuracy. According to [26, 27], a good overall classification map accuracy is accepted if and only it possesses an accuracy of at least 85%. The LU/LC was an input of the SWAT model to describe the HRU of the subwatershed. The SWAT Model has predetermined four-letter codes for each LU/LC category (Table 4). These codes were used to link the LU/LC map of the study area with the SWAT land database. To make it compatible, LU/LC with the SWAT model lookup table was prepared for each type of LU/LC.Table 4
LU/LC class SWAT code.
ValueLU/LC classSWAT databaseSWAT code1Forest landForest evergreenFRSE2WoodlandForest-mixedFRST3GrasslandRange-grassesRNGE4ShrublandRange-brushRNGB5SettlementResidential-medium densityURMD6Agricultural landAgricultural land-genericAGRL7Bare landBarrenBARR8WaterbodyWaterWATR


## 2.2.1. Meteorological Data

The SWAT model requires daily meteorological (Weather) data such as daily precipitation, minimum and maximum air temperature, relative humidity, wind speed, dew point, and daily sunshine. The standard deviation of each parameter was calculated. Thiessen method is used to adjust for nonuniform meteorological station distribution by proportion to the area which is closer. Eleven meteorological stations were selected (Table1) based on the quantity, quality, period, consistency, homogeneity, and their uniform distribution within and around the Upper Baro Basin. These data were obtained from the National Meteorological Agency of Ethiopia (NMAE) over the period (1987–2017). The missing values were estimated from other stations around the missed record station by using both arithmetic mean method and normal ratio method and also consistency test which shows double mass.Table 1
Summary of the meteorological stations.
StationCoordinateAltitude (m)Station classPeriod of dataLongitudeLatitudeGambella34.58338.25430C21983–2017Mettu35.566678.28331711C21981–2017Dembi Dolo34.88.51671850C21982–2017Ayira35.559.11555C21987–2017Gatira36.27.9832358C21987–2017Masha35.46677.752282C11981–2017Gore35.53338.13332033C11981–2017Seko Humbi34.98338.71661860C21987–2017Gimbi35.78339.16661970C21983–2017Bure35.18.23331750C21981–2017Alge35.66678.53331880C21987–2017The SWAT weather generator model (WGEN) was used to match the access table and lookup table. The daily precipitation and temperature of all gauging stations were prepared in comma-delimited (.csv) format. Solar radiation, relative humidity, and wind speed data were available only for principal Gore and Masha stations (synoptic stations).


## 2.2.2. Soil Data

Seventy soil types have been identified at Upper Baro Basin (Table2). The spatial soil data were obtained from the Ethiopian Ministry of Water, Irrigation, and Electricity (MWIE), which is prepared according to FAO soil classification. The most dominant soil in the basin is Dystric nitisols and Dystric gleysols.Table 2
Soil type’s area coverage and SWAT code.
ValueSoil typesArea (km2)Area covered (%)SWAT code1Dystric nitisols14151.7660.32DYNITISOLS2Dystric gleysols5423.5223.12DYGLEYSOLS3Orthic acrisols326.871.39ACRISOLS4Chromic luvisols120.450.51CHLUVISOLS5Orthic solonchaks778.063.32SOLONCHAKS6Calcic xerosols363.201.55CLXEROSOLS7Leptosols255.951.09LEPTOSOLS8Iatric cambisols1041.604.44EUCAMBISOLS9Eutric nitisols108.430.46EUNITISOLS10Eutric fluvisols276.991.18EUFLUVISOLS11Dystric cambisols58.400.25DYCAMBISOLS12Gypsic yermosols43.390.18GYYERMOSOLS13Dystric fluvisols231.440.99DYFLUVISOLS14Cambisols2.540.01CAMBISOLS15Calcic fluvisols92.400.39CLFLUVISOLS16Calcic cambisols167.830.72CLCAMBISOLS17Calcaric flubisols19.170.08CLFLUBISOLSSWAT model requires physical and chemical soil properties in the simulation of a hydrological component such as available water content, soil texture, hydraulic conductivity, bulk density, and organic carbon content for the different layers of each soil type where shown. The major soils of the Upper Baro Basin are as shown below table and a lookup table was prepared for each type of soil.


## 2.2.3. DEM

DEM with a resolution of 30 m was downloaded from United States Geological Surveyhttps://earthexplorer.usgs.gov/. The DEM is one of the essential inputs required by SWAT to delineate Upper Baro Basin (Figure 4), define the stream network, and determine sub-basin parameters such as slope area and slope length.Figure 4
DEM map (meters a.m.s.l.) of Upper Baro Basin.


## 2.2.4. Slope Classes in Subwatershed

To develop the hydrological response unit (HRU) in the SWAT model, the slope is essential. It generates from the resolution of 30 m ∗ 30 m DEM for the study area. Slope classification may be single class or multiclass. For this study, the slope option (an option for considering different slope classes for HRU definition) was selected. Therefore, the slope class in this study was classified into four classes. According to [22], slope classification was used to account for lower range in hydrological modelling. Depending on the slope of 0–5%, 5–10%, 10–15%, 15–20%, and >20% were selected for HRU creation of the study Upper Baro Basin. Finally, HRU definition analysis in SWAT helps to load LULC and soil type projects.


## 2.2.5. Land Use/Land Cover (LU/LC) Analysis

Landsat Satellite images of the years 1987 (Landsat 5), 2002 (Landsat 7) and 2017 (Landsat 8) were obtained from United States Geological Survey (USGS) (https://www.glovis.USGS.gov) (Table 3). The multidate multisensor satellite imagery was collected for successive cropping for each study year. Dataset selections were in dry season fixed with the lowest percent or zero monthly cloud cover. Therefore, all images are cloud-free (0%) in the study area.Table 3
Landsat data information for the study area.
Reference dataPath/rowAcquisition dataImageSensorResolution1987170/05422/01/1987Landsat 5TM30 m ∗ 30 m170/055171/0542002170/05423/01/2002Landsat 7ETM+30 m ∗ 30 m170/055171/0542017170/05427/01/2017Landsat 8OLF15 m ∗ 15 m170/055171/054The images were already terrain corrected projected to the Universal Transverse Mercator (UTM) of WGS84 zone 37 N. They were radiometrically calibrated to reflectance value and the Quick Atmospheric Correction algorithm was applied in ERDAS software. Radiometric and geometric corrections were applied to the removal of sensor noise, haze, correction data loss, and missing lines due to a solar position, satellite calibration [24, 25]. A supervised image classification method was employed using Maximum Likelihood Algorithm classification. For supervised classification, the Region of Interest (ROI) was made a signature for different LU/LC categories. The reference data collected using a GPS receiver were first converted to a vector file, then the region of interest (ROI) was made in ERDAS software. Using the ROI, the spectral signature of each LU/LC type has been extracted. Further investigations were carried out by taking these sample points in the land use types using a handheld GPS receiver. Based on this analysis and the existing land use distribution information, the area was categorized into eight broad LU/LC classes: cultivated land, woodland, forest land, grazing land, bare land, shrubland, water bodies, and built-up (Table 3).Overall accuracy and Kappa coefficient, the classify maps will check classification accuracy. According to [26, 27], a good overall classification map accuracy is accepted if and only it possesses an accuracy of at least 85%. The LU/LC was an input of the SWAT model to describe the HRU of the subwatershed. The SWAT Model has predetermined four-letter codes for each LU/LC category (Table 4). These codes were used to link the LU/LC map of the study area with the SWAT land database. To make it compatible, LU/LC with the SWAT model lookup table was prepared for each type of LU/LC.Table 4
LU/LC class SWAT code.
ValueLU/LC classSWAT databaseSWAT code1Forest landForest evergreenFRSE2WoodlandForest-mixedFRST3GrasslandRange-grassesRNGE4ShrublandRange-brushRNGB5SettlementResidential-medium densityURMD6Agricultural landAgricultural land-genericAGRL7Bare landBarrenBARR8WaterbodyWaterWATR


## 2.3. SWAT Model Description

SWAT is a physically based semidistributed continuous time-scale hydrological model, which works on a daily time step. This model can simulate runoff, sediment, nutrients, pesticide, and bacterial transport from agricultural watersheds [28]. It simulates the hydrological cycle parameters based on the water balance represented in equation (1) within the watershed [29]. More details on the information about the SWAT model can be found in the SWAT user manual [29].The hydrologic cycle simulated by the SWAT is based on the water balance equation [29] as follows:(1)SWt=SWo+∑i=1tRi−Qi−ETi−Pi−QRi,where SWt and SWo are the final and initial soil water content, respectively, (mm), Qi = daily surface runoff (mm), Ri = daily rainfall (mm), ETi = daily evapotranspiration (mm), Ri = daily lateral flow (mm), and Pi = daily percolation (mm).


## 2.4. Peak Runoff Rate

The peak runoff rate is a result of the erosive power of a storm and is used to predict sediment loss. SWAT calculates the peak runoff rate with a modified rational method for each HRU using the following equation:(2)Qpeak=αtc∗Qsurf∗Atconc,where Q peak is peak runoff rate (m3/s), αtc the fraction of daily rainfall that occurs during the time of concentration, Qsurf is the surface runoff (mm), A is the sub-basin area (km2), and tconc is the time of concentration (hr).


## 2.5. Surface Runoff Generation

The CN method was initially developed for small agricultural watersheds and the CN varies nonlinearly with the moisture content of the soil. It drops to zero as the soil approaches the wilting point and increases to be near 100 as the soil approaches saturation, with higher CNs associated with higher runoff potential watershed. This method is widely used [28].According to USDA, 1985 [30], this can be mathematically expressed using the following equation:(3)Qs=Ra−Ia2Ra−Ia+S,ForRa>Ia,Q=0forRa≤Ia,where Qs is the accumulated runoff (mm), Ra is the rainfall depth for the day (mm), and(4)S=25.4∗100CN−10,where CN is the curve number for the day and its value is the function of land use practice, soil permeability and soil hydrologic group. Conceptual framework showing the components and relationships have been used as a framework for the analysis of hydrology in research (Figure 5).Figure 5
General conceptual methodology used for modelling hydrology in the Upper Baro Basin.


## 2.6. SWAT Model Simulation, Sensitivity Analysis, Calibration, and Validation

The simulation result cannot be directly used for further analysis. Instead, the ability of the model to sufficiently predict the constituent streamflow should be evaluated through sensitivity analysis, model calibration, and model validation [31].


### 2.6.1. Sensitivity Analysis

Fifteen hydrological parameters related to streamflow were selected for sensitivity analysis in the study area. The sensitivity of each parameter was identified using theT test and P values. P values were used to determine the significance of the sensitivity. A value close to zero has more significance. T test provides a measure of sensitivity and the most sensitive parameters usually appear with P values less than the alpha level of 0.05 (larger absolute values are more sensitive) [32].


### 2.6.2. Model Calibration and Validation

The SUFI-2 algorithm was applied and frequently used in this study of calibration for the SWAT model at a large scale due to its easy application and the reduced number of model runs needed to achieve a decent prediction [33, 34]. In this study, automatic calibration was made every month from January 1, 1990, to December 31, 2002, until the average simulated value came closer to the measured value. Automatic calibration makes use of a numerical algorithm to increase the performance of the model and to optimize the numerical objective functions. After the simulation result for the calibration, the period had fulfilled the above statistical criteria, validation was performed for an independent period of records from January 1, 2003, to December 31, 2009, and the results were compared against an independent set of measured Upper Baro Basin.


### 2.6.3. Model Performance Evaluation

Model performance was carried out to verify the robustness of the model to simulate hydrological processes. A model performance framework proposed by [35] was used in this study.Nash and Sutcliffe simulation efficiency (NSE) indicates the degree of fitness of observed and simulated data, and it is calculated using the following equation:(5)NSE=1−∑i=0NQobs−QSim2∑i=0NQobs−Q¯obs2,where N = number of compared values, Qobs is the observed data, Q¯obs is the observed mean, and Qsim is the simulated data.The NSE indicates how well the plot of observed versus simulated value fits the 1 : 1 line. The closer the model efficiency is to 1, the more accurate the model, and if it is found between 0 and 1, it indicates deviations between measured and predicted values. If NSE is negative, predictions are poor, and the average value of output is a better estimate than the model prediction [36].Coefficient of determination (R2) measures the ability of a model to predict or explain an outcome in the linear regression setting (Equation (6)). R2 ranges from 0 to 1, with higher values indicating less error variance, and typically values greater than 0.5 are considered acceptable [35, 37, 38].(6)R2=∑Xi−XavYi−Yav∑Xi−Xav2∑Yi−Yav2,where Xi is measured value, Xav is average measured value, Yi is simulated value, and Yav is the average simulated value.Percent bias (PBIAS) measures the average difference between the simulated and measured values for a given quantity over a specified period (usually the entire calibration or validation period) and it is calculated using the following equation [39, 40]:(7)PBIAS=∑Yi−∑XiXi,where Xi is measured value and Yi is simulated value.Model evaluation classifies as follows:(i)
0.75 < NSE < 1.00, PBIAS < ±10—better predictor(ii)
0.65 < NSE ≤ 0.75, ±10 < PBIAS < ±15—good(iii)
0.50 < NSE ≤ 0.65, ±15 < PBIAS < ±25—satisfactory(iv)
(NSE ≤ 0.50, PBIAS >±25)—unsatisfactory


## 2.6.1. Sensitivity Analysis

Fifteen hydrological parameters related to streamflow were selected for sensitivity analysis in the study area. The sensitivity of each parameter was identified using theT test and P values. P values were used to determine the significance of the sensitivity. A value close to zero has more significance. T test provides a measure of sensitivity and the most sensitive parameters usually appear with P values less than the alpha level of 0.05 (larger absolute values are more sensitive) [32].


## 2.6.2. Model Calibration and Validation

The SUFI-2 algorithm was applied and frequently used in this study of calibration for the SWAT model at a large scale due to its easy application and the reduced number of model runs needed to achieve a decent prediction [33, 34]. In this study, automatic calibration was made every month from January 1, 1990, to December 31, 2002, until the average simulated value came closer to the measured value. Automatic calibration makes use of a numerical algorithm to increase the performance of the model and to optimize the numerical objective functions. After the simulation result for the calibration, the period had fulfilled the above statistical criteria, validation was performed for an independent period of records from January 1, 2003, to December 31, 2009, and the results were compared against an independent set of measured Upper Baro Basin.


## 2.6.3. Model Performance Evaluation

Model performance was carried out to verify the robustness of the model to simulate hydrological processes. A model performance framework proposed by [35] was used in this study.Nash and Sutcliffe simulation efficiency (NSE) indicates the degree of fitness of observed and simulated data, and it is calculated using the following equation:(5)NSE=1−∑i=0NQobs−QSim2∑i=0NQobs−Q¯obs2,where N = number of compared values, Qobs is the observed data, Q¯obs is the observed mean, and Qsim is the simulated data.The NSE indicates how well the plot of observed versus simulated value fits the 1 : 1 line. The closer the model efficiency is to 1, the more accurate the model, and if it is found between 0 and 1, it indicates deviations between measured and predicted values. If NSE is negative, predictions are poor, and the average value of output is a better estimate than the model prediction [36].Coefficient of determination (R2) measures the ability of a model to predict or explain an outcome in the linear regression setting (Equation (6)). R2 ranges from 0 to 1, with higher values indicating less error variance, and typically values greater than 0.5 are considered acceptable [35, 37, 38].(6)R2=∑Xi−XavYi−Yav∑Xi−Xav2∑Yi−Yav2,where Xi is measured value, Xav is average measured value, Yi is simulated value, and Yav is the average simulated value.Percent bias (PBIAS) measures the average difference between the simulated and measured values for a given quantity over a specified period (usually the entire calibration or validation period) and it is calculated using the following equation [39, 40]:(7)PBIAS=∑Yi−∑XiXi,where Xi is measured value and Yi is simulated value.Model evaluation classifies as follows:(i)
0.75 < NSE < 1.00, PBIAS < ±10—better predictor(ii)
0.65 < NSE ≤ 0.75, ±10 < PBIAS < ±15—good(iii)
0.50 < NSE ≤ 0.65, ±15 < PBIAS < ±25—satisfactory(iv)
(NSE ≤ 0.50, PBIAS >±25)—unsatisfactory


## 3. Results and Discussion

### 3.1. LU/LC Change

The spatial distribution and percentage of the total area of the eight LU/LC types between 1987–2002, 2002–2017, and 1987–2017 are shown in Figure6 (Table 5).Figure 6
Land use land cover map of 1987, 2002, and 2017.Table 5
The magnitude of land use/land cover type for 1987, 2002, and 2017 in the Upper Baro Basin.
Land use classSpatial convergeChange between years1987200220172002–19872017–20022017–1987(km2)%(km2)%(km2)%(km2)%(km2)%(km2)%Grass land4639.2119.773272.1513.95969.344.13−1367.05–5.83–2302.81–9.82–3669.87–15.64Agricultural2065.928.813659.0815.606290.5626.811593.156.792631.4811.224224.6318.01Bare land92.810.40166.240.71217.430.9373.430.3151.190.22124.620.53Urban471.302.012060.528.783638.1615.511589.236.771577.646.723166.8613.50Woodland2656.8111.322812.8011.992,447.7410.43155.980.66–365.06–1.56–209.08–0.89Forestland8619.9736.747819.9133.337323.3531.21–800.06–3.41-496.56–2.12–1296.62–5.53Schruble land4635.9419.763544.6315.112393.9510.20–1091.31–4.65–1150.69–4.90–2241.99–9.56Waterbody280.041.19126.660.5481.480.35–153.37–0.65–45.18–0.19–198.56–0.85The results obtained for 1987–2002 reveal thatforest land covers 36.74% and followed by grassland, 19.77% and shrubland and woodland also cover 19.76% and 11.32%, respectively. Between 2002 and 2017, the distribution of forest cover demonstrates highest in the catchment, which covered 33.33%, followed by agricultural land and shrubland, with the value of 15.6%, 15.1%, respectively. Forest land and shrubland decreased after 1987–2002, with −3.41% and −4.65%, respectively (Figure 7). The area of forest land decreased due to of change of land use (conversion) to the expansion of agricultural land in the southwestern part of the catchment. The urban area gradually expanded from 2.01% to 15.51% over the entire study period. These results indicate that government policy does not have a strong impact on land use.Figure 7
Comparison of LULC classes of the years 1987–2002, 2002–2017, and 1987–2017.In the last 30 years, the study displays an increase in agriculture, bare land, urban, and shrubland, whereas the forest, grassland, water, and wetland were noted declined. Agricultural land increased by 18.01% during the study period. This is due to the increase of population growth that causes the increase in demand for cultivation land for different agricultural products like cereal crops, field crops, industrial crops, whereas the forest area was diminished. This might be due to deforestation activities that have taken place for agriculture and urban expansion.Consequently, agricultural investment policies encouraged further deforestation that leads to land conversion [41] and various types of natural forests and woodlands in the Western part of Ethiopia, being cleared and replaced by commercial agriculture [42, 43]. In the study basin, cultivated land has gained/increased by 18.01%, mainly converted from the grazing land by 15.64%. Forest and shrubland conversion also contributed to the increasing cultivated land.In this study, the accuracy assessment was performed using the 340 referenced points collected through stratified sampling. Overall accuracies obtained were 92.14%, 94.63%, and 95.93% for the 1987, 2002, and 2017 images, respectively (Table6). The overall accuracy was computed by dividing the correct classification by the total number of reference pixels in the error matrix.Table 6
Confusion matrix for classified images.
LU/LC class198720022017producerUserProducerUserProducerUserForest land (FT)9795.5996.089810096.61Woodland (WD)9593.4497.4410096.15100Grass land (GR)81.0183.3393.5590.6392.1694Urban (UB)87.593.3395.7484.9110088Shrub land (SL)91.6793.6210093.88100100Agriculture (AG)92.1189.7486.5491.8491.0792.73Waterbody (WB)10010097.44100100100Bare land (BL)90.2486.0591.8410087.1100Overall accuracy92.14%94.6395.93Kappa coefficient0.90940.93850.9533


### 3.2. Stream Flow Modeling Using SWAT

#### 3.2.1. Sensitivity Analysis

Sensitivity analysis was carried out to identify which model parameter was the most important or sensitive. Flow sensitivity analysis was carried out for 13 years, the calibration period (from January 1, 1990, to December 31, 2002). The default and optimal parameters and their ranking in terms of the 10 highest relative sensitivity values used to develop, calibrate, and validate the model (Figure8).Figure 8
Global sensitivity analysis for streamflow.


#### 3.2.2. SWAT Model Calibration and Validation

The simulated monthly flows demonstrated a good agreement with the monthly observed discharge of Upper Baro Basin with a coefficient of determination (R2 = 0.86), Nash–Sutcliffe efficiency (NSE = 0.80), and percent bias (PBIAS = −17%) (Figure 9). However, the model overestimated the peak monthly flow for some years in the calibration period; it followed the trend of observed monthly Upper Baro basin discharge and presented a good response to extreme rainfall events, which resulted in high runoff volume.Figure 9
Calibration of average monthly observed and simulated flow (1990–2002).SWAT model also successfully validated streamflow for an independent period (2003–2010) withR2 = 0.84, NS = 0.76, and PBIAS = −4.5% during the validation period; the values fulfilled the statistical model performance criteria R2 > 0.6 and NSE > 0.5 recommended by SWAT developer [37]. The percent bias (PBAIS) values (−17% and 4.5%) during calibration and validation periods, respectively, also fulfilled the statistical model performance criteria PBAIS> ±15 [35] (Table 7).Table 7
Result of final calibrated flow parameters for Upper Baro Basin.
Parameter nameFitted valueRangeR_CN2−0.230±0.25V_ALPHA_BF0.0530-1V_GW_DELAY113.500–500V_GW_REVAP0.13430.02–0.2V_CH_N20.34360.01–1R_SOL_K0.54500-1R_SOL_AWC0.7550-1V_ESCO.hru0.25300-1V_REVAPMN.gw0.37000–10V_GWQMN.gw1.9820–2The percentage of observations bracketed by the 95 PPU was 75% for calibration and 72% for validation. The r-factor, which is the average thickness of the 95PPU band divided by the standard deviation of the measured data, was 1.08 for calibration and 0.84 for validation. According to [32], a p-factor greater than 0.70 and r-factor less than 1.5 display that the model is accurate in simulating streamflow. Thus, the results displayed that the model accurately simulated the biophysical processes influencing flow in the study area.Generally, the above information showed that the performance of the model was better during the calibration period more than the validation period (Table8) and also the calibration validation period (Figures 9–12) shows there was good relation with observed and simulated streamflow.Table 8
Summary of model performance evaluation for calibration and validation period on monthly time steps.
PeriodMean annual water yield (m3/s)Monthly model efficiency measuresObservedSimulatedR2NSEPBAIS (%)Calibration422.461535.960.860.81−18.7Validation374.42451.120.840.76−19.8Figure 10
Simulated versus observed monthly flow for calibration (1990–2002).Figure 11
Validation of average monthly observed and simulated flow (2003–2009).Figure 12
Simulated versus observed monthly flow for calibration (2003–2009).


### 3.3. Land Use/Land Cover Change Response to Hydrological Process

The simulated results of the Upper Baro Basin under different land use/land cover scenarios are shown in Table9.Table 9
Average annual hydrological summaries under different LULC.
ParameterLU/LC 1987 (mm)LU/LC 2002 (mm)LU/LC 2017 (mm)Change detection1987–20022002–20171987–2017SUR Q773.793.05816.5320.0523.4843.53GW Q270.8259.29243.22−11.08−16.07−27.58LAT Q87.3785.4981.74−1.88−3.75−5.63WYLD1131.711137.831141.496.123.669.78PERC302.62291.17275.00−11.45−16.17−27.62ET807.1800.4797−6.7−3.4−10.1SEDI62.72106.02124.4943.318.4761.77CN76.8977.8378.82


### 3.4. Mean Monthly Wet and Dry Month Streamflow Simulation and Their Variability

Months July, August, and September were considered as wet periods and Jan, Feb, and March were considered dry months; these seasonal variabilities of streamflow were evaluated.


### 3.5. Spatial Variability of Surface Runoff over Sub-Basin

Spatial variability of surface runoff shows that surface runoff is highly sensitive to the change of LU/LC over subwatersheds in the period of 1987, 2002, and 2017 (Figure13).Figure 13
Spatial distribution predicted annual surface runoff (SURQ) (mm) for different LULC conditions.


### 3.6. Comparison of Mean Monthly Surface Runoff and Groundwater

Decreases in base flow had a potential effect on the change in annual river flow (Figure14). Results from LU/LC change scenarios, the mean monthly (wet and dry) flow has increased from 82.9 mm and decreased 20.54 mm concerning the change of LU/LC from 1987 to 2017 (Figure 15).Figure 14
Mean monthly surface runoff (mm).Figure 15
Mean monthly groundwater recharge (mm).The results (Table10) demonstrate that the mean dry monthly flow for 1987–2017 land cover decreased by −20.54 (m3/s) when compared to the land cover in 1987–2017. Changes in land use occurred in all parts of the Basin but mainly in the southwestern part. Vegetation cover within the catchment had changed considerably due to human activities, mainly through the conversion of natural vegetation to agriculture, Urban, and largescale investment. These changes have altered water resources through biophysical and biogeochemical processes in the soil. Through these land-cover changes, the biophysical and morphological characteristics of the vegetation have been modified, and therefore, the changes have affected the exchange of water between the atmosphere and land surface. Similarly, the results indicate an increase in surface runoff increased from 773 mm/yr in 1987 to 793.05 mm/yr in 2002 and 816.53 mm/yr in 2017 due to a drastic increase in agriculture and urban activity at the expense of other land covers while groundwater was decreased from 283.21 mm to 257.13 mm in 1987 as compared to 2017 due to LULC change. Similarly, [44, 45] reported that an increase in annual surface runoff increased when forest lands were converted to the other classes especially cultivated land. The expansion of agricultural land by displacing forest and shrubland results in an increase in surface erosion generation following rainfall events and causes a reduction in soil moisture conditions and groundwater recharge to shallow aquifers. Similarly, studies suggested that grassland and lateral flow are positively correlated [46, 47].Table 10
Mean monthly wet and dry month streamflow simulation and their variability.
Mean monthly streamflow (m3/s)Mean monthly streamflow change1987200220172002–19872017–20022017–1987DryWetDryWetDryWetDryWetDryWetDryWet84.2871.272.1590763.66954.1−12.05+35.8−8.49+47.1−20.54+82.9


## 3.1. LU/LC Change

The spatial distribution and percentage of the total area of the eight LU/LC types between 1987–2002, 2002–2017, and 1987–2017 are shown in Figure6 (Table 5).Figure 6
Land use land cover map of 1987, 2002, and 2017.Table 5
The magnitude of land use/land cover type for 1987, 2002, and 2017 in the Upper Baro Basin.
Land use classSpatial convergeChange between years1987200220172002–19872017–20022017–1987(km2)%(km2)%(km2)%(km2)%(km2)%(km2)%Grass land4639.2119.773272.1513.95969.344.13−1367.05–5.83–2302.81–9.82–3669.87–15.64Agricultural2065.928.813659.0815.606290.5626.811593.156.792631.4811.224224.6318.01Bare land92.810.40166.240.71217.430.9373.430.3151.190.22124.620.53Urban471.302.012060.528.783638.1615.511589.236.771577.646.723166.8613.50Woodland2656.8111.322812.8011.992,447.7410.43155.980.66–365.06–1.56–209.08–0.89Forestland8619.9736.747819.9133.337323.3531.21–800.06–3.41-496.56–2.12–1296.62–5.53Schruble land4635.9419.763544.6315.112393.9510.20–1091.31–4.65–1150.69–4.90–2241.99–9.56Waterbody280.041.19126.660.5481.480.35–153.37–0.65–45.18–0.19–198.56–0.85The results obtained for 1987–2002 reveal thatforest land covers 36.74% and followed by grassland, 19.77% and shrubland and woodland also cover 19.76% and 11.32%, respectively. Between 2002 and 2017, the distribution of forest cover demonstrates highest in the catchment, which covered 33.33%, followed by agricultural land and shrubland, with the value of 15.6%, 15.1%, respectively. Forest land and shrubland decreased after 1987–2002, with −3.41% and −4.65%, respectively (Figure 7). The area of forest land decreased due to of change of land use (conversion) to the expansion of agricultural land in the southwestern part of the catchment. The urban area gradually expanded from 2.01% to 15.51% over the entire study period. These results indicate that government policy does not have a strong impact on land use.Figure 7
Comparison of LULC classes of the years 1987–2002, 2002–2017, and 1987–2017.In the last 30 years, the study displays an increase in agriculture, bare land, urban, and shrubland, whereas the forest, grassland, water, and wetland were noted declined. Agricultural land increased by 18.01% during the study period. This is due to the increase of population growth that causes the increase in demand for cultivation land for different agricultural products like cereal crops, field crops, industrial crops, whereas the forest area was diminished. This might be due to deforestation activities that have taken place for agriculture and urban expansion.Consequently, agricultural investment policies encouraged further deforestation that leads to land conversion [41] and various types of natural forests and woodlands in the Western part of Ethiopia, being cleared and replaced by commercial agriculture [42, 43]. In the study basin, cultivated land has gained/increased by 18.01%, mainly converted from the grazing land by 15.64%. Forest and shrubland conversion also contributed to the increasing cultivated land.In this study, the accuracy assessment was performed using the 340 referenced points collected through stratified sampling. Overall accuracies obtained were 92.14%, 94.63%, and 95.93% for the 1987, 2002, and 2017 images, respectively (Table6). The overall accuracy was computed by dividing the correct classification by the total number of reference pixels in the error matrix.Table 6
Confusion matrix for classified images.
LU/LC class198720022017producerUserProducerUserProducerUserForest land (FT)9795.5996.089810096.61Woodland (WD)9593.4497.4410096.15100Grass land (GR)81.0183.3393.5590.6392.1694Urban (UB)87.593.3395.7484.9110088Shrub land (SL)91.6793.6210093.88100100Agriculture (AG)92.1189.7486.5491.8491.0792.73Waterbody (WB)10010097.44100100100Bare land (BL)90.2486.0591.8410087.1100Overall accuracy92.14%94.6395.93Kappa coefficient0.90940.93850.9533


## 3.2. Stream Flow Modeling Using SWAT

### 3.2.1. Sensitivity Analysis

Sensitivity analysis was carried out to identify which model parameter was the most important or sensitive. Flow sensitivity analysis was carried out for 13 years, the calibration period (from January 1, 1990, to December 31, 2002). The default and optimal parameters and their ranking in terms of the 10 highest relative sensitivity values used to develop, calibrate, and validate the model (Figure8).Figure 8
Global sensitivity analysis for streamflow.


### 3.2.2. SWAT Model Calibration and Validation

The simulated monthly flows demonstrated a good agreement with the monthly observed discharge of Upper Baro Basin with a coefficient of determination (R2 = 0.86), Nash–Sutcliffe efficiency (NSE = 0.80), and percent bias (PBIAS = −17%) (Figure 9). However, the model overestimated the peak monthly flow for some years in the calibration period; it followed the trend of observed monthly Upper Baro basin discharge and presented a good response to extreme rainfall events, which resulted in high runoff volume.Figure 9
Calibration of average monthly observed and simulated flow (1990–2002).SWAT model also successfully validated streamflow for an independent period (2003–2010) withR2 = 0.84, NS = 0.76, and PBIAS = −4.5% during the validation period; the values fulfilled the statistical model performance criteria R2 > 0.6 and NSE > 0.5 recommended by SWAT developer [37]. The percent bias (PBAIS) values (−17% and 4.5%) during calibration and validation periods, respectively, also fulfilled the statistical model performance criteria PBAIS> ±15 [35] (Table 7).Table 7
Result of final calibrated flow parameters for Upper Baro Basin.
Parameter nameFitted valueRangeR_CN2−0.230±0.25V_ALPHA_BF0.0530-1V_GW_DELAY113.500–500V_GW_REVAP0.13430.02–0.2V_CH_N20.34360.01–1R_SOL_K0.54500-1R_SOL_AWC0.7550-1V_ESCO.hru0.25300-1V_REVAPMN.gw0.37000–10V_GWQMN.gw1.9820–2The percentage of observations bracketed by the 95 PPU was 75% for calibration and 72% for validation. The r-factor, which is the average thickness of the 95PPU band divided by the standard deviation of the measured data, was 1.08 for calibration and 0.84 for validation. According to [32], a p-factor greater than 0.70 and r-factor less than 1.5 display that the model is accurate in simulating streamflow. Thus, the results displayed that the model accurately simulated the biophysical processes influencing flow in the study area.Generally, the above information showed that the performance of the model was better during the calibration period more than the validation period (Table8) and also the calibration validation period (Figures 9–12) shows there was good relation with observed and simulated streamflow.Table 8
Summary of model performance evaluation for calibration and validation period on monthly time steps.
PeriodMean annual water yield (m3/s)Monthly model efficiency measuresObservedSimulatedR2NSEPBAIS (%)Calibration422.461535.960.860.81−18.7Validation374.42451.120.840.76−19.8Figure 10
Simulated versus observed monthly flow for calibration (1990–2002).Figure 11
Validation of average monthly observed and simulated flow (2003–2009).Figure 12
Simulated versus observed monthly flow for calibration (2003–2009).


## 3.2.1. Sensitivity Analysis

Sensitivity analysis was carried out to identify which model parameter was the most important or sensitive. Flow sensitivity analysis was carried out for 13 years, the calibration period (from January 1, 1990, to December 31, 2002). The default and optimal parameters and their ranking in terms of the 10 highest relative sensitivity values used to develop, calibrate, and validate the model (Figure8).Figure 8
Global sensitivity analysis for streamflow.


## 3.2.2. SWAT Model Calibration and Validation

The simulated monthly flows demonstrated a good agreement with the monthly observed discharge of Upper Baro Basin with a coefficient of determination (R2 = 0.86), Nash–Sutcliffe efficiency (NSE = 0.80), and percent bias (PBIAS = −17%) (Figure 9). However, the model overestimated the peak monthly flow for some years in the calibration period; it followed the trend of observed monthly Upper Baro basin discharge and presented a good response to extreme rainfall events, which resulted in high runoff volume.Figure 9
Calibration of average monthly observed and simulated flow (1990–2002).SWAT model also successfully validated streamflow for an independent period (2003–2010) withR2 = 0.84, NS = 0.76, and PBIAS = −4.5% during the validation period; the values fulfilled the statistical model performance criteria R2 > 0.6 and NSE > 0.5 recommended by SWAT developer [37]. The percent bias (PBAIS) values (−17% and 4.5%) during calibration and validation periods, respectively, also fulfilled the statistical model performance criteria PBAIS> ±15 [35] (Table 7).Table 7
Result of final calibrated flow parameters for Upper Baro Basin.
Parameter nameFitted valueRangeR_CN2−0.230±0.25V_ALPHA_BF0.0530-1V_GW_DELAY113.500–500V_GW_REVAP0.13430.02–0.2V_CH_N20.34360.01–1R_SOL_K0.54500-1R_SOL_AWC0.7550-1V_ESCO.hru0.25300-1V_REVAPMN.gw0.37000–10V_GWQMN.gw1.9820–2The percentage of observations bracketed by the 95 PPU was 75% for calibration and 72% for validation. The r-factor, which is the average thickness of the 95PPU band divided by the standard deviation of the measured data, was 1.08 for calibration and 0.84 for validation. According to [32], a p-factor greater than 0.70 and r-factor less than 1.5 display that the model is accurate in simulating streamflow. Thus, the results displayed that the model accurately simulated the biophysical processes influencing flow in the study area.Generally, the above information showed that the performance of the model was better during the calibration period more than the validation period (Table8) and also the calibration validation period (Figures 9–12) shows there was good relation with observed and simulated streamflow.Table 8
Summary of model performance evaluation for calibration and validation period on monthly time steps.
PeriodMean annual water yield (m3/s)Monthly model efficiency measuresObservedSimulatedR2NSEPBAIS (%)Calibration422.461535.960.860.81−18.7Validation374.42451.120.840.76−19.8Figure 10
Simulated versus observed monthly flow for calibration (1990–2002).Figure 11
Validation of average monthly observed and simulated flow (2003–2009).Figure 12
Simulated versus observed monthly flow for calibration (2003–2009).


## 3.3. Land Use/Land Cover Change Response to Hydrological Process

The simulated results of the Upper Baro Basin under different land use/land cover scenarios are shown in Table9.Table 9
Average annual hydrological summaries under different LULC.
ParameterLU/LC 1987 (mm)LU/LC 2002 (mm)LU/LC 2017 (mm)Change detection1987–20022002–20171987–2017SUR Q773.793.05816.5320.0523.4843.53GW Q270.8259.29243.22−11.08−16.07−27.58LAT Q87.3785.4981.74−1.88−3.75−5.63WYLD1131.711137.831141.496.123.669.78PERC302.62291.17275.00−11.45−16.17−27.62ET807.1800.4797−6.7−3.4−10.1SEDI62.72106.02124.4943.318.4761.77CN76.8977.8378.82


## 3.4. Mean Monthly Wet and Dry Month Streamflow Simulation and Their Variability

Months July, August, and September were considered as wet periods and Jan, Feb, and March were considered dry months; these seasonal variabilities of streamflow were evaluated.


## 3.5. Spatial Variability of Surface Runoff over Sub-Basin

Spatial variability of surface runoff shows that surface runoff is highly sensitive to the change of LU/LC over subwatersheds in the period of 1987, 2002, and 2017 (Figure13).Figure 13
Spatial distribution predicted annual surface runoff (SURQ) (mm) for different LULC conditions.


## 3.6. Comparison of Mean Monthly Surface Runoff and Groundwater

Decreases in base flow had a potential effect on the change in annual river flow (Figure14). Results from LU/LC change scenarios, the mean monthly (wet and dry) flow has increased from 82.9 mm and decreased 20.54 mm concerning the change of LU/LC from 1987 to 2017 (Figure 15).Figure 14
Mean monthly surface runoff (mm).Figure 15
Mean monthly groundwater recharge (mm).The results (Table10) demonstrate that the mean dry monthly flow for 1987–2017 land cover decreased by −20.54 (m3/s) when compared to the land cover in 1987–2017. Changes in land use occurred in all parts of the Basin but mainly in the southwestern part. Vegetation cover within the catchment had changed considerably due to human activities, mainly through the conversion of natural vegetation to agriculture, Urban, and largescale investment. These changes have altered water resources through biophysical and biogeochemical processes in the soil. Through these land-cover changes, the biophysical and morphological characteristics of the vegetation have been modified, and therefore, the changes have affected the exchange of water between the atmosphere and land surface. Similarly, the results indicate an increase in surface runoff increased from 773 mm/yr in 1987 to 793.05 mm/yr in 2002 and 816.53 mm/yr in 2017 due to a drastic increase in agriculture and urban activity at the expense of other land covers while groundwater was decreased from 283.21 mm to 257.13 mm in 1987 as compared to 2017 due to LULC change. Similarly, [44, 45] reported that an increase in annual surface runoff increased when forest lands were converted to the other classes especially cultivated land. The expansion of agricultural land by displacing forest and shrubland results in an increase in surface erosion generation following rainfall events and causes a reduction in soil moisture conditions and groundwater recharge to shallow aquifers. Similarly, studies suggested that grassland and lateral flow are positively correlated [46, 47].Table 10
Mean monthly wet and dry month streamflow simulation and their variability.
Mean monthly streamflow (m3/s)Mean monthly streamflow change1987200220172002–19872017–20022017–1987DryWetDryWetDryWetDryWetDryWetDryWet84.2871.272.1590763.66954.1−12.05+35.8−8.49+47.1−20.54+82.9


## 4. Conclusions

In the present study, the Upper Baro Basin has been simulated using the SWAT hydrological model for assessing the impact of LULC change on its hydrological regime successfully. Initially, the change in LULC between the years 1987, 2002, and 2017 has been assessed. It is evident that land cover changes had occurred between 1987 and 2017. Forest cover had decreased from 36.74% to 31.21%, attributable to logging for timber, firewood, and clearing for agricultural purposes. In contrast, the agricultural area increased over the years from 8.81% to 26.81%. Similarly, grasslands and bushes have also been converted to agricultural areas. Sensitivity analysis of the SWAT parameters indicates that runoff is most sensitive to curve number, an available water capacity of the soil layer, soil evaporation composition factor, manning’s “n” value for the main channel, and saturated hydraulic conductivity. Seasonal streamflow variability due to the LU/LC was assessed and comparisons were made on surface runoff, lateral flow, and groundwater flow contributions to streamflow based on the three simulation outputs. Surface runoff was increased from 713.85 mm to 760.49 mm, while groundwater was decreased from 283.21 mm to 257.13 mm in the years 1987 and 2017. Also, surface runoff was increased from 737.67 mm to 760.49 mm, while groundwater was decreased from 266.32 mm to 257.13 mm in the year 2002 and 2017. Land use/land cover change over the last 30 years also affect the soil and soil properties in the study area significantly. For example, forest land declined and agricultural land is increased. This dynamism would maximize soil erosion, which leads to the loss of important soil nutrients. Therefore, different soil and water conservation measures are required to reverse the challenges.



---

*Source: 6617541-2021-07-29.xml*