3.1 Weather and streamflow data

The daily weather variables used in this study for trend analysis and for driving the water-balance model are precipitation, minimum air temperature (T_min), maximum air temperature (T_max), relative humidity (RH), hours of sunshine (SH), and wind speed (WS). Weather data from 40 meteorological stations were obtained from the Ethiopian National Meteorological Services Agency (ENMSA) for the 1971–2010 period. Daily streamflow data for 25 gauging stations were collected from the Federal Ministry of Water, Irrigation and Electricity of Ethiopia for the same period 1971–2010. After screenings and rigorous analyses of the weather data, a considerable amount of time series data were found to be missing in most of the stations (see Table S1 in the Supplement). The occurrences of civil war and defective and outdated devices were the main causes for the missing data records. As a result, only the 15 stations (Fig. 1) in which rainfall data are relatively more complete proved to be suitable for trend analysis. Some 10 stations having complete climate variables, such as T_max, T_min, RH, WS, and SH, were used as input for the SWAT model (Fig. 1).

We resorted to spatial interpolation techniques, such as the inverse distance weighting (IDW) and linear regression (LR), to fill the gaps. Uhlenbrook et al. (2010) applied similar approaches or methods to the Gilgel Abbay sub-basin, which is the UBNRB's headwater. The selection and number of adjacent stations for interpolation are important for the accuracy of interpolated values. As mentioned by Woldesenbet et al. (2017a), different authors used different criteria to select neighbouring stations. Because of the relatively low number of network stations, a geographic distance of 100 km was considered for most stations when selecting neighbouring stations. If no station was located within 100 km of the target station, then the search distance was increased until at least one suitable station is reached. After the neighbouring stations were selected, the two methods (IDW and LR) were tested by means of cross-validation to fill in missing data sets. The candidate methods' performances were evaluated using the statistical metrics such as the root-mean-square error (RMSE), mean absolute error (MAE), correlation coefficient (R²), and percent bias (% bias) between observed and estimated values for the target stations. Equally weighted statistical metrics are applied to compare the performances of selected methods at target stations and to establish the ranking. A score was assigned to each candidate method according to the individual metrics. For example, the candidate achieving the smallest values of RMSE and MAE or the smallest percentage of bias received score 1, and score 2 was assigned to the one with the larger value. The final score is obtained by summing up the score pertaining to each candidate approach at each station. The method with the smallest score is the best. The monthly, seasonal, and annual weather data were aggregated from the daily time series data after filling the gaps. While filling in the missing data, uncertainty is expected due to low station density, poor correlations, and the considerable number of missing records. Similar techniques and approaches were used for the analysis and filling in of missing streamflow data records.

3.2 Spatial data

Spatially distributed data required for the SWAT model includes tabular and spatial soil data, tabular and spatial LULC information, and elevation data. A Shuttle Radar Topographic Mission digital elevation model (SRTM DEM) with a resolution of 90 m from the Consultative Group on International Agricultural Research – Consortium for Spatial Information (CGIAR-CSI; http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp, last access: 21 February 2015) – was used to represent land-surface drainage patterns. Terrain characteristics, such as slope gradient and the slope length of the terrain, and stream network characteristics, such as channel slope, length, and width, were derived from the digital elevation model (DEM).

The soil map (1:5 000 000) developed by the Food and Agriculture Organization of the United Nations (FAO-UNESCO) was downloaded from http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/faounesco-soil-map-ofthe-world/en/, last access: 4 April 2017. Soil information, such as the soil textural and physiochemical properties needed for the SWAT model, were extracted from the Harmonized World Soil Database v1.2, a database that combines existing regional and national soil information (http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-databasev12/en/, last access: 4 April 2017) with information provided by the FAO-UNESCO soil map (Polanco et al., 2017).

The LULC maps were produced from satellite-remote-sensing Landsat images for 1973, 1985, 1995, and 2010 at a scale of 30 m × 30 m resolution. Detailed descriptions on image processing and classification are available under Sect. 4.2.

4.1 Trend analysis

The non-parametric Mann-Kendall (MK; Kendall, 1975; Mann, 1945) statistic is chosen to detect trends for rainfall and streamflow time series data, as it is widely used for water resource planning, design, and management (Yue and Wang, 2004). Its advantage over parametric tests such as t test is that the MK test is more suitable for non-normally distributed and missing data, which are frequently encountered in hydrological time series (Yue et al., 2004). However, the existence of positive serial correlation in time series data affects the MK test results. If serial correlation exists in time series data, the MK test rejects the null hypothesis of no trend detection more often than specified by the significance level (Von Storch, 1995).

Von Storch (1995) proposed a prewhitening technique to limit the influence of serial correlation on the MK test. The effective or equivalent sample size (ESS) method developed by Hamed and Rao (1998) has also been proposed to modify the variance. However, the study by Yue et al. (2002) reported that von Storch's prewhitening is effective only when no trend exists, and the ESS approach's rejection rate after modifying the variance is much higher than the actual (Yue et al., 2004). Yue et al. (2002) then proposed trend-free prewhitening (TFPW) prior to applying the MK trend test in order to minimise its limitation. This study therefore employed TFPW to remove the serial correlation and to detect a trend in time data series with significant serial correlation. Further details can be found in (Yue et al., 2002). All the trend results in this study have been evaluated at the 5 % level of significance to ensure the effective exploration of the trend.

