Landslides and Gullies Interact as Sources of Lake Sediments in a Rifting Context: Insights from a Highly Degraded Mountain Environment

: Like many other lakes in the world, the interconnected Abaya and Chamo lakes in the Southern Main Ethiopian Rift are affected by rapid sediment accumulation. Although land degradation is a well-known issue in this part of the African continent, the main sediment sources, their spatial distribution and interaction in the Abaya–Chamo lakes’ basin have not yet been documented. Here, we present a systematic inventory, characterization, and spatial analysis of landslides and gullies as concentrated sediment sources, for four representative river catchments impacted by landscape rejuvenation. Using Google Earth imagery and ﬁeld surveys, we mapped with high accuracy a total of 7336 gullies and 430 landslides. Recent landslides observed during the last decade were few, small and shallow, and appear to have played a minor role in the current sediment dynamics. Large landslides are old and inactive. Although they do not contribute to the current sediment budget, they contribute indirectly to landscape dynamics by favoring the occurrence of gullies. Overall, large percentages of severe to extremely degraded areas of gully erosion are located in rejuvenated landscapes, especially at the level of the old landslides. Many active gullies are connected to the river network, as such acting as the source of sediment. Our analysis is a step towards understanding the nature and control of anthropic activities on sediment production in the region. We also highlight the importance of considering the interactions between sediment sources and the connectivity of the geomorphological system.


Introduction
Lakes are vital in preserving regional environmental and ecological functions and services [1][2][3]. However, these functions have been declining, especially in lower-middleincome and low-income countries due to rapid population growth that has accelerated land degradation, increasing sediment delivery to lakes [4,5]. Sedimentation caused by soil erosion significantly reduces the original storage capacity of many lakes. Each year throughout the world, 1-2% of the existing storage volume of reservoirs are lost due to The Abaya-Chamo lakes' basin, as part of the Main Ethiopian Rift, was formed by graben subsidence associated with tectonic extension and volcanic activity since the Miocene [28,38]. The geology of the study area is dominated by intrusive and extrusive igneous rocks, including basalts and trachyte lava flows (BTL), as well as lava flows and ignimbrite of trachytic to rhyolitic composition (ITR). Extended Holocene fluvial (fans) and lacustrine deposits cover the bottom of the rift valley [9]. Unconsolidated sedimentary material are also found along rivers. The major soil types identified in the area are Vertisols, Nitisols, Fluvisols, and Cambisols [41]. Nitisols and Cambisols are the dominant soil types [42].
The catchments are characterized by a tropical savanna climate [43] with two wet seasons end of March to mid-June, and mid-September to late November [44]. Meteorological conditions are highly dependent on elevation in the area, influenced by a combination of synoptic and local circulations that lead to high variability in the temperature, clouds, and precipitation patterns [40,45]. The average annual rainfall varies from~800 mm in the lowlands to~1000 mm in the highlands [46]. The main LULC types identified in the area are agricultural land, bare land, rangeland and forest land. LULC map of the year 2019 shows that, agriculture, being the main income-generating activity, dominates the landscape (i.e., 31%).

