Dangishta (57 km²) and Robit Bata (9 km²) watersheds are found in the northern Ethiopian Highlands, within the Lake Tana basin (Figure 1). Dangishta is located within in upper reaches of the Gilgel Abay that drains to Lake Tana. The climate is a moist sub-humid with an annual temperature variation between 10°C and 25°C and a yearly average rainfall of 1745 mm during the last 10 years at the National metrological Agency (NMA) weather station 5.6 km from the watershed outlet. The rainfall distribution is a unimodal type with 90% that occurs between May to October. Robit-Bata, at the south-eastern edge of Lake Tana, drains direct to Lake Tana. It has a sub-humid climate with temperatures between 11°C and 27°C. The annual rainfall is 1,420 mm, of which 85% falls between May and September.

The topographic characteristics of both watersheds vary. Dangishta’s elevations range from 2036 to 2,400 m (Figure 1). The stream length is 34 km. More than 78% of the area is flat to undulating, with slopes less than 8% (Table 1). Compared with Dangishta, Robit Bata is at a lower altitude with elevations that range from 1828 m near the outlet to 2030 m at the upper part of the watershed. The stream length is 6.1 km. Only 55% of the area in Robit Bata is flat to undulating.

The geology of both watersheds is characterized by massive basalt bedrock overlaid by different soil types (Walker and Yimam 2018; Tilahun et al., 2020). The primary soil types in Dangishta are Regosols, Alisols, Nitosols, and Vertisols (Yimam et al., 2018). The soil layer over bedrock thickness increases from the bottom to the upper part of the watershed (David Walker et al., 2018). In Robit Bata, Vertisols are predominantly found in the lowland areas, and Eutric Nitosols in the upper part of the watershed (Abidela Hussein et al., 2019). Deep soils are located around the middle slope and decline in depth downslope and upslope. The upper half of the watershed is severely degraded with shallow soils, and the hillsides are characterized by hillslope aquifers (Tilahun et al., 2020).

Both watersheds have highly conductive surface soils. The infiltration rate is only exceeded by rainfall intensity 2.5% of the time (Yimam et al., 2018; Tilahun et al., 2020). Otherwise, the watersheds are hydrologically very different. In Dangishta, after the rain infiltrates, excess water flows laterally through faults in the volcanic rock. It appears as springs downstream in the watershed where the flow is blocked by volcanic dikes, creating a water table close to the surface during the rain phase (Figures 2A,B). Like other seepage areas, the lateral conductivity is high (Mendoza and Steenhuis 2002; Zhao et al., 2020). Eucalyptus or grassland covering approximately 32% of the watershed (Table 2) are found on these seasonally saturated soils. The volcanic dikes (Figure 2A) which are across the stream also prevent gully formation.

In Robit Bata, the degraded upslope soils with thin soil (Figure 2C) become rapidly saturated during a rainfall event resulting in saturation excess runoff (Tilahun et al., 2020) and, therefore, a much quicker response time to runoff generation than Dangishta. Volcanic dikes are absent in Robit Bata. Hence an active gully is located upstream of the flume at the watershed outlet (Figures 2D,E). In Robit Bata, which has few seasonally saturated areas, the combined acreage of grassland and eucalyptus is under 10%.

In both watersheds, crop-livestock farming is dominant. Approximately 65% in Dangishta and 74% in Robit Bata of land is cultivated during the rain phase (Table 2; Figure 3). Rainfed cereal crops are “tef” (Eragrostis tef), maize (Zea mays), finger millet (Eleusine coracana), and barley (Hordeum vulgare). Soils are prepared at the onset of the rainy season in June and July using oxen-drawn plows. Fertilizers for rainfed crops are applied in the form of DAP (di-ammonium phosphate) and blended fertilizer NPSZnB (nitrogen, phosphorus, sulfur, zinc and boron), manure, and compost in split doses in July and mid-August. According to the local Bureau of Agriculture, the average applied P fertilizers in Robit Bata watershed was 21 kg P ha⁻¹, and in Dangishta, 20 kg P ha⁻¹.

In Robit Bata, small-scale irrigation in the dry season (October to December) using shallow groundwater is practiced on 117 ha or 13% of the land (Table 2; Figure 3). The main irrigated crop is khat, a mild narcotic perineal small bush, tomatoes, green pepper, and onions are grown as well with irrigation. The applied fertilizer dose for vegetables and khat varies among farmers, but the average application is 6 kg ha⁻¹ P. Irrigated acreage in Dangishta was negligible.