4.2 Change point test

The Pettitt test is used to identify whether or not there is a point change or jump in the data series (Pettitt, 1979). This method detects one unknown change point by considering a sequence of random variables $\left(X_t\right)=X_{\mathrm{1}},X_{\mathrm{2}},\mathrm{\dots },X_N,X_N+\mathrm{1},\mathrm{\dots },X_T$ that may have a change point at N, if the X_t variable for $t=\mathrm{1},\mathrm{2}$, …, N time step has a common distribution function, F1(x), and X_t for $t=N+\mathrm{1}$, …; T time step has a common distribution function, F2(x), where F1(x) ≠ F2(x).

4.3 Sen's slope estimator

The trend magnitude is estimated using a non-parametric median-based slope estimator proposed by (Sen, 1968), as it is not greatly affected by gross data errors or outliers and can be computed when data are missing. The slope estimation is given by

where x_j and x_k are the sequential data values, and n is the number of the recorded data. $\mathrm{1}<k<j<n$, and β is considered as the median of all possible combinations of pairs for the whole data set. A positive value of β indicates an upward (increasing) trend, and a negative value indicates a downward (decreasing) trend in the time series. All MK trend tests, Pettitt change-point detections, and Sen's slope analyses were conducted using the XLSTAT add-in tool from Excel (https://www.xlstat.com).

4.4 Remote-sensing land-use and land-cover map

4.5 SWAT hydrological model

The SWAT is an open-source-code, semi-distributed model with a large and growing number of model applications in a variety of studies, ranging from catchment to continental scales (Allen et al., 1998, 2012; Neitsch et al., 2002). It enables the impact of LULC change and climate change on water resources to be evaluated in a basin with varying soil, land-use, and management practices over a set period of time (Arnold et al., 2012).

In SWAT, the watershed is divided into multiple sub-basins, which are further subdivided into hydrological response units (HRUs) consisting of homogeneous land-use management, slope, and soil characteristics (Arnold et al., 1998; Arnold et al., 2012). HRUs are the smallest units of the watershed in which relevant hydrologic components, such as evapotranspiration, surface runoff and peak rate of runoff, groundwater flow, and sediment yield, can be estimated. Water balance is the driving force behind all of the processes in the SWAT calculated using Eq. (7),

where SW_t is the final soil-water content (mm H₂O), SW_o is the initial soil-water content on day i (mm H₂O), t is the time (days), R_day is the amount of precipitation on day i (mm H₂O), Q_s is the amount of surface runoff on day i (mm H₂O), Q_l is the amount of return flow on day i (mm H₂O), Q_b is the return flow from the shallow aquifer on day i (mm H₂O), E_a is the amount of evapotranspiration from the canopy and soil surface on day i (mm H₂O), Revap is the amount of water transferred from the underlying shallow aquifer reversely upward to the soil-moisture storage on day i (mm H₂O) in response to water demand for evapotranspiration, and DA_recharge is the amount of water recharge to the deep aquifer on day i (mm H₂O).

Figure 2Schematic representation of the SWAT model structure modified from Marhaento et al. (2017). P is precipitation; CS is canopy storage; TAE is total actual evapotranspiration; E_c is evaporation from the canopy surface; E_s is evaporation from the soil surface; E_t is transpiration from plants; Perc is percolation from the soil storage to shallow aquifer; SM is soil-moisture storage; SA is shallow aquifer; E_r = Revap is evaporation from the shallow aquifer; Q_t is total streamflow; DA is deep aquifer; HRU is hydrological response unit; Q_b is base flow; Q_l is lateral flow; and Q_s is surface runoff.

Download

Runoff is calculated separately for each HRU and routed to obtain the total streamflow for the watershed using either the soil-conservation-service (SCS) curve-number (CN) method (Mockus, 1964) or the Green–Ampt infiltration method (GAIM; Green and Ampt, 1911; see Fig. 2). However, spatial connectivity and interactions among HRUs are ignored. Instead, the cumulative output of each spatially discontinuous HRU at the sub-watershed outlet is directly routed to the channel (Pignotti et al., 2017). This lack of spatial connectivity among HRUs makes implementation and the impact analysis of spatially targeted management such as the soil-and-water conservation structure difficult to incorporate into the model. Different authors have made efforts to overcome this problem, for instance, a grid-based version of the SWAT model (Rathjens et al., 2015) or a landscape simulation on a regularised grid (Rathjens and Oppelt, 2012). Moreover, Arnold et al. (2010) and Bosch et al. (2010) further modified SWAT so that it allows landscapes to be subdivided into catenas comprising upland, hillslope, floodplain units, and flow to be routed through these catenas. However, the SWAT grid, developed to overcome this limitation, remains largely untested and computationally demanding (Rathjens et al., 2015).

