1. Introduction

Flood is excessive water availability, which is caused by above normal stream flow, leading to inundation of areas that are normally not covered by water. In flood generation, antecedent conditions which refer to the saturation of natural storage in the catchment are critical factors [1] [2]. It is the consequence of earlier precipitation or snowmelt in the catchment. If the saturation is at maximum capacity, a consecutive moderate amount of rain can also generate large floods [3]. Further, low permeability of the surface due to dry and crusted soil after a prolonged period without rain can also rapidly convert heavy rainfall to a runoff which usually results in a flash flood [4]. Intense and/or long-lasting rainfall is the main cause of large flood events in tropical regions [5] [6]. In these regions, more water vapor and heat in the atmosphere brings storms and consequently floods will be become more intense [7].

Ethiopia has experienced two major types of floods: flash and river (fluvial) floods. The occurrences and extents of river floods for the last six decades are increased from decade to decade [8] [9]. Spatial distribution of major flood events over major river basins shows: Awash Basin recorded 13 flood events and followed by Wabi Shebele, Rift Valley Lakes, and Genale Dawa with 12, 11, and 11 flood events, respectively [9]. In the eastern and southeastern part of Ethiopia, several areas have been afflicted by floods for decades, killing hundreds and displacing thousands of people. Akola et al., [8] reported that, in Diredawa city (Ethiopia), the frequency of floods has been increasing from 1945, 1977, 1981, 1997, 2001, 2004, 2005 to 2006. Similarly, the recent history of the Wabi Shebele River Basin is marked by frequent destructive floods [10] [11] [12] [13] (e.g., floods in 1996, 1999, 2000, 2003, 2005, 2006, 2008, 2010 and 2016).

Studying the changes in river discharge and precipitation patterns has critical importance as a climatic indicator for environmental risk problems such as droughts and floods. Among these analyzing a river discharge gives the entire picture of a catchment response for water resource planning and management [13] [14]. The seasonality of streamflow and cost of extreme weather events has been found a rapid upward trend in recent decades throughout the world including Ethiopia [11] [15] [16] [17] [18] [19]. This makes the search for trends in flood data series area of scientific interest. In addition, the study of trends in flood time series requires decision-makers to better understand the ongoing changes in hydrologic extremes to make preparations for the possibility of changing conditions. Flooding is the result of hydraulic basin condition, environmental susceptibility of areas to flooding, anthropic impacts and precipitation. A possible correlation with climate change needs to be considered a given scenario, which may be changed in rainfall regimes (probably the most common).

This paper aims to study possible flood changes through a flood simulation approach in watersheds, in conditions of data scarcity and identify correlations with potential driving forces. In Wabi Shebele River Basin obtaining historical river, records are hardly difficult. Some gauging stations that exist in the upper part of the basin have inconsistent data which needs extra data extension mechanisms. Thus, an alternative approach (hydrological modeling) to generate flow from weather and catchment variables is required in such areas.

2. Study Area and Data

2.1. Study Area

The Wabi Shebele Basin is a transboundary basin located at Horn of Africa, situated in between Ethiopia and Republic of Somalia. It originates from Bale highlands ranges of the Galama to Ahmar of Ethiopia, about 4000 m altitude and drains portion of Somalia before draining to Indi-an Ocean. More than 70% of the catchment (202,220 km²) is lying in Ethiopia. The Wabi Shebele Basin in this study is used to represent the catchment that is lying in Ethiopia within 4˚45'N to 9˚45'N latitude and 38˚45'E to 45˚45'E longitude [20]. The air temperature of the Wabi Shebele Basin varies with altitude [10]. The mean monthly temperature of the basin is 19.9˚C. The climate of the basin is dependent on the altitude and strong latitudinal movement of the intertropical convergence zone (ITCZ) [20]. The highlands are cool and densely populated while the lowlands are arid and sparsely populated with recorded annual rainfall of 1213 mm and 268 mm respectively [10] [20].

The basin is divided into four geographical areas based on morphometric characteristics and rainfall regimes [10] [21]: upper catchments which characterized by a mountainous area with abrupt valleys; middle catchments which is wider highland and rainy area; eastern catchments which characterized semi-arid areas and lower catchments which covers arid lowland area of the basin. The watershed characteristics analysis using Arc GIS indicate that flood characteristics of Wabi Shebele Basin is related to basin and relief morphometric characteristics. The mean peak flow (QMPF) in Wabi Shebele Basin has large positive association with linear morphometric parameters (like valley length, mean stream length) and with basin morphometric parameter (like drainage size, shape factor) and negative associations with relief morphometric characteristics (like basin elevation and valley slope). The basin has three climatological rainy seasons: spring (February-May), summer (June-September) and winter (October-January) [10] [18] [22]. While having the largest area coverage, the basin’s annual runoff is estimated to 3.4 BCM (Billion Meter Cube) (Figure 1).

2.2. Data