Discharge, suspended sediment concentration, total and bioavailable particulate and dissolved phosphorus were measured at the outlet of both watersheds in 2017–2018. In addition, the physical and chemical soil properties and precipitation were measured.

One-way ANOVA multiple mean comparison test was used to compare the mean scores of soil properties in two watersheds for different groups of land uses (eucalyptus trees rainfed cropland and grassland). The non-parametric Kruskal-Wallis test was employed to investigate the link between sediment-associated bioavailable phosphorus (BAPP) in the storm event samples and soil available P under three different land uses. R version 4.0.2 was used for data analysis and visualization, and SigmaPlot 12.5 for statistical analysis.

The soil texture of the topsoil in both watersheds was similar and could be classified texture-wise either as clay loam or loam. Soils in Dangishta had a clay content of 32%, which was slightly greater than Robit Bata at 25% (Table 3). Soils are of volcanic origin, with soils in Dangishta being more acidic than in Robit Bata. Cropland in Robit Bata had the lowest organic matter content. In Dangishta, the organic matter was also less than for eucalyptus trees and grasslands. The soil available phosphorus was significantly greater in Dangishta than in Robit Bata (Table 3). Due to fertilizer applications, the soil available P in the cropland was significantly greater than in either the grassland or eucalyptus plantations in both watersheds.

The daily rainfall of Dangishta and Robit Bata in the 2 years is shown in Figure 4. The highest recorded daily event was 123 mm on 7/21/2017 (Robit Bata) and 60.5 mm on 8/28/2017 (Dangishta). The intensity of the rainfall ranged between 0.02 and 1.4 mm min⁻¹ in Dangishta and between 0.05 and 1.8 mm min⁻¹ in Robit Bata (Supplementary Table S1). Generally, these rainfall intensities were less than the soil infiltration rates except in saturated areas (Tilahun et al., 2020).

The PED model was used for hydrograph separation by simulating the subsurface and direct runoff at the outlet. Five parameters were calibrated: the portion of the watershed with a shallow restrictive layer, periodically saturated and located near the stream, deep soils that contained the aquifer. The aquifer in Dangishta was simulated as a first-order reservoir with a halflife of 15 days and Robit Bata as a zero-order reservoir with a travel time of 80 days from the groundwater divide to the stream. The calibrated values of the parameters are discussed below and formally shown in Table S1.1 in Section 1 of the supplementary material.

In both watersheds, approximately half of the rainfall became recharge. However, the rainfall-runoff response was distinctly different. In Dangishta, the water balance closed, and the precipitation minus the evaporation was nearly equal to the discharge at the outlet over the 2 years. As expected, the greater rainfall in Dangishta resulted in greater discharge (Figure 4; Table 4). In Robit Bata, the water balance did not close. Analysis with the PED model showed that only 80% of the precipitation could be accounted for as discharge on an annual basis after subtracting the evaporation. In Dangishta, the subsurface flow was blocked by volcanic dikes and surfaced as springs before the outlet. In Robit Bata, flow through the faults was not blocked and therefore could pass by the watershed outlet.

The contribution of surface runoff and subsurface flow (both interflow and baseflow) to the total discharge is shown in Figure 4. The hydrograph separation with the PED model showed that almost all discharge at the outlet in the Dangishta was subsurface flow. Only a small portion (less than 5%) was direct runoff (Supplementary Table S1.1). It is not uncommon that most water leaves the watershed as interflow and baseflow in volcanic tropical watersheds (Toohey et al., 2018).

Although the subsurface flow contribution in Robit Bata was around 67%, it is much lower than in Dangishta (Table 4). A part of the simulated subsurface flow, the interflow originated from approximately 10% of the area near the river and had a short travel time as described in (Supplementary Table S1.1). The remaining subsurface flow from 50% of the watershed area had a long travel time. The river flow disappeared within about 4 months after the rain phase ended (Figure 4B). The direct runoff generated on the shallow soils over bedrock (Tilahun et al., 2020) was in addition to the subsurface flow. The direct runoff per unit area was 14–18%, much greater than in Dangishta. The greatest storm discharge was also recorded in Robit Bata. Approximately one-fifth of the surficial watershed did not contribute flow to the outlet. Fewer springs were observed in Robit Bata than in Dangishta and reflected less interflow contribution.