Hence, the standard SWAT CN method was chosen for this study, because it was applied in many Ethiopian watersheds (Gashaw et al., 2018; Gebremicael et al., 2013; Setegn et al., 2008; Woldesenbet et al., 2017b). Furthermore, its ability to use daily input data (Arnold et al., 1998; Neitsch et al., 2011; Setegn et al., 2008) as compared to the GAIM, which requires subdaily precipitation as a model input that can be difficult to obtain in data-scarce regions like the UBNRB. This study focused on the effects of LULC change and climate change on the basin's water-balance components, which include the components of inflows, outflows, evapotranspiration, losses, and the change in storage as shown in the general water balance in Eq. (8).

where $Q_{\mathrm{t}}=Q_{\mathrm{s}}+Q_{\mathrm{l}}+Q_{\mathrm{b}}$ and total actual evapotranspiration TAE = $E_{\mathrm{c}}+E_{\mathrm{s}}+E_{\mathrm{t}}+E_{\mathrm{r}}$, as shown in Fig. 2.

R is the amount of precipitation (mm d⁻¹) as the main inflow, Q_t is the total amount of streamflow (mm d⁻¹) as outflow, TAE is the total actual evapotranspiration (mm d⁻¹), E_c is evaporation from the canopy surface (mm d⁻¹), E_t is the amount of plant transpiration (mm d⁻¹), E_s is evaporation from the soil surface (mm d⁻¹), E_r or Revap is evaporation from the shallow aquifer (mm d⁻¹; Abiodun et al., 2018), losses are the amount of water lost from the system as a recharge to the deep aquifer (DA_recharge; mm d⁻¹), and ΔS is the change in soil-water storage (mm d⁻¹). SWAT has four storages: canopy storage (CS), soil moisture (SM), shallow aquifer (SA), and deep aquifer (DA). Water movement from the soil-moisture storage to the shallow aquifer is due to percolation, whereas water movement from the shallow aquifer reversely upward to the soil-moisture storage is Revap, and further water movement from the shallow aquifer to the deep aquifer is recharge. For a more detailed description of the SWAT model, refer to Neitsch et al. (2011).

The SWAT model set-up and data preparation were done using the ArcSWAT2012 tool in the ArcGIS environment, whereas parameter sensitivity analysis and model calibration and validation were performed using the SWAT-CUP (Calibration and Uncertainty Procedures) interface Sequential Uncertainty Fitting (SUFI-2) algorithm (Abbaspour, 2008). During model set-up, the observed daily weather and streamflow data from the given period were divided into three different periods: the first to warm up the model, the second to calibrate it, and the third to validate it. Determining the most sensitive parameters is the first step in the model calibration and validation process using the global sensitivity analysis option (Arnold et al., 2012). The second step is to complete the calibration process, making necessary adjustments for the model's input parameters to match model output with observed data, thereby reducing the prediction uncertainty. Initial parameter estimates were taken from the default lower and upper bound values of the SWAT model database and from earlier studies in the basin, such as that of Gebremicael et al. (2013). The final step, model validation, involves running a model using parameters that were determined during the calibration process and comparing the predictions to independently observed data not used in the calibration.

Table 2Data sets of the baseline and altered periods for the SWAT simulation used to analyse the combined and isolated effect of LULC and climate change on streamflow and water-balance components.

Download Print Version | Download XLSX

In this study, both manual and automatic calibration strategies were applied to attain the minimum differences between observed and simulated streamflows in terms of surface flow, peak flow, and total flow, following the steps recommended by Arnold et al. (2012). For the purpose of impact analysis, we divided the simulation period 1971–2010 into four decadal periods, hereafter referred as the 1970s (1971–1980), 1980s (1981–1990), 1990s (1991–2000), and 2000s (2001–2010), as shown in Table 2. The model's performance for the streamflow was then evaluated using statistical methods (Moriasi et al., 2007) such as the Nash–Sutcliffe coefficient of efficiency (NSE), the coefficient of determination (R²), and the relative volume error (RVE %), which are shown in Eq. (9)–(11). Furthermore, graphical comparisons of the simulated and observed data, as well as water-balance checks, were used to evaluate the model's performance. This evaluation may be represented by the following:

where Q_m,i is the measured streamflow in m³ s⁻¹, ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{m}}$ are the mean values of the measured streamflow (m³ s⁻¹), Q_s,i is the simulated streamflow in m³ s⁻¹, and ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{s}}$ are the mean values of simulated data in m³ s⁻¹.

4.6 SWAT simulations

Three different approaches were applied for assessing the effects of LULC change and climate change on streamflow and water-balance components. The first approach is to assess the response of streamflow to the combined effects of LULC change and climate change. We followed the approach in Marhaento et al. (2017) and divided the analysis period, 1971–2010, into four periods of similar length (four decades). These are periods when land-use changes are expected to change the hydrological regime within a catchment (Marhaento et al., 2017; Yin et al., 2017). The first period, the 1970s, was regarded as the baseline period. The other periods, the 1980s, 1990s, and 2000s, were regarded as altered periods. LULC maps of 1973, 1985, 1995, and 2010 were used to represent LULC patterns during the 1970s, 1980s, 1990s, and 2000s, respectively. For analyses, the SWAT model was calibrated and validated for each respective period using the respective LULC map and weather data (Table 2). The DEM and soil data sets remained unchanged. The differences between the simulation result of the baseline and altered periods represent the combined effects of LULC and climate change on streamflow and water-balance components.