We utilize both ground-based station observations and gridded analysis data in this paper. In Hydrologic model, digital elevation model (DEM), physiographic data of pedology, land use and occupation and classes of slopes and meteorological data were used to generate flow. Digital elevation model (DEM) with a spatial resolution of 90 meters obtained from SRTM GDEM official website.

The soil information was acquired from FAO Digital Soil Map of the World (DSMW) at a scale of 1:500,000 downloaded from FAO Soil Map and Database

website [23]. The soil data with the resolution of 1 km mainly include soil texture, soil depth, and soil drainage attributes needed for the SWAT model will be derived from Harmonized World Soil Databasev1.2, a database that combines existing regional and national soil information in combination with information provided by FAO-UNESCO soil map. Land-use and land cover data were obtained from Ministry of Water, Irrigation and Energy.

From the data described above, the hydrologic response units (HRU) were established. After HRU definition, the data from climatic stations located in the study basin were inserted on the SWAT model. These data refer to rainfall (mm), maximum and minimum air temperatures (˚C), relative humidity (%), solar radiation (KJ∙m²) and wind speed (m/s). These data were obtained from the National Meteorology Agency (NMA). The data sets provide daily observations for stations exist in the basin. Table 1 shows fourteen weather stations selected in the study, which has good quality and spatial distribution in the basin with a minimum record length of 10 years in between 1980 to 2010. Accordingly, five station from upper catchments; six stations from middle catchments; two stations from Eastern catchments and one station from lower catchments.

Table 2 shows measured discharge data used for model calibration and validation, collected from hydrology department of the Ministry of Water and Energy, Ethiopia (MoWRIE). According to the Ministry of Water and Energy, the measurements of river levels follow the guidelines of the World Metrological Organization (WMO) [10].

*UTM Adindan 1984, Zone 37N, NW = North western, N = North, M = Middle, E = East, S = South.

*UTM Adindan 1984, Zone 37N.

3. Methodology

The methods implemented in this paper comprises a semi distributed macroscale hydrological modeling for the simulation of daily runoff, implementation of Exploratory data analysis (EDA) to explore simulated data, data distributions and examine clusters in the data or relationships between variables and/or sample locations and a nonparametric trend test to detect trends in annual and seasonal maximum runoff.

3.1. Hydrological Model

Soil and Water Assessment Tool (SWAT) is a continuous time, spatially distributed model designed to simulate water, sediment, nutrient and pesticide transport at a catchment scale on a daily time step. In this study, the model was used to generate flows. The model is driven by metrological data like precipitation, temperature, relative humidity, solar radiation and wind speed and physiographic data of pedology, land use and occupation, and classes of slopes.

It uses hydrologic response units (HRUs) that consist of specific land use, soil and slope characteristics. The HRUs are used to describe the spatial heterogeneity in terms of land cover, soil type and slope class within a watershed [24]. The model estimates relevant hydrologic components such as evapotranspiration, surface runoff and peak rate of runoff, groundwater flow and sediment yield for each HRUs unit.

The water in each HRU in SWAT is stored in four storage volumes: snow, soil profile (0 - 2 m), shallow aquifer (typically 2 - 20 m), and deep aquifer. Surface runoff from daily rainfall is estimated using a modified SCS curve number method, which estimates the amount of runoff based on local land use, soil type, and antecedent moisture condition. Peak runoff predictions are based on a modification of the Rational Formula [25]. The watershed concentration time is estimated using Manning’s formula, considering both overland and channel flow. The SCS curve number is described by Equation (2).

$Q_{Surf}=\frac{{\left(P_i-I_a\right)}^2}{P_i-I_a+S}=\frac{{\left(P_i-0.2S\right)}^2}{P_i+0.8S}$ (2)

In which, Q_surf is the accumulated runoff or rainfall excess (mm/day), P_i is the rainfall depth for the day (mm), I_a is the initial abstraction lost from canopy interception, surface storage, and infiltration prior to runoff (mm H₂O; commonly approximated as 0.2S), S is the retention parameter (mm). The retention parameter is defined by Equation (3).

$S=25.4\left(\frac{1000}{\text{CN}}-1\right)$ (3)

The SCS curve number is a function of the soil’s permeability, land used and antecedent soil water conditions. Specifically, the CN values are based on the hydrologic soil group of the area, land use, management, and initial hydrologic condition; with the hydrologic soil group and land use being the most important variables.

For climate, SWAT uses the data from the station nearest to the centroid of each sub basin. Calculated flow, sediment yield, and nutrient loading obtained for each sub basin are then routed through the river system. Channel routing is simulated using the variable storage or Muskingum method. The soil percolation component of SWAT uses a water storage capacity technique to predict flow through each soil layer in the root zone. Downward flow occurs when field capacity of a soil layer is exceeded and the layer below is not saturated. Percolation from the bottom of the soil profile recharges the shallow aquifer. If the soil temperature in a particular layer reaches less than or equal 0˚C, no percolation is allowed from that layer. Groundwater flow contribution to total stream flow is simulated by routing a shallow aquifer storage component to the stream [26] [27]. The model computes evaporation from soils and plants separately. Potential evapotranspiration can be modelled with the Penman-Monteith [28], Priestley-Taylor [29], or Hargreaves methods [30], depending on data availability. In this study, the Penman-Monteith method was used to determine potential evapotranspiration.

