Evaluation of the impacts of land use/cover changes on water balance of Bilate watershed, Rift valley basin, Ethiopia

Land use/cover change is one of the responsible factors for changing the water balance of the watershed by altering the magnitude of surface runoff, inter ﬂ ow, base ﬂ ow, and evapotranspiration. This study was aimed at evaluating the impacts of land use/cover change on the water balance of Bilate watershed between 1989, 2002, and 2015. The water balance simulation model (WaSiM) was used to access the impacts of land use/cover change on water balance. The model was calibrated (1989 – 2003) and validated (2007 – 2015) using the stream ﬂ ow of at Bilate Tena gauging station. The result of land-use dynamics showed land use/cover change has a signi ﬁ cant impact on the water balance of the watershed like on runoff production, base ﬂ ow, inter ﬂ ow, evapotranspiration, and total simulation ﬂ ow. In the study watershed, the change in total simulated ﬂ ow increased by 77.83%, and surface runoff, inter ﬂ ow, and base ﬂ ow increased by 80.23%, 75.69%, and 87.79% respectively and evapotranspiration decreased by 6% throughout the study period (1989 – 2015). The results obtained from this study revealed that the watershed is under the land/cover change that shows its impacts on hydrological processes and water balance. Thus, effective information regarding the environmental response of land use/cover, change is important to hydrologists, land-use planners, watershed management, and decision-makers for sustainable water resource projects and ecosystem ser-vices. Therefore, the policy-makers, planners, and stakeholders should design strategies to ensure the sustainability of the watershed hydrology for the sake of protecting agricultural activities, and urban planning and management systems within the watershed area.