Data and Methods
Gully and landslide inventory maps are required to investigate the distribution of these sediment sources and the landscape characteristics associated with their occurrence. Different mapping methods, from field work to satellite images, can be used to create gully and landslide inventory maps [47][48][49][50]. Here, we used Google Earth imagery [24,51], which provides very-high (<1 m) spatial resolution images for our study area. We used satellite images acquired mainly in 2016-2019. Geomorphological field mapping was used extensively for validation.
The contribution of gullies to the sediment budget depends on their activity status, which determines how quickly their geometry evolves, the amount of sediment produced and delivered downstream and their associated impacts [52]. Gully activity was constrained from Google Earth imagery as active, partially active and inactive [52,53]. A gully is considered as active if it shows signs of active erosion, i.e., bare gully walls and bed, and/or fresh sediments deposited in the lower parts of the gully. An inactive gully is covered by dense vegetation and assumed not to produce much sediment currently. Partially active gully is an intermediate between active and inactive gully. Similarly, landslides are classified as active [54,55] when they present disturbed vegetation patterns and bare soil surfaces on an image. For the deep-seated landslides, the presence of fresh and sharply defined features such as steep headwall and flanks, and hummocky slopes are also indicators of active mass movements, independent of the vegetation disturbance. Inactive and old landsides are either fully vegetated or they can be partly (or heavily) dismantled by erosion. Finally, connection of the gullies and landslides to the river system is also an essential factor controlling the delivery of sediment to the lowlands [21,56]. A gully and/or landslide is considered to be connected to a river if its lower part is directly along a riverbed, and not connected if the sediments are deposited in fans away from the river network.
Remote sensing inventories need to be validated through ground truthing [48]. The mapped gullies and landslides were thus verified using handheld GPS along transects, over three months in 2018 and 2019. Transects were conducted in 30% of the sub-catchments over a total length of 308 km.
The gully and landslide inventories were then analyzed to understand factors controlling their occurrence and spatial distribution. Most of the gully analyses were conducted at sub-catchment scale, i.e., catchment of the first-order streams. The sub-catchments are hydrological regions, between drainage and divide lines. They were automatically derived from 1 arcsec (30 m) resolution Shuttle Radar Topography Mission digital elevation model (SRTM DEM) based on flow accumulation thresholds and minimum surface area [57,58]. A total of 277 sub-catchments ranging from one to 17 km 2 were delineated, dividing the four catchments into 55-81 geomorphological units. Gully length density (abbreviated as GLD), in km of gully length per km 2 , was computed for each subcatchment, excluding sub-catchments with no gullies. Using previous GLD classifications [59][60][61], we denote areas having GLD >5 km km −2 as extremely degraded, 2.5-5 (very severely degraded), 1-2.5 (severely degraded), 0.5-1 (moderately to severely degraded), 0.1-0.5 (moderately degraded), and <0.1 km km −2 as not degraded. The influence of potential controlling factors on the occurrence of gullies was analyzed by considering their spatial variation relative to the GLD at sub-catchment scale. The potential controlling factors considered for landslides and gullies are associated with:

•
The nature of the topography: knickpoints and landscape rejuvenation, elevation, slope angle, distance to river, and the presence of landslides; • The nature of the lithology and the regolith: obtaining information on the spatial distribution of the regolith, its nature and its depth can be challenging [35,62,63]. Here, we used information on the lithology and soil properties as proxy, together with the relative age of the topography inferred from presence of landscape rejuvenation.

•
The anthropogenic factors: LULC and its dynamics, population density, and road density.
To identify the locations of knickpoints, we utilized TopoToolbox's automatized algorithm [64] on SRTM DEM of the catchments from USGS Earth Explorer data, assuming a 1 km 2 drainage area threshold for channel initiation. Following [24], we distinguish knickpoints from background topographic variation using the maximum difference between the 10th and 90th percentile of relief in smoothed river profiles within each subcatchment. These knickpoints were used to separate the rejuvenated (i.e., zones situated at elevation between the main faults and the knickpoints that propagated upward from it) and relict (i.e., zones above the knickpoint or downward from the main faults) landscapes [24]. The SRTM data was also used for the extraction of elevation and slope angle. The lithologies of Shafe and Basso were digitized from the 1:250,000 geological map produced by the Geological Survey of Ethiopia [65], whereas lithology of Sile and Elgo were mapped using fieldwork. Soil texture was extracted from the Africa Soil Information Service (AfSIS) soil grids database with 250 m spatial resolution [66]. LULC maps for the four river catchments were produced using supervised spectral classification for the years 1986, 2003 and 2019 using data from Landsat 5, 7 and 8, respectively. Overall, the LULC classification accuracy (~90%) and kappa (~0.87) statistics showed excellent agreement between remote sensing imagery and the ground reference data [67]. Population data for the year 2020 from the 93 kebeles (the most local administrative level) that cover the study catchments were used to compute population density for each sub-catchment. Population density (# km −2 ) of each sub-catchment was computed by estimating the sub-catchment population based on an area-weighted sum of the population contribution of the different kebeles covering the sub-catchment, assuming homogeneous distribution of population, and dividing by total area (km 2 ) of the sub-catchment. The road network was mapped using Google Earth imagery and field experience, and then road length density computed. The gully number density within and outside the landslides were compared to infer on the interaction between landslides and gullies. As such, the presence of landslides was used as a potential controlling factor of the distribution of the gullies. This analysis was not carried out at a sub-catchment scale.

