Spatiotemporal Trends and Detection of Changes in Hydrological and Climatic Variables of Modjo River Watershed, Ethiopia

Understanding hydro-climatic trends in space and time is crucial for water resource planning and management, agricultural productivity and climate change mitigation of a region. This study examined the spatiotemporal variations in precipitation, reference evapotranspiration (ETo) and streamflow in a tropical watershed located in the central highlands of Ethiopia. Temporal trend implications were analyzed using the Mann-Kendall test, and Theil-Sen approach, whereas the inverse distance weighted interpolation method was applied for spatial trend variability analysis. The result showed that a significant decreasing trends in streamflow for the major rainy (Kiremt: Jun Sept) season and annual time scales. At the same time, the annual and monthly ETo followed significantly increasing trends, but there has been a trendless time series for most of the months and annual mean precipitation series for the period 1986 - 2015. The study indicated that the spatial variability of annual and seasonal precipitation series decreased from north to south and west to east, while this was increased for ETo both for annual and seasonal time series over the study watershed. The contribution of rainfall and mean temperature to streamflow decline was insignificant. It is pointed out that river flow regime is weakly affected by climate changes, hence human activities are stronger in explaining the river flow trends of the watershed. Therefore, urgent calls on the needs for reducing human-induced impacts, and implementing appropriate watershed management, conservation measures and an efficient use of water resources. and Sen’s slope tests to analyze trends in streamflow, rainfall, and temperature time series over the upper Awash sub basin, Ethiopia. The statistical test result revealed significant decreasing trends of annual rainfall and a significant increasing trend in annual temperature in some of the stations. They also indicated that a significant increasing trend in Bega season flow at some of the gauging stations. The other study recently conducted by Orke and Li (2021) investigated hydro climatic variability in the Bilate watershed using Mann –Kendell test, and the Sen’s slope estimator. They detected significant decreasing trends in annual mean rainfall and streamflow, while significant increasing trend in temperature for the period 1986 to 2015. In accordance with these views, therefore, the nonparametric tests are more suitable for presence of non-normally distributed data, outlier, and missing values in time series, which is a challenge and frequently encountered in hydrology (streamflow) and meteorological (rainfall, temperature, etc.) time series data. level for the studied period. The significant upward trend in ‘Tmax – Tmin’ and the numerical increase in Tmax have relevance in the increasing trends for the potential evapotranspiration of the watershed than the other terms in the Hargreaves equation. The findings of this study also suggest that the scarcity of available water is getting worse, along with rise in evapotranspiration. Availability of data and material The data used in the study were obtained from the Ethiopian Meteorological Service Agency (NMSA) and MoWIE. The data that support the findings of this study available upon request

study, June to September (summer rainfall) showed a significant decreasing trend while October to January (autumn rainfall) and February through May (winter rainfall) detected as trendless time series over the basin. Moreover, the very recent paper conducted by Mekonnen et al. (2020) using the Mann-Kendall and Sen's slope tests to analyze trends in streamflow, rainfall, and temperature time series revealed significant decreasing trends of annual rainfall and a significant increasing trend in annual temperature in some of the stations over the upper Awash sub basin. In view of the above research findings, therefore, there have been clear contradictions, and incompatibilities among trend results either seasonal or annual time scales that creates an illusion for planners, policy and decision makers considering the Awash River basin. Therefore, it is the believe of this paper that basin scale trend analysis can have contradictive results due to constitutes of months for seasons, length of records, number of stations considered and approaches followed for quality assessments.
Modjo river watershed falls within the upper reach of the Awash river basin. It is clearly observed that the watershed is affected by natural and human-induced factors with potentially adverse impacts for agricultural productivity, water resources and environment deterioration. Therefore, understanding the spatial and temporal variations of hydro-climate variables and their association is very crucial to the weather and climate forecasts for agricultural activity and water resources planning and development. However, in this watershed, few studies have been conducted in the domain of climate changes which mainly focused on temporal variations of precipitation, temperature and streamflow (Mahtsente et al. 2019;Mulugeta et al. 2019;Eshetu 2020). To our knowledge, no attention has been paid for the spatial and temporal trends in hydro-climate varables previously in the study area. In this study, the spatiotemporal trends and variability in precipitation, runoff and reference evapotranspiration (ETo) were examined in detail at the watershed scale. The specific objectives are (1) to identify the direction, magnitude, and rate of changes in precipitation, streamflow and evapotranspiration for the period 1986 -2015; (2) to indicate spatial variations in precipitation and ETo at seasonal and annual time scales, and (3) to evaluate the linkage of runoff variation with climate variables. The study results would be helpful for climate adaptation measures and to minimize the effects of vulnerable disasters of an agricultural dominated watershed.