The second approach included simulations that attribute the effects from LULC changes alone. It aimed to investigate whether LULC change is the main driver for changes in water-balance components. To identify the hydrological impacts caused solely by LULC, a “fixing–changing” method was used (Marhaento et al., 2017; Woldesenbet et al., 2017b; Yan et al., 2013; Yin et al., 2017). The calibrated and validated SWAT model and its parameter settings in the baseline period were forced by weather data from the baseline period, 1973–1980, while changing only the LULC maps from 1985, 1995, and 2010, keeping the DEM and soil data constant (Hassaballah et al., 2017; Marhaento et al., 2017; Woldesenbet et al., 2017b; Yin et al., 2017). We ran the calibrated SWAT model for the baseline period (1970s) four times, only changing the LULC map for the years 1973, 1985, 1995, and 2010 and retaining the constant weather data set from the 1970s (Table 2). The third approach is similar to the second, but the simulations are attributed only for climate change. The calibrated models for the baseline period were run again four times, corresponding to the LULC periods using a unique LULC map of the year 1973 but altering the four different periods of weather data sets for their respective periods.

5.1 Trend test

5.1.1 Rainfall

The summary of the MK trend test results for the rainfall recorded at the 15 selected stations located in and around the UBNRB revealed mixed trends (increasing, decreasing, and no changes). For daily time series, the computed probability values (p values) for seven stations were greater, although for eight stations, they were less than the selected significance level (α=5 %). This means that no statistically significant trend existed in seven stations, but a monotonic trend occurred in the remaining eight. Positive trends developed only at six stations, four of which were concentrated in the northern and central highlands (Bahir Dar, Dangila, Debre Markos, and Gimijabet Mariam). The other two stations, Assosa and Angergutten, are located in the south-western and southern lowlands (see Fig. 1). The other two stations, Alem Ketema and Nejo, which are located in the east and south-west of the UBNRB, respectively, showed a decreasing trend. On a monthly basis, the MK trend test result showed that no trend existed in 11 stations, while statistically non-significant increasing trends exist in 3 stations (Dangila, Gimijabet Mariam, and Shambu) and a decreasing trend exists in Alem Ketema station. On an annual timescale, MK trend test could not find any trend in 11 stations, but did exhibit a trend in 4 stations. The Debre Markos and Shambu stations showed statistically non-significant increasing trend, while Gimijabet Mariam and Alem Ketema showed statistically significant positive trends and non-significant decreasing trends, respectively. The trend analysis result for the annual rainfall time series agrees well with a previous study by Gebremicael et al. (2013), who reported no significant annual rainfall change at eight out of nine stations during the 1973–2005 period. Hence, it is interesting to note that the timescale of analysis is a critical factor in determining the given trends.

Table 3MK and Pettitt tests for the UBNRB's rainfall and streamflow after TFPW at different timescales.

^∗ Before and after TFPW; p probability at 5 % significance level.

Download Print Version | Download XLSX

Figure 3The Pettitt homogeneity test's (a) annual rainfall, (b) annual flow of the UBNRB, (c) linear trend of mean annual rainfall, and (d) linear trend of mean annual streamflow, where mu is the mean annual precipitation (mm), mu1 is the mean annual precipitation before change (mm), and mu2 is the mean annual precipitation after change (mm).

Download

The basin-wide areal UBNRB rainfall trend and change point analysis was again carried out on daily, monthly, seasonal, and annual timescales using the MK and Pettitt tests. We applied a widely used spatial interpolation technique, the Thiessen polygon method, to calculate basin-wide rainfall series from station data. A summary is provided in Table 3 and Fig. 3. The MK test showed increasing trends for annual, monthly, and long rainy-season rainfall series, whereas no trend for daily, short rainy, and dry season rainfall series appeared. The magnitude of trends for annual, monthly, and long rainy-season rainfall series are not statistically significant, as explained by the values of Sen's slope. However, the Pettitt test could not detect any jump point in basin-wide rainfall series except for daily time series rainfall (see Fig. S1 in the Supplement).

Previous studies (Conway, 2000; Gebremicael et al., 2013; Tesemma et al., 2010) conducted trend analyses of basin-wide rainfall and reported that no significant change in annual and seasonal rainfall series across the UBNRB exists, which contradicts the results of this study. This disagreement could be due to the number of stations and their spatial distribution across the basin, time period of the analysis, approach used to calculate basin-wide rainfall from gauging stations, and data sources. Tesemma et al. (2010) used monthly rainfall data downloaded from Global Historical Climatology Network (GHCN) database and 10 day rainfall data for the 10 selected stations obtained from the National Meteorological Service Agency of Ethiopia from 1963 to 2003. Conway (2000) also constructed basin-wide annual rainfall in the UBNRB for the 1900–1998 period from the mean of 11 gauges. Furthermore, (Conway, 2000) employed simple linear regressions over time to detect trends in annual rainfall series without removing the serial autocorrelation effects. Gebremicael et al. (2013), used only nine stations from the 1970–2005 period. However, in this study, we used daily observed rainfall data from 15 stations collected from Ethiopian Meteorological Agency from 1971 to 2010. The stations are more or less evenly distributed over the UBNRB.