Inventory, Spatial Distribution and Characteristics of Gullies
Ground control points were collected for 538 gullies using GPS. Gullies identified in the field showed no major deviation from the Google Earth inventory, with an accuracy of 96%. The few false negative and false positive were attributed to gullies of very limited size (≤1 m wide), with dense vegetation cover, and/or limited contrast between the gullies and the surrounding surface on some images available in Google Earth. In total, the inventory built from Google Earth imagery across the four studied catchments contains 7336 gullies ( Figure 2). The gullies are distributed across 215 of the 277 sub-catchments. Gullies are mainly found in the upper and central parts of the Shafe, in the midlands of the Sile and scattered around the central Dembele depression in the Elgo catchment. About 50% of all mapped gullies were identified in Shafe where they are highly concentrated (Table 1). In contrast, gullies are much less frequent in the neighboring Basso, except for one deeply   Of all the gullies mapped in the study area, 56% are considered as active, this proportion being larger in the Basso (>60%, Table 1), than in the other three catchments. Of the active gullies, 20% are directly connected to the river system. Connections of the gullies to the river network are more common in the Sile (Table 1, Appendix A Figure A1).
High contrasts in GLD are observed across the sub-catchments ( Figure 2): the highest GLD computed in the highland sub-catchments of Shafe for all gullies (10.6 km km −2 ) is 1.6 times higher than the highest GLD values recorded for Sile sub-catchments. The GLD of the study area is 1.6 km km −2 , indicating that a significant proportion of the area is severely degraded. According to the GLD classification [59][60][61], 54% of our study area can be categorized as severely to extremely degraded, 41% as moderately and moderately to severely degraded, and only 5% has a low degradation level.
Even when considering only active gullies, 29% of the study area is still characterized by severe to extreme degradation and 57% by moderate and moderate to severe degradation, highlighting that gully erosion is currently a significant and widespread erosion process in the study area. With its higher fraction of active gullies, the Shafe catchment (i.e., 53% of the area classified as severely to extremely degraded considering only active gullies) stands out even more, relative to the three other catchments, when considering only active gullies. As the average length of gullies is similar between catchments, a similar spatial contrast is obtained when considering the number of gullies per km 2 instead of their length. The average number of gullies per km 2 in the area is 10.

