Impact of Land Use/Land Cover Change on Hydrology of the Catchment: The Case of Upper Ribb Catchment, Lake Tana Sub Basin, Ethiopia

Land use/land cover dynamics has been recognized as one of the factors responsible for altering the hydrologic response of watersheds. Various extent of water resources projects planning and implementation will require knowledge of these changes on catchment hydrology. This study assesses the impact of land use and land cover change effects on stream flow using SWAT model in Upper Ribb Catchment. The analysis of stream flow was done for three different years of 1973, 1995 and 2016 using SWAT model. During the study period of over 43 years (1973 – 2016) stream flow was increased to 6.143m/s due to an increase of cultivated land by 29.947%. Model calibration and validation for stream flow were done at Upper Ribb gauging station. The performance of the model was also checked at this station. Both the monthly calibration and validation results showed good match between measured and simulated stream flow data with the coefficient of determination (R) of 0.85, NashSutcliffe efficiency (NSE) of 0.824 for the calibration, and R of 0.829 and NSE of 0.814 of the validation period.


INTRODUCTION
The relationship between land use/cover change and watershed hydrology is very complicated, with linkages existing at a wide variety of spatial and temporal scales; but, land use change has a strong influence on global water yield. Land cover and use directly impact the amount of evaporation, groundwater infiltration and overland runoff that occurs during and after precipitation events. These factors control the water yields of surface streams and groundwater aquifers and thus the amount of water available for both ecosystem function and human use (Mustard et al., 2004).
The total environmental effects such as change in vegetation cover, soil characteristics, flora and fauna population and hydrological cycle have been strongly influenced by the conversion of land and forest resources (Hurni et al., 2010).
Changes in stream flow records were a result of changes in land cover observed in the area (Rientjes et al., 2011).
Land use/cover changes are highly pronounced in the developing countries that are characterized by agriculture based economies and rapidly increasing human populations. It is caused by a number of natural and human driving forces (Meyer and BL Turner, 1994). Natural effects such as climate change are only over a long period of time, whereas the effects of human activities are immediate and often direct. From the human factors population growth is the most important in Ethiopia (Tekle and Hedlund, 2000) as it is common in developing countries. Population growth has significant effect on land degradation, poverty and food insecurity in the northern Ethiopian highlands (Pender et al., 2001). Some 85% of the population lives in rural areas and directly depend on the land for its livelihood. This means the demands of lands are increasing as population increases.
Population growth causes degradation of resources particularly forests that rely on the available land. This increases the run off volume by decreasing infiltration.
Land use planning and management are closely related to the sustainability of water resources as changes of land use are linked with amount of water through relevant hydrological processes (Guo et al., 2008). To maintain water sustainability, effective methods and mechanisms should be used. Nowadays, the hydrological models are good to represent the hydrological characteristics (Surur, 2010).
There are two basic advantages using hydrological models instead of relying only on collected data (Droogers study area (Garede and Minale, 2014) and (Ephrem , 2011), seven different types of land use and land cover has been identified for the Upper Ribb Catchment. These are described as follows: Cultivated land; Areas used for crop cultivation, both annually and perennially. It is the area with standing crop, tree crops, and crop lands where the crops were harvested. Bush and Shrub land; Areas with bushes, shrubs and small trees, with little wood and mixed with some grasses. It includes plantation trees and scrub vegetation at the fringes of forest cover and areas dominated by scattered trees. Grazing land; Areas covered with grass used for grazing, as well as bare lands that have little grass or no grass cover. It also includes other small sized plant species. Water bodies; Areas with surface water in the form of ponds, reservoirs, lakes, streams, rivers and its main tributaries. Forest land; Land covered with dense trees which are mainly ever green forest land Urban and Settlement Area; Areas with low density to high density residential. It comprises both dispersed rural and urban settlement areas. Woody Savanna Land; Areas cover with wood mixed with higher grass cover including seasonal as well as permanent wetlands.
In this particular study the image was classified based on the pixel based Maximum Likelihood supervised classification method. This is due to the fact that unlike other classifiers it considers the spectral variation with in each category and the overlap that may occur among different classes (Campbell and Mortenson, 1989).
Accuracy assessment must be done to determine how well the classification process accomplished the task. The most widely used classification accuracy is in the form of error matrix which can be used to derive a series of descriptive and analytical statistics (Manandhar et al., 2009). The columns of the matrix depict the number of pixels per class for the reference data, and the rows show the number of pixels per class for the classified image. From this error matrix, a number of accuracy measures such as overall accuracy, user's accuracy, producer's accuracy and Kappa statistics are determined. 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 actually represents that class on the ground or reference data, whereas product's accuracy represents the probability that a pixel on reference data has been correctly classified. Kappa statics is a type of technique used in accuracy assessment. It expresses the agreement between two categorical data sets.
In this study, accuracy assessment was performed using the available and the Google Earth Image together with previous knowledge of the area which used as reference data to generate testing data set by generating certain random testing points. A total of 200 testing sample points were selected randomly for the recent year 2016 and accuracy assessment was done.

Stream flow data
Daily stream flow data is required for SWAT simulation result calibration and validation. The daily stream flow data  was found downstream of the Ribb dam outlet at Upper Ribb gauging station (844km 2 ) near Debretabor road and were collected from the Ministry of Water, Irrigation and Electricity of Ethiopia (MoWIE). Depending on the extent of calibration and validation, flow data was collected and arranged as per the requirement of SWAT model. Checking consistency and homogeneity of data A time series of hydrological data may exhibit jumps and trends called inconsistency and non-homogeneity (Yevjevich and Jeng, 1969). Inconsistency is a change in the amount of systematic error associated with the recording of data. It can arise from the use of different instruments and methods of observation. Non-homogeneity is a change in the statistical properties of the time series. Its causes can be either natural or man-made. These include alterations to land use, relocation of the observation station and implementation of flow diversions.
A consistent record is one where the characteristics of the record have not changed with time. Adjusting for gage consistency involves the estimation of an effect rather than a missing value.
Double Mass Curve (DMC) analysis was the method that was used to check homogeneity and consistency of rainfall as well for adjustment of inconsistent data. The curve is determined by plotting the cumulative values of observed time series of station for which consistency and homogeneity need to be checked on y-coordinate on versus cumulative value of observed time series of group of stations on x-axis and the station affected by trend, a break in slop of the curve would indicate that conditions have changed that location. The data series, which was inconsistency and non-homogenous, must adjusted to consistent and homogenous values by proportionality. This proportionality for the stations to be adjusted for consistency and homogeneity was using equation. The double mass curve plot below shows (Figure 3.2) five of station that found in and around the Upper Ribb catchment has better correlation because plot of cumulative annual rainfall of neighboring versus each station are aliened on single straight line so that correction for consistency and homogeneity will not be done i.e the observed precipitation data of the stations shows consistence and homogeneous except Addis Zemen station. However, the graph of Addis Zemen station was somewhat deviated from the original slope. Therefore, equation 3.1 was used for correcting Addis Zemen station to be homogeneous with other surrounding stations.

FIGURE 3. 2 DOUBLE MASS CURVE PLOT MADE FOR FIVE STATIONS Rainbow Homogeneity Checking
Rainbow is a software package for hydro meteorological frequency analysis and testing the homogeneity of historical data sets. It was used for checking the homogeneity of rainfall and flow data. RAINBOW homogeneity test is based on the cumulative deviations from the mean. The homogeneity of the data of a time series is tested by evaluating the maximum and the range of the cumulative deviations from the mean. In this study the flow and rainfall stations used were homogeneous at 90%, 95% and 99% probability as shown the following figure.

Addis zemen
Bahir Dar Debre Tabor Ebinat Gassay Before the consistency and homogeneity of data tests filling missing must be done. The missing of rainfall values required for SWAT input was filled by weather generator. However, filling the stream flow missing data is very important for calibration and validation of the SWAT model. Whenever data missing and insufficiency exists, some information transfer techniques that can be appropriate for filling the missed observations and extending records, such as area ratio, average value, and regression methods, were used. In this study there was small no of missing data between 1980 and 2004. Thus, filling of daily missed flow data were conducted by time series average value method.

SWAT Model Setup 3.2.1 Watershed Delineation
Delineation of the watershed and sub watershed from 30 m resolution DEM using Arc SWAT model watershed delineation tool was the first step in creating SWAT model input. First the soil map, LULC map and the DEM were projected into the same projection called UTM Zone 37N, which is projection parameter for Ethiopia. The watershed and sub watershed delineation process in Arc SWAT consists of five major steps, DEM setup, stream definition, outlet and inlet definition, watershed outlets selection and definition and calculation of sub basin parameters. After the DEM setup was completed and the mask data was provided on the DEM, the model automatically calculates the flow direction and flow accumulation. Consequently, stream networks, sub basin outlet, whole and sub watersheds and topographic parameters were calculated using the respective tools. The threshold based stream definition option in the stream definition was used to define the minimum size of the subbasins. As suggested by the Arc SWAT interface (the smaller the specified number of hectares, the more detailed the drainage network delineated by the interface) a minimum threshold area in hectares were selected. The watershed was sub-divided into 27 sub-watersheds using the default (minimum) threshold area which was 1356.29ha.

Hydrological Response Units (HRUs) Analysis
SWAT model require land use/cover, soil and slope data in order to determine the area and the hydrologic parameters of each land use and soil category simulated in each sub watershed. The land use/cover, slop and soil map were imported into the interface and reclassified. The SWAT database has only values of hydrological property parameters of the most common type of land use/cover classes. But, some of the land use/cover classes and their parametric values did not exist in SWAT default data base. Therefore, it was necessary to replace these classes with land use/cover classes of the SWAT database which have similar hydrological properties (Table 3.1). The soil map of the study area was reclassified according to Arc SWAT requirements. During the reclassifying process there was a problem of obtaining the values of soil parameters that represent physical and chemical properties of each soil class which were used as SWAT input data. Food and Agriculture Organization of the United Nations (FAO) soil classification system which was supported by other additional method was used to determine soil types and properties of each soil class. Partly the values of the parameters of hydrological properties have been determined by studying typical textural characteristics of an existing soil material and estimating their values by referring other similar previous works.
Slope classification was carried out using a multiple slope option into five number of slope classes. Next, all the reclassified three maps were overlaid. Then, the sub basins were divided into Hydrologic Response Units (HRUs) by assigning the threshold values of land use/cover, soil and slope percentage. When multiple HRUs was assigned to each sub basin the thresh hold level should be defined in which the user can specify sensitivities for land use/cover, soil and slope data that will be used to determine the number and kind of HRUs in each watershed. According to (Setegn et al., 2008), HRU definition with multiple options that account for 10% land use, 20% soil and 10% slope threshold combination gives a better estimation of runoff and sediment components. Therefore, for this study, HRU definition with multiple options that accounts for 10% land use, 20% soil and 10% slope threshold combination was used. Hence, there were created 76, 81 and 79 HRUs for land uses of 1973, 1995 and 2016 respectively.

Sensitivity Analysis, Calibration and Validation of the Model 3.3.1 Sensitivity Analysis
The results from SWAT simulation cannot be directly used for further analysis but instead used for further analysis to sufficiently predict the constituent stream flow should be evaluated through sensitivity analysis, calibration and validation of the model (White and Chaubey, 2005). After feeding the required input data for SWAT model, flow simulation was performed for 25 years of recording periods from 1980 to 2004. The first year flow record was used as a warm up period and the simulation then used for sensitivity analysis of hydrologic parameters and for calibration and validation of the model. When a SWAT simulation was taken place there was a discrepancy between measured data and simulated results. So, to minimize this discrepancy, it was necessary to determine the parameters which are affecting the results and the extent of variation. Hence, to check this, sensitivity analysis is one of SWAT model tool to show the rank and the mean relative sensitivity (MRS) of parameters identification and this step ordered to analysis. Sensitivity analysis helps to determine the sensitivity of parameters by comparing the output variance due to input variability. It also facilitates selecting important and influential parameters for a model calibration by indicating the parameters that shows higher sensitivity to the output due to the input variability. Therefore, the number parameters that can be involved for calibration will be less in number and influential. It also evaluates the model capacity and helps to understand the behavior of the system being modeled.
Sensitivity analysis was performed to determine the influence a set of parameters had on predicting total flow. It was performed on Twenty-six different SWAT flow parameters.
By applying default lower and upper boundary parameter values, the parameters were tested for sensitivity analysis for the simulation of the stream flow. The sensitivity analyses were run for flow parameters using Upper Ribb gauging station measured flow. In the analysis, the sensitive parameters of the stream flow of the watershed were identified. The parameters, which resulted from the analysis, were ranked with their category of classification according to the magnitudes of the mean relative sensitivity (MRS) values (Lenhart et al., 2002). Very high Therefore, based on the above classification parameters producing very high, high and medium MRS values gave high and prior attention for calibration process.

Model Calibration
After the sensitivity analysis result, model calibration was done to obtain related values of predicted output of  Vol.9, No.6, 2019 20 interest with that of measured data. Model calibration is a means of adjusting or fine tuning model parameters to match with the observed data as much as possible, with limited range of deviation accepted (Neitsch et al., 2002). It is the procedure of adjustment of parameter values of a model to reproduce the response of reality within the range of accuracy specified in the performance criteria. Parameters for adjustment are selected from those identified by sensitivity analysis. Additional parameters other than those identified in the sensitivity analysis are used primarily for calibration due to the hydrological processes naturally occurring in the watershed. However, sometimes it is necessary to change parameters in the calibration process other than those identified during sensitivity analysis because of the type of miss match of the observed variables and predicted variables (White and Chaubey, 2005). The process of adjustment can be done manually or automatic methods. The manual method is the most common, and especially recommended for the application of more complicated models in which a good graphical representation is a prerequisite (Refsgaard and Storm, 1990). A calibration of stream flow was carried out at Upper Ribb gauging station. This site was selected due to the availability of measured flow.
In this study manual and automatic calibration methods were applied. First the parameters were automatically calibrated by using automatic calibration tool built in Arc SWAT and by using SWAT CUP until the model simulation result becomes acceptable as per the model performance measures. Then, the final parameter values that were calibrated using previous two calibration methods were used as the initial values for the manual calibration procedure.
The statistical and graphical approaches were also be used to evaluate the SWAT model performance a number of times until the acceptable values were obtained for surface runoff independently. SWAT developers in (Santhi et al., 2001) assumed an acceptable calibration for hydrology at R² >0.6 and NSE > 0.5 and these values were considered in this study as a reference.

Model validation
Model validation is testing of calibrated model results with independent data set without any further adjustment at different spatial and temporal scales (Neitsch et al., 2002). In order to utilize the calibrated model for estimating the effectiveness of future potential management practices, the model tested against an independent set of measured data. Model calibration determines the best or at least a reasonable, parameter set while validation ensures that the calibrated parameters set performs reasonably well under an independent data set. Stream Flow validation was carried out at a station similar to the calibration for 7 years. The three statistical (R 2 , NSE & RSR) and graphical model performance measures used in calibration procedure were also used in validating stream flow.

Model Performance Evaluation
The evaluation of hydrologic model behavior and performance is commonly made and reported through comparisons of simulated and observed variables. Frequently, comparisons are made between simulated and measured stream flow at the catchment outlet. There are various methods to evaluate the model performance during the calibration and validation periods. Among those, the following performance evaluation criteria's were used in this study:

Nash-Sutcliffe Efficiency
The Nash-Sutcliffe efficiency (NSE) is a normalized statistic that determines the relative magnitude of the residual variance ("noise") compared to the measured data variance ("information"). NSE indicates how well the plot of observed versus simulated data fits the 1:1 line. NSE is computed using the following equation: Where Si = model simulated output; Oi = observed hydrologic variable; Omean = mean of the observations that the NSE uses as a benchmark against which performance of the hydrologic model is compared; and N = total number of observations. The value of NSE ranges from negative infinity to 1 (best. NSE value < 0 indicates the mean observed value is better predictor than the simulated value, which indicates unacceptable performance. While NSE values greater than 0.5, the simulated value is better predictor than mean measured value and generally viewed as acceptable performance (Santhi et al., 2001).

Coefficient of Determination
Coefficient of determination (R 2 ) is an indicator of the extent to which the model explains the total variance in the observed data. The R 2 value is an indicator of strength of linear relationship between the observed and simulated values. R 2 ranges from 0 (which indicates the model is poor) to 1 (which indicates the model is good), with higher values indicating less error variance, and typical values greater than 0.6 are considered acceptable (Santhi et al., 2001). The R 2 is calculated using the following equation: . 9.; . / Where, Smean = mean of the model simulations A major limitation of R 2 is that it describes the linear relationship between the two data sets, and one may obtain  (Table 4.1) shows that there was a dramatic increase of cultivated and grazing lands for the first period (1973 -1995) with +32.412 and +21.554 respectively. Conversely bush/shrub lands were decreased by 51.985% for this period. However, the bush/shrub lands show appreciable increase during the second period (1995 -2016) with +17.790%. On the contrary, the grazing land showed a significant decrease in the second period (1995 -2016) with 18.107%.
Previous similar studies in this watershed and other parts of the country also reflect similar results. For instance, (Garede and Minale, 2014) showed the cultivated & settlement land, shrub land and grassland coverage during 2011 were 70.43%, 14% and 7.58% respectively for the whole Ribb catchment in the north western part of Ethiopia. (Yeshaneh et al., 2013) stated that crop field coverage in Koga watershed in 2010 were 76.83%. (Geremew, 2013) shows that the cultivated area was increased by 45%, while forest, grassland, shrub land and water was decreased by 2%, 34%, 5.7% and 4.9% respectively from 1986 to 2001. (Hadgu, 2008) indicating a sharp reduction of natural habitats and an increase in agricultural land in the highlands of Tigray, northern Ethiopia over a period of 41 years (1964 -2005). He reported that shrub land was dominant in 1964 covering 46% of the area followed by woodland with coverage of 28% of the area. However, agricultural land was dominant in both 1994 and 2005 covering 34% and 40% respectively. The next dominant LULC types in 1994 and 2005 were shrub land with coverage of 21% and 39%. (Andualem and Gebremariam, 2015) reported that there was an increase of cultivated lands and a decrease of forest cover by 33.79 and 1.4 percent respectively in Gilgel Abbay watershed, north western Ethiopia from the periods 1986 -2011. (Bewket, 2003) identifies agricultural conversion of 79 % of the Riverine forests of the Chemoga watershed within the Blue Nile basin from 1957 to 1998. (Rientjes et al., 2011) also presented the agricultural land, shrub land and grassland coverage during 2001 were 62.7%, 8.9% and 8.8% respectively for Upper Gilgel Abay catchment.

Accuracy Assessment
An accuracy assessment was made by using a confusion matrix with 200 randomly selected points (Table 4.2) by using land use maps, ground truth points and Google Earth. Great importance was given to the representation of different LU/LC classes by these randomly chosen points.  Vol.9, No.6, 2019 24 The 2016 land use and land cover classification has showed, user's accuracy and producer's accuracy are greater than 85%, as well the overall accuracy of 96% (Table 4.2). These values indicate the land sat and the methodologies used were so accurate. The Kappa coefficient also calculated, with a value of K= 0.92 which indicated the classification is almost perfect since it is between 0.81 and 1.00 (Landis and Koch, 1977).

Stream Flow Modeling
For the Upper Ribb catchment sensitivity analysis was carried out using the daily observed flow at Upper Ribb gauging station. For this analysis, 26 parameters were considered and only 9 parameters were identified to have significant influence in controlling the stream flow in the watershed. The sensitive parameters identified were presented in the table below. After sensitivity analysis has been done, the calibration of SWAT model simulated stream flow was carried out with observed average monthly stream flow data. Manual Stream flow calibration was performed for the simulated results based on the sensitive parameters and calibration results of using Sequential Uncertainty Fitting program (SUFI).    FROM 1983-1996FOR 1995 The following figure showed that the values of the scatter plots of the measured and simulated monthly flows data for the calibration periods. There is a good linear correlation between the two datasets (measured and simulated). The measured and simulated average monthly flow for Upper Ribb Catchment was obtained. During the calibration period, they were 24.11and 26.42m 3 /s respectively. The measured and simulated average monthly flow for the validation period was 23.85 and 26.13m 3 /s respectively. These indicate that there is a reasonable agreement between the measured and the simulated values in both calibration and validation periods (Table 4.5). Therefore, these results of estimated stream flows indicate that SWAT model is good predictor of stream flow of Upper Ribb Catchment.
However, the calibrated and validated stream flow shows that simulated flow over predicts peak flows but under predicts all other times due to very high surface runoff in the watershed. Such too high surface runoff was due to low amount of soil water evaporation, which in turn controlled by soil evaporation compensation factor (ESCO).

Evaluation of stream flow due to land use and land cover change
One of the main objectives of this study was to estimate effect of the land use and land cover change on the stream flow at the study outlet. Therefore, by using the calibrated flow parameters (Table 4.4 ), which was obtained by using the monthly flow data at Upper Ribb gauging station(844km 2 ), and the land use map upstream of the reservoir (678km 2 ), calibrated and validated results of simulated annual average stream flow for the 1973, 1995 and 2016 land use are presented in Table 4.6. The results showed that an increased stream flow of the first simulation period. But, there is a decrease of stream flow of the second simulation periods.  (Table 4.6). Stream flows showed a higher increase in the first period with mean annual stream flow value of 8.25m 3 /s. This effect was due to the fact that cultivated lands and grazing lands were increased by 32.412% and 21.554% respectively. However, there was a decrease in the second period with mean annual stream flow value of 2.104m 3 /s. This was due to the cause that bush/shrub lands were increased by 17.790% and grazing lands were decreased by 18.107% though cultivated lands showed a slight decrease by 2.464%.
In general stream flows have increased throughout the study period for over 43 years period with a magnitude of 6.143m 3 /s due to an increase of cultivated land by 29.947%.

Conclusions
This study has focused the impact of land use and land cover changes on the hydrology of the Upper Ribb Catchment using SWAT model.
The result shows that Upper Ribb Catchment has experienced a significant response in stream flow due to land use and land cover over the past 43 years. The sensitivity analysis using SWAT model has resulted nine important parameters that control the stream flow of the studied catchment.
From those nine parameters mostly sensitive to change for calibration of stream flow were the soil evaporation compensation factor (ESCO), curve number (CN2), soil available water capacity (SOL_AWC), threshold depth of water in the shallow aquifer required for return flow (GWQMN) and maximum canopy storage (Canmx). ESCO ranks the first level in stream flow sensitive parameters due to the reason that soil water evaporation in the watershed controls other hydrological components like infiltration, lateral subsurface flow and surface runoff. Model calibration and validation have showed that the SWAT model simulated the flow quit satisfactorily. Performance efficiency indicators show that the SWAT model simulates well for stream flow in the Upper Ribb Catchment. Thus, the Nash-Sutcliffe model efficiency (NSE) and coefficient of determination (R²) are found to be 0.824 and 0.85 in calibration and 0.814 and 0.829 in validation for stream flow analysis.
After calibration and validation of the model, impacts of the land use and land cover change on stream flow was carried out. The changes in land use has resulted changes in stream flow, in which the expansion of cultivated lands results an increase of surface runoff. As a result the stream flow was increased from year to year during the 43 years period due to the conversion of bush/shrub lands to cultivation and grasslands. Over 43 years period (1973 -2016) an increase of cultivated land by 29.947% resulted in an increase of stream flow by 6.143m 3 /s.
To sum up stream flow has showed a direct relationship with cultivated land as a result they increased from year to year.
In this study I have attempted to show only the effect of land use and land cover change on stream flow in the Upper Ribb Catchment using the calibrated and having good performance model ''SWAT''. However, I would like to suggest for other future researchers, planners and policy makers of water resource projects in this watershed to consider the effect of climate change as well as other different management practices on stream flow. The gauging station used in this catchment has stream flow records of only up to 2004. Therefore, the time series flow records must be continued to be used for future calibration and any other analysis which will done on the Upper Ribb Catchment if that needs current flow data.