5.1.2 Streamflow

The MK test's result for daily, monthly, annual, and seasonal (long and short rainy season and dry season) streamflow time series showed a positive trend, the magnitude of which is statistically significant, as summarised in Table 3. The Pettitt test also detected change points for daily, annual, and short rainy-season streamflow but did not detect change points for monthly, long, and dry season streamflow (see Figs. 3 and S2). The change point detected by the Pettitt test for annual streamflow is in 1995, whereas for daily and short rainy seasons, it is in 1985 and 1987, respectively. The result obtained from the MK test agrees well with the findings in Gebremicael et al. (2013), who reported an increasing trend in the observed annual, short, and long rain seasons' streamflow at the El Diem gauging station, but our results disagree with findings for dry-season streamflow. Furthermore, the increasing trend of long rainy-season streamflow agrees well with the result of Tesemma et al. (2010), but it disagrees with the results of short rainy season and annual flows. Tesemma et al. (2010) reported that the short rainy season and the annual flows are constant for the 1964–2003 period. This disagreement is likely attributable to the difference in analysis period, as can be seen from Fig. 3. The last 7 years, 2004–2010, had relatively higher streamflow records.

Although, the results of the MK test for annual and long rainy-season rainfall and streamflow show an increasing trend for the last 40 years in the UBNRB, the magnitude of Sen's slope for streamflow is much greater than it is for rainfall (Table 3). Moreover, short rainy-season streamflow shows a statistically significant positive increase, whereas the rainfall shows no change. The mismatch between the rainfall and the streamflow trend magnitude could be associated with evapotranspiration and could be attributable to the combined effect of LULC change, climate change, the infiltration rate due to changing soil properties, rainfall intensity, and extreme events.

Table 4Confusion (error) matrix for the 2010 land-use and land-cover classification map.

Download Print Version | Download XLSX

5.2 LULC change analysis

According to the confusion-matrix report, with an overall accuracy of 80 %, a producer's accuracy values for all classes ranged from 75.4 % to 100 %, user's accuracy values ranged from 83.7 % to 91.7 %, and a kappa coefficient (k) of 0.77 was attained for the 2010 classified image (Table 5). Monserud (1990) suggested a kappa value of <40 % as poor, 40 %–55 % as fair, 55 %–70 % as good, 70 %–85 % as very good, and >85 % as excellent. According to these ranges, the classification in this study has very good agreement with the validation data set and meets the minimum accuracy requirements to be used for further change detection and impact analysis.

Figure 4Land-cover map of UBNRB derived from Landsat images from (a) 1973, (b) 1985, (c) 1995, and (d) 2010a (b).

Download

Figure 5(a) LULC composition and (b) LULC change in the UBNRB during the period from 1973 to 2010.

Download

The classified images of the basin (Fig. 4) have shown different LULC proportions at four distinct time periods, as shown in Fig. 5. In 1973, cultivated land dominantly covered (62.9 %) the UBNRB, followed by bushes and shrubs (18 %), forest (17.4 %), and water (1.74 %). In 1985, cultivated land area increased to 65.6 %, followed by bushes and shrubs (18.3 %), while forest decreased to 14.4 %, and water remained unchanged at 1.7 %. In 1995, cultivated land area further increased to 67.5 %, followed by bushes and shrubs (18.5 %). Forest further decreased to 12.2 %, and water remained unchanged at 1.7 %. In 2010, cultivated land decreased to 63.9 %, bushes and shrubs increased to 18.8 %, forest increased to 15.6 %, and water remained unchanged at 1.7 %. During the entire 1973–2010 period, cultivated land, along with bushes and shrubs, remained the major proportions compared to the other LULC classes. The highest increase (2.7 %) and the largest decrease (−3.6 %) in cultivated land occurred during the 1973–1985 and 1995–2010 periods, respectively. The largest increase in bushes and shrubs was 0.3 % from 1973 to 1985, whereas the largest increase in forest coverage (3.4 %) was recorded during the 1995–2010 period. Water coverage remained unchanged from 1973 to 2010.

Although, the image classification results show very good accuracy, uncertainties in classification could be expected. First, as elsewhere in Ethiopia, LULC may change rapidly over the land surface of the basin, and image reflectance may be confusing due to the topography and variation in the image acquisition date. Landsat images were not all available for a particular year or season (as described under Sect. 4.2.1); images from different years and different seasons might harbour errors. Secondly, the workflow associated with LULC classification involves many steps and can be a source of uncertainty. The errors are observed in the classified LULC map, as shown in Fig. 4. On the western side of the map in Fig. 4a, a rectangular section with forest appears, which completely disappears in Fig. 4b. Rectangular forest cover appears in the northern part of the country in Fig. 4b, which again disappears completely in Fig. 4c. In Fig. 4d, forest cover with linear edges (north–south) appears on the map's eastern side. That being recognised, the land-cover mapping is reasonably accurate overall, providing a good base for land-cover estimation and for providing basic information for the hydrological impact analysis.

