Climate and human activities impact on runoff and sediment yield in the central rift valley of Ethiopia

Abstract

This study analyzed the effects of climate variability and human activities on runoff and sediment yield of a watershed in the Central Rift Valley of Ethiopia since the 1980s. In the study, innovative trend analysis, the Pettitt test, and sequential regime shift detector methods were used to analyze trends and abrupt changes during the period 1981–2020. The Soil and Water Assessment Tool model was then integrated to disaggregate the contributions of climate variability, land-use change, and direct human activities on runoff and sediment yield at the watershed scale. The result indicates that the runoff and sediment yield decreased significantly with a sudden change in 2001. The study found that direct human activities played a dominant role, accounting for 57.9%, and 53.4% of annual runoff and sediment yield reduction, respectively. On the contrary, changes in climate and land use have led to a slight increase in runoff and sediment yield after 2001. The greater impacts of human activities can be attributed to the growth rate of water abstraction, and riverbed sand excavation during the past four decades. The findings of the result could be useful in supporting decision-makers in achieving watershed sustainability.

Keywords:

• Central rift valley
• climate change
• human activity
• ITA method
• SWAT model

REVIEWING EDITOR:

• Sanjay Kumar Shukla, Edith Cowan University, Australia

SUBJECTS:

• Agriculture & Environmental Sciences
• Earth Sciences
• Civil, Environmental and Geotechnical Engineering

1. Introduction

Climate variability and anthropogenic activities are the two non-independent factors that affect the spatiotemporal distribution of river discharge and the corresponding transporting sediments. It is increasingly reported that changes in runoff and sediment yield result from the impacts of climate variability, land use change, and human activities in different basins (Dai et al., 2023; Luan et al., 2021; Zhao et al., 2015). A wide range of studies noted that the intensification of climate change and human activities is undergoing significant changes in sediment generation and runoff processes (Guo et al., 2022; Poerbandono et al., 2014; Zhang et al., 2020). Shifts in climate characteristics influence the amount and distribution of precipitation and temperature, which in turn affects the spatiotemporal patterns of surface runoff (Wang et al., 2008) and sediment transport (Zhang et al., 2020). Human-induced land use dynamics can affect sediment transport from catchments (Zhao et al., 2012). Therefore, quantitative estimation regarding the hydrological impacts of climate and human activities greatly helps managers and decision-makers for watershed management and sustainability in the future.

Understanding how climate and human activities affect a watershed’s water-related changes is complex because of the intricate connections between elements like the landscape, climate, and water-related factors. Therefore, it is very important for land use planners and water resource managers to understand and quantify these impacts at the catchment scale. Various techniques, including statistics (like flow duration curve, double mass curve, multi-statistical regression model, redundancy analysis), conceptual models, and physical hydrologic models, are commonly employed to measure how these factors influence runoff and sediment in rivers (Chawla & Mujumdar, 2015; Guo et al., 2022; Pirnia et al., 2019; Wu et al., 2017; Xue et al., 2017; Zhang et al., 2020). The statistical methods are simple to apply and require less data but basically lack a physical basis. Hence, such methods are impractical to use for spatiotemporal analysis of hydrological components under climate and environmental changes. The conceptual models are mainly important for analysis on long-time scales and do not reflect changes in the flow process (Wu et al., 2017). On the other hand, the physically based hydrologic models (comprising the semi and fully-distributed models) are by far the most commonly used and dominant tools to quantify the impacts of climate variability and human interventions on hydrologic processes in different regions (Guo et al., 2016; Kumar et al., 2022; Mekonnen et al., 2018; Shang et al., 2019; Teklay et al., 2021; Zhao et al., 2018). However, owing to their extensive parameterization and intensive data requirements, the fully distributed models are difficult to employ in a data-scarce watershed. Among the physically-based and semi-distributed models, the Soil and Water Assessment Tool (SWAT) model is one of the widely used models to quantify the impacts of climate change, and human activities on water and sediment yield in river basins (Arnold et al., 1998).

Streamflow and sediment loads are essential elements of watershed ecology and environmental dynamics. The analysis of changes in these parameters, along with the factors influencing them, is vital for sustainable management and has garnered significant attention among researchers worldwide (Gao et al., 2011; Yao & Liu, 2018; Zhang et al., 2020; Zuo et al., 2016). It is also widely observed that the impacts of climate variability and human activities on river discharge and sediment load are significant in several river basins worldwide (Pirnia et al., 2019; Zhao et al., 2015; Zuo et al., 2016). In this regard, previous researchers have reported different findings on the attribution to river flow and sediment yields. For example, the studies by Ma et al. (2009), and Karlsson et al. (2016) reported LULC changes have a larger impact than climate change in runoff. Zhao et al. (2015) found that climate variability and human activities reduced contributions by about 72 and 28% as well as 14 and 86% of the river discharge and sediment load, respectively. Ahmad et al. (2022) indicate that the main driver of runoff change is climate change. They also add that the individual impacts are more important than the combined effects for the design of irrigation facilities, the planning and management of water resources, and the regional decision-making policies. These divergent outcomes suggest that distinct river basins possess varying climate change patterns, land cover conditions, anthropogenic influences, surface and topographical characteristics, all of which culminate in diverse effects on river flow and sediment transport across different scales.

Central Rift Valley of Ethiopia exhibits the highest level of variability in terms of socioeconomic activities, land use changes, climate patterns, and environmental conditions. The combination of a growing population and environmental changes is causing a decline in natural resources in this area. Additionally, the region is experiencing increased sediment accumulation in rivers due to both natural processes and environmental changes. This buildup of sediment raises the concern that river channels could become filled over time and low water availability (Yao & Liu, 2018). In this region, deforestation and direct human activities (e.g. multipurpose water abstraction) have already disturbed the hydro-environmental of the river systems during the past decades. Most studies in central, Ethiopia and its surrounding areas have emphasized the hydrological impacts of climate variably (Bambrick et al., 2015; Daniel & Abate, 2022; Gebrechorkos et al., 2019; Godebo et al., 2021), or land use land/cover change changes (Ayalew et al., 2022; Belihu et al., 2020; Desta & Lemma, 2017; Elias et al., 2019).

However, the combined effects of climate and human activities on streamflow and transporting sediments have received less attention in this part of the country. As a result, the lack of complete information on water availability, and sediment transport poses serious challenges and massive problems for sustainability. In addition, in recent years, river flow has experienced dramatic changes in the central rift valley of Ethiopia (Alemayehu et al., 2006; Belete et al., 2015) that constrained the sustainable development of the society economy, and hydro-ecology of the area. However, there are no scientific investigations on the streamflow and sediment yield changing characteristics in the area that support decision-makers in achieving watershed sustainability. For this, quantitative estimation and attribution analysis of sediment load and surface runoff changes can help to illustrate reasons for the change and understand the complex hydrological response in watersheds situated in this region of the country.

The current study is designed to enhance the understanding of how climate variability, changes in LULC, and direct human activities impact river discharge and the corresponding sediment dynamics at the watershed scale. To achieve this, a combination of statistical methods, geospatial techniques, and hydrologic modeling approaches were employed within the context of the Modjo watershed. The specific objectives of the study include: (i) assessing annual and seasonal trends in hydro-climatic time series; (ii) estimating the quantity of surface runoff and sediment yield under various scenarios; and (iii) quantifying the contributions of influencing factors behind changes in runoff and sediment yield during the period from 1981 to 2020. The originality of this research resides in the exploration of the intricate relationships among climatic fluctuations, anthropogenic interventions, and their collective impacts on the dynamics of runoff and sediment yield within the central rift valley of Ethiopia. The hypothesis asserts that alterations in both climate conditions and human activities exert substantial influence over the patterns of runoff and sediment yield within the central rift valley of Ethiopia, resulting in quantifiable alterations to water flow dynamics and sediment accumulation processes. The study findings can help in providing detailed information regarding water resource management, utilization as well as environmental sustainability in the future.

2. Materials and methods

2.1. Characterization of the study area

This study was carried out in the Modjo watershed, which is located in the Central Rift Valley of Ethiopia within the upper Awash basin (Figure 1). The watershed extended between 8°35′00″N to 9°05′11″N and 38°54′35″E to 39°15′30″E (Zena et al., 2022). The river watershed has a high topographic variation that ranges between 1579 and 3083 m above sea level. The Modjo River is the main river of this watershed that originates from the northwest mountainous area of the watershed and travels a total distance of about 128.4 km in the south, and eventually meets the Awash River, upstream of the Koka hydropower dam station.

Figure 1. Location of the study watershed.