Inventory, Spatial Distribution and Characteristics of Landslides
Ground validation was conducted for 99 landslides. The landslides identified in the field showed no major deviation from Google Earth imagery, with an accuracy of about 90% (2% false positive and 8% false negative). The completeness of the Google Earth inventory was affected by small and shallow landslides noticed in the field, but not identified on the images, as well as other features that were miscataloged as landslides, based on images. In total, the inventory built from Google Earth imagery across the four catchments contained 430 both active and inactive landslides ( Figure 3). The landslides cover 131 km 2 (~15%) of the study area. Active landslides (44%) are commonly small and shallow active processes; reactivation of large landslides being only a few. The active landslides represent 1.3% of the total area affected by slope instabilities, highlighting the contrast in size with the large inactive landslides (Table 2). In Shafe and Basso, >70% of the area is basaltic and rhyolitic lava flows (BRL), and~10% is rhyolites and trachytic lava (RTL). In Sile and Elgo, more than two thirds of the area is dominated by ignimbrite of trachytic to rhyolitic (ITR) composition. The hillshade is from the 30 m SRTM DEM [68]. Background imagery is from Google Earth (15 June 2020). Landslides active and connected to the river (%) 32 52 50 32 NB. N represent the absence of landslide in relict landscape.
Inactive landslides are commonly larger and older, their age being undetermined. These landslides are partially or fully covered by vegetation and are, to some extent, reworked by human activity. The largest landslides (>0.1 km 2 ) clearly have a long-term morphologic legacy on the landscape, and are located in the midlands of the Shafe and Basso, as well as the upper part of Sile. These are deep-seated landslides, with some earthflows also being present. About 50% of the active landslides are developed within older inactive landslides. Of the total landslides mapped, 89% have their toe in direct contact to the rivers (out of the connected landslides, 45% are active), showing the potential capacity of landslides to directly deliver sediment to the hydrological network.
Active landslides are commonly observed along gully banks ( Figure 4). These slope instabilities are usually very small and are difficult to distinguish from the images in Google Earth; larger landslides are much less frequent ( Figure 4A,B). Gully erosion is commonly present in large landslides ( Figure 4C,D, Table 3, see Section 3.3.4 for further details).
The spatial distribution of knickpoints along the channels is illustrated in Figure 3C,D. Several of these knickpoints were confirmed during field work through the presence of waterfalls or steep river stretches. Shafe, Basso and Elgo are characterized by four landscape zones with three levels of knickpoints: a relict landscape upstream from the highest knickpoints, rejuvenated phase I, rejuvenated phase II and a transition zone between phase I and phase II. The Sile catchment is characterized by relict, rejuvenated phase I and rejuvenated phase II landscapes, the latter occupying the largest area.
The spatial connection between catchment slopes, mapped faults, and the river profiles suggests that many of these knickpoints likely originate from faults ( Figure 5) as already highlighted in the Central Main Ethiopian Rift [39]. In the Shafe ( Figure 5A), two sets of knickpoints at~1700 and 1900 m elevations are likely associated with ruptures of the faults downstream of each cluster. Furthermore, a single knickpoint in the upper reach of the Shafe defines the boundary between the relict and rejuvenated landscape with higher mean slopes and more landslides downstream of the knickpoint. This knickpoint occurs within the southern branch of the upper Shafe, but we note that an associated knickpoint is missing from the northern branch, where the boundary between the relict and rejuvenated landscape is further upslope from the channel network, suggesting this transient signal has already propagated through the northern catchment. Within the mid and upper reaches of the Basso ( Figure 5B), a clearly defined cluster of knickpoints at~2400 m elevation is likely related to mapped faults downstream, also suggesting recent fault activity. Such associations are less evident for the Sile and Elgo catchments ( Figure 5C,D). This, in part, is due to the lack of topographic signatures for faulting within this area, as well as a lack of correlation between knickpoints and regional slope angle distributions. Despite this, some knickpoints are distinguishable from the river profiles. Within the Sile, a group of knickpoints that span two branches of the river network at~2000 m elevation do not coincide with slope angle increases. Similar sets of knickpoints in the Elgo at~1800 m and 2200 m elevation show no relationship with slope angle, and suggest less tectonic activity and/or at the incipient stage of rifting compared to the Basso and Shafe [39]. This interpretation is also supported by the more dendritic topology and amount of channel branching within the Sile and Elgo catchments compared to Shafe and Basso. . Gullies connected to the rivers can be observed (C,D). Most of these gullies are active, they can transport significant amount of soil to the lowlands during intense rainfalls. Landslide head scarps (blue arrows); active gullies (in red); partially active gullies (in yellow); landslides (in violet). GD in and GD out are gully number density inside and outside the landslides (# km −2 ), whereas GLD in and GLD out are gully length density inside and outside the landslides, respectively. All gullies include active, partially active and inactive gullies. N either means the landscape does not exist in the catchment and/or landslides do not host gullies in that landscape.
Most of large deep-seated landslides are located downslope of the knickpoints in the rejuvenated landscapes. Within two rejuvenated landscapes, 21% of the hillslope areas are affected by landslides. In the Shafe, a spatial relationship is noticed-a series of landslides occur downstream of the western-most knickpoint, while the sub-catchment to the north has the largest amount of small active and inactive landslides towards the top of the rejuvenated zone. The lithological control on the landslides is less evident to demonstrate here since most of the rejuvenated landscapes are characterized by one dominant group of rock types, namely BRL for the Shafe and Basso catchments and ITR composition for the Sile and Elgo catchments (Figure 3).

Controlling Factors of Gully Erosion
The potential factors controlling the spatial variation in GLD were tested for the whole gully inventory, and for active and inactive gullies separately. For facilitating the comparison between catchments, the GLD value of each sub-catchment was standardized by considering the deviation relative to the average of each catchment [69]. A negative GLD value indicates a value less than the average of each catchment, whereas a positive value indicates the opposite. This analysis was limited to sub-catchments with observed gullies.

Topographic Controls
We observed a non-linear relationship between GLD, elevation and mean slope angle of each sub-catchment ( Figure 6). GLD was highest for elevation between 1700-2600 m and for slope angle between 10-20 • in the Shafe, while no clear association was observed in the Basso. In the Sile, the density was highest at an elevation of 1400-1600 m, with one exception of a catchment at 2500 m, and no clear association with mean slope angle was observed. For the Elgo, the highest GLD was concentrated at an elevation of 1400-1900 m and for a slope angle between 15-35 • .  In all catchments, the highest GLD and percentage of severely to extremely degraded areas were in the rejuvenated landscapes (i.e., rejuvenated phase I, II and transitional), showing a positive relationship between these landscapes and GLD (Figures 3 and 7). In the Shafe, high GLDs were dominant in the rejuvenated phase I, whereas the rejuvenated landscape II was characterized by higher GLD in the Sile and Elgo. Note that in Basso, the highest GLD was observed in relict landscape at the level of one sub-catchment.

Lithologic and Soil Controls
In the Shafe, a large range in the GLD was noted for the most common lithology of basalt and rhyolitic lava (BRL) (Figure 8). Rhyolite and trachytic lava (RTL) is associated with the highest GLD median. Although the GLD values are generally lower for the Basso, we observed a broader range with some high density values for BRL lithology, including two outlier sub-catchments. In the Sile and Elgo, there was little contrast in the range of GLD among the different lithologies, although higher values were reached in ITR, the dominant lithology in the rejuvenated landscapes.
The clay ratio of soil, calculated as the percentage of sand and silt fractions over the percentage of clay, is used as an index of soil erodibility [70], with low clay-ratio soils considered to have a low erodibility and higher-ratio soils to be more erodible [70]. Except in the Sile, where the clay ratio is increasing as the GLDs decrease, no clear relation between clay ratio and density was observed in the other three catchments (Appendix A Figure A3).

Anthropogenic Controls
Over the 1986-2019 period, we observed an increase in agricultural land at the expense of forest and rangeland ( Figure 9). However, the pattern is contrasting between the different catchments, with the largest increase in agricultural land in the Sile and Elgo (~32% each), and more forests in the Shafe and Basso. In the Shafe, Basso and Sile, the largest percentages of severely to extremely degraded active gullies occur in stable agricultural land, mostly within rejuvenated landscapes (Figure 9). In the Elgo, the largest percentage is in a zone of agriculture expansion, which is also located in the rejuvenated landscape. For inactive gullies, the largest percentages of severely to extremely degraded areas occur in stable agricultural land of the Shafe and Elgo, and in the zone of rangeland expansion for Sile. A total of~340,000 people were living in the area in 2020 [71], resulting in a high average population density (310 inhabitants km −2 ) for a rural landscape. Higher population densities are observed on flatter slopes <10 • , at elevations <1300 m and >1800 m; however, agricultural activities are common in the steeper midlands of the study catchments. Overall, there is an inverse relationship between population density and GLD; although in the Shafe, several populated sub-catchments are more prone to gullying (Appendix A Figure A4).
Overall, road density and GLD show an inverse relationship for the Shafe and Sile, whereas no clear relationship is observed for the Basso and Elgo. Road density, similar to the population density, is higher in the lowlands and upper plateau where fewer gullies are observed. In the different catchments, 6-25% of the gullies are connected to the roads ( Table 1) and half of the gullies connected to roads are active; i.e., a proportion of activity similar to that observed for the whole gully inventory.

Landslides
Seventeen percent of landslides are characterized by gully incision. These landslides host more than 1200 gullies (Table 3) within the older and larger landslides. In the Basso and Sile, two of the largest landslides host over 100 gullies each. Some examples of gullies observed within the landslide are shown in Figure 4, including active and partially active gullies. As exemplified in Figure 4, we commonly observe that gullies affect the "head" of the slide where steep slopes are formed by depletion, but also the toe (convex steep slope formed by accumulation of landslide debris) or the lateral edge of the slide (flow convergence). The proportions of gullies within and outside the landslides vary among the different landscapes (Table 3). In the rejuvenated landscapes of all catchments, the proportions of active gullies (both number and length) are much higher within the landslides than outside. A similar trend was also observed for inactive gullies, except in the Shafe.

Gully and Landslide Inventories, Data Reliability and Age of the Processes
We mapped a total of 7336 gullies and 430 landslides in an area of 871 km 2 -excluding the lowlands, flat zones and depression where no gullies and landslides developed. Google Earth imagery proved to be an effective mapping tool with errors of only 4% for gullies and 10% for landslides. This level of accuracy is similar or better than that reported in other studies using Google Earth imagery, for example [24,[72][73][74].
Despite the high accuracy of the inventories, care must be taken with their reliability. Errors are mainly associated with small features or high vegetation cover preventing correct identification using Google Earth imagery. With respect to landslide activity, we identified two populations: active landslides that are relatively small and shallow, and large deep-seated landslides that are inactive. Small active landslides are slope-related processes whose signature in the landscape can disappear quickly, sometimes over a few years due to vegetation regeneration and land reclamation [75,76]. The oldest image in Google Earth used for the mapping being from 2016, we thus argue that the inventory of these small active failures provides information on the state of the landscape over a period of~10 years. Our inventory is therefore insufficient to consider the role of extreme events with decadal return periods or human-induced landscape changes, for example [23,24,77]. Nevertheless, small landslides are clearly contemporary of the current natural and humaninduced environmental conditions. The ages of large landslides are undetermined and clearly much older, potentially over 1 kyr [25], and have natural origins associated with environmental factors that have not been recorded during the last decade [78]. Although some of these large landslides show altered morphologies from erosion, we argue with confidence on the completeness of this inventory [13].
Overall, our landslide inventory provides information that complements two recent mapping efforts that only partially cover the study area [13,79]. The authors of [13] looked at a few rather large-scale and old landslides, whereas [79] investigated recent small-scale landslides. Note that the latter study provides an inventory that includes, under the name "landslide", a lot of mass-movements that are distributed along gully banks and therefore directly associated with the formation of these erosion processes. Such small-scale slope instabilities that we also observed in the field are not considered in our inventory as landslides since they are above all the direct consequence of soil erosion processes and respond to different drivers.
Our inventory of more than 7000 gullies is a reliable dataset to highlight spatial contrast in gully occurrence across the study catchments and identify potential factors controlling their spatial distribution. The gully inventory discriminates between (partially) active and inactive processes. When a gully is initiated, its extension can be very quick, reaching sometimes more than 100 m per year [17]. In some cases, gully extensions of more than one kilometre in a few decades have been reported by [80]. When a gully forms, it can remain active over several decades [76,81]. For active gullies, we therefore interpret most of them to correspond with current natural environmental conditions, although the timing of their initiation and rate of expansion remains unconstrained. Once a gully has stabilized, it can remain visible in the landscape for several thousands of years [36,82].

Interactions between River Incision, Landslides and Gullies
The role of faults in explaining the presence of large landslides has been highlighted in previous studies within the region [13,79]. Here, we further consider the seismotectonic control on the occurrence of old bedrock landslides through the analysis of the rejuvenated landscapes created by the upstream migration of fault-related knickpoints. Faults, fractures, knickpoint migration and landscape rejuvenation are parameters that act as environmental controls on bedrock landslides in other regions [31,34,83], including in the context of a rifting system such as the central portion of the western branch of the East African Rift [24,76].
Our research is one of the very few that clearly demonstrate the role of landslides as processes favoring the occurrence of gullies [20,27,76], especially over such a large region. Of course, we cannot ignore the role of landscape rejuvenation within the scheme of landslide-gully interaction. Indeed, regolith is needed for the development of gullies [33,36], with regolith thickness often thinner and less continuous on steeper slopes [37,63]. The reduction of regolith availability would be further exacerbated through the rejuvenation of the landscape [32,35] and, to some extent, the presence of bedrock landslides [34]. Based on this model, we would thus expect that GLD would be lower in rejuvenated landscapes. However, the opposite was observed, suggesting that an optimum balance between regolith availability and slope steepness in and outside the landslides is reached. In other words, the recent landscape is old enough for providing regolith conditions favorable to gullying.
The higher GLD in the rejuvenated landscapes could suggest that the hillslopes are not at threshold angles and that rates of regolith production outpace or have reached an equilibrium with the rates of erosion, despite the presence of the large bedrock landslides. In the Shafe, Sile and Elgo, we observed that the highest GLD occur overall in sub-catchments with slope values 15-25 • (Figure 6), which are typically below the threshold angles [24,31,34,83]. Nevertheless, our analysis of slope-specific GLD showed that a clear relationship with slopes was not met in the study area ( Figure 6); the GLD patterns did not show any clear tendency towards a spatial trend. Various reasons may exist for this non-linear behavior and require the disentangling of the characteristics of regolith, weathering and rock properties [34,63] that are unknown in the region. Furthermore, the role of human activities should not be ignored (see Section 4.3). Nevertheless, we argue that one reason for the presence of gullies on steep slopes in a rejuvenated context is that retreating knickpoints are developing in areas where the lithology is more prone to weathering and regolith formation. Drainage networks may develop in places where the lithology is more prone to incision [84] and, for example, the basalts and trachybasalts may display profiles with a rather higher degree of weathering [13]. The rejuvenation is also associated with the presence of active faults. More intense fracturing could not only favor bedrock river incision and landslide, but also lead to higher weathering [34,75,84].

The Role of Humans in the Occurrence of Gullies and Landslides
Due to population increase in the area, pressure on natural resources is increasing. For example, agricultural land is expanding at the expense of forest and rangeland. Large areas of severely to extremely degraded gullies were also observed in agricultural land, in the period 1986-2019. Their occurrence could be associated with the practice of annual crops and livestock breeding in the steeper midland zones, where the land remains bare during off seasons. This can result in soil crusting and sealing, which can contribute to high runoff generation [85]. Human activities might, therefore, contribute to the occurrence of active gullies, with the highest densities in agricultural land. Even though the ages of inactive and partially active gullies are undetermined, the role of human-induced environmental conditions on their initiation might be minimal; the landscape rejuvenation being a more likely factor controlling their distribution.
The occurrence of gullies and landslides along roads might be due to the slope modification during cut and fill or concentration of runoff water by the establishment of artificial drains and increased catchment area [27,86]. However, road density and GLD have an inverse relationship for the Shafe and Sile, and no clear association for the Basso and Elgo catchments. Despite this, the role of road construction should not be ignored, as new and large gullies are observed, mainly along the main road in the Elgo catchment.

Gullies and Landslides as Sediment Source
The GLD of 1.6 km km −2 in our study area highlights significant land degradation. This density is significantly higher than those reported in previous studies in Ethiopia and other parts of the world with similar climate, and LULC, for example [60,87,88] (see Table A1 in Appendix A). Although gully incision is quite severe in the area, the majority of GLD values do not reach the threshold value proposed by [89] for badland topography (i.e., gully number density >70 km −2 and GLD > 10 km km −2 ).
An average of 10 gullies per km 2 was measured in our study area, which is more than most areas mapped across the Horn of Africa [51], with the notable exception of the Afar region of north eastern Ethiopia and Eritrea where more than 50 gullies per km 2 are concentrated. When only active gullies are considered, the area is characterized by moderate and moderate to severe degradation. This suggests that a significant amount of soil may be transported from uplands to lowlands during intense rainfalls, as a high GLD has been shown to be associated to higher soil loss [90]. With 15-30% of active gullies connected to the river network, i.e., proportions that are higher than that typically found in other regions [52,53], the potential of gullies acting as direct sediment sources is high. Partially active gullies connected to a river have the potential to deliver eroded sediments directly to the river as well.
Compared to gullies, landslides are less frequent in the area in terms of number and current activity. Evidence of landslides affecting 15% of the area, however, highlights that this geomorphological process has played a significant role in shaping the landscape over geological timescales. Although in some cases, the legacy of large landslides on the sediment patterns may persist over thousands of years [91], the current contribution to the sediment budget by these landslides here is limited as they are all old features that look stabilized. This does not mean that new large landslides cannot occur, as has been reported in other regions [23,92]. Since our period of observation was short, we did not observe clusters of small landslides. In the Kivu rift, the authors of [25] observed clusters of more than 500 mostly shallow and highly mobile landslides associated with convective rainfall. When this happens, these shallow landslides can have a strong impact on sediment delivery over a period of years to decades [21]. In our study area, active landslides are limited to small scale instabilities along the rivers, representing direct sediment sources to the hydrological network. These are commonly debris avalanches that might be triggered by river erosion. In the absence of extreme landslide events associated with intense rainfall during the short period of observation, recent landslides appear to play a minor role in the current sediment dynamics, representing a limited volume of mobilized sediment. Sediment produced by gullies and landslides disconnected from the channel network is locally stored, and might later be actively reworked through surface runoff, creep and debris flows [75,93]. For the current sediment deposition problem of both the Abaya and Chamo lakes, sediment contributed by gullies and landslides may be significant, as many active gullies and landslides connect to the rivers that feed these lakes. Further research is required to determine the amount of sediment that gully erosion and landslides contribute to the two lakes, and to quantify soil loss and the temporal dynamics of gully initiation, evolution and stabilization.

Conclusions
The present study demonstrates that satellite image analysis through Google Earth coupled with targeted field validation is an effective method for deriving an exhaustive inventory of gullies and landslides at regional scale. Identifying the spatial distribution of gullies and landslides, their level of activity and connectivity with rivers is critical for identifying the areas contributing most to sediment production. The proportion of active gullies in the study area is three times higher than inactive gullies, suggesting that land degradation through gullying is an active process in the studied catchments and may have expanded in recent years, although this needs to be confirmed by landscape reconstruction over the last decades. This is not the case for landslides, where only small features along rivers seem to be actively contributing to the sediment budget, whereas large, old landslides that show very little sign of activity appear to be preferential zones for gully development. These large, older landslides are natural in origin, and mainly controlled by landscape rejuvenation.
Our analysis across four catchments of similar size highlights large contrasts in GLD, reaching up to 10.6 km km −2 in the most affected places. High GLD is mainly noted in sub-catchments with average slopes steeper than 10 • and a landscape rejuvenated by the upslope propagation of knickpoints. These gullies affect the soils formed in volcanic rock, mainly BTL and ITR, but further documentation of the lithological and pedological variations is required to further characterize their control on gully occurrence. At a regional scale, the occurrence of gullies does not seem to correlate with human factors such as population and road density. The recent construction of a tarmac road across one of the studied catchments, however, has caused a series of new and highly dynamic gullies to form. High GLDs are observed in agricultural land; however, further work is needed to better constrain gully initiation and evolution in these areas. To quantify the contribution of gully erosion and landslides in the overall sediment budget of the catchments of lakes Abaya and Chamo, an in-depth understanding of the factors controlling the distribution and level of activity of gullies, as well as quantification of the total amount of sediment delivered by each river catchment, is essential. This is also required for identifying suitable mitigation and rehabilitation measures of the gullies and landslides, and to prevent increasing soil loss and its consequent sedimentation in the Abaya and Chamo Lakes.
Our study further demonstrates the important interactions between gullies and landslides. Although large landslides do not contribute to the current sediment budget, the high GD and GLD of active gullies within landslides show that they favor the development of gullies. Although less frequent, the presence of active landslides along the gully banks also shows that gully formation locally favors the initiation of shallow landslides.
Overall, the role of the large landslides, together with landscape rejuvenation, highlights the dominance of natural, long-term, geomorphological factors on land degradation processes. Our detailed inventory and characterization of these land degradation processes in the catchments of lakes Abaya and Chamo are a step towards understanding the natural and anthropic controls of sediment production, and a contribution to the long-term objective of more sustainable land management practices preserving soils and reducing sediment delivery to the lakes.
Author Contributions: L.B. conceived and design the study, carried out the data collection, data analysis, writing-original draft preparation, and review and editing. M.K. participated in conceptualization, methodology, writing-review and editing, supervision, and project administration. O.D. participated in conceptualization, data collection and validation, methodology, writing-review and editing, and supervision. J.P. and G.G. participated in writing-review and editing, and supervision. D.O., T.E. and A.K. participated in review and editing. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
We are very grateful to VLIR-UOS for supporting the PhD scholarship of Liuelsegad Belayneh and the field research in the framework of the Inter-University Cooperation with Arba Minch University. Arba Minch University is acknowledged for its facilitation and logistic arrangement for this research work. The authors would like to acknowledge those who helped during fieldwork.

Conflicts of Interest:
The authors declare no conflict of interest.     Figure A3. Relationships between GLD with clay ratio of the soil and slope angle for each subcatchment of the four study catchments. Clay ratio is % (sand + silt)/% (clay). Slope represents the mean slope angle (degree) of each sub-catchment and n represents the number of sub-catchments for each slope-angle class. Figure A4. Relationships between population density for the year 2020 and GLD, for active gullies of each river catchment. Slope represents the mean slope angle (degree) of each sub-catchment and n represents the number of sub-catchments of each slope-angle class. Figure A5. Relationships between road density and GLD, for active gullies of each river catchment. Slope represents the mean slope angle (degree) of each sub-catchment and n represents the number of sub-catchments of each slope-angle class. Figure A4. Relationships between population density for the year 2020 and GLD, for active gullies of each river catchment. Slope represents the mean slope angle (degree) of each sub-catchment and n represents the number of sub-catchments of each slope-angle class. Figure A5. Relationships between road density and GLD, for active gullies of each river catchment. Slope represents the mean slope angle (degree) of each sub-catchment and n represents the number of sub-catchments of each slope-angle class.