The rate of expansion of cultivated land before 1995 was higher than after 1995. Conversely, the area of the forest land decreased in 1985 and 1995 with reference to the 1973 baseline. However, after 1995, the forest's size increased again, whereas cultivated land decreased. The increased forest coverage and the decrease in cultivated land over the period 1995 to 2010 showed that the environment was recovering from the devastating drought, and forest clearing for firewood and cultivation due to population growth has been minimised. This could be due to the afforestation programme, which the Ethiopian government initiated, and to the extensive soil-and-water conservation measures carried out by the community. Since 1995, eucalyptus tree plantations expanded significantly across the country at the homestead level for firewood, construction materials, charcoal production, and income generation (Woldesenbet et al., 2017b). In summary, forest coverage decreased by 1.8 %, while both bushes and shrubs as well as cultivated land increased by 0.8 % and 1 %, respectively, from the original 1973 level to the 2010 period. This result agrees well with other studies (Gebremicael et al., 2013; Rientjes et al., 2011; Teferi et al., 2013; Woldesenbet et al., 2017b), which reported a significant conversion of natural vegetation cover into agricultural land.

5.3 SWAT model calibration and validation

The SWAT model's most sensitive parameters for simulating streamflow were identified using a global sensitivity analysis of SWAT-CUP. Their optimised values were determined by the calibration process that Arnold et al. (2012) recommended. Parameters such as the SCS curve number (CN2), base-flow alpha factor (ALPHA_BF), soil evaporation compensation factor (ESCO), threshold water depth in the shallow aquifer required for return flow to occur (GWQMN), groundwater “Revap” coefficient (GW_REVAP), and available water capacity (SOL_AWC) were found to be the most sensitive parameters for the streamflow predictions.

Table 5The SWAT model's statistical performance measure values.

Download Print Version | Download XLSX

Figure 6Calibration and validation of the SWAT hydrological model (left and right, respectively) at monthly timescale: (a) 1970s, (b) 1980s, (c) 1990s, and (d) 2000s.

Download

Figure 6 shows the calibration and the validation results for monthly streamflow hydrographs for each model. These results revealed that the model represents the monthly hydrographs well, as also indicated by R², NSE, and RVE (%) statistical performance measures (Table 6). For the calibration period, the values of R², NSE, and RVE (%) range from 0.79 to 0.91, 0.74 to 0.91, and −3.4 % to 4 %, respectively. For the validation period, they ranged from 0.84 to 0.94, 0.82 to 0.92, and −7.5 % to 7.2 %, respectively. According to the rating of Moriasi et al. (2007), the SWAT model's performance in the UBNRB can be categorised as very good, although underestimation was observed in the base-flow simulation. The optimal parameter values of the four calibrated-model runs are shown in Table 7. A change was obtained for CN2 parameter values, which can be attributed to the catchment's response behaviour. For instance, an increase in the absolute average (basin-wide) CN2 value from the 1970s to the 1980s and 1990s, from 72.9 to 74.7 and 75.6, respectively, indicates a reduction in forest coverage and an expansion of cultivated land. By contrast, a decrease in the CN2 value was attained during the period 1990s to 2000s, from 75.6 to 73.6, attributed to the increase in forest coverage and reduction in cultivated land.

5.4 Combined effects of LULC change and climate change on streamflow and water-balance components

The simulation results of the four independent, decadal-time-scale calibrated and validated SWAT model runs reflect the combined effect of both LULC and climate change during the past 40 years (Table 8). From the simulation result, mean annual streamflow increased by 16.9 % between the 1970s and the 2000s, while the observed mean annual streamflow increased by 15.3 % for the same period. However, the rate of change is different in different decades. For example, it increased by 3.4 % and 9.9 % during the 1980s and 1990s, respectively, from the baseline 1970s period.

The ratio of mean annual streamflow to mean annual precipitation (Q_t∕P) increased from 19.4 % to 22.1 %, while the actual evaporation to precipitation (E_a∕P) ratio decreased from 61.1 % to 60.5 % from the 1970s to 2000s. Moreover, the ratio of surface runoff to streamflow (Q_s∕Q_t) increased notably, from 40.7 % in the 1970s to 50.1 % and 55.4 % in the 1980s and 1990s, respectively, and it decreased to 43.7 % in the 2000s. In contrast, the base-flow-to-streamflow ratio (Q_b∕Q_t) notably decreased from 17.1 % in the 1970s to 10.3 % and 3.2 % during the 1980s and 1990s, respectively, but it increased to 20 % in the 2000s. The result for surface runoff agrees with findings in Gebremicael et al. (2013), but disagreement is observed for base flow. The study reported that the surface runoff (Q_s) contribution to the total river discharge increased by 75 %, while the base flow (Q_b) decreased by 50 % from the 1970s to 2000s.

Table 6Summary of the UBNRB's precipitation indices at decadal time series.

SDII is the ratio of total precipitation (mm) to R1 mm (days).

Download Print Version | Download XLSX