water retention (Tufa et al. 2014). The land use/cover plays a fundamental role in driving hydrological processes within a sub-basin (Gwate et al. 2015). These include changes in water demands such as irrigation, changes in groundwater recharge, and runoff, and changes in water quality from agricultural runoff (Guo et al. 2008).
Therefore, a far better understanding of land use/cover change, its effect, and interaction to the hydrology of a basin are highly essential in water supply from altered Hydrological processes of infiltration, groundwater recharge, and runoff, and changes in water quality from agricultural runoff (Guo et al. 2008).
Ethiopia is one of the most populous countries in Africa with over a population of 94 million people (CSA 2013). 85% of the population lives in rural areas and directly depends on the land for its livelihood (Asmamaw 2013). This means the demand for lands is increasing as the population increases. Agriculture, which depends on the availability of seasonal rainfall, is the main economy of the country. People need land for food production and housing and it is common practice to clear the forest for farming and housing activities. Therefore, the result of these activities is the land use/cover changes due to daily human intervention. Hence, understanding how the land cover changes influence on hydrology of the watershed enables planners to formulate policies to minimize the undesirable effects of future land cover changes. Small-scale sub-basinbased hydrological information considering land use/cover change is crucial for hydrological processes assessment for irrigated agriculture or any use of water. Water availability is becoming a critical factor in so many sectors so that assessing the anticipated impacts of land use/cover change on hydrology is unquestionable (Tubiello 2007).
Bilate watershed which is one of the sub-watersheds of the Rift Valley basin is facing the above types of challenges. Deforestation is a day-to-day activity for people living in and around the watershed due to increasing demand for charcoal, construction, domestic needs, expansion of arable land, and grazing areas (Degelo 2007). This continuous change in land use/cover is expected to impact the water balance of the watershed by changing the magnitude and pattern of the components of hydrological processes which are surface runoff, baseflow, interflow, and evapotranspiration, which results in increasing the extent of the water management problem.
So, studying about impacts of land use/cover change impacts on water balance for Bilate watershed was crucial to solving a wide variety of water resources problems, including design of hydraulic structures; urban and highway drainage; planning of flood-control works; source pollution; disposal of waste material; evaluation of environmental impacts of land use and management practices; planning of soil conservation works and agricultural products. This enables the local governments and policymakers to formulate and implement effective and appropriate response strategies to minimize the undesirable effect of future land use/cover change (PHE; Ethiopia Consortium 2011). Similar research was done by Schulla & Jasper (2012), Kebede et al. (2013), Hagos (2014, and Maria (2015) using the WaSiM model. Therefore, this research aims to determine the land use/cover change on the water balance of the Bilate Watershed using the WaSiM model.

Description of the study area
Bilate River is one of the inland Rivers of Ethiopia whose source is located at the Gurage Mountains in central Ethiopia. The Bilate Watershed is a part of the Main Ethiopian Rift valley basin which is part of the Great Rift Valley. The Bilate River watershed (BRW) covers an area of about 5,625 km 2 and is located in the Southern Ethiopian rift valley and partly in the western Ethiopian Highlands ( Figure 1) and its elevation is about 1,175 meters.

Input data's for WaSiM model
WaSiM allows various model configurations depending on the aim of the application and on the amount and quality of input data. It is possible to combine various sub-model components and to run the model in various spatial and temporal discretization.
(B) Soil data FAO/UNESCO-Soil Map of East Africa (2012), available in Arc/Info format with a scale of 1: 1,000,000, and were obtained from GIS and Remote Sensing Department, Ministry of Water, Irrigation, and Electricity of Ethiopia (MOWIE). These data were used as input for the WaSiM model.

(C) land use/cover data
The land use/cover image with three years of spatial resolutions of 30 meters (1989, 2002, and 2015) (https:// www.usgs.gov/core science systems/science analytics and synthesis/gap/science/land-cover-data) or https:// earthexplorer.usgs.govwere downloaded from USGS Earth Explorer and prepared using ERDAS software and ArcGIS. The WaSiM model was run for three different years' intervals of Landsat data (1989, 2002, and 2015) of land use such as Grassland/Pasture, Barren lands, Rangeland Scrublands, Cultivation/Agriculture, and Mixed Forest.
The following parameter values used as input in the WaSiM model were generated for those land uses within study period intervals: albedo, leaf area index (LAI), leaf surface resistance (Rsc), Intercept Cap, rs_evaporation, aerodynamic roughness length (Z0), vegetation cover fraction (VCF) and root depth within those different land use/cover types were obtained from different works of literature and website; http://www.unigiessen.de/ gh1461/plapada/php/list/contentor http://www.unigiessen.de/~gh1461/plapada/plapada.html.

(D) Meteorological data
The meteorological data required for this study was collected from the National Meteorological Agency of Ethiopia. The daily meteorological data that was collected for this study includes precipitation, Maximum and minimum temperature, relative humidity and wind speed, and sunshine hours from the year 1987-2015 in and around the watershed area for 4 stations as shown in Table 1.
The missing data of meteorological daily data were filled by using the Arithmetic mean values method hence; the total missed values were counted and compared with the data for each year, the percentage of missed data of all stations was less than 10%. Discharges of two gauging stations, Alaba Kulito and Tena (on Bilate River) are found in the watershed and daily flow data were collected from the Ministry of Water, Irrigation, and Electricity of Ethiopia for both gauging stations.
From the two gauging stations, the Bilate Tena gauging station was selected. Because BilateTen a gauging station was found in the same location as the outlet of the Bilate watershed. Flow data was required for performing sensitivity analysis, calibration, and validation of the model from 1989 to 2015 for the period of 27 years ( Table 2).

Preprocessing of data
2.2.2.1. Missing meteorological data estimation. The missing data of meteorological daily data were filled by using the Arithmetic mean values method hence; the total missed values were counted and compared with the data for each year, the percentage of missed data of all stations was less than 10%.
2.2.2.2. Consistency test for meteorological data. The data consists of the given meteorological stations was checked with the help of a double mass-curve method within reference to their neighborhood stations. It was tested using the XLSTAT 2017Software SNHT test (Amiri and von Rosen 2011).
In the case of the R statistic (R stands for Range), the null and alternative hypotheses are given by H 0 : The T variables are not homogeneous for what concerns their mean. Two-sided test: Ha the T variables follow one or more distributions that have the same mean. The double mass curve was used to check the consistency of the rainfall stations in the study area, and the analysis shows that the stations were consistent over the considered period.
2.2.2.3. Filling missed hydrological data. The daily flow data are archived based on m 3 /s and transformed into mm/time step before implementation into WaSiM Since the available meteorological and hydrological data cover the same period from 1989 to 2015 used.
The missed hydrological data were filled using Regression Equation (1) with a correlation coefficient of 0.82. A regression equation was used to fill the missed hydrological data using Equation (1) from Alaba Kulito (nearby  -1987-2015 gauging station) and with a correlation coefficient of 0.82.
where Y dependent variables X explanatory/independent variables 2.2.2.4. Areal rainfall. Areal rainfall is the average rainfall over an area, and is referred to as the areal rainfall distribution and is restricted to a long-term average value. It is expressed as a mean depth (mm) over the catchment area and used to know the distribution of rainfall for the calibration and validation period (Elizabeth 1994). Figure 2 shows the Thiessen polygon and areal proportion of each four selected stations in the sub-basins.

Land Use/Cover data preparation and processing (A) Landsat image processing
After delineating the watershed of the study area, land use/cover data preparation and processing is very crucial to have land cover data for the watershed. Landsat ETM þ was selected for the period of 1989, 2002, and 2015 respectively. To avoid a seasonal variation in vegetation pattern and distribution throughout a year, the selection Uncorrected Proof of dates of the acquired data was made as much as possible in the same annual season of the acquired years. The images used in this study area were orthorectified to a Universal Transverse Mercator projection using datum WGS (World Geodetic System) 84 zone 37N. The acquisition dates, sensor, path/row, resolution, and the producers of the satellite images used in this study are summarized in Table 3 below.

(B) Land Use/Cover Classifications
The Land use/cover change studies were differentiated using the available data source such as remote sensing, any other relevant information, and the previous local knowledge. Hence, based on the prior knowledge of the study area, ERDAS Imagine, and additional information from previous research in the study area (Degelo 2007;Wagesho 2014;Getahun 2017), the types of land use and land cover have been identified for the Bilate watershed. The descriptions of these land use and land covers are given as follows in Table 4.
All the three raster land use of the watershed were classified into six major types (Grassland/Pasture, Range and Shrub lands, cultivation/agriculture, Mixed Forest, Settlements, and Barren land). To differentiate the cultivation land from barren lands, the season of the land use /land cover downloading was selected in December and November months. To parameterize the land use in a distributed way, a land-use grid was required. The land use grid was parameterized with a land-use table that describes each grid cell with a parameter data set according to the grid classification. A specific value, which refers to a land-use type in the control-file, is assigned to each cell of this grid. The characteristics (e.g., root depth, resistance, LAI, VCF) of these types are declared in the land use table in the control-file. Most of the parameters describe a seasonal cycle with maximum (e.g., leaf area index -

Agricultural/cultivationlands
These include a diverse class of cultivated land, plots covered by food and commercial crops (croplands) and land units covered by residuals after immediate harvest.
Mixed-forest/Forest-lands Forestlands have usually tree crown areal density capable of modulating the micro climate and water holding capacity of the watershed. They range from densely populated tall trees of tropical rain forest used for timbering to moderately grown green forest. Forestlands could be evergreen, deciduous or mixed forestland.
Range and Shrub-lands Range lands are typical to arid and semiarid lands characterized by xerophytic vegetation and transition zones from forest land to sparse woodlands whereas Shrub lands are a plant community characterized by vegetation dominated by shrubs, often also including grasses, herbs, and geophytes.
Grass-lands/Pasture Grasslands are land units where the potential natural vegetation is predominantly grasses and grass like plants. It is dominated by naturally occurring grasses as well as those areas of actual range land that have been modified to include grasses whereas pasture land is an area covered with grass or other plants suitable for the grazing of livestock.
Water and Marshy land Area that remains water logged and swampy throughout the year, and rivers. But water or marshy land (Boyo Lake) was not considered for this study because there is no full data to construct a lake model module in WaSiM control file.

Barren land
Land of limited ability to support life and in which less than one-third of the area has vegetation or another cover. LAI) or minimum (e.g., stomata resistancersc) values during the vegetation period. It is easily shown that the increase of cultivation/agriculture, barren lands, and settlement area causes the decrease of mixed forest area, grasslands/pasture, and range and shrub land over the last 27 years for three LULC maps.

(C) Accuracy Assessment
The overall accuracy is used to indicate the accuracy of the whole classification (i.e. number of correctly classified pixels divided by the total number of pixels in the error matrix), whereas the other two measures indicate the accuracy of individual classes. User's accuracy is regarded as the probability that a pixel classified on the map represents that class on the ground or reference data, whereas the product's accuracy represents the probability that a pixel on reference data has been correctly classified.
In this study, the assessment was carried out using the original image for 1989 maps and the Google Earth Image for 2002, and 2015 together with previous knowledge of the area was used as reference data to generate testing data set. A total of 83, 85, and 85 testing sample points were selected randomly for the years 1989, 2002, and 2015 respectively. After completing the accuracy process as indicated in Table 6 the overall accuracy estimated as 87% is acceptable. The land cover vector data were converted into an appropriate raster format, grid size, and different land covers. The raster format of the land use map is converted to vector, to ASCII, and then to grid format which is required as input for the WaSiM hydrological model.
2.2.2.6. Soil data preparation. The parameterization of the soil's physical properties is crucial for any hydrological model application. The soil hydraulic properties, which are saturated and unsaturated hydraulic conductivity and water retention, control the main hydrological processes (Elsen 2001).
The watershed was discretized into five different soil types presented in Table 5. For the soil parameterization, the method of multiple soil horizons was used, where each soil type may have a different number of horizons. Each soil horizon has different hydraulic properties and may consist of a different number of layers of various thicknesses Figure 3.
The parameterization of the soil physical properties for each horizon was based on the Van Genuchten parameters after Wösten et al. (1999) and HWSD (Harmonized World Soil Data) Viewer was used to determine percentages of silt, clay, and sand in each layer of the soil.
2.2.2.7. WaSiM model setup. The control file of the model was adjusted as per the watershed characteristics and available input data. Meteorological input data of the model were interpolated for each grid cell in the watershed and are followed by the simulation of the main hydrological processes like evapotranspiration, interception, infiltration, and the separation of discharge into the direct flow, interflow, base flow, and total simulated flow. These calculations are modularly built and can be adapted to the physical characteristic of the watershed. All spatial data were prepared in a raster data set (grid) with a resolution of 30 m by 30 m.
2.2.2.8. Sensitivity analysis. The sensitivity analysis includes test runs in which the value of only one coefficient or parameter is changed, while the values of the others remain constant. The WaSiM model sensitive parameters for this study were selected from different findings which are done using the WaSiM model in different watersheds in Ethiopia and other basins out of our country. The sensitivity analysis was checked using the manual method by setting the values of the sensitive parameters on the WaSiM control file one after the other. From the model runs, sensitivity analysis mainly focused on the unsaturated zone model parameters, land use model (rsc or leaf surface resistance), and soil model (K recession) in the WaSiM-control file but the most sensitive parameters were found, unsaturated zone model as shown in Table 6.

Determining the water balance of the watershed
The water budget simulation section of WaSiM comprises a chain of modules that combine both the physical and empirical descriptions of water flow. In this study, constant climatic conditions under changing land use or land cover considered. The following components or model modules were used to calculate the water balance of watersheds.

Soil model. WaSiM versions with a physically based soil model use the RICHARDS-equation for
modeling the fluxes within the unsaturated soil zone. The modeling is done one-dimensional in the vertical direction using soil with several numeric layers. The continuity equation for this type of problem is given by: The dependencies of the hydraulic properties on the water content of the soil are considered discretely. The flux q between two layers with indices u (upper) and l (lower) is then given by: where: q flux between two discrete layers m/s;k eff: effective hydraulic conductivity m=s; h h hydraulic head, dependent on the water content and given as difference of geodetic altitude hgeo[m] and suction ψ(Θ) after equation; d thickness of the layers under consideration [m]; The basin under this study was subdivided into hydrological -sub-basins. The discharge at the outlet of each sub-basin was calculated by the -sub-models mentioned above then, the discharge at the outlet of the entire watershed area was calculated by routing the discharge of the individual -sub-basins through the interconnecting rivers and channels. The simulated water balance components which were generated from LULC_1989, LULC_2002, and LULC_2015, total simulated flow, interflow, base flow, potential evapotranspiration, real evapotranspiration, and infiltration were analyzed for those land use maps. Thus, the general methodology was continued by applying the three LULC map data (Landsat images) and analyzing the impacts of land use land cover change on water balance for three LULC maps. Finally, the water balance of the watershed area was determined using Equation (7) for Bilate Watershed for those land use land cover maps, and change in storage was also computed using Equation (8).
where, Q is the runoff, P the precipitation, ETR the evapotranspiration, and ΔS the change in soil storage. All variables were represented in mm per time step for the whole watershed area.

RESULTS AND DISCUSSION
3.1. Land Use/cover analysis

Overall accuracy, producer's accuracy, and user's accuracy
The accuracy assessment is used to determine the correctness of the classified image. It was performed using a confusion matrix. The overall accuracy gives the overall results of the confusion matrix. It is calculated by dividing the total number of correct pixels (diagonals) by the total number of pixels in the confusion matrix. The overall accuracy for the maps of 1989, 2002, and 2015 were 87, 80, and 91% respectively hence, they fulfill the minimum requirements. The producer's accuracy tells us how well a certain area can be classified. It is obtained by dividing the number of correctly classified pixels in the category by the total number of pixels of the category in the reference data. The producer's accuracy is also known as an Omission Error, which is the probability of reference pixels being classified correctly. It gives only the proportion of correctly classified pixels. The overall result of the producer's accuracy ranges from 69% to 93% as indicated in Table 7.
User's accuracy is the ratio between the total number of pixels correctly belonging to a class (diagonal elements) and the total number of pixels assigned to the same class by the classification procedure (row total). This quantity explains the probability that a pixel of the classified image truly corresponds to the class to which it has been assigned. In this study, the user's accuracy ranges from 80% to 93%.

Land use/ cover maps
The land use /land cover map of 1989 in Tables 8 and 9 shows that the total cultivated land/agriculture coverage class was about 28.92% of the total area of the watershed. It increased rapidly and became 34.12% and 43.3 of the watershed in 2002 and 2015 land use land cover maps respectively. This is mainly because of the population growth that caused the increase in demand for new cultivation land and settlement which in turn resulted shrinking of other types of land use land cover of the area. The forest coverage in 1989 was about 21.79% of the total area of the watershed. However, in the year the 2002 and 2015 was reduced to almost 12.79% and 6.4% of the total area respectively. This is most probably because of the deforestation activities that have taken place for agriculture and the expansion of settlements.
The individual class areas and change statistics for the three periods are summarized in Table 9. The results of previous studies showed that the same fact in the Bilate watershed. For example, Wagesho (2014) reported that cultivation and settlements of Bilate Watershed were increased by 61.6% and the mixed forest decreased to 67.7% from 1976 to 2000. Hence, the impacts of land use/cover change of the Bilate watershed were indicated in Figures 4 and 5.

. Sensitivity analysis
The executed sensitivity analysis mainly focused on the unsaturated zone model parameters. Table 10 shows the results of the sensitivity analysis. Both Q 0 and k B are quite sensitive to the base flow, as expected. Furthermore, k D is considerably sensitive to peak flow. Finally, the drainage parameter dr seems to be quite sensitive to the base flow. An increase in the k B value always results in an increased base flow value for low flow conditions at the beginning of the calibration period. Finally, the higher the k D value, the lower is the direct flow and when the value of K i increased, the value of the interflow becomes lower. A similar analysis was done by Maria (2015).

Model calibration and validation
The WaSiM sensitive parameters were identified/ selected from other authors' findings Schulla & Jasper (2012), Kebede et al. (2013), Hagos (2014), and Maria (2015. Finally, the sensitive parameters for this study were listed in Table 10 and the maximum, minimum, and optimum values of the sensitive parameters of this study were identified, as indicated in Table. As reported in Kebede et al. (2013), the parameters such as rs_evaporation (soil surface resistance for evaporation only) and rsc (leaf surface resistance) were also calibrated manually. The other parameters in the soil water dynamics of the WaSiM-ETH were K D, and K I (Table 10), which control the surface runoff and interflow storage effects in the Richards equation which was used for this study. Finally, the parameters Krec, dr, K D, and K I were found to be very sensitive. From all of these sensitive parameters, d r was the most sensitive parameter.
As the model analysis of this research indicated that the hydrographs were good at simulating the daily, weekly, and monthly scales. The monthly simulation indicates better than daily and weekly simulation for this study.       The analysis of the LULCC contribution was made on surface runoff; interflow, base flow, total simulated discharge, and evapotranspiration as characteristics of the hydrological process of the watershed. The contribution of surface runoff, total simulated flow, and interflow have increased from 1989 to 2015 indicated in Figure 11. This was related to the surface cover of the watershed since the forest land of the watershed was changed to agricultural land accelerates the runoff rate and reduced the infiltration of soil moisture content (Table 8). From the result of the land use land cover map, areas of forest have been decreased from 1989 to 2015 which has contributed to the increasing surface runoff contribution. In the same manner, Wagesho (2014) reported that the simulated surface runoff component increased progressively since the 1970s in Bilate Watershed. As Bahati et al. (2021), About the historical/current minimum, maximum, and mean annual flow of Muziz river, future minimum, maximum, and mean annual flow will increase respectively.
On the other hand, the rate of evapotranspiration has decreased from 1989 to 2015 indicating losses are mainly through evapotranspiration. These result revealed that the land use land cover change has significant impacts on infiltration rates, on runoff production, total simulation flow, interflow, base flow, evapotranspiration, and water retention capacity of the soil or change in storage of the soil, hence, it affects the water balance of watershed. This is because; the land cover under little vegetation is subjected to high surface runoff, low water retention, and low evapotranspiration (Tufa 2014). The changes in water balance (hydrological process) under the land use land cover changes are summarized in Figure 10.
As shown in Tables 12 and 13, the simulated water balance for the Bilate watershed using the WaSiM-ETH reveals the interflow component of the water balance takes a higher fraction of simulated discharge. The change in soil water storage ΔS is the result of the balance, being positive when the profile has a net gain of water, and negative for a net loss (Reichardt et al. 1995). From the results, the mean annual stream flows were evaluated due to land use/land cover change in the Bilate watershed shown in Figure 11.
The annual simulation of hydrological processes was analyzed for LULC_1989, LULC_2002, and LULC_2015 data. The result indicated that; there was a change in total simulation flow, evapotranspiration, surface runoff, base flow and inters flow in each land use land cover data (Table 12). In the study year intervals  change in total simulated flow increased by 77.83%, and direct runoff, interflow, and base flow increased by 80.23%, 75.69%, and 87.79% respectively. Hence, evapotranspiration is decreased by 6% throughout the study time.
Hydrological cycling in a watershed can be characterized and quantified by a water balance, which is the computation of all water fluxes at the boundaries of the system under consideration. It is an itemized statement of all  gains, losses, and changes of water storage within a specified elementary volume of soil. From this study, rainfall was considered gains where evapotranspiration, total simulated flow (runoff, interflow, and base flow) were considered as the losses.

CONCLUSION
In this study, we analyzed the impact of land use land cover change on the water balance of the Bilate watershed. As part of our analysis, we considered six dominant land use land covers including mixed forest, cultivation/agricultural land, barren land, grassland/pasture, range and shrub land, and settlement on the Bilate watershed for LULC_1989, LULC_2002, and LULC_ 2015. Like many before us, we found ArcGIS to be a very important tool for the preparation of input data for analyses and WaSiM model to be important for considering land use land cover data, soil, and DEM data. The advancement of computational power and the availability of spatial and temporal data have made hydrological models attractive tools to examine and analyze the characteristics of watersheds and how the hydrological process of the catchment functions under varying land-use dynamics. Particularly, in this study hydrological modeling is a useful tool for investigating interactions among the watershed components and hydrologic response analysis to LULCC at various spatial and temporal scales.
The simulated water balance for the Bilate watershed using the WaSiM-ETH showed the interflow component of the water balance takes a higher fraction of simulated discharge and also surface runoff and total simulation flow increased through the study period. Additionally, the ETR, ETP, and soil storage capacity decreased throughout the study periods. Overall, the hydrological model describes changes in the water balance from 1989 to 2015, which indicate that the change in total simulated flow increased by 77.83%, and direct runoff, interflow, and base flow increased by 80.23%, 75.69%, and 87.79% respectively. Additionally, evapotranspiration decreased by 6% throughout the study time.
The future sustainable land and water resources in the Bilate Watershed depend on the spatial planning of land use to optimize environmental benefits. Factors that must be considered include managing surface runoff control,  Uncorrected Proof erosion protection, flood protection, and water availability. Finally, educating the community on the effect of the unplanned land-use practices on the environment, natural resources, and ecosystem are of paramount importance for the future sustainability of the watershed.'