Representative examples of concentration of the suspended sediment and P in the discharge at the Dangishta and Robit Bata outlet during individual storms are shown in Figure 5. Differences in flow paths between the catchments resulted in different sediment and P removal rates. In Dangishta, except for the two events (Figures 5D,E), sharply increasing peak concentrations of particulate-associated P fractions were absent because interflow and baseflow were dominant at the outlet. The suspended sediment and associated particulate P are negligible compared with the concentration in the surface runoff (Steenhuis et al., 2009). The event hydrographs (grey shaded) in Robit Bata are characterized by surface runoff and resulted in shorter duration, higher and sharper discharge peaks than Dangishta, mainly due to subsurface flow (Figure 5). In Robit Bata, the sediment and total P (TP and BAPP) similar to the discharge exhibited high peaks. The peaks occurred predominantly before the peak discharge (Figures 5H,J–L).

In Robit Bata, most sediment at the outlet originated from the channel section with unstable banks slipping in the channel (Figures 2D,E). The sediment peak concentrations in Robit Bata ranged from 20 g L⁻¹ or 2% (Figure 5L) to 60 g L⁻¹ or 6% (Figures 5H,J). These high sediment concentrations have been observed in watersheds with actively eroding gullies in Ethiopia (Zegeye et al., 2018; Akale et al., 2017).

The DP concentration in the discharge at the outlet behaved differently compared to the particulate P and sediment concentrations (Figure 5). Overall, DP concentrations varied little throughout the events in Dangishta. In Robit Bata, the greatest DP concentrations were generally observed during the recession flow (i.e., the falling limb). In the discussion, we will examine the differences in sediment and P concentrations and loads in the two watersheds.

Over the entire 2017–2018 study period, on average, the concentrations in the discharge at the outlet in Robit Bata compared to Dangishta were 11 times greater for sediment, four times for TP, a factor of two for BAPP, and in the same order for DP (Figure 6). Precipitation (and therefore discharge) starts in Dangishta earlier than in Robit Bata, and thus Dangishta reports a concentration in May, while in Robit Bata, the creek was still dry.

Despite the higher recorded discharge in Dangishta, sediment and associated total P loss were much lower than in Robit Bata because the sediment concentrations were much smaller (Table 5). The Dangishta on average annual sediment load at the outlet is fifteen times lower than in Robit Bata (Table 5). The TP load was three times lower.

The BAPP flux was found to be over the 2 years in Robit Bata was 1.4 times that observed in Dangishta (Table 5). However, the proportion of potential bioavailability of particulate P of all P (sum of TP and DP) was greater for Dangishta. The BAPP made up 30% of the P lost in Dangishta but only 15% in Robit Bata (Table 5). The DP contributed only 14 and 4% (quotient of the export of DP and all P lost). The DP lost in the discharge in Robit Bata was less than half of Dangishta (Table 5).

The biological available particulate phosphorus (BAPP) is an important factor in the eutrophication of the shallow Lake Tana, where the wave action suspends the bottom sediment. Consequently, it is of interest to compare available particulate P in the soil (P_ca, Table 3) in the watershed with that at the outlet. The biologically available phosphate concentration on the sediment can be calculated as the quotient of the BAPP concentration and the sediment concentration in the water (BAPP/SS).

Comparing watersheds, Figure 7 shows that despite the two times greater BAPP concentration recorded in Robit Bata than in Dangishta (Figure 6), the sediment biological available P (BAPP/SS) was lower than in Dangishta. Averaged over the 2 years, the BAPP/SS was 34 mg kg⁻¹ in Dangishta and significantly greater than the 4 mg kg⁻¹ in Robit Bata.

The available soil P concentration (P_ca, Table 3) is compared with the sediment biological available P (BAPP/SS) in Figure 7A). Available soil P was measured with the Bray method and BAPP with 0.1 M NaOH during storm events. In the USA, the soil data of Kleinman et al. (2001) shows, after removing the outliers, that the 0.1 M NaOH extracted 1.6 times as much P as the Bray method. Within the scatter of the data and taking the conversion factor account of the two extraction methods, the sediment available P concentration (BAPP/SS) in Dangishta was less on the average but in the same order as the grassland and the eucalyptus trees in the valley bottom (Figure 7). It is consistent with the subsurface flow path of the water, which does not transport the fertilized upland soils to the outlet. Hence the sediment and soil concentrations are similar (Huijun et al., 2010). Although in the same order, the sediment concentrations are less than the valley bottom soils because of phosphorus desorption in the interflow with a low P (Weihrauch and Weber 2021).