In general, 1.8 % forest cover loss and 1 % increased cultivated land combined with 2.2 % increased rainfall from the 1970s to the 2000s led to a 16.9 % increase in simulated streamflow. The 1990s was the period during which the greatest deforestation and expansion of cultivated land was reported; meanwhile, it was the time when the rainfall intensity and the number of rainfall events significantly increased compared to the 1970s and 1980s, as shown in Table 4. Hence, the increased mean annual streamflow could be ascribed to the combined effects of LULC and climate change. In the case of (Q_s∕Q_t), the increasing pattern could be ascribed to increasing rainfall intensities, the expansion of cultivated land, and the diminution of forest coverage, which might adversely affect soil and water storage and decrease rainfall infiltration, thereby increasing water yield or streamflow. In contrast, the decreasing Q_b∕Q_t is positively related to the increasing evapotranspiration linked to both LULC and climate factors (Table 8). This hypothesis can be explained with the change in CN2 parameter values obtained during calibration of the four SWAT model runs. The CN2 parameter value, which is a function of evapotranspiration derived from LULC, soil type, and slope, increased in the 1980s and 1990s, relative to the 1970s, and could be associated with the expansion of cultivated land and shrinkage of forest land. The increasing CN2 results reflect more surface runoff and less base flow being generated.

Table 7SWAT sensitive model parameters and their (final) calibrated values for the four model runs.

R value from the SWAT database is multiplied by a given value; V replace the initial parameter by the given value; a adding the given value to initial parameter value.

Download Print Version | Download XLSX

Table 8Mean annual water-balance-component analysis in the upper Blue Nile River basin, considering LULC and climate change over their respective periods. All streamflow estimates are for the El Diem station, where water yield (Q_t) = $Q_{\text{s}}+Q_{\text{l}}+Q_{\text{b}}$ and change in soil storage = $P-Q_{\text{s}}-Q_{\text{l}}-Q_{\text{b}}-E_{\text{a}}$-Revap-Recharge.

Download Print Version | Download XLSX

Another important factor contributing to decreasing surface runoff and an increasing base-flow ratio from the 1990s to the 2000s could be the establishment of soil-and-water conservation (SWC) measures. According to Haregeweyn et al. (2015), various nationwide SWC initiatives, such as Food for Work (FFW), the Managing Environmental Resources to Enable Transition (MERET) programme for more sustainable livelihoods, the Productive Safety Net Programme (PSNP), community mobilisation through free-labour days, and the National Sustainable Land Management Project (SLMP), have been undertaken since the 1980s. Haregeweyn et al. (2015) evaluated these initiatives' effectiveness and concluded that community labour mobilisation seems to be the best approach. This can reduce mean seasonal surface runoff by 40 %, with broad spatial variability ranging from 4 % in Andit Tid (north-western Ethiopia) to 62 % in Gununo (southern Ethiopia).

Figure 7Ratio of water-balance component analysis at the El Diem station using an isolated effect (LULC/climate change).

Download

5.5 Effects of an isolated LULC change on streamflow and water-balance components

Yan et al. (2013) used a fixing–changing method, which was also applied in this study to identify the hydrological impacts of LULC change alone. The calibrated and validated SWAT model and its parameter settings in the baseline period was forced by weather data from the baseline 1973–1980 period, while changing only the LULC maps from 1985, 1995, and 2010, keeping the DEM and soil data constant, as suggested by Hassaballah et al. (2017). The results from Fig. 7 indicated that Q_s∕Q_t ratio changed from 40.7 % to 41.2 %, 41.9 %, and 40.9 % by using the LULC maps from 1973, 1985, 1995, and 2010, respectively, whereas the Q_b∕Q_t ratio changed from 17.1 % to 16.8 %, 16.5 %, and 16.9 %, respectively. The largest Q_s∕Q_t ratio (41.9 %) and the smallest Q_b∕Q_t ratio (16.5 %) were recorded with the 1995 LULC map. This could be attributed to the 5.1 % reduction in forest coverage and the 4.6 % increase in cultivated land with the 1995 LULC map, relative to the 1973 LULC map.

On a basin scale over a decadal time period, water gains are mainly from precipitation. The losses are mainly due to runoff and evapotranspiration (Oki et al., 2006), as the losses due to the deep percolation over the whole UBNRB are negligible (Steenhuis et al., 2009). The long-term mean annual deep percolation simulated in this study is about 16.7 mm and is constant in four decadal periods, which is about 6 % of the total water yield. With the fixing–changing approach, the change in streamflow attributable to LULC change was essentially the change in evapotranspiration between the two periods, as the amount of precipitation was constant (1970s), and the change in water storage during the two periods was similar (Yan et al., 2013). Annual E_a losses from seasonal crops are smaller than those from forests, because seasonal crops transpire during a relatively shorter time interval than perennial trees do (Yan et al., 2013). As a result, the actual mean annual E_a simulated by the SWAT model was 871.6 mm at the baseline. It decreased to 871.4 and 871 mm in 1985 and 1995, respectively, and increased to 872.1 mm in 2010. This could be due to the simultaneous expansion of cultivated land and the shrinkage in forest coverage in the 1985 and 1995 LULC maps, relative to the 1973 baseline. Furthermore, this deforestation may reduce the canopy's interception of the rainfall, decrease soil infiltration by increasing raindrop impacts, and reducing plant transpiration, which can significantly increase surface runoff, reducing base flow (Huang et al., 2013). Here, the evapotranspiration change caused by the LULC change is minimal. As a result, the change for surface runoff and base flow is not significant.