3.2. Defining an Extreme Event

In this paper six extreme hydrologic indices: Annual maximum discharge (AMAX), Peak over threshold (3rd quartile) frequency (POTF), Peak over threshold (3rd quartile) magnitude, seasonal maximum discharge for winter (SMW), Seasonal maximum discharge for spring (SMSp) and Seasonal maximum discharge for summer (SMSu) are used to define extreme high discharges.

In extreme value analysis ensuring independence of samples is initial task. In this study the time interval approach is used to ensure the independence of flow discharges. Time intervals between 5 to 14 days between successive peaks; 5 days for catchments < 10,000 km² and 14days for catchments ≥ 10,000, is used in this study. This approach is reported, a strong flood-frequency estimations approach e.g., [31] [32].

3.3. Flood Change Detection and Attribution

3.3.1. Flood Change Detection

Two distribution-free (nonparametric, e.g., rank-based) tests and exploratory data analysis (EDA) are used to detect changes in flood discharge. Exploratory data analysis (EDA) is used as preliminary analysis to explore initial hypotheses on changes in data time series to be confirmed by statistical analysis, nonparametric Mann-Kendall trend test to detect trends in flood discharges, and Quantile perturbation method (QPM) approach to see temporal variabilities in extreme discharges. In practice there is a continuum between “trend” and “change”.

Exploratory data analysis (EDA): involves mainly plotting graphs, and allowed to explore some features in data and assess the first hypotheses to be confirmed by the statistical analysis. In addition, the linear regression gradient plot in the EDA allowed testing of potential trends data time series.

Mann-Kendall (MK) test: The MK test statistic, S is defined as [33] [34] and the test is conducted by applying Equation (4) to Equation (9):

$S={\displaystyle \underset{i=1}{\overset{n-1}{\sum }}{\displaystyle \underset{j=i+1}{\overset{n}{\sum }}\mathrm{sgn}\left(X_j-X_i\right)}}$ (4)

where X_j are the sequential data values, n is the length of the data set, and:

$\mathrm{sgn}\left(\theta \right)=\{\begin{array}{l}1\text{for}\theta >0\\ 0\text{for}\theta =0\\ -1\text{for}\theta <0\end{array}$ (5)

When n ≥ 8, the statistic S is approximately normally distributed with mean and variance given by Mann [34] and Kendall [35]:

$E\left[S\right]=0$ (6)

$V\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)-{\displaystyle {\sum }_{i=1}^n{t}_1\left(l-1\right)\left(2l+5\right)}}{18}$ (7)

where t_l is the number of ties of extent l. the standardized test statistic Z is computed by:

$Z=\{\begin{array}{l}\frac{S-1}{\sqrt{Var\left(S\right)}}\text{for}S>0\\ 0\text{for}S=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}}\text{for}S<0\end{array}$ (8)

The standardized MK statistic Z follows the standard normal distribution with a mean of zero and variance of one under the null hypothesis of no trend. A positive Z value indicates an upward trend, while a negative one indicates a downward trend.

The Kendall rank coefficient is often used as a test statistic to establish whether two variables may be regarded as statistically dependent. Under the null-hypothesis of independence of X_i and X_j, the sampling distribution of “tau” has an expected value of zero. The value of “tau” ranges from −1 (100% negative association, or perfect inversion) to +1 (100% positive association or perfect agreement). A value of zero indicates the absence of association.

The p-value of the MK statistic S of sample data can be estimated using the normal cumulative distribution function:

$p=0.5-\varnothing \left(\left|Z\right|\right)\text{where}\varnothing \left(\left|Z\right|\right)=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_0^{\left|Z\right|}{\text{e}}^{\frac{r^2}{2}}\text{d}t}$ (9)

If the p-value is small enough, the trend is quite unlikely to be caused by random sampling. At the significance level of 0.05, if p ≤ 0.050, then the existing trend is assessed to be statistically significant.

Quantile Perturbation method (QPM): is a statistical analysis specially designed for analyzing extreme conditions, to study trends and multi time period oscillation patterns in hydro-climatic extremes [14] [35] [36] [37]. The method has two concepts: 1) the frequency aspect which focuses on how often an extreme event (quantile) may occur and, 2) the perturbation aspect which determines the changes in the extremes for a particular return period. The method compares the long-term baseline period extreme value quantiles with that of a selected sub-period quantile. The selected sub-period (block period) is a subseries taken from long-term baseline time series representing the period of interest. To select appropriate value of block length in between 5- and 15-year intervals, Tabari et al., [37] recommends applying QPM to extreme time series to different block of years and select the one which shows a high variability at a given time interval. To check the significance of perturbation factor in extreme quantiles 95% CI (Confidence Interval) is computed using non-parametric bootstrapping method is performed. The perturbation values fall outside of confidence intervals are identified as statistically significant perturbation value.