In Robit Bata, in contrast to Dangishta, the biological available particulate P in the sediment (BAPP/SS) was many times lower than soil available P, P_ca, measured for all three land uses (Figure 7B). It is caused by pickup of subsoil low in P by the stream in the gully, offsetting the sediment high in P from the upper part of the watershed. The subsoil in the highlands is still low in P, unlike in developed countries where soil eroded from the bank might contain legacy phosphorus (Dolph et al., 2019).

The difference in source and transport pathways for the sediment and the phosphorus during runoff events can be understood by studying the hysteresis relationship between discharge and concentration of SS, TP, BAPP and DP (Tananaev 2015; Thomas et al., 2016; Moatar et al., 2017; Mhiret et al., 2019; Heathwaite et al., 2021). For the two watersheds, 45 runoff events with a minimum of two observations discharge and concentrations pairs of SS, TP BAPP, DP on the rising limb and two on the falling limb were analyzed. Three types of hysteresis loops were identified: clockwise, counterclockwise, and mixed. Examples are given in Figure 8, and all loops are summarized in the excel files in the supplementary material. The clockwise loops occur when peak concentration arrives before the peak discharge. The counterclockwise arise when the concentration peak occurs after the discharge peak. Mixed loops occur when the concentration peaks are found both before and after the discharge peak (Bowes et al., 2005; Lefrançois et al., 2007; Zegeye et al., 2018; Vale and Dymond 2020).

There was a significant difference among the watersheds. In Dangishta, 18 loops were analyzed, 75% of the sediment SS-Q loops were clockwise, and 20–30% were counterclockwise (Figure 9A). Although the stream channel is stable, sediment is deposited between storms and settled in the stream. When the discharge increased, the sediment was remobilized, but since sediment was limited, the concentrations were lower during the falling limb as the sediment source has run out. Mixed loops were not observed in Dangishta. In humid-topics of Brazil, Gellis (2013) reported a similar finding that the previous storm events supply sediment to the channel. In-channel sediment became the source of sediment on the rising limb of the next event.

In Robit Bata, 27 storms were analysed. In 2017, 60% of the suspended sediment concentrations and discharge (SS-Q) loops were clockwise loops, and 40% were mixed (Figure 8B). In 2018, 65% had clockwise loops, and the remainder was mixed with one counterclockwise (Figure 9D). The clockwise loops are expected when the sediment source is close to the outlet (Bender et al., 2018). Zegeye et al. (2016) observed the same pattern of clockwise and mixed hysteresis loops at the outlet in a gully with unstable and slipping banks. Since most of the sediment originates from unconsolidated sediment of the slipping bank in the gully, any sediment signal from the upper part is obscured.

In Robit Bata, 70% of the TP-Q loops over the 2 years were clockwise and thus behaved similarly to the sediment SS-Q loops, although mixed SS-Q loops became clockwise TP-Q loops (Figure 9E). The phosphorus (TP-Q) loops were clockwise, such as on August 23, 2017 (Figure 8B). The concentration of P is greater before the peak, and the sediment concentration was both after and before. It means that TP concentration per Gram of sediment increases before the peak of the discharge. During falling limb in mixed loops, the sediment originated from gully with unstable and slipping banks and has a lower TP concentration per Gram of sediment (Figure 7). In Dangishta, the type of phosphorus loops diverged from the kind of sediment loops in several cases. Figure 9B shows in 2017, 60% of TP loops are counterclockwise, and in 2018, 50% of the loops are either counterclockwise or mixed. Less than half of the loops are clockwise. In the cases where the SS-Q loops were clockwise and the phosphorus (TP-Q) were counterclockwise. It indicates that TP concentration per Gram of sediment increases after the peak of the discharge. The sediment is initially high because it is the sediment that was deposited in the steam between storms when the DP concentration was low (Figure 5). The sediment after the peak originates from the land surface with a greater TP concentration per Gram of sediment.

The dissolved P (DP-Q) loops behaved the same in both watersheds but differed from the SS, TP and BAPP (Figure 8). In Dangishta, 65% of the dissolved phosphorus (DP-Q loops) are counterclockwise in the 2 years (Figure 9C). In Robit Bata, counterclockwise loops increased above that for sediment and TP (Figure 9F). The DP is a mobile fraction and can be transported via surface runoff, interflow, and baseflow (Evans and Davies 1998). At the onset of the rain phase, the DP followed a clockwise response predominantly similar to the other P fractions. However, as the season progressed, a perched groundwater table developed and the upland became connected with the downhill area. The contribution of interflow and baseflow to streamflow became the major portion of the discharge at the outlet (Figure 4). Higher DP contributions that entered streamflow via sub-surface processes resulted in counterclockwise loops. Fortesa et al. (2020) noted that when the pathways in the watershed become connected, the hysteresis loops become counterclockwise.