In the region, rainfall occurs during two wet seasons, the summer (June to September) and spring (March to May), which depends on the seasonal fluctuation of the Inter-Tropical Convergence Zone (ITCZ). Precipitation is maximum during July and August, while minimum during the month of December (Figure 2). Based on the data from 1981 to 2020 (https://www.ethiomet.gov.et), the average annual precipitation varies from 846.7 to 1121.2 mm. The annual mean minimum and mean maximum temperatures vary between 10.6–11.6 and 21–27 °C, respectively. Meanwhile, the annual average temperature fluctuates between 16.9 and 19.3 °C (Gessesse et al., 2015). The maximum mean temperature occurs during the month of May and the minimum is in December.

Figure 2. Monthly mean precipitation and temperature in the Modjo watershed.

The dominant land use type is agricultural land, covering over 80% of the total watershed area. Additionally, there are smaller areas with shrubs, woodland, and grassland. The primary agricultural system in the study area is mixed farming, which integrates both crop and livestock production (Zelalem & Kumar, 2018). Crop productions are entirely based on rainfall but some root crops (carrot, potato, onion, etc.), and vegetables (cabbage, tomato, chili pepper, etc.) are produced using small-scale irrigation, along the Modjo River. In general, there are eleven major soil types present. Among these, Pellic Vertisols, Vertic Cambisols, and Luvic Phaeozems are the most extensively distributed, encompassing 42, 33.56, and 12.62% of the watershed area, respectively.

2.2. Data collection

The Digital Elevation Model (DEM) of the study watershed has been downloaded from the Alaska satellite facility website (freely available at https://vertex.daac.asf.alaska.edu/). The DEM has a resolution of 12.5 by 12.5 m, and it was utilized for tasks, such as watershed delineation, slope classification, and generating hydrological response units (HRUs). The stream network dataset provided by the Ministry of Water and Energy (https://www.mowe.gov.et) was superimposed onto the DEM to improve hydrographic segmentation and sub-basin boundary delineation.

A digital soil map of 1 km resolution was obtained from the Ministry of Water and Energy (MoWE), and Harmonized World Soil Database (HWSD_v1.2). After extracting the soil map of the watershed using ArcGIS; the FAO (Food and Agriculture Organization) soil classification system was used (Xijun & Ying, 2019). For each soil category, the percentage of clay, sand, and silt is estimated by averaging the associated SoilGrids250m (https://soilgrids.org) data. SPAW (the Soil-Plant-Atmosphere-Water) model (Available at https://hrsl.ba.ars.usda.gov/soilwater/Index.htm) is also used to determine soil water content based on clay, sand, and silt fractions. Land use data for the years 1990, 2000, and 2020 at moderate resolution (30 by 30 m) were downloaded from the United States Geological Survey (USGS) database (https://earthexplorer.usgs.gov). The first two images were from Landsat 5 Thematic Mapper (TM), and the third image was from Landsat 8 Operational Land Imager (OLI). Using the hydro-climate trend results, two LULC raster maps representing the baseline and changed periods were employed for scenario development to quantify the consequences of influencing factors between the two periods.

Daily meteorological data (1981–2020), including precipitation and minimum/maximum air temperature, were collected from the National Meteorological Service Agency of Ethiopia (https://www.ethiomet.gov.et). The meteorological datasets underwent checks for missing records, homogeneity, correctness, and sufficiency. With missing values amounting to <3%, linear regression was employed for the rainfall data, while the simple average method was applied to fill missing air temperature data.

Climatic data homogeneity was assessed using the Standard Normal Homogeneity Test (Alexandersson, 1984). The Double Mass Curve (DMC) graphical analysis method was used to address inconsistencies in station data by comparing its time trend with that of adjacent stations. The results indicated consistency at a 5% significance level for all stations except the Modjo station. Inconsistencies at the Modjo station were adjusted using the DMC method, and the corrected data were then utilized as inputs for hydrologic modeling and trend analysis. Meteorological and hydrological station details are presented in Table 1.

Table 1. Ground stations for hydro-climatic datasets of the study.

Data type	Station	Elevation (m)	Data used information
Meteorological stations	Aleltu	2663	Daily rainfall and min/max Temp (1981–2020)
	Bishoftu	1900	Daily rainfall and min/max Temp (1981–2020)
	Chefe donsa	2392	Daily rainfall and min/max Temp (1981–2020)
	Edjere	2254	Daily rainfall and min/max Temp (1981–2020)
	Koka dam	1618	Daily rainfall and min/max Temp (1981–2020)
	Modjo	1763	Daily rainfall and min/max Temp (1981–2020)
Hydrological station	Modjo river	1593	Daily river discharge (1983–2020)

Similarly, daily discharge and monthly sediment yield data were obtained from the MoWE at https://www.mowe.gov.et. These datasets were measured at the Modjo River gauging station (Elevation = 2180 m, Latitude = 8.40°N, Longitude = 39.023°E) and exhibited missing values ranging from days to months. To address this, the datasets underwent quality control, involving the adjustment of outliers and filling in missing records.

The Markov Chain Monte Carlo method, utilizing Multiple Imputation Algorithms, available in the XLSTAT software package (accessible at https://www.xlstat.com), was employed to fill in the missing records. Due to the low frequency of sediment data collection in the country, the relationship between discharge and sediment transport (i.e. the assessment curve) was utilized to estimate sediment yield data for in-situ calibration and validation of discharge. Sediment data, originally obtained based on concentration, were converted into sediment yield.

2.3. Methodology

In this study, we assessed how surface runoff and sediment yield change over time and space due to climate, LULC, and human activities, utilizing hydrological modeling and statistical methods. First, we examined changes in hydro-meteorological factors between 1981 and 2020. Subsequently, different scenarios were developed to understand the factors driving these changes. Finally, a calibrated SWAT model was employed to quantify the impact of climate, LULC changes, and human activities on the spatiotemporal variations of streamflow and sediment yield, both individually and in combination.

2.3.1. Temporal trend analysis

In this study, various statistical methods were employed to analyze time series data, including trend examination and change point detection. The Innovative Trend Analysis (ITA) method, a non-parametric technique introduced by Sen (2012), was applied to explore temporal trends within the hydro-meteorological datasets. ITA has gained popularity in various global regions (Li et al., 2019; Shah et al., 2022; Zena et al., 2022). This method was chosen for its ability to detect both monotonic and non-monotonic trends without imposing conditions of serial autocorrelation, presence of outliers, or specific record length in the time series, which are requirements for other similar methods (Shah et al., 2022). Moreover, ITA is considered more sophisticated compared to traditional methods, such as Mann–Kendall and linear regression (Phuong et al., 2020; Zena et al., 2022).

The ITA method starts by dividing the time series data into two equal parts and arranged in ascending order. The first part of the time series (X_i: i = 1, 2, 3, …, n/2) is placed on the horizontal x-axis, and the other time series (X_j: j = n/2 + 1, n/2 + 2, …, n) is on the vertical y-axis. The slope parameter of the ITA method was considered as the indicator of changes in the trend value of time series data. As per Şen (2017), the straight-line slope of both monotonic and non-monotonic trends can be estimated using Equation (1)(1) $$S_{\mathrm{\text{ITA}}}=\frac{2x(y_{\mathrm{j}}‐\mathrm{ }{x}_i)}{n}$$(1) : (1) $$S_{\mathrm{\text{ITA}}}=\frac{2x(y_{\mathrm{j}}‐\mathrm{ }{x}_i)}{n}$$(1) where S_ITA is the ITA slope value; x_i and y_j are the arithmetic mean of the first and the second halves of the sub-series, respectively, and n is the total number of observations. The positive and negative value of S_ITA indicates an increasing and decreasing trend in time-series data in that order. Equation (2)(2) $${\sigma }_{S_{\mathit{\text{ITA}}}}^2=\frac{8{\sigma }^2}{n^3}(1-{\rho }_{\overline{y1}}{\rho }_{\overline{y2}})$$(2) is applied to estimate the standard deviation of the sampling slope. (2) $${\sigma }_{S_{\mathit{\text{ITA}}}}^2=\frac{8{\sigma }^2}{n^3}(1-{\rho }_{\overline{y1}}{\rho }_{\overline{y2}})$$(2)

The term ${\rho }_{{\overline{y}}_1}{\rho }_{{\overline{y}}_2}$ is the correlation coefficient between the ascendingly/discerningly sorted two values, and can be estimated using Equation (3)(3) $${\rho }_{{\overline{y}}_1}{\rho }_{{\overline{y}}_2}=\frac{E\left({\overline{y}}_1{\overline{y}}_2\right)-E\left({\overline{y}}_1\right)E\left({\overline{y}}_2\right)}{{\sigma }_{{\overline{y}}_1}{\sigma }_{{\overline{y}}_2}}$$(3) : (3) $${\rho }_{{\overline{y}}_1}{\rho }_{{\overline{y}}_2}=\frac{E\left({\overline{y}}_1{\overline{y}}_2\right)-E\left({\overline{y}}_1\right)E\left({\overline{y}}_2\right)}{{\sigma }_{{\overline{y}}_1}{\sigma }_{{\overline{y}}_2}}$$(3)

Ultimately, the upper and lower confidence limits (CL) of the Innovative Trend Analysis (ITA) slope at the significance level α can be computed utilizing Equation (4)(4) $${\mathit{\text{CL}}}_{(1-\alpha )}=0\pm S_{\mathit{\text{crt}}}x{\sigma }_{S_{\mathit{\text{ITA}}}}$$(4) , as detailed by Zena et al. (2022). (4) $${\mathit{\text{CL}}}_{(1-\alpha )}=0\pm S_{\mathit{\text{crt}}}x{\sigma }_{S_{\mathit{\text{ITA}}}}$$(4)

Here S_crt represents the critical slope for standardized time series. The ITA method has increasingly been used in different basins of Ethiopia (Gedefaw et al., 2018; Girma et al., 2020; Zena et al., 2022), among others.

2.3.2. Change point detection

In this study, a hydrological model and statistical techniques were combined to figure out how climate change and human activities affect river flow and sediment at the watershed scale. Before using the model, the Pettitt test and sequential regime shift detector (SRSD) methods were used to find breakpoints in hydro-meteorological data and divide the period into smaller parts for comparison. The non-parametric Pettitt test is one of the commonly used methods to detect a single change point in hydro-climate datasets (Pettitt, 1979). Considering X₁, X₂, X₃, …, X_t as a set of random variables, the test statistic (KT) that signifies the location of the change point is computed using Equation (5)(5) $$K_T\hspace{0.17em}=\hspace{0.17em}\text{max|}{U}_t{,}_T|$$(5) , as proposed by Pettitt (1979): (5) $$K_T\hspace{0.17em}=\hspace{0.17em}\text{max|}{U}_t{,}_T|$$(5)

Here, K_T is the point of change year occurrence, and the non-parametric test statistics U_t,T is described in Equation (6)(6) $$U_{t,T}=\sum _{i=1}^t\sum _{j=t+1}^T\mathit{\text{sgn}}(X_i-X_j)$$(6) : (6) $$U_{t,T}=\sum _{i=1}^t\sum _{j=t+1}^T\mathit{\text{sgn}}(X_i-X_j)$$(6)

Similarly, the SRSD is a non-parametric trend change detection method (Rodionov, 2004) that was used to identify abrupt change points in time series data. Unlike the Pettitt test, it has the advantage of detecting multiple change points and measuring changing confidence levels in time series data. Additionally, the method does not require an initial visual inspection of a time series and a priori hypothesis about the timing of the shift of the former (Liu et al., 2013). The details of the SRSD method are as follows. For time series {x_k, k = 1, 2, …, n}, the mean of the first regime (R1), is determined first using Equation (7)(7) $${\overline{x}}_{R1}=\sum _{k=1}^l{x}_k(1\le k\le l)$$(7) . (7) $${\overline{x}}_{R1}=\sum _{k=1}^l{x}_k(1\le k\le l)$$(7)

Where; l is the cut-off length of the regimes to be determined. Then, the difference from the mean of the second regime (R²) that would be statistically significant (Student’s t-test) can be estimated as per Equation (8)(8) $$\mathrm{\text{diff}}={\overline{X}}_{\mathrm{\text{new}}}-{\overline{X}}_{\mathrm{\text{curr}}}=t\sqrt{2{\sigma }_l^2/l}$$(8) . (8) $$\mathrm{\text{diff}}={\overline{X}}_{\mathrm{\text{new}}}-{\overline{X}}_{\mathrm{\text{curr}}}=t\sqrt{2{\sigma }_l^2/l}$$(8)

Where t is the t-distribution value with 2l—2 degrees of freedom at the given probability level p, ${\sigma }_l^2$ is the average variance for l-year intervals in the time series.

In the SRSD method, each x_k (k ≥ l + 1) is evaluated in sequential order (Rodionov, 2004). If x_k is within the range of [${\overline{x}}_{R1}‐\mathrm{\text{diff}},$ ${\overline{x}}_{R1}+\mathrm{\text{diff}}$], the ${\overline{x}}_{R1}$ is recalculated with the x_k value and (l − 1) previous x_k values. If x_k is out of [${\overline{x}}_{R1}‐\mathrm{\text{diff}},$ ${\overline{x}}_{R1}+\mathrm{\text{diff}}$], it is considered as a possible start point (k = j) of the new (second) regime (R²). In this case, the regime shifts index (RSI) is determined to confirm or reject the null hypothesis of the abrupt change starting at year j is estimated using Equation (9)(9) $${\mathrm{\text{RSI}}}_{i,j}=\sum _{i=j}^{j+m}\frac{x_{\mathrm{i}}-{\overline{x}}_{R1}+\mathrm{\text{diff}}}{l_{{\sigma }_l}}(m=0,1,\dots ,l-1)$$(9) : (9) $${\mathrm{\text{RSI}}}_{i,j}=\sum _{i=j}^{j+m}\frac{x_{\mathrm{i}}-{\overline{x}}_{R1}+\mathrm{\text{diff}}}{l_{{\sigma }_l}}(m=0,1,\dots ,l-1)$$(9)

The performance of the SRSD method can be controlled by three parameters (Rodionov, 2015): (1) the significance level (p); (2) the cut-off length (l), and (3) Huber’s tuning constant (h). The significance level p (also called the target probability level) controls how sensitive SRSD is to shifts in the variance. Since the critical value t_crt decreases with the increase of p, SRSD becomes sensitive to smaller and smaller shifts. The cut-off length (l) controls the time scale of the detected regime shift. Huber’s tuning constant is important only in the presence of outliers. For this particular study, a cut-off length of (l) = 20 years, Huber’s tuning constant (h) = 1, and a probability level at p = 0.05 were considered. If there are multiple change points, the change point with the lower p-value is chosen as the final one (Rodionov, 2015).

2.3.3. Land use classification and change detection

Preprocessing of satellite images before image classification and change detection has been done by considering a series of sequential operations, including atmospheric correction, geometric correction, image registration, and masking for clouds, water, and irrelevant features. The geometric correction was performed based on control points and topographic maps obtained from the Ministry of Water, Resource, and Energy of Ethiopia (https://www.mowe.gov.et). The images used had already been geo-referenced and corrected for all sensor irregularities. The image is re-projected into the local system (WGS_84_UTM_Zone_37N).

After image preprocessing, a supervised classification using the maximum likelihood decision rule was employed to identify the primary land use and land cover (LULC) categories within the study watershed (Kumar et al., 2022). Field surveys conducted between December and March 2020, along with existing aerial photographs from MoWE, were used to gather classification training samples. Changes in land use and land cover (LULC) were detected using a LULC map from 1990 as a base. The primary LULC categories included barren land, built-up areas, eucalyptus trees, rain-fed agricultural land, small-scale irrigation, natural forests, grasslands, scattered shrubs, and water bodies.

For simplicity, dense and moderate natural forests, open shrub land, sparse natural forests, and plantation forests were grouped together as forest land. Similarly, closed, scattered, and open shrub land cover types were combined as shrub land. Additionally, the scattered irrigation agriculture and rain-fed agricultural lands were merged into cropland due to their small size and distribution. As a result, the main LULC types within the studied watershed were barren land, built-up areas, cropland, forest, grassland, shrub land, and water bodies.

In the remote sensing literature, numerous precision measurements have been developed to summarize information obtained from an error matrix (Liu et al., 2007). In this study, three LULC maps from the years 1990, 2000, and 2020 were utilized. The accuracy of the 2020 LULC categories was tested and verified based on 245 real ground control points recorded by GPS during a field survey conducted from December to March 2020. The ground reference data were collected only after confirmation by local elders in the area. Similarly, 185 random points were created for the years 1990 and 2000. For ‘n’ number of categories (Table 2), ‘A_i+’ and ‘A_+i’ can be used to denote, respectively, the sum of rows and columns for category ‘i’. Accordingly, for the total categories, the following identities are established (Equation 10(10) $$\sum _{i=1}^n{A}_{i+}=1\text{and}\sum _{i=1}^n{A}_{+i}=1$$(10) ): (10) $$\sum _{i=1}^n{A}_{i+}=1\text{and}\sum _{i=1}^n{A}_{+i}=1$$(10)

Table 2. An n by n transition matrix for classification accuracy assessment.

Classified	References
	1	2	…	…	n	Row totals	Ua
1	A11	A12	…	…	A1n	A1+	A11/A1+
2	A21	A22			A2n	A2+	A22/A2+
…	…	…	…	…	…	…	…
…	…	…	…	…	…	…	…
n	An1	An2	…	…	Ann	An+	Ann/An+
Column totals	A+1	A+2	…	…	A+n	1
Pa	A11/A+1	A22/A+2	…	…	Ann/A+n

A_ij represents the proportion of area in the mapped LULC category i and the reference land cover category j.

User and producer accuracies, common measures for assessing individual land use categories (Liu et al., 2007; Story & Congalton, 1986), are part of various approaches to evaluate thematic map accuracy based on error matrices. Table 3 outlines category and map-level accuracy measures, considering Overall Accuracy (OA), Aii (total reference value of category ‘i’), and ‘n’ (number of categories) from the error matrix in Table 2. There is no consensus on a standard approach for accuracy assessment in thematic maps within the remote sensing community (Foody, 2002; Liu et al., 2007). While the Kappa coefficient is widely used, Foody (1992) suggests the modified Kappa (Brennan & Prediger, 1981) to address the overestimation of chance agreement. This study incorporates not only the Kappa coefficient but also quantitative disagreement and allocation disagreement, as proposed by Pontius and Millones (2011).

Table 3. Category and map-level accuracy measures used in the study.

No	Formula	Name	References
1	Ua = Aii/Ai+	User’s accuracy	Story & Congalton, 1986
2	Pa = Aii/A+i	Producer’s accuracy	Story & Congalton, 1986
3	OA = ${\sum }_{i=1}^n{A}_{\mathit{\text{ii}}}$	Overall accuracy	Story & Congalton, 1986
4	K = (OA − Ea)/(1 − Ea) Ea = ${\sum }_{i=1}^n{(A}_{i+}{A}_{+i})$	Kappa coefficient	Cohen, 1960; Rosenfield & Fitzpatrick-Lins, 1986
5	mK = (OA − 1/n)(1 − 1/n)	Modified kappa	Aickin, 1990; Foody, 1992
6	DQ = 1/2$\times {\sum }_{i=1}^n\left|A_{i+}-A_{+i}\right|$	Quantitative disagreement	Pontius & Millones, 2011
7	DA = 2 min (Ai+; A+i) − 2Aii	Allocation disagreement	Pontius & Millones, 2011

2.3.4. Hydrologic modeling

The study utilized the Soil and Water Assessment Tool (SWAT) model (Arnold et al., 1998) to simulate surface runoff and sediment yield in the watershed. SWAT is a physically-based semi-distributed model with a continuous daily time-step, adaptable to various climate and land-use change scenarios (Omer et al., 2017). The model divides the basin into sub-basins connected to the river network based on a digital elevation model, further classifying them into hydrological response units (HRUs) based on a unique combination of land use, soil type, and slope. SWAT then simulates hydrological components from each HRU to the sub-basin level, aggregating them to the basin outlet through the stream network.

SWAT (Soil and Water Assessment Tool) simulates watershed hydrology, as outlined by Neitsch et al. (2011), through two key phases: (i) the land phase, which governs the water and sediment loading to the main channel, and (ii) the routing phase, responsible for estimating the movement of water and sediments through the channel network to the outlet. The hydrological component of the SWAT model utilizes the water balance equation (Equation 11(11) $${\mathrm{\text{SW}}}_i={\mathrm{\text{SW}}}_o+\sum _{i=1}^t(R_i-{\mathrm{\text{SQ}}}_i-{\mathrm{\text{ET}}}_i-W_i-{\mathit{\text{GW}}}_i)$$(11) ) within the soil profile to simulate all watershed processes, as outlined by Arnold et al. (2012). (11) $${\mathrm{\text{SW}}}_i={\mathrm{\text{SW}}}_o+\sum _{i=1}^t(R_i-{\mathrm{\text{SQ}}}_i-{\mathrm{\text{ET}}}_i-W_i-{\mathit{\text{GW}}}_i)$$(11) where SW_t = the total soil water content after t days (mm); SW_o = the initial soil water content on a given day (mm), R_i = the amount of precipitation depth on a given day (mm), SQ_i = the amount of surface runoff on a given day (mm), ET_i = the amount of evapotranspiration on a given day (mm), W_i = the amount of water entering the vadose zone from the soil profile on a given day i (mm), GW_i = the amount of base flow on a given day (mm), and t = time (in days).

2.3.5. Simulation of runoff and sediment yield

SWAT simulates the hydrology of a watershed (Neitsch et al., 2011): (i) by using the land phase (which regulates the amount of water and sediment loading to the main channel), and (ii) by the routing phase (that estimates the movement of water and sediment discharge through the channel network to the outlet). The hydrological part of the SWAT model uses the water balance equation in the soil profile to simulate all the watershed processes at the sub-basins (Arnold et al., 2012). In the SWAT model, two methods are available for surface runoff calculation: one is the Soil Conservation Service (SCS) runoff curve number method and the other is the Green–Ampt infiltration method (Wang et al., 2008). The SCS curve number (CN) method was considered to estimate (Equations 12(12) $$S=\frac{25.4\left(1000\right)}{\mathrm{\text{CN}}-10}$$(12) and 13(13) $$Q_{\mathrm{\text{surf}}}=\frac{{(R_{\mathrm{\text{day}}}-I_a)}^2}{(R_{\mathrm{\text{day}}}-I_a+S)}$$(13) ) the total surface runoff obtained at the outlet of the watershed (Neitsch et al., 2011): (12) $$S=\frac{25.4\left(1000\right)}{\mathrm{\text{CN}}-10}$$(12) (13) $$Q_{\mathrm{\text{surf}}}=\frac{{(R_{\mathrm{\text{day}}}-I_a)}^2}{(R_{\mathrm{\text{day}}}-I_a+S)}$$(13) where S is the retention parameter; CN is the runoff curve number, which would vary across all HRUs as each one has a special combination of land use, soil, and slope; I_a is the initial abstractions that include interception, surface storage, and infiltration before runoff. The study used the variable storage method, which is based on a simple continuity equation to route the flow of water in the channel. Based on the availability of data, this study also selected the Hargreaves method (Hargreaves & Samni, 1985) which uses the minimum and maximum temperature to calculate evapotranspiration.

Sediment transport and erosion are simulated by using the Modified Universal Soil Loss Equation (MUSLE) (Neitsch et al., 2011) in the sub-basins as in Equation (14)(14) $$\mathrm{\text{Sed}}=11.8{(Q_{\mathrm{\text{surf}}}\times q_{\mathrm{\text{peak}}}\times {\mathrm{\text{area}}}_{\mathrm{\text{HRU}}})}^{0.56}\times K_{\mathrm{\text{USLE}}}\times C_{\mathrm{\text{USLE}}}\times P_{\mathrm{\text{USLE}}}\times {\mathit{\text{LS}}}_{\mathrm{\text{USLE}}}\times \mathrm{\text{CFRG}}$$(14) : (14) $$\mathrm{\text{Sed}}=11.8{(Q_{\mathrm{\text{surf}}}\times q_{\mathrm{\text{peak}}}\times {\mathrm{\text{area}}}_{\mathrm{\text{HRU}}})}^{0.56}\times K_{\mathrm{\text{USLE}}}\times C_{\mathrm{\text{USLE}}}\times P_{\mathrm{\text{USLE}}}\times {\mathit{\text{LS}}}_{\mathrm{\text{USLE}}}\times \mathrm{\text{CFRG}}$$(14)

Where: Q_surf = the surface runoff volume (mm/ha), q_peak = peak runoff rate (m³/s), area_hru = area of HRU (ha), K = Universal Soil Loss Equation (USLE) soil erodibility factor, C = USLE cover and management factor, P = USLE support practice factor. Most parts of the watershed have no erosion control practices, as a result, we considered a value of 1 for the P factor. LS = USLE topographic factor and CFRG = course fragment factor. In this work, the values of the C_USLE factor were adopted from the SWAT manual (Agriculture = 0.2, Barren land = 0.2, Grass = 0.003, Forest = 0.001, water = 0, and shrubs = 0.003).

The sediment transport in the channel acts simultaneously which comprises degradation and deposit. The occurrence of either of these processes is determined by comparing the concentration of the sediments in the reach at the beginning of the time step (C_s,i) with the maximum sediment concentration (C_s,max). Sediment degradation and deposits in stream channels are calculated using the simplified Bagnold Equations (15)(15) $${\text{Sed}}_{\text{deg}}=\left({\text{conc}}_{\text{sed},\text{ch},\text{mx}}-{\text{conc}}_{\text{sed},\text{ch},i}\right)\times V.K_{\text{CH}},C_{\text{CH}};{\text{if conc}}_{\text{sed};\text{ch},i}>{\text{conx}}_{\text{sed};\text{ch};\text{mx}}$$(15) and (16)(16) $${\text{Sed}}_{\text{dep}}=\left({\text{conc}}_{\text{sed};\text{ch};i}-{\text{conc}}_{\text{sed};\text{ch};\text{mx}}\right).V;{\text{if conx}}_{\text{sed};\text{ch};\text{mx}}>{\text{conc}}_{\text{sed};\text{ch};i}$$(16) (Bagnold, 1977; Neitsch et al., 2011): (15) $${\text{Sed}}_{\text{deg}}=\left({\text{conc}}_{\text{sed},\text{ch},\text{mx}}-{\text{conc}}_{\text{sed},\text{ch},i}\right)\times V.K_{\text{CH}},C_{\text{CH}};{\text{if conc}}_{\text{sed};\text{ch},i}>{\text{conx}}_{\text{sed};\text{ch};\text{mx}}$$(15) (16) $${\text{Sed}}_{\text{dep}}=\left({\text{conc}}_{\text{sed};\text{ch};i}-{\text{conc}}_{\text{sed};\text{ch};\text{mx}}\right).V;{\text{if conx}}_{\text{sed};\text{ch};\text{mx}}>{\text{conc}}_{\text{sed};\text{ch};i}$$(16)

Where: Sed_deg is the amount of sediment re-entrained in the reach segment; Sed_dep is the amount of sediment deposited in the reach segment; conc_{sed, ch, i} is the initial sediment concentration in the reach (kg/L); conc_{sed, ch, mx} represents the maximum concentration of sediment that can be transported by the water (kg/L); K_CH is the channel erodibility factor; V is the volume of water in the reach segment (m³), and C_CH is the channel cover factor.

2.3.6. Model calibration and evaluation

Parameter sensitivity, calibration, and validation were conducted within the Sequential Uncertainty Fitting-2 (SUFI-2) optimization framework, using the SWAT-CUP tool (Abbaspour et al., 2007). Sensitivity analysis is the first step in the calibration process, which is needed to identify the parameters that significantly influence the hydrological processes of the watershed. In the present study, 32 parameters related to runoff and 20 parameters related to sediment transport were used to execute the sensitivity analyses. The most sensitive model parameters were selected by using the p-value and t-statistics, which identify the relative significance of each parameter (Arnold et al., 2012). The p-value determines the significance of the sensitivity, in which a value close to zero has greater significance and the t-test provides a measure of sensitivity, in which larger absolute values are more sensitive (Abbaspour et al., 2007).

Monthly streamflow and sediment concentration data from the baseline period were used to calibrate and validate the SWAT model. For the calibration purpose, simulation was performed by considering 60% of the baseline period, with the corresponding land use map. The SWAT model is then validated using an independent dataset without changing the calibrated parameters to ensure that the calibrated parameter set performs quite well in an independent dataset.

The performance of the SWAT model to predict the observed values was then evaluated based on statistical evaluators and graphical methods. Statistically, the Nash–Sutcliffe Efficiency (NSE) (Nash & Sutcliffe, 1970), coefficient of determination (R²) (Idrissou et al., 2022), as well as Percent Bias (PBIAS) (Gupta et al., 2009) were applied in the study. NSE measures how much of the variability in observations is explained by the simulations, and can be estimated using Equation (17)(17) $$\mathrm{\text{NSE}}=\frac{({\sum }_{i=1}^N{(O_i-O_{\mathrm{\text{avg}}})}^2)-(\sum _{i=1}^N{(S_i-S_{\mathrm{\text{avg}}})}^2)}{\sum _{i=1}^N{(O_i-O_{\mathrm{\text{avg}}})}^2}$$(17) : (17) $$\mathrm{\text{NSE}}=\frac{({\sum }_{i=1}^N{(O_i-O_{\mathrm{\text{avg}}})}^2)-(\sum _{i=1}^N{(S_i-S_{\mathrm{\text{avg}}})}^2)}{\sum _{i=1}^N{(O_i-O_{\mathrm{\text{avg}}})}^2}$$(17)

Statistically, R² ranges from 0 to 1, with the best model being that with R² approaching one. The R² value can be estimated using Equation (18)(18) $$R^2=\frac{\sum _{i=1}^N\left(O_i-O_{\mathrm{\text{avg}}}\right)(S_i-S_{\mathrm{\text{avg}}})}{\sum _{i=1}^N{(O_i-O_{\mathrm{\text{avg}}})}^2(S_i-S_{\mathrm{\text{avg}}}{)}^2}$$(18) : (18) $$R^2=\frac{\sum _{i=1}^N\left(O_i-O_{\mathrm{\text{avg}}}\right)(S_i-S_{\mathrm{\text{avg}}})}{\sum _{i=1}^N{(O_i-O_{\mathrm{\text{avg}}})}^2(S_i-S_{\mathrm{\text{avg}}}{)}^2}$$(18)

The positive and negative values of PBIAS are equivalent to overestimation and underestimation, respectively of the model results. The value for the PBIAS statistical criteria can be estimated by employing Equation (19)(19) $$\mathrm{\text{PBIAS}}=\frac{\sum _{i=1}^n(O_i-S_i)}{\sum _{i=1}^n{O}_i}\times 100\mathrm{\%}$$(19) : (19) $$\mathrm{\text{PBIAS}}=\frac{\sum _{i=1}^n(O_i-S_i)}{\sum _{i=1}^n{O}_i}\times 100\mathrm{\%}$$(19)

where O_i is the observed value; S_i is the simulated value; O_avg is the mean of measured data; S_avg is the mean of simulated value for the evaluated range of data, and N is the number of observations under consideration.

2.3.7. Strategy for attribution analysis

Climate change and human activities, both direct and indirect, are anticipated to influence streamflow and sediment transport in heavily human-altered watersheds (Omer et al., 2017). Notably, activities like agricultural expansion, urbanization, and water abstraction are key contributors to altering river flow and sediment load in the studied watershed. Direct human activities encompass soil conservation technologies and water resource abstraction for various purposes, while climate variability includes changes in precipitation and temperature, and land-use change encompasses agricultural expansion, urbanization, and removal of natural vegetation.

To disentangle the hydrological impacts of these factors, the study employed a one-factor-at-a-time approach, altering one factor while keeping the others constant. The specific procedures followed in this work include:

1. Determination of the baseline period. The entire period was divided into baseline and changed periods to quantitatively estimate the contribution of climate change and human activities to runoff and sediment yield variations. For this purpose, Rodionov’s sequential approach (Rodionov, 2004) and cumulative anomaly methods were used in the study. However, the separation of the baseline and changed periods determined by these methods may be characterized by frequent human activities during the baseline periods. For the betterment of the results, it is better to use both abrupt change detection and human-designed analysis approaches. However, there was no known and substantial human-designed analysis to determine the baseline period. Hence, the study applied the statistical methods alone for mutation analysis.

2. Developing scenarios to separate influencing factors. In the study, five simulations (S1, S2, S3, S4, and S5) consisting of two LULC maps (LULC maps before and after a change period) and two climate data periods (records before and after a change period) (Table 4). Scenario 1 (S1) is treated as a baseline simulation affected by the combination of LULC change, direct human activities, and climate variability. Scenario 2 (S2) considers the effects of clime change alone that developed using climate data from the changed period and LULC map from baseline. Scenario 3 (S3) is a land use change scenario that is developed based on climatic data during the baseline period and the LULC map from the changed period. Scenario 4 (S4) represents combined climate and LULC change impacts, which is based on climate records & LULC map of the changed period. Scenario 5 (S5) is used to assess the impacts of direct human activities on the changes in streamflow and sediment yield.

3. Estimating the contributions of climate variability, LULC change, and direct human activities. The total change in runoff (ΔQ) and sediment yield (ΔS) can be regarded as a function of climatic variables and the integrated effects of topography, soil, and human activities over watersheds. Under the assumption that the topography and soil of the studied watershed did not vary after the mutation year, the runoff and sediment changes between the changed (i.e. Q₂ and S₂) and baseline (i.e. Q₁ and S₁) periods were assumed to be the result of climate variability and human activities (Zhang et al., 2012). The effects of human activities on the variations of runoff and sediment load can be categorized as direct or indirect. In this study, the direct human activities assumed are human-related innervations other than the LULC changes over the watershed.

Table 4. Simulation scenarios of the SWAT model to separate the impacts of climate and human activities.

Scenario	Scenario description	Climate	LULC
S1	Baseline period	1981 to a change year	Before a change year
S2	Isolated impact of climate change	A change year to 2020	Before a change year
S3	Individual impact of LULC changes	A change year to 2020	After a change year
S4	Combined impacts of climate and LULC changes	1981 to a change year	After a change year
S5	Impacts from direct human activities	1981 to a change year	–

Referring to Table 4, the simulated values between S2 and S1 correspond to the distinct impacts of climate change on runoff and sediment yield. The variation in simulation values between S3 and S1 is attributed solely to changes in land use and land cover (LULC). In contrast, direct human activities encompass actions that directly influence the hydrological cycle, excluding LULC alterations. This includes activities, such as extracting water from river channels, constructing dams, managing reservoirs, implementing soil and water conservation projects, and withdrawing water from surface sources and river channels for agricultural, industrial, or domestic purposes. In this study, the hydrological response of the watershed was assessed by quantifying the impact of direct human activities. This was achieved by comparing the reconstructed values under S4 with the coincident observed sediment yield and runoff values during the changed period (2001–2020), and the resulting impact is represented in scenario 5 (S5).

In a similar vein, the difference between S2 and S4 can exacerbate the impact of LULC (Δ$Q_L^C$) and Δ$S_L^C$) caused by climate variability (Zang et al., 2013). The interactive effect of LULC by climate (ΔS_INT or ΔQ_INT) can also be estimated from the difference between the individual and amplified effect of LULC change (i.e. ΔQ_L − Δ$Q_L^C$ or ΔS_L − Δ$S_L^C$).

4. Estimating contribution rates of driving factors. The quantification of the relative and integrated contributions of influencing factors was conducted using expressions in Equations (20a)–(20c), as specified by Omer et al. (2017). (20a) $${ɳ}_C=\frac{Q_2-S1}{\Delta }\times 100\mathrm{\%}$$(20a) (20b) $${ɳ}_L=\frac{S3-S1}{\Delta }\times 100\mathrm{\%}$$(20b) (20c) $${ɳ}_{\mathrm{\text{DHA}}}=\frac{Q_2-S4}{\Delta }\times 100\mathrm{\%}$$(20c) where S1–S4 are simulated values under the baseline, climate, LULC, and combined climate and LULC scenarios in that order; Q₂ represents the recorded value during the changed period; ɳ_C, ɳ_L, and ɳ_DHA represent the relative contributions of climate variability, land use change, and direct human activities, respectively. A similar approach was followed for the quantification of sediment load under climate and human activities.

3. Results and discussion

3.1. Trends of hydro-climate variables

The hydro-meteorological records from 1981 to 2020 showed that the mean annual precipitation, runoff, sediment yield, and the runoff coefficient (the ratio of runoff to rainfall) of the watershed have experienced a significant downward trend, with S_ITA < 0, S_crt > S_ITA (S_ITA = −2.50 mm, 8.4 mm, 4.8 × 10⁴ ton, and 0.01 mm, respectively) at the 5% significant level (Table 5). The precipitation record from 1990 to 2015 showed a significant increasing trend, with S_ITA > 0, S_crt > S_ITA. Whereas, no significant trend for the mean annual evapotranspiration was detected in the same period (S_ITA = 3.30 > S_crt = 0.1062).

Table 5. Trends and abrupt changes for hydro-meteorological and sediment variables (1981–2020).

Variables	Time frame	ITA method				SRSD method
		SITA	CL(1−α)	Scrt	Trend	Cj	Conf	RSI
Precipitation	Annual	0.52	±0.376	0.859	Increasing	–	–	–
	Summer	2.39	±0.444	0.3751	Increasing	2018	6.53 × 10^−6	1.75
	Spring	−0.44	±0.0409	−0.2833	Decreasing	–	–	–
Runoff depth	Annual	−8.40	±0.1834	−0.600	Decreasing	2001	1.7 × 10^−6	−0.717
	Summer	−8.32	±0.2003	−0.627	Decreasing	2001	1.8 × 10^−6	−1.30
	Spring	−0.136	±0.0004	−0.028	Decreasing	2018	5.54 × 10^−3	0.457
Evapotranspiration	Annual	3.30	±0.00576	0.1062	Increasing	1999	1.12 × 10^−5	0.8680
	Summer	1.782	±0.00258	0.0715	Increasing	2011	6.83 × 10^−5	15.87
	Spring	0.3273	±0.00081	0.0397	Increasing	1999	3.57 × 10^−6	0.5097
Runoff coefficient	Annual	−0.010	±3.4 × 10^−7	−0.0008	Decreasing	2001	9.3 × 10^−7	−0.801
	Summer	−0.009	±4.1 × 10^−7	−0.0009	Decreasing	2001	1.1 × 10^−6	−0.841
	Spring	−0.0002	±1.5 × 10^−9	−2.78 × 10^−5	Decreasing	2018	3.375 × 10^−3	1.102
Sediment yield	Annual	−47749	±7,921,540	−3940.3	Decreasing	2001	6.9 × 10^−7	−0.891
	Summer	−47654.2	±9,008,839	−4202.1	Decreasing	2001	6.9 × 10^−6	−0.890
	Spring	−115.8	±58.86	−10.74	Decreasing	1999	0.016	−0.226

S_ITA: ITA trend slope value; CL: upper/lower confidence limit; S_crt: critical slope (at α = 5%); RSI: regime shift index; C_j: change points; Conf: confidence level.

Bold and italic numerical values indicate significant trend changes.

Statistical findings revealed decreasing trends in parameters, such as runoff depth, sediment load, and runoff coefficient for both the summer (June to September) and spring (March to May) seasons as indicated by S_ITA < 0, S_crt > S_ITA. While the summer season’s precipitation exhibited no significant trend, the minor rainy season displayed a notable decreasing trend, supported by S_ITA < 0, S_crt > S_ITA. Conversely, no significant trends were detected in the seasonal evapotranspiration of the watershed during the same period. The ITA method verified the dependability of seasonal and annual trends in hydro-meteorological parameters (including precipitation, evapotranspiration, sediment yield, and river runoff) within the study area.

The results of trend analysis showcased that negative trends predominantly influenced the flow and sediment load of the analyzed watershed. Precipitation serves as a fundamental input in watershed modeling for surface runoff and sediment yield outputs. This study concludes that the mean annual and summer season precipitation exhibited no statistically significant trends throughout the specified period. This implies that activities influenced by humans might play a more significant role in shaping runoff and sediment parameters in the studied watershed compared to precipitation. A study conducted by Qin et al. (2019) and Tian et al. (2009) elsewhere also demonstrated that alterations in the distribution of annual runoff are primarily attributed to human activities.

3.2. Abrupt change points

A year is considered a significant change point when both methods identify the same breakpoint in that year. The results indicate that abrupt changes in annual runoff, sediment yield, and runoff coefficient occurred in 2001 (Table 5). Figure 3 provides graphical representations of the abrupt changes based on the Pettitt test. The graph highlights significant alterations for the mean annual, summer, and spring seasons, while no trend changes were observed for dry season flow between 1981 and 2020.

Figure 3. Change points for (a) annual, (b) summer, and (c) spring seasons flow during the period 1981–2020.

Similarly, for the same period, the annual evapotranspiration exhibited a change point in the year 1998 at the 0.05 significance level (Table 5). The variation in abrupt change years between climatic and hydrologic variables may be attributed to the influence of climate variability and human activities during the analyzed period.

The streamflow and sediment yield underwent a notable decrease around 2001, contrary to the precipitation trend during the analysis period. A comparison of average values between the baseline (1981–2000) and changed periods (2001–2020) reveals a substantial reduction in runoff depth and sediment yield following increased human activities. Specifically, the annual runoff depth, sediment yield, and runoff coefficient exhibited significant declines after 2001, indicating negative impacts of human activities on surface runoff and sediment yield. For instance, the mean annual runoff depth decreased by approximately 73.7% (p-value: 1.2 × 10⁻⁶) (Table 5), and the sediment yield experienced a substantial 94.9% decrease compared to the baseline period. Concurrently, the runoff coefficient decreased by 74.03%, and runoff depth decreased by 73.70% (Table 6). Despite small increases in annual precipitation, evapotranspiration, and drought index (1.13, 6.56, and 5.1%, respectively) during the changed period, intensive human activities, such as land conversion, water abstraction, riverbed excavation, and conservation measures likely contributed significantly to the observed declines in runoff and sediment yields over the past decades.

Table 6. Annual runoff, precipitation, sediment yield, and evapotranspiration change between periods.

Variables	1981–2000	2001–2020	Change	Rate of change (%)
Runoff, R (mm)	422.94	111.10	−311.84	−73.7
Precipitation, P (mm)	915.5	925.8	10.3	1.13
R/P	0.462	0.120	−0.342	−74.03
Evapotranspiration, ETo (mm)	1006.1	1072.1	66.0	6.56
ETo/P	1.099	1.155	0.056	5.10
Sediment yield (10^5 t)	10.1	0.518	−9.58	−94.9

In recent times, several soil and water conservation measures, such as terraces, and check dams, were implemented in a scattered manner, which could be considered as the major factors to decrease surface flow in the watershed. These have substantial trapping effects on runoff generation, and sediment transports from the upstream parts of the watershed. Therefore, the impact of direct human activities (e.g. water resource abstractions, riverbed sand excavation, and soil-water conservation measures) had more significant impacts on decreasing river flow and sediment yield than the insignificantly increasing tendencies of climate and LULC change, particularly built-up area and barren land effects in the Modjo watershed.

3.3. Land use change detection

The Modjo Watershed has experienced significant land use and land cover (LULC) changes over the past decades due to factors, such as population increase, agricultural land expansion, and settlement areas. This study focuses on three periods of LULC data (1990, 2000, and 2020) for change detection, as presented in Table 7. The data reveals that a substantial portion of the watershed remained covered by cropland/farmland across all periods.

Table 7. Absolute and relative changes in LULC of the Modjo Watershed.

LULC types	LULC (ha)			LULC changes (in %)
	1990	2000	2020	1990–2000	2000–2020	1990–2020
Barren land	2134.71	2436.84	3181.32	14.2	30.6	49.0
Built-up	2071.62	2312.1	6621.03	11.6	186.4	219.6
Cropland	174752.91	181607	193894	3.9	6.8	11.0
Forest	8896.68	15051.4	8883.36	69.2	−41.0	−0.1
Grassland	6293.16	1390.59	880.47	−77.9	−36.7	−86.0
Shrub land	24394.23	15541.8	5109.66	−36.3	−67.1	−79.1
Water bodies	964.62	1168.29	938.07	21.1	−19.7	−2.8

The analysis further highlights a notable shift in the percentage of agricultural land, increasing from 79.61% in 1990 to 88.33% in 2020. Forested areas in the watershed, mainly consisting of plantation trees in high topography and rural settlement areas, showed no significant change between 1990 and 2020, indicating a dominance of plantation-type forest. In contrast, the percentage of land covered by woody/shrubs decreased by 79.1%, declining from 11.11% in 1990 to 2.33% in 2020.

The accuracy assessment matrices for the LULC maps of 2000 and 2020 are detailed in Tables 8 and 9. These matrices, commonly known as confusion matrices, illustrate the assignment of pixels to classes during the image classification process. In these tables, diagonal values indicate correctly classified pixels. The overall classification accuracy, calculated as the ratio of correct points to the total number of points, is 93.18% for the 2000 land use map and 96.2% for the 2020 map.

Table 8. A classification accuracy matrix for the LULC map of 2000.

LULC classes	Barren	Built-up	Cropland	Forest	Grassland	Shrubs	Water	Totals	Correct (C)
Barren	100	9	5	0	0	0	0	114	100
Built-up	17	199	0	0	0	0	0	216	199
Cropland	4	0	143	0	23	0	0	170	143
Forest	0	0	0	105	0	12	0	117	105
Grassland	0	0	10	0	100	0	0	110	100
Shrubs	0	0	0	10	0	38	0	48	38
Water	1	0	0	0	0	0	26	27	26
Totals	122	208	158	115	123	50	26	802	711
K	0.917
mK	0.956

Table 9. A classification accuracy matrix for the LULC map of 2020.

LULC classes	Barren	Built-up	Cropland	Forest	Grassland	Shrubs	Water	Totals	Correct (C)
Barren	73	1	1	0	0	0	0	75	73
Built-up	3	31	0	0	0	0	0	34	31
Cropland	1	0	170	0	1	0	0	172	170
Forest	0	0	0	97	0	2	0	99	97
Grassland	0	0	2	0	26	0	0	28	26
Shrubs	0	0	0	9	0	82	0	91	82
Water	0	0	0	0	0	0	34	34	34
Totals	77	32	173	106	27	84	34	533	513
K	0.953
mK	0.956

Additionally, Kappa (K) and modified Kappa (mK) coefficients were determined for both LULC maps. The coefficients for LULC 2000 were 0.917 (K) and 0.920 (mK), while for LULC 2020, they were 0.953 (K) and 0.956 (mK). The overall disagreement was found to be 3.8%, affirming the validity of the supervised image classification.

It is noteworthy that, according to Table 8, the K and the mK values are statistically the same. In contrast to Foody’s (1992) suggestion that the degree of random agreement might be overestimated with the kappa coefficient, we disagree. Our rationale is that there is no definitive way to determine the correct model of random agreement, rendering it impractical to label a measure as overestimating or underestimating random agreement without a standardized criterion for comparison. Consequently, based on our findings, no single measure, including K and mk, can be deemed the standard or default measure for randomly adjusted accuracy.

An error matrix for the accuracy assessment of LULC map 2020 is presented in Table 10. The proportion of correctly classified units is widely used in map-level image classification accuracy assessment. Based on Table 5, the correctly classified units, C = 0.962. As a result, the map-level classification disagreement is 0.038, which is an offset of the proportion correct (Pontius & Millones, 2011). Pontius and Millones (2011) decomposed the overall disagreement for image classification into the quantity disagreement (D_Q), and the allocation disagreement (D_A). The quantity disagreement reflects the amount of difference between the reference and comparison maps that is due to the imperfect matching in the total number of categories (Pontius & Millones, 2011). Summing all the category totals (i.e. A_i+ and A_+i) from Table 10 and dividing the result by 2 gives a D_Q value of 0.04.

Table 10. Summary of error matrix for accuracy assessment of LULC map 2020.

LULC classes	Barren	Built-up	Cropland	Forest	Grassland	Shrubs	Water	Ai+	C	Ua (%)
Barren	0.14	0.002	0.002	0.0	0.0	0.0	0.0	0.14	0.14	97.3
Built-up	0.01	0.06	0.0	0.0	0.0	0.0	0.0	0.06	0.06	91.2
Cropland	0.002	0.0	0.32	0.0	0.002	0.0	0.0	0.33	0.32	98.8
Forest	0.0	0.0	0.0	0.18	0.0	0.004	0.0	0.19	0.18	98.0
Grassland	0.0	0.0	0.004	0.0	0.049	0.0	0.0	0.05	0.05	92.9
Shrubs	0.0	0.0	0.0	0.02	0.0	0.15	0.0	0.17	0.15	90.1
Water	0.0	0.0	0.0	0.0	0.0	0.0	0.06	0.06	0.06	100
A+i	0.15	0.06	0.33	0.20	0.05	0.16	0.06	1.0
Pa (%)	93.3	96.8	98.8	91.0	96.3	97.6	100
DQ	0.04
DA	0.04
D	0.08

U_a: user’s accuracy; P_a: producer’s accuracy; C: correctly classified units; D_Q: quantity disagreement; D_A: allocation disagreement.

Similarly, the allocation disagreement shows a degree of disparity between the reference and the comparison maps because the correspondence is less than the maximum in the spatial distribution of the categories (Pontius & Millones, 2011). In Table 10, D_A = 0.04, which is the same value as that of the D_Q. Since the A_i+ = A_+i, any disagreement is an allocation disagreement. According to Pontius and Millones (2011), the total disagreement is 0.08 (i.e. 0.04 + 0.04). In general, the widely used map-level accuracy measure is the value of proportion correct, and its complement is the overall disagreement. In this study, the value of the disagreement measure is low, indicating a high degree of agreement between the reference and the comparison maps, and the distinction between the quantity disagreement and allocation disagreement suggested by Pontius and Millones (2011) is less important. The value of a disagreement measure can be considered substantial if it exceeds 0.1 (Warrens, 2015). If any of the map-level disagreement measures are high, the category-level measures can be used to identify which categories are causing the disagreement.

To evaluate the effect of climate variability, and human activities on runoff and sediment yield, the LULC of 2000 and 2020 were assumed to represent the periods 1981–2000 and 2001–2020, respectively (Table 11). There was a considerable increase in the areas of built-up, barren, and cropland, with respective proportions of 6.47, 47.54, and 184.58% during the changed period from 2001–2020. In contrast, remarkable downward trends were detected in the categories of shrub, grassland, forest, and water body, by 75.1, 39.62, 31.97, and 5.99%, respectively. The observed changes in LULC patterns (e.g. deforestation, afforestation, urbanization, and agricultural land expansion) indirectly and slowly affected the water cycle process by changing the earth cover and surface conditions. In recent years, population growth, urbanization, and unsustainable industrial development have accelerated the imbalance between the supply and demand of natural resources, such as land and water in the studied watershed. The expansion of urban land has caused an increase in impervious areas and the number of populations, resulting in greater water demand and reduced runoff depth monitored at the watershed’s outlet. Other studies in different parts of the country have also reported similar results. For example, Arragaw and Bewket (2017) obtained settlement area, and agricultural land increased while forest, water bodies, and woodland significantly decreased in the central rift valley of Ethiopia. Population dynamics is one of the important reasons for the changes in LULC across the study area. Population increase could be ascribed to the natural migration of people for the need for jobs and better living conditions, especially over the last three decades. The increase in population prompted deforestation and the clearing of shrubs along river banks for the need for additional land mainly for cultivation and settlement.

Table 11. Land use/land cover (LULC) change comparisons before and after the change year.

LULC type	1981–2000				2001–2020
	Change (ha)	% change	ha/yr	% rate	Change (ha)	% change	ha/yr	% rate
Barren	302.13	14.15	27.47	128.7	744.48	30.03	37.2	1.50
Built up	240.48	11.61	21.86	105.5	4308.93	105.75	215.4	5.29
Cropland	6854.09	3.92	623.10	35.7	12287	6.46	614.4	0.32
Forest	6154.72	69.18	559.52	628.9	−6168.04	−48.81	−308.4	−2.44
Grassland	−4902.57	−77.90	−445.69	−708.2	−510.12	−50.37	−25.5	−2.52
Shrub land	−8852.43	−36.29	−804.77	−329.9	−10432.14	−131.29	−521.6	−6.56
Water bodies	203.67	21.11	18.52	191.9	−230.22	−22.30	−11.5	−1.11

3.4. Parameter sensitivity analysis

Summarized in Table 12 are the most sensitive parameters to streamflow and sediment yield. Out of the 32 flow-related parameters, runoff is highly sensitive to the SCS runoff curve number (CN2), followed by TRNSRCH, CH_K2, CH_N2, OV_N, and other parameters. The parameters representing the surface runoff, channel properties, and groundwater flow are more sensitive in the study watershed. This indicates accurate estimation of these parameters is critically important for the hydrologic simulation in the considered watershed. Surface runoff is extremely sensitive to the SCS runoff curve number for soil moisture condition II (CN2). The CN2 is a function of watershed properties that include soil type, land use and treatment, ground surface condition, and antecedent moisture condition which has a great role in surface runoff generation and controls the fraction of water to infiltrate into the soil layer (Neitsch et al., 2011). Decreasing in the CN2 values results in decreased runoff, and an increase in soil infiltration, base flow, and groundwater recharge. In general, the parameters representing the surface runoff, and channel properties are more sensitive in the studied watershed. This indicates accurate estimation of these parameters is critically important for the hydrologic simulation of the watershed.

Table 12. Sensitivity analysis results for flow and sediment parameters.

Process	Parameter	Description	Range	t-Value	p-Value	Rank
Streamflow
	CN2^a	SCS runoff curve number II value	−0.4–0.4	15.50	0.000	1
	TRNSRCH^b	Fraction of transmission losses from main channel that enter deep aquifer	0.0–1	−12.40	0.000	2
	CH_K2^c	Effective hydraulic conductivity in main channel [mm/h]	5–130	−10.80	0.000	3
	CH_N2^c	Manning’s n value for main channel [m/m]	0–0.3	−7.78	0.002	4
	REVAPMN^c	Threshold water level in shallow aquifer for water movement [mmH2O]	0–100	6.79	0.007	5
	SOL_AWC^a	Available water capacity of the soil layer [mm/mm]	−0.8–0.3	−5.77	0.031	6
	GW_REVAP^b	Groundwater ‘revap’ coefficient	−0.02–0.2	4.28	0.043	7
Sediment
	CH_EQN	Sediment routing method	0.001–1	−16.00	0.000	1
	CH_D^b	Average depth of main channel [m]	0.0–25	−14.85	0.000	2
	CN2^a	Runoff curve number II	−0.4–0.4	−14.67	0.000	3
	USLE_K^a	USLE equation soil erodibility (K) factor	0.0–0.6	−12.63	0.009	4
	CH_COV2^c	Channel cover factor	0.0–0.4	10.76	0.019	5
	HRU_SLP^a	Average slope steepness [m/m]	0.0–0.2	−9.72	0.037	6
	SPCON^b	Linear parameter for calculating the maximum amount of sediment during channel sediment routing	0.0–0.001	−7.59	0.041	7
	SPEXP^b	Exponent parameter for calculating sediment re-entrained in channel routing	1.0–1.4	5.50	0.046	8

^aMultiplying of existing value by (1 + a given value).

^bReplacing of existing values.

^cAdding of a given value on existing value.

Similarly, changes in the ranking were observed for the most sensitive sediment-related parameters, such as the SCS runoff curve number (CN2), Effective hydraulic conductivity in the main channel (CH_K2), Maximum canopy storage (CANMX), and channel erodibility factor (CH_COV1) (Table 12). The CN2 parameter is found to be highly sensitive both during streamflow and sediment yield simulations. These results corroborated with the findings of Zelalem and Kumar (2018), which indicated that the CN2, CH_K2, and others are sensitive parameters in sediment tests in a similar watershed. These parameters were then considered for model calibration. To avoid interactions among parameters, the streamflow was calibrated first, and then, the calibration for sediment load was followed.

3.5. Performance evaluation

To effectively explain changes caused by climate and environmental factors, it is imperative to correctly calibrate and validate the underlying model. calibration and validation of the SWAT model were carried out after parameter sensitivity analysis. The years 1981–1993 and 1994–2000 were treated as the calibration and validation periods, respectively. The calibration and validation of the baseline models were performed using monthly streamflow and sediment yield records using the land use/land cover of 2000. The SWAT model performance in simulating streamflow, as well as sediment load, was evaluated based on R², NSE, and PBIAS statistics. Figures 4, 5, and Table 12 present the results of the graphical and statistical performance in simulations of these variables during the calibration (1984–1993) and validation (1994–2000) periods. Model reliability was evaluated based on R², NSE, and PBIAS statistics. In the process, the Nash–Sutcliffe coefficient (NSE) was considered as the objective function with threshold values of 0.6 and 0.5 for streamflow and sediment load calibrations, respectively. The obtained R² and NSE values during calibration and validation periods were higher than the acceptable value suggested by Moriasi et al. (2015).

Figure 4. Monthly streamflow results during calibration (1984–1993) and validation (1994–2000) periods.

Figure 5. Monthly sediment yield results during calibration (1984–1993) and validation (1994–2000) periods.

Based on Table 13, NSE = 0.82, R² = 0.83, and PBIAS = −6.5% during the calibration phase, whereas these are 0.87, 0.82, and −1.3% during validation, respectively. Moriasi et al. (2015) suggested that 0.5 < NSE ≤ 0.70 and 0.60 < R² ≤ 0.75 for streamflow simulation as satisfactory model performance values for the SWAT model. The NSE and R² values for sediment simulation were 0.66 and 0.64 during the calibration and 0.57 and 0.55 during the validation period. The result shows a satisfactory performance for sediment simulation according to Moriasi et al. (2015) suggestions. Moriasi et al. (2015) suggested 0.45 < NSE ≤ 0.70, and 0.40 < R² ≤ 0.65 for sediment yield simulation as satisfactory model performance values for the SWAT model.

Table 13. Performance statistics summary for model parameter test.

	Period	Observed mean	Simulated mean	p-Factor	r-Factor	R^2	NSE	BIAS (%)
Streamflow
Calibration	1984–1993	18.36 mm	19.57 mm	0.87	1.01	0.83	0.82	−6.5
Validation	1994–2000	22.97 mm	23.26 mm	0.88	0.92	0.82	0.87	−1.3
Sediment
Calibration	1984–1993	41706.6 t	41311.1 t	0.53	1.68	0.68	0.64	0.90
Validation	1994–2000	48973.5 t	52675.9 t	0.61	1.09	0.57	0.55	−7.60

R²: coefficient of determination; NSE: Nash–Sutclife efficiency coefficient; PBIAS: percent bias.

The runoff hydrographs of the calibration and validation periods are shown in Figure 4. The model overestimated the dry period flow and underestimated the wet season flows during calibration. Figure 5 shows the simulation results of monthly sediment load during the calibration (1984–1993) and validation (1994–2000) periods, respectively. In general, the result indicated that SWAT simulated streamflow and sediment yield of the Modjo watershed in very good and good performance ranges (Figures 4 and 5), respectively according to Moriasi et al. (2015). Therefore, the statistical and visual inspection results showed that the SWAT model is capable of simulating and reproducing the hydrological processes in the watershed. The performance SWAT model was also evaluated in a similar watershed in simulating stream flow (Zelalem & Kumar, 2018), and in the upper Awash basin (Shawul et al., 2019), in which our study watershed is situated. Their finding showed that SWAT performs very well in predicting hydrologic components at the monthly time scales. Other studies in the Ethiopian highland, such as Gashaw et al. (2019) in the Andassa watershed, Gessesse et al. (2019) in the Choke watershed, and Sime et al. (2020) in the Ketar watershed, among others that reported the performance of the SWAT model in simulating in their respective study areas.

3.6. Spatial variability of runoff and sediment under different scenarios

3.6.1. Spatial variability of runoff at the watershed scale

The calibrated and validated SWAT model was run separately for the two periods of 1981–2000 and 2001–2020, using the corresponding LULC maps of 2000 and 2020, respectively. Figures 6 and 7 display the spatial distribution of annual runoff and sediment load variations of the Modjo watershed. Regarding the baseline condition (S1), the highest runoff amount occurred in the southeastern and central parts of the watershed with a maximum value of 28.4 mm/year (Figure 6), this is due to relatively high urban and built-up areas compared with the other parts of the watershed. Whereas the lowest values occurred in the central and southern parts with a value of <2.0 mm/year. The spatial pattern of the runoff under climate variability (scenario 2) was generally consistent with that of the baseline scenario (S1) due to the insignificant changes in climate from S1 to S2. However, the maximum value contributed by some parts in the central, southern, and eastern regions of the watershed increased by about 0.103 mm/year as compared with S1. Moreover, some parts in the southern parts of the watershed, increased by 2.32–10.06 mm/year (Figure 6). The average precipitation in the changed period was 915.5 mm, which was greater than the value during the reference period (925.8 mm).

Figure 6. Spatial distributions of annual runoff under different scenarios.

Figure 7. Spatial distributions of the annual sediment load under different scenarios.

The spatial distributions of surface runoff under scenarios 3 and 4 were basically similar. However, the maximum flow value has a significant difference between the two scenarios. In scenario 3, a decrease in runoff value in a few parts of the middle regions was detected compared with scenario 1. However, an increase in surface runoff was detected in the southern and northeastern parts of the watershed during scenario 3 (Figure 6). This is attributed to land use changes, particularly the removal of natural woody shrubs and an increase in agricultural land area in parts of the sub-basins located in the southern parts of the watershed. For the combined condition (scenario 4), the minimum value was 0.58 mm/year that was observed in the lower and some parts of the central regions, and the maximum surface runoff of 21.4 mm/year contributed by the southeastern, and northwestern parts of the watershed (Figure 6). Compared with the baseline (S1) scenario, the simulated maximum flow value showed a declining trend from 28.4 to 21.4 mm/year. This may be attributable to human activities, such as exploitation and utilization of water resources as well as climate change after the change year. As both climate change and human activities intensify, surface runoff of the watershed is undergoing significant spatiotemporal variability.

3.6.2. Spatial variability of sediment at the watershed scale

Figure 7 shows the spatial patterns of sediment yields over the Modjo watershed. The spatial distribution of sediment yield has been found to be higher in the lower altitude than in the higher altitude parts of the watershed (Figure 7). This is because of the high sediment deposition in the lower altitude of the watershed. During the baseline period, the southern and most parts of the middle regions were known with the lowest sediment load distribution (Figure 7). For scenario 2, the western, parts of the middle and eastern regions have shown increased sediment distribution due to climate changes in the sub-basins of the watershed (Figure 7).

In scenario 3, an increase in the spatial distribution of sediment load was observed in the western, southeastern, and parts of middle regions compared with scenario 1. This is attributed to land use changes, particularly the removal of natural woody shrubs and the increase in agricultural lands in parts of the sub-basins located in the southeastern and middle parts. Additionally, the sediment distribution attributed to the combined impacts of climate and LULC changes decreased compared with the baseline scenario (Figure 7). The minimum and the maximum sediment load values in the middle, and eastern parts of the watershed decreased from 350,000 to 38,000 and 5,410,000 to 1,140,000 t, respectively. In general, the combined impacts of climate and LULC changes decreased the spatial distribution of sediment load in the different sub-basins of the Modjo watershed, due to land use change, river sand excavation, and other causes.

3.6.3. Spatial distributions of runoff depth and sediment load among scenarios

The spatial impacts of climate variability, LULC change, and combined changes were assessed by comparing results between S2 and S1, S3 and S1, and S4 and S1, respectively. The spatial distributions of changes in the runoff depth and sediment load among scenarios are shown in Figures 8 and 9, respectively. The spatial patterns of changes in surface runoff caused by different effects were generally consistent, which is characterized by the more significant decrease in surface runoff in the upper region than that in the middle and downstream parts of the catchment (Figure 8). The increase in surface runoff due to LULC change (S3 − S1) with the largest increase value of 0.61 mm/year, was less significant than that caused by climate change (S2 − S1), which reached the maximum increase of 2.1 mm/year. The runoff decreased particularly in the lower as well as parts of the central region of the catchment (Figure 8), which is mainly attributable to surface and river channel water abstraction for industrial development and irrigation of agricultural lands.

Figure 8. Spatial variability of annual surface runoff at the sub-basin scale among scenarios.

Figure 9. Spatial variability of annual sediment load at the sub-basin scale among scenarios.

Similar to surface runoff, sediment yield in the upper parts of the watershed decreased more significantly than that in the lower parts of the studied watershed. As can be seen in Figure 9, the spatial variability of sediment load increased in the lower part of the watershed both under climate and LULC change impacts. The maximum decrease in sediment yield caused by LULC change was 206 t/ha/year, and the decrease in sediment yield caused by climate change reached the largest value of 365 t/ha/year.

3.7. Temporal impacts of climate and human activities

3.7.1. Impacts of climate and human activities on surface runoff change

Table 14 shows the SWAT model simulated annual runoff under different scenarios. Between the periods 1981 to 2000, and 2001 to 2020, the simulated precipitation of the watershed showed a small increase of 10.3 mm/year (Table 6). The simulated runoff between scenarios 2 and 1 as well as scenarios 3 and 1 represent the individual effect of climate and LULC change, respectively. Simulations S1 and S2 showed that the streamflow increased on average by 30.3 mm/year. This indicates the individual effect of climate change to increase streamflow of the studied watershed was 9.72%. The annual average runoff simulated under S1 and S3 indicated that the streamflow increased on average by 12.4 mm/year, indicating the LULC change between 2000 and 2020 accounted for 4% of total runoff increases in the watershed. Streamflow is the interface of LULC and climate changes of a watershed, which mainly affects the temporal changes and spatial variability runoff, respectively. Quantitatively estimating the simulation values between S4 and S1 yields a 100.11 mm decrease in the annual total runoff of the watershed. The statistical results show that the combined effects of climate and LULC changes decreased 32.1% of the annual runoff in the Modjo watershed. However, the combined effects of LULC change and climate variability have far greater contributions than their individual contribution. This phenomenon indicates that the interactions between climate and LULC changes were high in the Modjo watershed.

Table 14. Annual surface runoff depth under different scenarios.

Scenario	Qo (mm/yr)	QS (mm/yr)	Comparative changes in annual runoff depth
			Climate variability		Human activity		LULC change		Climate and LULC
			mm	in %	mm	in %	mm	in %	mm	in %
S1	422.94	392.67	/	/	/	/	/	/	/	/
S2	111.9	422.97	30.30	9.72	/	/	/	/	/	/
S3	422.94	405.08	/	/	/	/	12.41	4.0	/	/
S4	111.09	292.56	/	/	/	/	/	/	−100.11	−32.10
S5	111.09	–	/	/	−181.47	−58.19	/	/	/	/

Q_o: observed runoff value; Q_S: simulated runoff value.

As presented in Table 14, the simulated annual runoff for the baseline period (392.67 mm) was lower than the observed runoff during the same period (422.94 mm), which indicates that there might be some direct human innervations in the considered watershed that impacts the hydrological components, including surface runoff. As explained earlier, the impacts of direct human activities on runoff variation can be estimated from the difference between the observed runoff after the change period (Q₂) and simulated runoff from the combined change scenario (Omer et al., 2017). The result showed that about 57.93% of the total annual runoff reduction was attributed to direct human activities (Table 14). Although the increase in barren and built-up areas (Table 11) is expected to increase the surface water of the area, however, due to the small magnitude and spatial coverage both did not significantly affect the runoff in the watershed. Therefore, the runoff decline in the Modjo watershed has been influenced more by direct human activities (such as water abstraction from surface sources and river channels for irrigation, industry, domestic supply, and others) compared to the single and combined effects of climate and LULC changes. Hence, the impact of human activities (for example, agricultural intensification, abstraction, and utilization of water resources) could be considered the main contributor to the runoff decline while implementing river management plans in the future. Because, this has a direct effect on the availability, exploitation, allocation, and utilization of water resources as well as on the physiochemical and biological processes of river ecosystems.

Although no similar kinds of literature have been found in the study area for the sake of comparison, there have been several studies elsewhere that demonstrated similar findings in their respective study areas. For example, Rientjies et al. (2011) reported a significant decreasing trend in annual discharge with losses in forest cover in the Gilgel Abbay watershed. The study conducted elsewhere by Pirnia et al. (2019) found that human activities are the dominant driving force of annual runoff changes. A study by Gao et al. (2011) in the middle reaches of the Yellow River found that human activities were the dominant driving factors for runoff reduction. Huang et al. (2016) analyzed the contribution rate of climate and human activities to runoff variations in the Wei River Basin. They showed that human activities are the main driving factors of runoff reduction. A study by Tian et al. (2022) analyzed the variation of sediment yield in response to the controlling factors. The study found a significant decrease in sediment load, which is attributable to rainfall, vegetation, and runoff. Similarly, the characteristics and attribution analysis of runoff and sediment evolution were conducted by Guo et al. (2023). In this study, the contribution of human activities to runoff change is large, which exceeds 90%.

3.7.2. Evaluating climate and human activity impacts on sediment yield

The simulated sediment values under the different scenarios are presented in Table 15. Between the baseline and changed periods, the total sediment load decreased by about 4.93 × 10⁶ t, this was due to the impacts of climate change and human activities after the abrupt change year. As determined by the SWAT model, the individual contributions of climate and LULC changes on sediment yield accounted for about 16.23 and 4.1%, of the total increase, respectively. The simulation result between S4 and S1 found that about 26.4% of the total sediment load decline was attributed to the combined impacts of climate and LULC changes. Therefore, the integrated impacts of natural and human-derived LULC change is another main factor that causes the sediment yield reduction of the study watershed. On the other hand, the effect of direct human activities on annual sediment load variation was determined from the difference between the observed value after the change period (2001–2020) and the simulated value from the combined change scenario (S4). The result showed that about 53.4% of the total annual sediment load reduction of the study watershed was attributed to direct human activities.

Table 15. Contributions of influencing factors on annual sediment yield change.

Scenario	So (t/yr)	SS (t/yr)	Relative changes in annual sediment yield
			Climate variability		Human activity		LULC change		Climate and LULC
			ton	%	ton	%	ton	%	ton	%
S1	10.1 × 10^6	9.1 × 10^6	/	/	/	/	/	/	/	/
S2	5.17 × 10^6	9.9 × 10^6	0.8 × 10^6	16.23	/	/	/	/	/	/
S3	10.1 × 10^6	9.3 × 10^6	/	/	/	/	0.2 × 10^6	4.06	/	/
S4	5.17 × 10^6	7.8 × 10^6	/	/	/	/	/	/	−1.3 × 10^6	−26.37
S5	5.17 × 10^6	–	/		−2.63 × 10^6	−53.35	/	/	/	/

S_o: observed sediment load value; S_S: simulated sediment load value.

3.8. Interaction of climate and land use changes on runoff and sediment changes

The difference between climate variability (S2) and the combined effect scenario (S4) can exacerbate the impact of LULC (Δ$Q_L^C$) caused by climate variability (Zang et al., 2013). The quantitative result estimated by comparing the simulation flow results of S2 with S4 accounted for 41.81% of the total runoff decline (Table 16). This indicates the aggravated impacts of LULC change by climate variability in affecting runoff are significant and clear, which cannot be ignored in the study area.

Table 16. Interactive and effects of LULC change climate variability on runoff and sediment load.

Runoff	Simulated runoff		Observed runoff		% of comparative change in annual runoff due to
Period	Qs (mm)	ΔQs (mm)	Qo (mm)	ΔQo (mm)	ΔQC	ΔQL	Δ$Q_L^C$	ΔQINT
S1	392.67	/	422.94	/	/	/	/	/
S4	292.56	−100.11	111.09	311.85	9.72	4.0	−41.81	47.8
	Simulated sediment		Observed sediment		% of comparative change in annual runoff due to
Period	Ss (ton)	ΔSs (mm)	So (ton)	ΔSo (ton)	ΔSC	ΔSL	Δ$S_L^C$	ΔSINT
S1	9.1 × 10^6	/	10.1 × 10^6	/	/	/	/	/
S4	7.8 × 10^6	−1.3 × 10^6	5.17 × 10^6	4.93 × 10^6	16.23	4.1	−30.43	46.65

ΔQ_s/ΔS_s: changes in annual simulated flow/sediment; ΔQ_o/ΔS_o: changes in annually observed flow/sediment load.

Complex interactions occur regarding the impacts of climate change and human activities on the hydrology and sediment yield over the Modjo watershed. As presented in Tables 14 and 15, the combined effect of LULC and climate change leads to 32.1 and 26.4% of runoff and sediment decrease which is higher than the sum of their individual impacts (i.e. 13.72 and 20.3%, respectively). This is an indication of the interaction between LULC change and clime variability in affecting watershed values, such as surface runoff and sediment. In this work, the LULC change and climate variability (LULC-climate) interactive effect (ΔQ_INT) was estimated from the difference between the individual and amplified effect of LULC change (ΔQ_L − Δ$Q_L^C$ and ΔS_L − Δ$S_L^C$). Accordingly, the LULC-climate interactive effect on surface runoff and sediment load accounted for 47.8 and 46.65% of surface runoff and sediment load changes in the period 2001–2020, respectively (Table 16). A similar approach was applied by Omer et al. (2017) who reported that the exacerbated impact of LULC by climate variability should not be ignored.

Based on the results, climate, and LULC changes manifested in increasing trends both for the runoff and sediment transport of the studied watershed during the past four decades. In comparison, the effects of climate and LULC change had a greater impact on annual sediment yield than that of surface runoff. However, the combined effects of climate and LULC changes were greater in decreasing surface runoff (32.1%) than decreasing the sediment load (26.4%) of the watershed. Similar to our findings, using a watershed in the Loess Plateau of China Zuo et al. (2016) found that the combined effects of climate and LULC changes decreased the sediment yield.

4. Conclusions

This study has contributed valuable insights into the complex interplay of climate variability, land use/ land cover (LULC) changes, and direct human activities in shaping the hydrology of the central rift valley of Ethiopia, particularly the Modjo watershed. The SWAT model served as a robust tool, calibrated and validated for the baseline period (1981–2000), incorporating statistical tests like ITA, Pettitt test, and SRSD to discern significant changes.

The findings reveal a noteworthy decreasing trend in surface runoff (8.4 mm/year) and sediment yield (4.8 × 10⁴ t/year), juxtaposed with statistically insignificant increasing trends in precipitation and evapotranspiration over the period 1981–2020. The pivotal year 2001 delineated a distinct reference phase (1981–2000) and a changed phase (2001–2020) for further analysis. Direct human activities emerged as the primary driver of change, demonstrating substantial reductions of 53.35 and 57.93% in annual sediment yield and surface runoff, respectively.

The integrated impacts of climate variability and LULC changes exhibited a synergistic effect, surpassing their individual contributions. Simulations indicated a yearly decrease of 32.1% in runoff and 26.4% in sediment load due to these combined factors. Climate variability exacerbated the influence of LULC changes, resulting in a reduction of 47.8% in annual runoff and 46.65% in sediment yield between the baseline and change periods. The intricate relationship between climate and LULC underscored significant reductions of 47.8 and 46.65% in annual runoff and sediment yield after 2001, respectively.

Furthermore, the study elucidated spatial variations in surface runoff and sediment yield, emphasizing more pronounced reductions in upper catchment areas due to the compounded impacts of climate and LULC changes. This enhanced understanding of the Modjo watershed’s hydrology highlights the nuanced dynamics among climate, land use, and human activities. The study underscores the crucial role of human-driven factors in influencing streamflow and sediment load, thereby impacting the hydrological and eco-environmental balance of the watershed. These insights are of paramount importance for policymakers and decision-makers involved in ensuring the sustainability of the watershed and informing resource management strategies.

Ethical approval

NA.

Consent to participate

No human or animal participants were involved in this research.

Consent to publish

All authors have read and approved the final manuscript for publication.

Author contributions

The authorship of this work is attributed to the collaborative efforts of the research team. The conceptualization and design of the study were a joint undertaking. The analysis and interpretation of the data were conducted collectively, with each author contributing expertise in their respective areas. The composition of the conclusions and the articulation of their implications were achieved through collaborative deliberation. All authors have reviewed and approved the final version of this manuscript, ensuring its accuracy and scholarly integrity.

Acknowledgments

The authors thank the Ethiopian Ministry of Water and Energy, USGS, the editors, and the anonymous reviewers for their invaluable comments and constructive suggestions that improve the quality of this manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

All relevant data are included in the paper.

Additional information

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Notes on contributors

Kokeb Zena

Kokeb Zena is an accomplished instructor actively pursuing his doctoral studies at the Jimma Institute of Technology, where he is under the esteemed guidance of Prof. Dr. Eng. Tamene Adugna and Dr. Eng. Fekadu Fufa. His academic journey has been marked by substantial contributions to the realms of hydraulic and water resources engineering, both in terms of teaching and scholarly endeavors. With an unwavering commitment to the pedagogical landscape, Kokeb Zena has significantly enriched the teaching and learning processes within his field of expertise. Kokeb's multifaceted engagement underscores his profound commitment to the scholarly community and the ongoing evolution of hydraulic and water resources engineering.

Tamene Adugna

Prof. Dr. Eng. Tamene Adugan is a highly regarded scholar and esteemed professor of Hydraulic and Water Resources Engineering at Jimma Institute of Technology, Jimma University in Oromia, Ethiopia. Prof. Tamene’s extensive knowledge and expertise in the fields of environment, water resources, and climate change are reflected in his numerous publications in internationally renowned peer-reviewed journals. His exceptional contributions and research have earned him global recognition as a leading authority in these areas.

Fekadu Fufa

Dr. Eng. Fekadu Fufa is a distinguished scholar and associate professor of Water and Environmental Engineering at Jimma Institute of Technology, Jimma University in Oromia, Ethiopia. With his extensive expertise in the field, Dr. Eng. Fekadu has made a significant impact through his research and contributions to a range of areas including water resources, environmental engineering, and climate change. His exceptional work has been recognized by the international community, with numerous publications featured in prestigious peer-reviewed journals.