3.3.2. Flood Change Attribution

Studying trend and variability in hydroclimatic time series alone has no value unless the cause or driving factors of changes are followed. Merz et al., [38] reported in their research that the main goal of trend and variability studies in hydroclimatic time series is to test hypothesis about the influence of driving forces of floods. The detected trends were explained by the correlation with changes in another variables; climate variables, environmental background condition variables and external factors. In addition, changes in catchment and environmental background condition factors, deforestation and reforestation and alteration in agricultural management practices are analyzed.

Six watersheds representing mountainous and plain regions of Wabi Shebele River Basin were selected to identify most powerful variables in estimating the mean peak flood (QMPF) of the basin from climate, watershed’s physical characteristics and external factors. The mean peak flow (QMPF) was expressed in this study as the arithmetic mean value of peak over threshold (3^rd quartile) flows for the period of record. The environmental background factors are extracted from delineated watershed using Arc Hydro interface in Arc GIS. Population density and land use and land cover data were collected from MoWIE. Daily precipitation data from the meteorological stations in Wabi Shebele Basin and the surrounding areas were obtained from the National Meteorology Service Agency (NMA). Detail of data used and sources are described in Table 3. To highlight the most powerful parameters in flood discharge characteristics PCA is used.

Principial Component Analysis (PCA): is used in this study to see multivariate relationships between potential driving factors and mean peak flow discharge (QMPF). If the number of predictor variables increases and they are highly correlated, Multiple linear regression (MLR) models gets more unstable. PCA is one of the multivariate statistical techniques that can be used to deal with highly correlated variables in regression [39] [40]. In PCA, different types of variables: hydrologic and watershed variables are treated together. The original dataset of n variables, which are correlated to various degrees are transformed to n numbers of uncorrelated PCs. The PCs are linear transformation of the original variables in such a way that the original and the new variables have equal sums of the variances. Although the number of PCs and original variables are equal, the first few PCs explain the majority of the variance in the data set, reducing the dimensionality of the original data set [39].

The PCs are sequenced from the highest to the lowest variance, i.e., the first PC describes the data’s highest variance proportion. The next highest variance is explained by the second PC and so on. The values of PCs can be obtained from Equations (10) and (11):

$\text{PC}1=a_{11}{x}_1+a_{12}{x}_1+\cdots +a_{1n}{x}_n={\displaystyle {\sum }_{J=1}^n{a}_{1j}{x}_j}$ (10)

$\text{PC}2=a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n={\displaystyle {\sum }_{J=1}^n{a}_{2j}{x}_j}$ (11)

where $x_1,x_2,\cdots ,x_n$ are the original variables and a_jj are the eigenvectors. The eigenvalues are the variances of the PCs. The covariance or correlation matrix of the data set is used to derive the coefficients a_jj, which are the eigenvectors. The eigenvalues of the data matrix can be calculated by Equation (12):

$\left|C-\lambda I\right|=0$ (12)

where C is the correlation/covariance matrix, λ is the eigenvalue, and I is the identity matrix. The PC coefficients or the weights of the variables in the PCs are then calculated by Equation (13):

$\left|C-\lambda I\right|a_{jj}=0$ (13)

3.4. Flood Frequency Analysis

Parallel to changes in magnitudes, changes in the frequency of extreme events, which can be described by the occurrence rate has significant role in flood events. Changes in the occurrence rate reflect the clustering properties of flood events, caused by natural variability, regime shifts in the atmosphere [40] [41] or land use changes. For each sample watersheds, peak over threshold (3^rd quartile) time series was extracted to see trends in magnitude and number of occurrences in annual and seasonal time level.

4. Results and Discussion

4.1. Calibration and Validation of the Hydrological Model

For model calibration and validation, the observed daily and monthly streamflow data were used from 1988 to 2000 with three years warming period. To evaluate the model performance, three parameters have been used, namely R² and NSE and P-bias. NSE is a normalized statistic, ranges from −∞ to 1, used to indicate the relative value of residual variance compared to the variance of the observed data and values close to one shows a perfect match of the modeled with the observed data [42], Equation (10). R² is the proportion of the total variance in the observed data that can be explained by the model, Equation (11).