Similar to an earlier study (Guzman et al., 2013), the rainy monsoon phase was divided into intervals of cumulative effective rainfall to assess the transport of sediment, particulate and dissolved fractions of P as a function of discharge (Figure 10). Daily concentrations of sediment and P were plotted as a function of average daily discharge. Although some trends were the same (i.e., high sediment concentration at low flows at the beginning of the rainy period), the regression coefficients were very low, unlike Guzman et al. (2013) (Figures 10A,E). In the watersheds analyzed by Guzman et al. (2013), the runoff was generated by overland flow, and gully erosion was minor. In this study, the decreasing rate of rill erosion during the rain phase was responsible for the decreasing sediment concentrations at the outlet. In Dangishta, surface runoff was minimal, and the daily concentration decreased with the increase in interflow and baseflow as the rain phase progressed (Figure 4). In Robit Bata, dilution by interflow occurred for the smaller storms. Still, the maximum sediment concentration remained the same during the growing season due to the pickup of sediment from the gully (Figure 10E).

The daily total P concentration behaved similarly to the sediment concentration (Figures 10B,F). With the progression of the rain phase, the TP concentrations decreased for the same amount of discharge. However, the R² coefficients for the relationship between the total P concentration and discharge were poor similar to the sediment. Curious was that the BAPP concentration and discharge had a decent R² for cumulative effective rainfall intervals of 450–800 mm in the Robit Bata watershed (Figure 10G).

Remarkably, the dissolved P concentration and discharge were linearly related in the Robit Bata watershed with R² > 0.55 for 300 m or more cumulative effective rainfall (Figure 10G). In Dangishta, the relationship was linear for the cumulative rainfall interval between 450 and 800 mm (Figure 10D). Theoretically, this linear relationship can be obtained by assuming that the interflow concentration is small compared to the storm flow and the interflow is constant during each cumulative rainfall interval. These are reasonable assumptions since the DP concentrations are the same as the groundwater in the valley bottom, and as Figures 10D,G show, the intercept with the x-axis is small.

Appropriate landscape measures need to be implemented to reduce P export in these watersheds to improve downstream water quality, such as Lake Tana, where water hyacinths have rapidly expanded. Measures are proposed to decrease sediment loss to improve the water quality of the lake (Ayele et al., 2016; Bartley et al., 2020). However, this research shows that reducing sediment loss, actually might decrease lake water quality. The reasoning is as follows. Recent investigations in the relatively shallow Lake Tana found a positive correlation between the BAPP sediment concentration and the dissolved P concentration in the lake water (Kebedew et al., 2020). This study shows that in watersheds with gullies (e.g., Robit Bata), the BAPP concentration per unit sediment exported is much lower than watersheds without gullies (e.g., Dangishta). Thus, stopping the gully sediment from entering the lake will increase the BAPP sediment concentration of the deposited soil on the lake bottom. Hence, the mixing increases the dissolved phosphorus concentration in the lake water, promoting water hyacinth growth.

Another management practice advocated in the USA and other developed countries consists of buffer strips near the stream to remove P from the runoff water before entering the stream (Georgakakos et al., 2018; Young et al., 2019). In the (sub) humid Ethiopian highlands, the near stream area is saturated after the middle of the rain phase. They are, therefore, unsuitable for growing crops and are used for grazing animals and eucalyptus and thereby act as a natural buffer. So near stream buffer strips exist already. Conversely, drainage of these periodically saturated near stream areas for crop cultivation is expected to enhance the transport of P to the stream significantly.

Finally, a practice implemented in the Ethiopian Highlands to stop soil loss and associated particulate P are 50–60 cm deep infiltration furrows with bunds (Dagnew et al., 2017). It has not been successful in (sub) humid highlands due to the lack of maintenance (Gebregergs et al., 2021; Mhiret et al., 2021). The current trend of replacing cropland with eucalyptus trees might effectively decrease P loads since increased acreage of eucalyptus trees reduced the runoff volume by delaying the start of runoff after the dry phase (Mhiret et al., 2019).

Consequently, designing landscape interventions in the humid tropics will require further study since the results of our research lead to the conclusion that commonly implemented practices need further adaptation for the (sub)humid highlands.