Study area
This study is carried out for the Modjo river watershed located in the central Ethiopia ( Fig. 1), within the upper Awash river basin. The geographical location of the watershed is between latitude of 8°35′00″N to 9°05′11″N and longitude of 38°54′35″E to 39°15′30″E). It is important to note out that the watershed has natural unregulated flow regimes that cannot impact on annual and seasonal streamflow discharge, which is the focus of this study.
The region experiences a bi-modal type of rainfall regime, with two maxima mostly in April and July/August months. On the basis of recorded data of the period 1986 -2015, the annual rainfall varies between 806 mm and 1175 mm (https://www.ethiomet.gov.et). Based on recent studies (Eshetu 2020), the annual mean minimum and mean maximum temperatures vary between 10.6 -11.6 o C and 21 -27 o C, respectively, whereas the annual average temperature of the watershed varies between 16.9 o C and 19.3 o C (Gessesse et al. 2014). Agricultural cultivated land use is the dominant land cover type (Gessesse et al. 2014). The main crops in the area includes Teff, Lentils, Haricot beans and Chickpea, with the highest dominancy by Teff. In the area, natural forest cover is characterized by sparse woody vegetation, with Eucalyptus trees close to the towns and around the rural villages. In general, it is observed that the watershed has a highly deforested land covers. According to the soil data obtained from the Harmonized World Soil Database (HWSD) v1.2 (Nachtergaele et al. 2008), the major dominant soil types are Pellic Vertisol and Vertic Cambisols.

Fig. 1 The study area 2.2 Data and quality assessment
The daily precipitation, and temperature data from 1986 to 2015 at seven stations were used in this study (Table 1). The observed climate data (rainfall, min and max temperature) data was collected from the National Meteorological Service Agency (https://www.ethiomet.gov.et) of Ethiopia. As for many developing countries in the world, the density of climate stations here with long-term records is low, and the available gauges do not cover entire drainage area (Fig. 1). Therefore, in this study, the three rain-gauge (Aleltu, Koka dam and Zequala stations) and two temperature stations (Aleltu and Koka dam) were considered from the neighboring watersheds. In the agricultural activities times, three seasons are experienced in the area, dry season (local name is Bega), is defined from Oct to Jan, short rainy season (locally known as Belg) from Feb to May, and long rainy season (locally called as Kiremt), from Jun to Sept (Weldegerima et al. 2015). The daily streamflow data for the period 1986 to 2015 were obtained from the Ministry of Water, Irrigation & Electricity (https://www.mowie.gov.et) of Ethiopia. To facilitate comparison with precipitation depth, the streamflow values were standardized in terms of runoff depth. Moreover, a 12.5 m by 12.5 m spatial resolution digital elevation model (DEM) was downloaded from https://vertex.daac.asf.alaska.edu/. This data was implemented for study area boundary delineating and for geostatistical interpolation of precipitation and evapotranspiration from gauges.
In the present study, 7 rain-gauge stations, 4 temperature stations and one hydrometric station have been used for analysis. In order to ensure that the data accurately evaluated in analysis, pre-examination in data series should be applied. In this context, the daily data obtained from the stated sources were subjected to different quality control procedures such as checking missing records, break points, inappropriate values, etc. In the first place, the preliminary quality of the daily climate and hydrology data was confirmed by the National Meteorological Service Agency (NMSA), and Ministry of Water, Irrigation and Electricity (MoWIE) of Ethiopia, respectively. In the process, considering data from 1986, we tried to eliminate the majority of missing values in the datasets. In the study, filling of missing records in climate data was done using Multivariate Imputation by Chained Equations (MICE) algorithm using R statistical software (R core Team 2020). The MICE creates multiple predictions for each missing value and considers uncertainty in the calculations and provides standard errors. One important advantage of the MICE algorithm over other approaches is that it involves filling in the missing values multiple times, and creating a complete dataset. Unlike the climate (rainfall and temperature) data records, the observed streamflow record had no gaps that needed to be filled. This is due to the consideration of the time series data from 1986 onwards. Table 1 Lists and locations of climate stations whose data records were analyzed in this study The other pre-processing analysis carried out was detection of outliers. An outlier is extreme value whose presence in time series can affect the reliability of estimation outcomes. In the present study, it is observed that the hydro-climate time series data have suffered from large outliers. For this purpose, the Box-Whisker Plot approach, was used in detecting outliers in time series data. It is a plot of presenting datasets by aid of central and distribution criteria (Sajad et al. 2014). In the plot, the interquartile range (IQR) is the difference between the third and first quartiles (i.e. IQR=Q3-Q1) which represents length of the rectangle. Any data smaller than Q1-1.5IQR and larger than Q1+1.5IQR considered as weak outlier, but data placed out of the external boundary (i.e. lower than Q1-3IQR and greater than Q1+3IQR) considered as strong outlier. Unlike missing data records, however, the streamflow data contain outliers (extreme values), which is objectively identified, removed and replaced by more plausible values to improve the information content of the data for trend analysis.

Estimation of evapotranspiration
The evapotranspiration rate from a reference surface, not short of water, is called the reference crop evapotranspiration or reference evapotranspiration (ETo). According to Allen et al. (1998), ETo expresses the evaporating power of the atmosphere at a specific location and time of the year and does not consider the crop characteristics and soil factors. This is to mean that the only factors affecting ETo are climatic parameters. Among the large number of empirical or semi-empirical equations, the FAO Penman-Monteith method is recommended as the sole method for determining ETo because it is physically based, and explicitly incorporates both physiological and aerodynamic parameters. Moreover, this method requires the availability of air temperature, radiation, air humidity and wind speed data. In our study area, due to the unavailability of all these data, the temperature based, Hargreaves method (Hargreaves and Samni (1985) was used to estimate the daily reference evapotranspiration (ETo). In the studies by Seong et al. (2018) Mubialiwo et al. (2020), the Hargreaves method is selected because it yielded comparable results with the FAO Penman-Monteith method. Overall, the equation used for ETo estimation: Where, ETo is reference evapotranspiration (mm/day), Tmin and Tmax are the minimum and maximum air temperature (℃), respectively and Ra is the extraterrestrial solar radiation (W/m 2 ) estimated based on the location's latitude and the calendar day of the year. This equation has a link to solar radiation through Ra and takes into account the impact of radiation warming the surface to the ground by the term of change between maximum and minimum temperature.

Application of spatial interpolation method
To discover the connections among the recorded values at different sites in the area, the trend test statistic values from measured locations are weighted spatially using deterministic interpolation techniques. In the present work, the inverse distance weighted interpolation (IDW) method was applied for spatially interpolate the meteorological (i.e., rainfall and potential evapotranspiration) time series data at the annual and seasonal scales over the Modjo watershed. this method is deterministic interpolation technique because it is directly based on the surrounding measured values, which works based on the assumption that the interpolating surface should be influenced most by nearby points and less by more distant points (Gemmer et al. 2004). It weights the points closer to the prediction location greater than those do farther away. By the IDW method, each station has a local influence, which decreases with distance by means of the use of a power parameter (Pingale et al. 2014). The following formula (Yavuz and Erdogan 2012) were applied in using the IDW technique: Where, Z (Lo) is the prediction value for location Lo, N is the number of measured sample points surrounding the prediction location, λi is the weight assigned to each measured point, Z(Li) is the observed value. The power parameter (p) can be any real number greater than 0. A higher power results in less influence from distant points. A large p results in nearby points wielding a much greater influence on the unconsidered location than a point further away resulting in an interpolated output (Hussain et al. 2021). To minimize the interpolation errors, the power parameter (p) =1 was adopted at the annual and seasonal time scales to have more detail output surface and to ensure a high degree of local influence.

Trend detection and serial correlation effects
In order to understand the direction and magnitude of streamflow and climate changes, the nonparametric Mann-Kendall test (hereafter M-K test) trend test (Mann 1945;Kendall 1975) have been applied in our study area (Amo-Boateng et al. 2014;Daniel et al. 2017;Pirnia et al. 2019b;Worku et al. 2019;Yimer et al. 2020;Latif et al. 2021). However, to apply this test, the time series data should be serially independent, otherwise it may lead to underestimation or overestimation of trend results (Kulkarni andVon Storch 1995, Zhang et al. 2000). Therefore, to reduce the influence of serial correlation on the M-K trend results, the procedure of trend free pre-whitening (TFPW) suggested by Yue et al. (2002) was considered and sufficient in the present study. The steps used in TFPW start with de-trending the original series by removing linear trend in the data using: Yt= Xt x βt (4) Next, checking the presence of auto-correlation with lag-1 serial correlation coefficient (r1) in the de-trended series using the following equation: Where Yt is the de-trended series, Xt is the original data series value at time t, ym is mean of the trended series data, β is Sen's slope and n is the data series. In this study, the Z-statistics was tested for the significance of trend by comparing it with two-tailed threshold levels of 1.96, 1.65 and 2.58 for 5%, 1% and for 10 % levels of significance, respectively. Following that, the condition was checked using: If a satisfied condition is reached, then it can be considered that the data series is independent at the 5%, 1% and 10% levels of significance, respectively, and the M-K test needs no removal of the significant autocorrelation (i.e. trend analysis can be applied to the original data series). Otherwise, the autocorrelation part will be removed using Eq. (9) and linear trend is added to the new series using Eq. (10) Y1= Yt-r1Yt-1 (9) Y2= Y1-βt Here, Y1 is a series without auto-regressive part, and Y2 is new series without auto-regressive and linear trend in the original data series. In this study, refer to the results obtained with the pre-whitening approach, all the time series are independent at the considered levels of significance, except annual rainfall data in Modjo station, which is dependent at 5% and 1% significant levels and the significant autocorrelation was removed as per the above correction procedure for this particular station before applying the M-K test. In addition, using the Kolmogorov-Smirnov test, normality test was conducted at α= 0.05 significant level, and none of the time series are following normal distribution. Hence, there is no ground for parametric approach for trend detection in this study.
Following the autocorrelation and normality tests, the temporal trends on the climate and hydrology data was assessed through the nonparametric M-K test. The null hypothesis (Ho) of randomness states that the data (x1, . . ., xn) are a sample of n independent and identically distributed random variables (there is no trend). The alternative hypothesis (Ha) is that there is a monotonic (not necessarily linear) trend. The M-K test, operation is based on the following equation: Here, S denotes the test statistics, j and k are the years, and n is the length of data (years) for a time series. In general, values of S indicate the direction of decreasing or increasing trend.
This is to mean that the difference (xj − xk), when j>k, represents a function (indicator) which exhibits 1, 0, or −1 (values). The standard normal test statistic Z of the trend is computed using: For presence of a statistically significant trend, the calculated standard Z value is compared with standard normal distribution with two-tailed significance levels. In this study, trends were checked at α = 5%, 1% and 10% significant levels, unless stated otherwise. In the trend analysis, the null hypothesis (Ho) of no trend is rejected if |Z| > 1.65 at 10% significant level (90% confidence interval); if |Z| > 1.96 at 5% significant level (95% confidence interval); and if |Z| > 2.58 at the 1% significant level (99% confidence interval). Moreover, an increasing (upward) trend should be interpreted if Z > 0 or the time series experiences a decreasing (downward) trend if Z < 0, where as if Z= 0, a trend due to random fluctuation, entirely there is no trend (both increase and decrease trends have similar impacts) in the time series.

Estimation of trend magnitude
Basically, the M-K test is incapable in determining the magnitude of trend changes in time series. For this purpose, the nonparametric Theil-Sen approach (hereafter TSA) (Theil 1950;Sen 1968) was applied to estimate the magnitude of trend changes in the study. Like the M-K test, the TSA is also insensitive to normality and outliers which are common in hydro-climatological datasets. In this approach, the magnitude of the slope line (djk) year, the Theil-Sen Slope (βi) and percentage of trend per year (Tp) were estimated as follows, respectively.
Where βi is the slope of trend line, x denotes the data point, n is the number of data, Tp is percentage trend per year, πave is the average value of the precipitation, streamflow, and evapotranspiration over the period of record, j and k are the corresponding times. If the βi value is positive, indicates an upward or increasing trend, when the estimated βi value is negative, a downward or decreasing trend.

Correlation between runoff and climate
In the study, the dependency of runoff changes on climate variables (precipitation, reference evapotranspiration and temperature) was assessed in terms of temporal correlation analysis. Initially, it is important to understand the relationship between independent (precipitation, evapotranspiration and temperature) and dependent variables (runoff) may be linear or non-linear. In this purpose, the first step would be to determine which model best describes the relationship between the two variables significantly. For that matter, linear model, was selected to analyze the variable most affects the streamflow of the catchment for the period 1986 -2015 by considering climate data near the streamflow gaging station. For this purpose, Pearson correlation coefficients (r) (Pearson 1904;McCuen 2003) was estimated, using the equation below: Where Q and C represent time series of runoff and climate (rainfall, temperature and evaporation) variables, respectively and n is the length of the data series. The r value varies between −1 and 1, the closer an r value is to 1 or −1, the stronger the correlation; the closer the r value is to 0, there less to correlation (close to uncorrelated situation). Following the estimation of the coefficient of correlation (r), coefficient of determination is determined by squaring the r value (denoted as r 2 ). The r 2 value expresses proportion of the total variation in the dependent variable.

Runoff characteristics
The preliminary analyses of hydrological time series data are important to understand if the data possess the appropriate characteristics for trend analysis. In the study, different types of statistical analysis were applied for the stated purpose. The statistical parameters of monthly, seasonal and annual runoff series are presented in Table 2. The annual average runoff of Modjo watershed is about 211.5 m 3 for the period from 1986-2015. Annual runoff depth varies between 33.67 m 3 (1987) to 1298.2 m 3 (1995). The coefficient of variation (CV) was 122.4% (Table 2), indicating that the runoff was highly instable and varied greatly from year to year, and difficult to predict. August is the month with the highest, while February is the month with lowest streamflow in the watershed. In dry months, it is also observed that precipitation fall, however, base flow from the highlands and mountainous regions in the northwestern part of the watershed is source of water for these months. Similar to the annual case, monthly flows are highly variable. No month has been experienced with moderate flow variation over the watershed. However, it can be seen that high flow months are relatively less variable than low flow months in the watershed. In general, due to the highly flow variations, monthly streamflow patterns are difficult to predict in the study area.
The important task in hydrology and climate distribution is to characterize the variability and location of dataset. For this purpose, skewness and kurtosis statistical parameters are fundamentally important. Skewness is the measure of relative symmetry. A dataset is symmetric if it looks the same to the right and left of the center point. The runoff frequency distribution showed positive skewness (Sk > +1) at the monthly, annual and seasonal time scales (Table 2), indicating runoff in the watershed is asymmetric and it lies to the right of mean. In general, the runoff time series data distribution is highly skewed. Kurtosis is measure of the relative tail (potential outlier) character of the distribution. A dataset with large kurtosis tend to have outliers or heavy tails. In the present study, positive kurtosis was found for all the months except for August and September (Table 2). August and September are with negative kurtosis (kurtosis <3), that is less extreme events and less likely in severe events than the normal distribution. In the rest of the months, a distribution with larger kurtosis (kurtosis >3) was obtained. This is leptokurtic distribution (i.e. with longer tails than normal distribution), which have a high likelihood of extreme events and unpredictable incidents. According to Westfall (2014), because of higher kurtosis values in the runoff distribution, more of the variability is resulted from infrequent extreme deviations.

Temporal trends in streamflow
Streamflow/runoff is output of catchment interplay between natural and human over a drainage basin. Therefore, trend analysis of streamflow provides useful information for understanding changes associated with these factors. For this purpose, the runoff trend of Modjo watershed was detected at monthly, annual and seasonal time scales over the period 1986 to 2015 (Table 3). Based on the M-K and TSA tests, trends with increasing slope were determined from October to January, with a weak variation in magnitude of +0.026, +0.035, +0.024 and +0.017 mm every year. Whereas, the trend values in the months of Feb, Apr, June, July, August and September were declined, with a weak magnitude in the month of June (β= -0.357 mm per year), while strong significant reduction in July (β= -2.57 mm), in August (β= −4.03 mm), and in September (β= -3.25 mm). The temporal trends of seasonal runoff were also analyzed for the same period. Table 3 shows the seasonal trend characteristics of runoff. Accordingly, the long rainy (Kiremt; Jun-Sept) and short rainy (Belg: Feb-May) seasons experienced downward trends, with Z = −3.68, and Z=-0.464, respectively. However, the M-K test result indicated that the runoff series in Belg season did not reach the considered significant levels, which means that the declined trend was not significant. In contrast, runoff in Bega (Oct-Jan) season was significantly increased, (Z value = 2.78) at 5% significant level over the studied period (Table  3). This may be due to base flow contribution during the season. Similar finding was reported by Mekonnen et al.
. From the temporal analysis result (Table 3), it can be seen that the annual flow discharge exhibited a highly significant decreasing trend (Z>2.58, p<0.001), with a strong trend magnitude of -1.61 mm every year (or 0.46 % decline per annum). It can be seen that trends of the annual precipitation, reference evapotranspiration, and runoff are not consistent temporally. During the period, precipitation shows a slight decreasing trend (insignificantly) whereas the annual ETo experienced a highly significant increasing trend for the whole study. Regarding trend direction, however, the declining trend in annual streamflow of Modjo watershed was consistent with the slightly decreased mean annual rainfall magnitude over the watershed (see Table 9). This incompatibility in value may be due to overexploitation of the water resources for domestic and industrial activities in the region. The result obtained here indicates that due to characterization of a negative trend in river flow series, a possible future water resource scarcity can be expected in the watershed and its surroundings. Moreover, if the dry years are successive, there will be drying tendency and insufficient water availability. Therefore, there is a need for enhancing water use efficiency in different sectors, and appropriate watershed management as well as conservation practices in this watershed. The findings of this study is consistent with the study conducted by Mahtsente et al. (2019) and Eshetu (2020) who documented significantly decreasing trend for annual and seasonal flow over the periods they considered.

Trends in precipitation
Precipitation is the major input in the rainfall-runoff relation over a natural basin. Any change in this variable can significantly affect the final output of a basin. Mean annual precipitation of Modjo watershed was estimated using the Thiessen polygon method (Thiessen 1911). Thiessen weighted factor (ρ=Ai/A) is areal weight assigned to each gauge by the polygon to identify their contributions. Mean annual precipitation and the weights of each station are given in Table 4. It has been observed that out of the seven analyzed gauging stations six of them had areas of influence inside the watershed (i.e., ρ > 0). The mean annual precipitation estimated for the period 1986 -2015 was 889.2 mm ( Table 4). Coefficient of variation (CV) is a statistical measure of dispersion around the mean, which is the ratio of standard deviation to the mean. Based on the results, the variability of annual precipitation was from less to moderate ranges, except in Koka that revealed high variability (39.3%) during the studied period. The finding of this study is consistent with other studies recently carried out by Daniel et al (2017) over the Awash basin, and Eshetu (2020) in Modjo watershed, who reported a moderate variability in annual rainfall records for some of the stations. Skewness is the measure of symmetry or asymmetry of data. In this paper, the data is negatively skewed (Sk= -0.33) indicating annual precipitation dataset is negatively skewed or skewed left, meaning that the left tail is longer. In general, the annual precipitation dataset is approximately symmetric (i.e. -0.5<Sk<0.5). The other statistical parameter estimated was kurtosis. It is a shape parameter that characterize the relative tail (potential outlier) character of the distribution. The precipitation of the watershed characterized by a Ku value of -0.15 (Table 4), compared to the normal distribution, it is less sharply peaked with the same mean and standard deviation (Westfall 2014). The results obtained here are compatible with the finding reports by Mahtsente et al (2019) for the stations in common.
It can be seen that watershed's precipitation peaks in August, with a well-defined beginning of the major rainy season period in June (Fig. 2). The mean annual precipitation of the watershed varies between 677.5mm (1986) and 1178.5mm (1993), with standard deviation of 121 mm/year (Table 4). The long rainy (kiremt: Jun to Sept) season is the major rainy period, which is characterized by heavy rainfall accompanied by a substantial surface runoff in the study region. On average, rainfall during this period subsidizes about 75% of the watershed's total annual precipitation showing the existence of high intensity of rainfall. The short rainy period (Belg: Feb to May), contributes to a light amount of precipitation (20% in average) in the area. The rainfall precipitated during this season is highly valuable, especially for cropping of short growing plants. The major problem as far as rainfall distribution is concerned in the study area is more of the inconsistency as well as a change in onset and cessation periods than the total amount. The rain onsets late and ends up very early, which makes the cropping time being shorter than before. In the area, the dry (Bega: Oct -Jan) season is the period, which contributes a very low amount of precipitation. Trend analysis of rainfall was examined using the M-K test and TSA. For monthly precipitation, six stations showed a significant declined trends in February. The precipitation in Aleltu station, was significant decreasing trend for some of the months (February, March and April), while significantly upward trend for the month of November (Tables 5). Statistically insignificant trend was exhibited in monthly precipitation for all months in Bishoftu station (Table 5). Overall, most of the monthly precipitation trends for this station are in the upward direction. Monthly trend was statistically significant in Chafe donsa station only for February and May months at α=0.05 and 0.001 significant levels, respectively (Table 5). The Theil-Sen's test also confirmed the downward trend (DT) for February and upward trend (UT) for May with a weak value of β= -0.82 and strong value of β= +1.31 mm per year, respectively. Edjere station is located in eastern part has experienced a significant trend only in two months (February and July) at the considered significant levels. The rainfall distribution for May and July of this station showed strong trend magnitudes of β= -1.2 and β= 2.89 mm per year, respectively (Table 5). Overall, monthly precipitations at Edjere station are in upward trends for most of the months. Modjo station revealed statistically significant trend for February, May and July months at 10% significant levels. Based on Theil-Sen (Table 5), Modjo showed a significantly downward changes for February (β= -0.26), but significantly upward trend for May and July months, with strong magnitudes of β= 1.11, and β= 3.82 mm per year, in that order. Zequala station located west of the watershed, showed a significant trend only for February and April months, having a strong declining magnitudes of -1.88 and -4.3 mm every year, respectively (Table 5). In general, monthly precipitation rate in this station show downward trends, except in November which experienced numerically increasing trend for the period 1986 -2015.
Trend analysis was also conducted on annual precipitation of the seven gauging stations. Three stations (Aleltu, Modjo and Zequala) experienced significant trends over the period 1986 -2015 (Table 5). Among the three stations,   (Table 5). Generally speaking, it was observed that the annual rainfall in high altitude stations is decreasing while this was increasing trend in low altitude stations. The maximum decline in annual precipitation was observed at Zequala station, at 0.1% significant level. Whereas, the maximum increase in annual precipitation trend was observed at Modjo station at the 5% significant level (Table 5). In general, the majority of the stations show a non-significant trends. Similar to this study, Daniel et al. (2017) reported that a non-significant decreasing tendency in annual rainfall at Debrezeit (Bishoftu) station and an increasing trend for Koka dam station. Mahtsente et al, (2019) found a declined trend in annual rainfall for Aleltu, and Bishoftu stations, whereas an inclined annual rainfall changes for Koka dam and Modjo stations with difference in significance. On the contrary, some findings are not compatible with our results for different ranges of reasons. For example, Eshetu (2020) reported an insignificant declined trend in annual rainfall at Aleltu and Edjere stations, while increasing insignificantly in annual rainfall series was detected for Modjo station for the study period she considered. The possible reasons for such difference goes to the data period and approach followed for quality control assessment.  Table 7 presents the summary of results of the M-K and TSA statistical methods applied in ETo trend analysis. As can be seen, positive trend types were observed in the majority of the stations in Modjo watershed. In the study, three out of the four stations exhibited statistically significant increasing trends (p < 0.10) for the annual time series.

Trends in evapotranspiration
According the M-K test, the magnitude of trend in the annual ETo of Aleltu station located in the far north (outside the watershed) revealed insignificant and decreasing trend slope with a value of -0.215 mm per year. The annual trend for Koka, Bishoft and Modjo stations revealed highly significant increasing trend (Table 7). Regarding the monthly PET time series, the majority of months have a positive trends and trend magnitude ( Table 7). As can be seen, Koka dam station showed upward trends of ETo in all months. For this station the minimum and maximum values of trend magnitude are estimated in January (0.11 mm/year) and in June (0.703 mm/year), respectively. The monthly ETo trend over Modjo station showed statistically significant trend except for May. In this station 11 out of 12 months revealed significantly increasing trend over the study period. In Modjo station, the maximum trend magnitude was determined for the month of September (20.3% per period), and the minimum estimate was for January (0.063 mm/year). In Bishoftu station, statistically significant and increasing trends were detected in all months, except for June and December that showed insignificantly increasing and decreasing trends, respectively. In Bishoftu, the maximum value of slopes in ETo (0.317 mm/year) was observed in the month of October, whereas a minimum value (0.04 mm/year) was found in the month of August (Table 7).

Precipitation trend variability
The spatiotemporal trends of rainfall time series were examined for Modjo watershed. As discussed before, the nonparametric statistical methods of M-K statistics and TSA trend magnitude estimator were applied to check the presence of temporal trends in seasonal and annual precipitation time series. These temporal trends are then interpolated using the inverse distance weighted (IDW) interpolation method to show their spatial distributions over the study watershed (Fig. 3). Accordingly, the spatial distribution of annual rainfall indicated that the rainfall pattern of Modjo watershed had great variation aerially from high to low altitude areas. As can be seen from Fig. 3a, the annual total rainfall distribution of the watershed declined north to south and west to east. The northern and part of the western of the watershed received the highest annual rainfall of the period that ranges from 1105 mm to 1145 mm (Fig. 3a). In the period, the lowest annual rainfall was observed in the southwestern and eastern parties of the watershed, which had annual rainfall of the time between 821 mm and 861 mm. The spatial distribution of seasonal rainfall is shown below (Fig. 3b to 3d). Fig. 3b presents the spatial distribution of Kiremt (major rainy: June -September) season rainfall over the study area. In the period, Kiremt season rainfall was almost followed the same spatial distribution as that of the annual rainfall, except some changes in the central part of the watershed (Fig. 3b). During Kiremt, the highest rainfall values were recorded in the western, central, and north-eastern parts while the lowest rainfall values recorded in the southeastern and northwestern parts of the study area. The Belg season (short rainy: February -May) is the second contributor of rainfall in the area. Fig. 3c displays Belg season rainfall trends over the period 1986 -2015. During this season, the highest rainfall values were observed in the western (which has the highest topography in the watershed), followed by the southern parts, whereas the lowest rainfall values were recorded in most parts of the northern, and southwestern parts of the watershed (Fig. 3c). This may have an impact for the agricultural productivity of the area as the area is one of the main farming of cereal crops. The dry (Bega: Oct-Jan) season is the period with insignificant precipitation contribution to area. Fig. 3d shows the spatial distribution of rainfall during the dry (Bega) season. During this season, the western, parts of eastern, and southeastern of the Modjo watershed received relatively maximum rainfall value, while the northern and southwestern parts received the lowest rainfall amounts over the studied period.   4 demonstrates interpolated maps of temporal trends in annual, and seasonal precipitation series of the study area. As the area is dominated by rain-fed agricultural activities, this is crucial to understand the spatial variabilities of rainfall patterns over the area. Moreover, the results would be useful for planning and development measures of future extreme cases such as flooding, drought and early warning systems. Accordingly, the results indicate that the highland rainfall stations showed decreasing trends on annual basis, whereas this was detected as decreasing trend for the lowland stations except for Bishoftu station. Regarding the annual rainfall trend, the result reveals negative trends were followed in the western, north and northwestern parts of the watershed (Fig. 4a). In the same period, the annual precipitation of the watershed was revealed insignificant positive trends in eastern, and southeastern areas of the watershed.
Considering seasonal precipitation trends, the trend was not statistically significant at the considered confidence intervals for Kiremt, and Bega seasons during the period 1986-2015. However, the trend detected for Belg season was significant at the watershed scale. Kiremt rainfall (Jun -Sept) declined insignificantly over the study area (Fig.  4b). During this season, negative trends in rainfalls were detected in the southern, northern, southwestern and northwestern parts of the watershed. The spatial distribution of Belg (Feb -May) season rainfall is shown in Fig. 4c. In this season, decline in rainfall in most parts of the watershed, were detected. Among all areas, the highest decreasing trend in rainfall was observed in the central and southwestern locations. The Belg season is an important period in the study area in which some short period crops are cultivated. Based on the obtained results, the effect of rainfall variability on crop production and food security was high over the study region, but this is severe in the central and southwestern parts of the study area. Bega (Oct -Jan) season is the time with insignificant rainfall amount over the study area. It can be seen that a negative trend was found in the western part of the watershed (Fig.  4d). Fig. 4 Spatial distribution of temporal trends in annual and seasonal rainfall as measured by the M-K test.

Evapotranspiration trend variability
The spatial distribution in annual and seasonal reference evapotranspiration (ETo) was done using the IDW interpolation method (Fig. 5a to 5d). Accordingly, the spatial distributions of the annual ETo decreased from north to south and from west to east (Fig. 5a), and the annual ETo value was below 800 mm in north and western parts of the watershed. As demonstrated in Fig. 5, the annual mean ETo in the southeast was larger than that in the southwest. The spatial variation of ETo indicates the maximum annual ETo (>845 mm) was occurred at the southeastern and the minimum were from the western and northern parts of the region (Fig. 5a). However, the rage of variation is not significantly high, indicates low spatial variability in annul ETo was experienced over the watershed for the period 1986 -2015.
It is analyzed that the seasonal variability in ETo followed similar distribution as that of the annual spatial distribution of ETo. For example, the long rainy (Kiremt) season spatial variability of ETo was the lowest in north and western parts of the watershed, while the south and southeastern regions do have the highest spatial distributions (Fig. 5b). The areal averaged ETo for Belg season is shown in Fig. 5c. During this season, the south and central parts experienced with high spatial distribution in ETo. It is also observed that the north and western parts of the region were with the lowest distribution in ETo for the analysis period, which are relatively with high altitude areas of the watershed. A bit difference in areal coverage of spatial distributions of ETo was observed during the dry (Bega) season of the area. As can be seen in Fig. 5d, the spatial distributions ETo for Bega was greater than 315 mm in most areas of the watershed, and the maximum was at the low altitude areas of the watershed. Overall, low values of ETo were observed for the high-altitude stations, and high values were characterized for the low altitude areas in the watershed.   6 demonstrates the spatial distributions of temporal trends in seasonal and annual ETo across the study area. On the annual scale, ETo trends exhibited that higher value in the lowland areas of the basin mainly in the southern, southeastern and south western parts while the lower values were found in the far northern and eastern portions of the watershed (Fig. 6a). The spatial distribution of annual ETo shows upward trends in most parts of the watershed. In sum, the lowest spatial and decreasing trend values were detected in the far northern part (outside the study watershed).
Unlike the seasonal precipitation trends, the seasonal ETo trends were changed significantly over the studied period. The spatial variation of trends in ETo during long rainy (Kiremt: Jun -Sept) season is demonstrated in Fig.  6b. It was observed that similar trends with the annual ETo, except some variations in the central region of the watershed. In this season, the northern highlands were detected the lowest trend value than the southern lowland areas. While, the rest of the catchment had the highest trends in the Kiremt season. The ETo of the short rainy (Belg: Oct-Jan) season (Fig. 6c) demonstrates some spatial variations of trends in the central region of the catchment. For example, in the eastern, and northeastern areas of the watershed, the ETo trend showed declined trend in addition to the far northern parts of the regions during the Belg season. For the same period, a bit difference in areal coverage of spatial distributions of ETo was observed during the dry (Bega: Oct-Jan) season of the area. Unlike spatial distribution in the wet seasons, in Bega season increasing trend value covers large areas in the western parts, with Z value ranges between 4.45 and 7.95. In general, ETo in the Bega season exhibited an increasing trend when going south to north and west to east parts of the study area.

Examining relation between runoff and climatic variables
Areal averaged trend analysis in precipitation and evapotranspiration time series data was detected for the period 1986 -2015 ( Fig. 7 and Table 8). In view of that significantly decreasing trend in precipitation was found in February, and April months, numerically declined trend for January, March, June, and August months, whereas trendless time series was detected for the rest of the months (Table 8). The Theil-Sen slope test also confirmed that a strong decreasing trend for February (β= -1.13 mm per year, or 143% of decline) and April (β= -1.83 mm each year, or 89.7%). February is a critical month, the time in which Belg rain starts in the area. Therefore, this result indicated late beginning of precipitation during this season over the study area. On the other hand, a statistically significant increasing trend was observed for November (Z= +3.13, p=0.004), the month used for crop harvesting in the area. Therefore, this situation has an adverse impact on harvesting time of crops in the region.
Annual and seasonal trend analysis were conducted at the watershed scale considering precipitation data of 1986 -2015. On annual basis, watershed's mean precipitation decreased (Z= -0.71) numerically, with a magnitude of 0.71 mm/year ( Fig. 7a and Table 8). Even though annual precipitation experienced a downward trend, it was statistically insignificant (p=0.097) at the 95% confidence interval. Similarly, the wet seasons (i.e. Kiremt and Belg) experienced downward trend, but it was significant only for Belg season (Z<1.96, p<0.05). The decrease in wet seasons rainfall over the watershed could not be compensated by an increase either monthly or annual time scales. For instance, the decrease in Kiremt and Belg seasons rainfall could not be compensated by an increase in the constitute months or annual time scale (Table 8). Our result agrees with the findings of previous papers in similar and related areas. For example, Mulugeta et al. (2019) reported a non-significant decreasing trend of annual rainfall on the entire Awash river basin (in which our study area is situated). Among others, Cheung et al. (2008) and Viste et al. (2013) found insignificant trends in annual, and summer rainfall for the entire Awash River Basin considering different periods. However, some paper reported a contradicting results. For example, Wegesho et al. (2013) reported a significant decreasing trend for the period 1951-2000, over most parts of the Awash River Basin. Cheung et al. (2008) and Viste et al. (2013) documented a no trend condition in spring (February to May) rainfall over the entire Awash basin. This incompatibility may be due to record length and the months constituting in each season. Similar to areal precipitation, trend analysis of reference evapotranspiration (ETo) was also conducted at monthly, seasonal, and annual data in Modjo watershed for the period 1986-2015 (Table 8). For monthly time series, the M-K test showed an upward trend direction for all months over the watershed (Table 8). Out of 12 months, 9 of them were detected to had significant increasing trends at the considered significant levels. The minimum and maximum values of trend magnitude are in January (0.033 mm/year, or 5.8% of total change) and in September (0.231mm/year, or 6% of total change), respectively. The ETo trend in the watershed is significantly upward at the annual basis, with a weak magnitude of β= 0.147 mm every year (Table 8). Overall, the annual ETo of Modjo watershed increased by 0.01% per year between 1986 and 2015.
Regarding seasonal trend, Kiremt (Jun-Sept), Belg (Feb-May) and Bega (Oct-Jan) seasons, the ETo trend increased by 5.51, 7.9 and 4.27% per period, at significant levels of 5%, 1% and 10%, respectively. Analyzing the ETo trend in Modjo watershed is an innovative, for this reason the result obtained here cannot be compared with other studies within the study area. Trend tests for the terms in the Hargreaves equation including Tmin, Tmax, Tmax − Tmin and (Tmax + Tmin)/2 was also performed to identify the factorial influences on PET trend at the watershed scale. According to the M-K test, Tmin decreased significantly at the annual and seasonal time scales. In this case, the annual, and seasonal (i.e. Kiremt, Belg and Bega seasons) Tmin trend values declined at a rate of -0.0025, -0.002, -0.0024 and -0.0029 o C/year, respectively, which were significant at 1%, 1%, 10% and 1% levels. Regarding Tmax, trendless time series were detected at annual, and seasonal time scales. However, the maximum temperature of the watershed increased insignificantly during the studied period. In the same period, the annual, and seasonal Tmax -Tmin term increased significantly (Z= +3.996, +3.64, +4.103 and +3.211, respectively) at the considered confidence interval, with magnitude of +0.0031, +0.0025, +0.0037 and +0.0032 o C/year in that order. On the other hand, the term '(Tmax + Tmin)/2' decreased insignificantly for the annual, Kiremt, and Belg time series, with magnitude of +0.0007, +0.0007, and +0.0005 o C/year, respectively, but significantly declined for the dry (Bega) season at 10% significant level for the studied period. The significant upward trend in 'Tmax -Tmin' and the numerical increase in Tmax have relevance in the increasing trends for the potential evapotranspiration of the watershed than the other terms in the Hargreaves equation. The findings of this study also suggest that the scarcity of available water is getting worse, along with rise in evapotranspiration. The contribution of climate (rainfall, temperature and evapotranspiration) fluctuations over time to the variability in the annual mean streamflow was done using the correlation approach. The result indicated that strong dependence of monthly river discharge on precipitation, with correlation coefficient (r) of 0.91 at α=0.05 significance level. However, the correlational relationships of runoff with ETo and temperature were not significant (with r values of 0.22 and 0.13, respectively). The results suggested that precipitation was the main climatic factor, whereas ETo and temperature have limited contribution for the observed monthly flow changes in the Modjo watershed. Similarly, multiple linear regression of precipitation and runoff as well as temperature and runoff showed weaker correlation with insignificance level (p>0.05) at the annual scale (Fig. 8). However, the coefficient of correlation between annual runoff and ETo is stronger than annual runoff with other climate parameters. Therefore, this suggests that annual rainfall and surface air temperature made insignificant contributions (r= 0.14 and 0.31, respectively) to runoff changes between 1986 and 2015. The portion of the annual runoff values explained by precipitation and temperature fluctuations is about 1.8% and 9.7%, respectively (Fig. 8). More than 90% of the annual runoff declined in the study watershed, therefore, influenced by the temperature derived evapotranspiration, land use changes, progressive human-induced factors (such as excessive water withdrawal for municipal, industries, etc.) and other factors in the near past decades. 0 500 1000 1500 1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 2006 2008 1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 2006 2008 1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 2006 2008

Conclusion
In the present study, using Modjo watershed as a concentration area, the spatiotemporal variability in precipitation, reference evapotranspiration (ETo), and streamflow was investigated over the period 1986 -2015. Theil-Sen approach and the Mann-Kendall test were used for identifying the temporal trends whereas, the spatial variability of trends in mean precipitation and ETo time series was analyzed using the inverse distance weighted (IDW) interpolation technique in the ArcGIS environment. We observed a declined trend of runoff in some of the months, annually and wet season time scales according to the calculated Z and β values obtained for the period 1986 -2015. However, this is not true for the areal averaged precipitation and ETo of the watershed, which revealed an insignificant decreasing and a significant increasing trends, respectively. In this point, analysis of monthly and seasonal runoff showed that the significant decreased of annual flow is the result of decreasing trends in some of the months and wet seasons during the studied period. Moreover, the runoff characteristics was highly variable and more difficult to interpret in the study area.
Precipitation trend revealed both downward and upward trends but not significant in most of the months except for February, April and December. Analysis of seasonal precipitation data demonstrated that substantial declined trend during Belg season was resulted from the declined trends during the constitute months of February, March and April. In this time, nearly 143% and 89.7% of significant declined trends were detected in precipitation for February and April, respectively. The spatial interpolation results showed that the western and northern parts of the watershed were with high rainfall distribution. However, these parts of the watershed were affected by decreasing precipitation and increasing ETo trends during the studied period. We observed an increasing trend in reference evapotranspiration (ETo) for most stations and at catchment scale. A significant incline in ETo during 1986-2015, with a consistent increase in the amount of ETo in most data series of the stations and the catchment. However, the increased trend in the watershed scale was generally less than 10% for all time scales, except for the month of March (i.e. 10.14% per period). From the correlation analysis, we observed that the coefficient of determination (r 2 ) between annual runoff and climate variables (precipitation, temperature & ETo) was small (<50%), this indicates influence of these parameters on annual flow is weak for the considered period. This shows that other factors such as land use, soil infiltration, water extraction, deforestation, etc. seems to be more important than the temporal climate changes in explaining the observed changes in flow series. However, the vulnerability is growing in terms of significantly increased in ETo, declined in precipitation, and in the end decreased river flows during the recent period, particularly since the 1990s. The vulnerability can be stabilized by implementing sustainable watershed managements, appropriate conservation measurements and efficient use of the water resources.