$\text{NSE}=1-\frac{{\displaystyle {\sum }_{i=1}^N{\left(X_i-Y_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^N{\left(X_i-\bar{X}\right)}^2}}$ (14)

${\text{R}}^2={\left[\frac{{\displaystyle {\sum }_{i=1}^N\left(X_i-\bar{X}\right)\left(Y_i-\bar{Y}\right)}}{\sqrt{{\displaystyle {\sum }_{i=1}^N{\left(X_i-\bar{X}\right)}^2}}\sqrt{{\displaystyle {\sum }_{i=1}^N{\left(Y_i-\bar{Y}\right)}^2}}}\right]}^2$ (15)

$\text{Pbias}=\frac{{\displaystyle {\sum }_{i=1}^N\left(Y_i-X_i\right)}}{{\displaystyle {\sum }_{i=1}^N{X}_i}}\ast 100$ (16)

where x and y are observed and modelled streamflow, respectively, N is the number of data pairs. Table 4 shows statistical evaluation of model performance. Based on the model the physiographic characteristics of watersheds are presented in Table 5.

4.2. Annual Maximum Discharge

The presence of seasonal cycles and clusters in annual maximum discharges (AMAX) are examined using exploratory data analysis (EDA) (Figure 2). The magnitude of annual maximum discharge oscillates at 5 - 10-year intervals in

Relief is calculated as the difference between the highest and the lowest elevation, Aspect is the averaged aspect of the basin (N North; NE North-East; S South; SW South-West).

most of sample stations. Annual maximum discharges indicate less than mean annual maximum values in 1990s for the middle and eastern catchments of Wabi Shebele Basin. However, annual maximum discharge in 2000s is observed above mean maximum annual discharge in all gauging stations entire the basin.

The MK test applied at each site for the period 1981-2010 showed less decreasing trend in upper and lower catchments. For the longest period, 55% of the stations indicate weak to significant decreasing trends in annual maximum discharge. However, the watersheds in middle basin indicate a majority of significant increasing trend annual maximum discharge, p < 0.05. Too see multi-temporal changes in annual maximum discharge Mann-Kendall trend tests are analyzed at 5, 10, 15, 25 and 30 years intervals.

From both Figure 2 & Figure 3, some conclusions are extracted on multi-temporal trend analysis. First, blue colors are more frequent than red colors,

meaning that positive trends are more frequent than negative, especially for the most recent period since 1996 in all stations (Figure 3). Second, negative trends (red colors) appear more significantly in eastern Fafen catchments and lower stations of Wabi Shebele River Basin before 2000 (as shown on the Fafen @ Jijiga, Fafen @ Kebridehar, Wabi @ Gode and Wabi @ Burkur). Most of gauging stations in upper and middle catchments indicates weak to significant increasing trends and while stations in lower and eastern catchments showed weak to strong decreasing tendencies in annual maximum discharge during the past 30 years in Wabi Sheble River Basin. The years, 1980s and 2000s are the decades were weak to significant increasing flood discharges are occurred. For 1990s, most stations indicate decreasing trends in annual maximum discharge in the study area.

4.3. Flood Frequency and Seasonality

In Wabi Shebele Basin there are nine mean annual flood events, discharges greater than third quartile (3^rd quartile). Some years (i.e., 1981, 1982, 1983, 1990, 1991, 2006, 2007 and 2010) were particularly rich in flood events, with more than 14 events on average. In contrast, in some other years (i.e., 1997, 2001, 2002 and 2004) only two to four events (POTF), was observed in the basin. The decades 1981-1990 and 2001-2010 are noticed as the decades of the richest flood events.

The averaged POTF in upper and middle catchments reveals increasing trends of flood events for longest time period, 1981-2010 (Figure 4). However average POTF in lower stations of Wabi Shebele Basin and Fafen catchments indicates decreasing trends in number of events. At sample gauging stations the result of MK trend test for the long-term period (1981-2010) reveals positive trends of POTF for 55% of gauging stations. Among these significant increasing trends in

POTF is observed over middle catchments (Erer watershed at Hamaro, and Golocha) and the rest 45% of stations indicates negative trends, mainly present in Fafen watershed and Wabi Shebele River at Gode and Burkur.

From Figure 5, the QPM analysis using peak over threshold (3rd quartile) indicates that most of extreme flows varies within a confidence with high oscillation patterns in the entire the basin. Extreme flows vary above reference line (above mean) in upper, middle and lower valley of Wabi Shebele River stations, whereas it varies below reference line in eastern catchments.

Figure 6, illustrates the importance of time windows in analysis of seasonal flood trends. The EDA of spring and winter discharge does not present trends in the longest period (1981-2010) in most catchments of the basin. Almost in all sample catchments, seasonal maximum discharge indicates oscillation pattern at a decade interval. Like annual maximum discharge, maximum discharge in all seasons also indicates less than mean seasonal maximum discharge in 1990s.

From Figure 7, the multi-temporal analysis in seasonal flood discharge indicates similar patterns in summer and winter and different in spring are observed in throughout the basin. In Fafen catchments spring flood discharge indicates weak increasing trend while indicates decreasing trend in summer and winter season for the last 30 years. Darker colors show statistically significant trends, and they are more frequent in summer and less in winter in eastern catchments. For 1990s, upper and lower Wabi Shebele Basin flood discharge indicates significant decreasing trends in all seasons and weak decreasing trend in eastern catchments. In all catchments flood discharge indicates increasing tendency in 2000s. This pattern is stronger in summer and winter when 72% of stations showed significant increasing trends for the period 2001-2010. The same result was previously observed in the Wabi Shebele River Basin, where a significant increase in spring, summer and winter floods was identified [22]. Table 6, illustrates seasonal significant anomalies in extreme discharges investigated using QPM.

4.4. Summary of Flood Changes

The presented results revealed that there has been an increasing trend in the flood magnitude and frequency over early 21^st century entire the basin. Similar result is reported in literatures that the magnitude and frequency of floods indicates increasing trend in Wabi Shebele River Basin since 2000 [10] [11] [22]. The positive Kendall’s Z values indicates increasing trend within analysis period (Table 7). For the period 1981 to 2010 the annual maximum flood discharge shows upward trends in upper and middle catchments while downward trends in eastern and lower catchments of Wabi Shebele Basin. The annual maximum stream flow for middle catchments (i.e., Erer at Hamaro and Gololcha at Wabi junction) shows a positive significant trend because the computed p-value is in both watersheds are lower than the significance level (α = 0.05) in the region. However, significant decreasing trends in annual maxima are observed in Fafen watersheds at Jijiga and Kebridehar gauging stations. Seasonal trend analysis reveals similar trends and patterns with annual maximum stream flow almost in all stations during past 30 years.

Extreme discharge variability analysis using peak over threshold (3^rd quartile) based on QPM showed significant increasing trend in early 1990s & 2000s and

decreasing trends in 1980s particularly in upper and middle catchments (Table 8). Over eastern catchments 1980s is the decade in which significant increasing trends observed and decreasing trends in 1990s and 2000s. The lower Wabi Shebele river stations indicate general decreasing trend in analyses period, 1980-2010.

In seasonal extreme variability analysis, significant anomaly occurrence season varies with catchments as presented in Table 6. In upper and middle catchments (Wabi at Dodola, Maribo, robe Wabi at Legehida and Gololcha watersheds), spring season is the season in which highest extreme variabilities are occurred. Similarly, in eastern catchments (Erer watersheds) highest extreme discharge anomalies are occurred in winter season and in lower Wabi Shebele catchments (at Gode and Burkur stations) during summer season. The years; 2^nd half of 1980s, 1^st half of 1990s and 2^nd half of 2000s are the years of positive anomalies (significant increasing trends) in flood discharges. Whereas, the years; 1^st half of 1980s, 2^nd half of 1990s and 1^st half of 2000s are the years of significant negative anomalies occurrence in all season. It is known that, the Wabi Shebele Basin is characterized by two rainfall regimes [43]: the area characterized by a quasi-double maximum rainfalls pattern with a small peak in April and maximum peak in August which covers the west-east highland of the basin (bimodal type I); and the area dominated by double maximum rainfall pattern with peaks during April and October covers the south-eastern low-lying areas of the basin (bimodal type II).

4.5. Potential Drivers of Flood Changes

4.5.1. Principal Component Analysis (PCA)

PCA of flood peak discharge with major driving forces like climate drivers, environmental background conditions and external factors is performed to highlight most important variables in Wabi Shebele Basin. Therefore, PCA is applied to eleven selected predictors including those which have no significant correlation with each other (less collinearity), i.e., Drainage area (DA), Basin elevation (BE), Basin slope (BS), Shape factor (SF), valley slope (VS), clay, sand, loam, Mean Annual Rainfall (MAR), forest, Agricultural land (AGR) and Population Density (PD) to achieve uncorrelated six PCs. The eigenvalues represent the quantity of variability in the data and they are presented in Table 9. Table 9 confirms that the first two PCs explain the maximum degree of variability of the data set with the proportion of 45.6% and 28.7%, respectively. The proportions of other PCs (PC3, PC4, PC5 and PC6) range 0.0% - 17.3%. Both PC1 and PC2 accounts for 74.3% variance, meaning more than 2/3 of variability in dataset is explained in the first two PCs. To explain at least 90% of variation in the data the first three components are used.

The coefficients in Table 9 show the linear combinations variables that make each principal component. Absolute values near zero indicate that a variable

contributes little to the component, whereas larger absolute values indicate variables that contribute more to the component. In the analysis, first principal component has large negative associations with BE, VS, MAR and forest, so this component primarily measures the basin altitude difference and land cover. The second component has large positive associations with BS, SF and loam, so this component primarily measures the slope and shape of the catchment. The third component has large positive association with sand, AGR and PD, so this component primarily measures the basin farm land and population density.

The loading plot in Figure 8, visually shows the results for the first two components. From the graph, DA and sand indicates small angle (<90˚) from QMPF line, meaning the variables positively correlated to QMPF. The variables: forest, PD, BE and VS indicates angles related to 180˚, meaning they are significantly negatively correlated to QMPF. Whereas the variables: BS, SF and loam has no significant correlation with QMPF in Wabi Shebele River Basin.

4.5.2. Climate Drivers

Most of the rivers in Ethiopia exhibit typical characteristics of tropical rainfall-dependent flow regimes. According to previous studies, the spatial and temporal distribution of rainfall governs amount and intra- and inter-annual variability of water availability in Ethiopia [6]. General understanding of hydroclimatic variables such as precipitation, temperature and discharge is important for water resource planning and management. Among this temperature is a factor indirectly influencing streamflow as it is responsible for the evaporation and moisture, while precipitation is the major driving factor changes in streamflow especially in tropical river basins like Wabi Shebele River Basin. Table 10 shows that the correlation matrix in between peak flood discharge and driving forces. PCA conducted in this study reveals that rainfall factor is less related than other driving forces (Table 11). To see relationships of rainfall and floods in the basin,

QMPF = Mean peak flow (m³/s), DA = Drainage area (Km²), BE = Basin elevation (m), BS = Basin slope (m/km), BP = Basin perimeter (km), VL = Valley length (km), SF = Shape factor, DD = Drainage density, VS = Valley slope (m/km), clay = fraction of clay, sand = fraction of sand, loam = fraction of loam, MAR = Mean annual rainfall (mm), AGR = fraction of agricultural land, forest = fraction of forest coverage, PD = Population density (pop/km²).

Pearson’s correlation test is performed in between ground based gauged precipitation data and flood discharges at different catchments of Wabi Shebele River Basin. Maximum flow discharge at annual and seasonal aggregation levels are used to correlate with total rainfall at annual and seasonal levels. The result indicates that flood discharge in the basin is positively correlation values in most watersheds both at annual and seasonal aggregation levels.

4.5.3. Environmental Background Conditions

The environmental background condition factors, i.e., drainage area (DA), elevation, slope, and soil type are analyzed in the study. The elevation varies greatly in Wabi Shebele River Basin, with moderate to low elevations in eastern and lower catchments (i.e., Fafen watershed and lowlands of Wabi Shebele Basin), moderate and high elevations in the upper and middle catchments (i.e., middle, northern and north west of the basin). Elevational differences of watersheds reflect differences in land cover type and associated hydrologic runoff response and geomorphic disturbance [44]. Less land cover results rapid response to rainfall, increased flood magnitudes, increased potential for debris flows, abundant erosion and consequently increase slope failures. Among the environmental background condition factors elevation and drainage size are among the most important factors in formation of high floods in terms of magnitude and frequency in Wabi Shebele River Basin (Table 10).

Watershed slope (WS; %) is the mean watershed slope, measured by calculating the maximum rate of change between each cell. It provides an indication of the steepness of the drainage area. As the slope decreases, catchment soils become more permeable and thus the effect of infiltration becomes more significant [45]. In Wabi Shebele Basin the east and downstream part of the basin is characterized by low slopes, while the west and upstream part of the basin is characterized by high slopes.

In surface runoff generation, soil infiltration rate is another sensitive variable. Course textured soils have large well-connected spaces and allow more water to infiltrate through it quite rapidly while fine grained soils dominated by clay have low infiltration rates due to their smaller sized pore spaces [46]. Soil containing large amount of the sand and silt tend to form crust and become compacted, which significantly reduces the infiltration rate. The amount of organic matter on soil surface can enhance infiltration because organic matter has more porous than mineral soil particles and it can hold much greater quantity of water. In Wabi Shebele Basin soil distributions vary spatially; loamy sand in the east and downstream part of the basin, clay in middle part of the basin, sandy loam in the center of the province, and silty clay in north west (Figure 9). Sand soil is identified as the most powerful variable among soil types in flood formation of Wabi Shebele Basin (Figure 8 and Table 10).

4.5.4. External Factors

The results presented under Sections 4.5.2 and 4.5.3 indicates the relationships

between flood trends and changes in some atmospheric variables (i.e., meteorological) and environmental background conditions, but these attributes explain variability of flood discharge only at partial level. Therefore, non-climatic changes in catchment and river parameters must also be taken into account.

During the last four decades some environmental changes occurred in the studied catchments that influence the conditions of flood runoff. Land use and land cover has been changing in the northern part of Wabi Shebele River Basin (Table 12). In the basin the extent of shrublands indicates significant increasing trends while grassland and cultivated area showed decreasing trend from 1984 to 2004 [47]. Similarly, the extent of riparian woodland in the basin indicates decrement in this period interval. Shrubland class is the areas with extensive physical limitation; like very steep slopes, shallow soils, rock outcrops, series of deeply dissected gorges, dry and rugged areas. Due to human activities and pressures in large semi-arid areas (middle and eastern upper catchments) of the basin, most coverage of Riparian Woodland, Grassland and Perennial and seasonal Swamp and Marshland covers are changed to Shrubland in the basin.

Another potential driver of floods in watersheds is population density. The population growth and cultivated land density has strong correlation [48]. Increment of cultivated land results in accelerated runoff process, especially as a consequence of rapid development of gullies [17]. In Wabi Shebele Basin, human activities are concentrated in the west and eastern upper highland areas of the Basin. The cultivation land density has strong correlation value with population density with correlation value (r) of 0.657, while negatively correlated to flood discharge in catchments with correlation value of −0.395 (Table 10).

Source: [10].

4.6. Discussion

The peak flood discharge shows good agreement with some variables from climate, environmental background conditions and external factors in Wabi Shebele Basin. Climate factors, specifically precipitation indicates positive correlation with flood events in the basin. The correlation analysis indicates that annual maximum flow discharge has positively correlated to total rainfall in 55% of sample watersheds. Seasonal flood discharges also positively correlated to seasonal precipitations in most of sample watersheds (>78%). Therefore, our study identified the precipitation as one of driving forces of flood events in Wabi Shebele Basin, which confirms the findings of Ethiopian Ministry of Water, Irrigation and Energy (MoWIE) [10] during basin master plan study. This study reveals that severe hydrologic extremes in the basin especially in 1973, 1979, 1984-1985 is caused by natural atmospheric variability.

Watershed characteristics like drainage size, elevation, slope and soil types are identified as the most powerful variables in flood events of Wabi Shebele Basin. The mean peak flow (QMPF) is positively correlated with variables: drainage area (DA), basin perimeter (BP), valley length (VL), shape factors (SF), fraction of sand and loamy, where maximum correlation with DA, BP, VL and sand with correlation coefficient of 0.92, 0.93, 0.96 and 0.54 respectively. This means that larger watersheds are expected to have a higher mean peak flow. There are different studies confirms this result that the drainage area is a significant factor that positively affect peak discharge e.g., [45] [49]. Drainage size affects not only the flow collecting ability but also the time to peak discharge [45].

External factors like land use change and population density are also exhibit differences in hydrologic runoff response, which can be directly linked to flood events. Fraction of forest coverage and population density is found negatively correlated to flood discharges in Wabi Shebele Basin. Loss of land cover, thinner forest canopies, grass lands and reduced infiltration of rainfall result in rapid hydrologic response, increased flood magnitudes and frequency [49]. In Wabi Shebele River Basin, the magnitudes and frequency of flood events particularly in middle watersheds indicates increasing trend in recent decades (Sections 4.2 and 4.3). In mountainous zones, increased potential for intense convective storms, increase cultivated land and more highly confined river valleys results more rapid runoff response to precipitation and cause for extreme floods. In the basin, human activities are concentrated in the west and eastern upper highland areas of the Basin [10].

5. Conclusions

In this study variabilities in flood and relationships with major driving factors are investigated. Both exploratory data analysis (EDA) and non-parametric tests (i.e., Mann-Kendal trend test and quantile perturbation (QPM) methods) are used to see temporal variabilities in flood discharge. The risk level used was 5 %. The p-value of the MK statistic S of sample data is used to measure the significance of trend value; if p ≤ 0.050 (significance level), then the existing trend is assessed to be statistically significant. Finally, Principal Correlation Analysis (PCA) is performed in between flood discharge and potential driving factors to identify the most powerful drivers on flood changes. A larger absolute value of the Eigen coefficient indicates variables that contribute more to the component.

The multi-temporal trend analysis at 5, 10, 15, 15, 25 and 30-year intervals starting from 1980 in annual maximum discharge showed increasing trends for most recent periods (since 1996) in all sample stations while decreasing tendency of flood discharges is observed before the 2000s particularly in eastern and lower catchments in the basin. Peak over threshold (3^rd quartile) values of discharge analysis using QPM show that there are significant positive anomalies (outside confidence interval) in the early 1990s & 2000s and negative anomalies in the 1980s, particularly in upper and middle catchments. However, in eastern catchments, the 1980s is the decade in which significant positive anomalies were observed and negative anomalies in the 1990s and 2000s.

The PCA reveals that drainage area (DA), watershed mean elevation (BE), valley slope (VS), fraction of sand coverage (sand), population density (PD), forest fraction (forest) and mean annual rainfall (MAR) were initially found to be the best predictors of flood peak discharge. Among these, drainage area (DA), mean annual rainfall (MAR) and the fraction of forest coverage (forest) were discovered as the principal driving factors for flood peak discharge in Wabi Shebele River Basin. However, fraction of clay soil, fraction of loam soil, basin shape factor (SF), basin slope (BS) and fraction of agricultural land (AGR) were found to be less important in flood peak discharge prediction of Wabi Shebele Basin watersheds.

Generally, the study tried to answer, the trends and variability status of extreme flood discharge and its relations ships with potential driving factors in the Wabi Shebele Basin. The study can provide information on which driving factors should be prioritized in mitigation measures to decrease extreme weather disasters in the basin. Accordingly, basin slope and drainage size from environmental background conditions and forest coverage and population density from external factors are the major driving factors one can improve to minimize the hydrological disasters in the area.

Acknowledgements

We would like to thank the Ethiopian National Meteorology Service Agency and Ministry of Water, Irrigation, and Electricity for providing the necessary data.