5.6 Effects of isolated climate change on streamflow and water-balance components

The impacts of climate change are analysed by running the four models using a unique LULC map from 1973 with its model parameters, while changing only the weather data sets from the 1970s, 1980s, 1990s, and 2000s. The simulated water-balance components shown in Fig. 7 indicate that the Q_s∕Q_t ratio increased from 40.7 % to 45.2 %, 45.6 %, and 46.2 % during the 1970s, 1980s, 1990s, and 2000s, respectively, while the Q_b∕Q_t ratio changed from 17.1 % to 13.5 %, 14.9 %, and 12.7 % for the same simulation periods. The decreasing Q_b∕Q_t ratio for the altered periods, compared to the baseline period, could be attributed to evapotranspiration changing from 872 to 854 mm, 906 mm, and 884 mm in 1970s, 1980s, 1990s, and 2000s, respectively, which can be linked to temperature and amount of rainfall. However, it is important to know the dominant rainfall–runoff process in the study area to fully understand the effect of climate change on the water-balance components.

Although no detailed research has been conducted on the upper Blue Nile River basin to investigate the runoff-generation processes, Liu et al. (2008) investigated the rainfall–runoff processes at three small watersheds located inside and around the upper Blue Nile River basin, namely Mayber, Andit Tid, and Anjeni. Their analysis showed that, unlike in temperate watersheds, in monsoonal climates, a given rainfall volume at the onset of the monsoon produces a different runoff volume than the same rainfall at the end of the monsoon. Liu et al. (2008) and Steenhuis et al. (2009) showed that the ratio of discharge to precipitation minus evapotranspiration, Q∕(P–ET), increases with cumulative precipitation from the onset of a monsoon. This suggests that saturation excess processes play an important role in watershed response.

Furthermore, the infiltration rates that Engda (2009) measured in 2008 were compared with rainfall intensities in the Maybar and Andit Tid watersheds located inside and around the UBNRB. In the Andit Tid watershed, which has an area of less than 500 ha, the measured infiltration rates at 10 locations were compared with rainfall intensities considered from the 1986–2004 period. The analysis showed that only 7.8 % of rainfall intensities were found to be higher than the lowest soil infiltration rate of 25 mm h⁻¹. Derib (2005) performed a similar analysis in the Maybar watershed (with a catchment area of 113 ha). The infiltration rates measured from 16 measurements ranged from 19 to 600 mm h⁻¹, with a 240 mm h⁻¹ average and 180 mm h⁻¹ median, whereas the average daily rainfall intensity from 1996 to 2004 was 8.5 mm h⁻¹. Hence, he suggested from these infiltration measurements that infiltration excess runoff is not a common feature in these watersheds.

From the above discussion points, it is to be noted that surface runoff could increase with increasing total rainfall amount, regardless of rainfall intensity. However, the mean annual rainfall amount in this study decreased from the 1970s to the 1980s (1428 and 1397 mm, respectively), while the (Q_s∕Q_t) ratio increased from 40.7 % to 45.2 %. Similarly, the mean annual rainfall amount in the 1990s (1522 mm) was greater than the mean annual rainfall amount in the 2000s (1462 mm) while the (Q_s∕Q_t) increased from 45.6 % to 46.2 %. In contrast, in climate indices such as 99 % rainfall and SDII (ratio of total precipitation amount to number of days when rainfall >1 mm; R1 mm), the number of days when rainfall >20 mm (R20 mm) increases consistently from 1970 to the 2000s, as shown in Table 4. This indicates that the increasing surface runoff might be due to an increasing of number of extreme rainfall events and rainfall intensity. Although we did not use hourly rainfall data for the SWAT model, this study suggested that the excess infiltration of overland flow dominates the rainfall–runoff processes in the UBNRB, not the saturation excess of overland flow. The contradiction from the previous studies might be due either to the limitation of the SWAT CN method when applied in monsoonal climates or the overlooked tillage activities, which significantly impact the soil infiltration rate. Extensive tillage activities are carried out across the basin at the beginning of the rainy season. Soils get disturbed as a result, which can increase the infiltration rate and ultimately decrease the amount of rainfall converted to runoff.

Although the CN method is easy to use and provides acceptable results for discharge at the watershed outlet in many cases, researchers have concerns about its use in watershed models (Steenhuis et al., 1995; White et al., 2011). The SWAT CN model relies on a statistical relationship between soil-moisture condition and CN values obtained from plot data in the United States, with a temperate climate that was never tested in a monsoonal climate, exhibiting two extreme soil-moisture conditions. In monsoonal climates, long periods of rain can lead to prolonged soil saturation, whereas during the dry period, the soil dries out completely, which may not happen in temperate climates (Steenhuis et al., 2009). Hence, further research that considers biophysical activities such as tillage and the seasonal effects on soil moisture at representative watersheds of the basin is necessary to properly assess the rainfall–runoff processes.