Morphological dynamics of gully systems in the subhumid Ethiopian Highlands : the Debre Mawi watershed

Gully expansion in the Ethiopian Highlands dissects vital agricultural lands with the eroded materials adversely impacting downstream resources, for example as they accumulate in reservoirs. While gully expansion and rehabilitation have been more extensively researched in the semiarid region of Ethiopia, few studies have been conducted in the (sub)humid region. For that reason, we assessed the severity of gully erosion by measuring the expansion of 13 selected permanent gullies in the subhumid Debre Mawi watershed, 30 km south of Lake Tana, Ethiopia. In addition, the rate of expansion of the entire drainage network in the watershed was determined using 0.5 m resolution aerial imagery from flights in 2005 and 2013. About 0.6 Mt (or 127 t ha−1 yr−1) of soil was lost during this period due to actively expanding gullies. The net gully area in the entire watershed increased more than 4-fold from 4.5 ha in 2005 to 20.4 ha in 2013 (> 3 % of the watershed area), indicating the growing severity of gully erosion and hence land degradation in the watershed. Soil losses were caused by upslope migrating gully heads through a combination of gully head collapse and removal of the failed material by runoff. Collapse of gully banks and retreat of headcuts was most severe in locations where elevated groundwater tables saturated gully heads and banks, destabilizing the soils by decreasing the shear strength. Elevated groundwater tables were therefore the most important cause of gully expansion. Additional factors that strongly relate to bank collapse were the height of the gully head and the size of the drainage area. Soil physical properties (e.g., texture and bulk density) only had minor effects. Conservation practices that address factors controlling erosion are the most effective in protecting gully expansion. These consist of lowering water table and regrading the gully head and sidewalls to reduce the occurrence of gravity-induced mass failures. Planting suitable vegetation on the regraded gully slopes will in addition decrease the risk of bank failure by reducing pore-water pressures and reinforcing the soil. Finally, best management practices that decrease runoff from the catchment will reduce the amount of gully-related sediment loss. Published by Copernicus Publications on behalf of the European Geosciences Union. 444 A. D. Zegeye et al.: Morphological dynamics of gully systems


Introduction
Gully erosion is likely the most serious form of land degradation (Lamb et al., 2005;Daba et al., 2003). Gullies form because they are an energy-efficient way for runoff to travel from uplands to valley bottoms and permanent channels (Gyssels and Poesen, 2003;Simon et al., 2011). They can contribute more than 90% of catchment sediment yield (Tebebu et al., Zegeye et al., 2014). Gullies have also been found to damage structures and transport routes (Nyssen et al., 2004b;Valentin et al., 2005).
Gullying is a threshold-dependent process controlled by a wide range of factors (Valentin et al., 2005), such as mechanic action of water, soil properties, and drainage area. The mechanic actions of water at the gully bottom can result in a rapid mass movement on the gully sides (Lanckriet et al., 2015). When these mechanic actions at the gully head exceed the cohesive strength of soil, erosion proceeds upslope through a headward cutting gully (Munoz-Robels et al., 2010). Interactions between such processes are important as hydraulic erosion promotes bank collapse, which then modifies subsequent hydraulic erosion (Thorne, 1990;Avni, 2005). Similarly, gully formation is initiated with the occurrence of saturation-induced high pore-water pressures, which lead to seepage erosion that ultimately gives rise to gully formation (Daba et al., 2003;Tilahun et al., 2013a).
Active gully networks are predominantly found in the saturated valley-bottomlands (Tebebu et al., 2010;Steenhuis et al., 2014), and the deepest and the most spectacular gullies occur in the bottom of the watershed where in subhumid monsoonal and wetter climates, the soil becomes saturated starting around the middle of the rainy phase and then remain saturated until the end of the rainy season Tebebu et al., 2014).
Soil properties and soil types also play a role in gully formation and expansion. For example, V ertisols, heavy clay soils with a high proportion of swelling clays (IUSS Working Group WRB, 2015), form deep wide cracks from the surface downward when they dry out and are prone to the development of pipes that can collapse and thereby turn into rills or gullies (Valentin et al., 2005).This may be one of the reasons that most severe gully areas are often associated with Vertisols (Valentin et al., 2005;Tebebu et al., 2014). Similarly, in pasture bottom lands, piping often leads to development of permanent gullies (Jones, 1987;Zegeye et al., 2014). These pipes are part of gully networks and during the rainy season, the infiltrating rainfall discharges through the pipes, which increases the lower soil horizon's vulnerability to erosion. After a while, all the pipes become connected, increasing the size of the eventual gully.
The drainage area at the gully head is the most important parameter explaining linear, areal and volumetric gully headcut retreat (Vandekerckhove et al., 2003;Frankl et al., 2012Frankl et al., , 2013. Runoff-contributing drainage area can be used as a surrogate for runoff, especially if it is assumed that the rainfall amount is equal for all drainage areas and that surface conditions and land use are also very similar (Oostwoud Wijdenes and Bryan, 2001). Frankl et al. (2012) reported that among all environmental characteristics in the catchment, only the drainage area had a strong positive association with gully headcut retreat. Similarly, when data from stable and unstable sub-catchments were combined, the main factors related to gully volume were drainage area of gully and stream heads (Muñoz-Robles et al., 2010). Nyssen et al. (2002) claimed that increasing the drainage area of the gully head enhances gully development. Additionally, long-term retreat rates often show negative-exponential trends with runoffcontributing area of the gully head, which could be explained by the declining runoff-contributing area of the gully head as it moves upslope (Begin et al., 1980).
Most gully erosion studies in Ethiopia have been carried out in the semi-arid Tigray region, northern Ethiopia (e.g., Billi and Dramis, 2003;Nyssen et al., 2006;Tamene et al., 2006;Nyssen et al., 2008;Frankl et al., 2012). In this region, rehabilitation of gully erosion has been relatively successful (soil loss has been decreased to a range of 1-6 t ha -1 ) by using check dams and upland soil and water conservation (SWC) measures (Nyssen et al., 2004a(Nyssen et al., , 2006(Nyssen et al., , 2009. However, such physical structures have been ineffective to control gully erosion in the (sub) humid Ethiopian highlands, where gullies are formed in saturated Vertisol areas and where water flow often bypasses the check dams Dagnew et al., 2015). This is because the runoff mechanism in the humid region is different from that in the arid and semi-arid regions (Bayabil et al, 2010;Engida et al, 2011;Tilahun et al., 2013a,b;Steenhuis et al., 2014). Conservation structures that are effective in preventing gullying by overland flow in semi-arid regions (Nyssen et al., 2004a(Nyssen et al., , 2006, may not be effective on the saturated areas that promote gully formation and expansion in the humid Ethiopian highlands (Tebebu et al., 2010).
Thus, there is a clear need for a better understanding of gully erosion processes and factors governing these processes in the subhumid Ethiopian highlands before effective gully control or rehabilitation strategies can be designed. The objectives of this study were therefore: (1) to understand gully erosion processes in sub-humid Ethiopian highlands and identify factors controlling these processes for effective conservation practices, and (2) to provide quantitative estimates of the severity of gully erosion in the sub- Rainfall is unimodal with an average value of 1240 mm yr -1 . The majority of the annual rainfall (>70%) falls between June and the beginning of September. The dry season lasts between 8-9 months. The minimum and maximum monthly temperatures are 9°C and 26 °C, respectively.
The geomorphology of the Debre Mawi watershed is strongly controlled by geological settings. The lava dykes in the watershed affect the hydrology upslope, forcing subsurface flow to the surface resulting in saturated source areas for surface runoff (Abiy, 2009 Here, we focus on the gully processes in the above mentioned bottomlands, and study both the medium-term (8 years, from 2005 to 2013) and short term (2 years, from 2013 to 2014) gully advancement in the watershed. Specifically, we conducted a comprehensive study of the dynamics of 13 gully headcuts (hereafter referred to as G1 through G13), and factors controlling these dynamics. Gullies G1 and G9 through G13 are located in the central part of the watershed, and gullies G2 through G8 are located at the bottom flat area of the watershed (Fig. 1). All gullies except G11 and part of G6 are situated on communal grazing lands that were enclosed recently.

Gully widening and headcut retreat during the 2013 and 2014 rainy seasons
Short-term monitoring of gully headcut retreat and expansion was conducted for the 13 gullies during the 2013 and 2014 rainy seasons (June to early September). Headcut retreat and expansion are characterized by the change in gully length (L), width (W), and volume (V). The headcut retreat and gully width were measured repeatedly (about 5-8 times or mostly following rainstorms, otherwise on average every 2 weeks) using a tape meter and benchmark pins installed 5 to 10 m from the gully edges. However, a few gullies expanded more than this distance and the affected pins were reinstalled 5 to 10 m upslope of the newly formed gully edges. Gully cross-sectional geometry was surveyed by dividing each cross section into several trapezoidal segments at abrupt changes in profile, and measuring the width and depth of each segment (Fig. 2). Cross-sectional area (A) and surface area (S) were then calculated using Eqs. 1 and 2 respectively.
where i is segment index, w is width and h with overbar is average height of a segment, the subscript j is cross section index, N is number of cross sections, and Lj is channel length of the subsection between cross sections j and j+1. The total soil loss volume over the monitoring period was then obtained by taking the difference in VT after and before the 2013 and 2014 rainy seasons.
The mass of the soil loss was calculated by multiplying the soil loss volume for each subsection by the measured average bulk density (see Sect. 3) of the soils.
The relationships between the change in gully headcut retreat dimensions (lateral, headward and volumetric erosion) and the controlling factors (daily rainfall, cumulative rainfall, water table, drainage area, headcut height, and soil physical properties such as bulk density and soil texture) were analysed. Additionally, empirical relationships between the volumetric retreat (V) and the lateral (W) and longitudinal (L) retreat were developed.

Gully erosion dynamics from 2005 to 2013
To place the short-term (2-year) gully dynamics assessment described above into a broader context, we also monitored gully dynamics over a longer, 8-year period. Following the approach of Frankl et al. (2013), we determined the surficial land loss area due to gullying for the Debre Mawi watershed by digitizing all gully edges in Google Earth on aerial imagery flown on 6 Mar 2005 and 23 Mar 2013. This was not only done for the 13 gullies discussed above, but for all gullies found in the imagery ( Table   1). The horizontal resolution of the imagery was 50 cm.
Gullies were digitized by determining the location of each gully in the watershed using a hand-held GPS with a horizontal accuracy of about 3 m, after which its coordinates were imported into Google Earth to situate all gullies on the aerial imagery.
The gully edges were then digitized using Google Earth's polygon mapping tool. Finally, the digitized polygons were converted to shape-file format using ESRI's ArcGIS software, which was also used to calculate the surface area and the length of each gully. Since gully volume could not be obtained from aerial measurements, it was derived from the digitized gully surface area through a surface area to volume relationship using Eq. (4) that was obtained from the detailed measurements of the 13 gullies in 2013. Since the last aerial imagery was flown on 23 Mar 2013 (before the start of the short-term monitoring study discussed above), the surface area of the 13 gullies after the 2013 and 2014 rainy seasons was calculated using Eq. (2).

Additional measurements to determine factors controlling gully expansion
Ground water elevation is believed to be one of the most important factors for gully formation and bank instability (Tebebu et al., 2010). Therefore, ground water depths were measured near the gully heads in the 2013 rainy season with piezometers that were installed with an auger to a maximum depth of 3 m, as the water table is usually less than 3 m from the surface during the monsoon season (Tebebu et al., 2010). Intrusion of silt and sand to the piezometer was prevented by wrapping filter fabric around the 40 cm-long screened bottom end. All piezometers were capped to prevent rainwater entry and were set in concrete to prevent any physical damage. Groundwater table elevations were read using a measuring tape twice a day: in the morning and in the evening.
Daily precipitation was recorded using an automatic rain gauge installed in the northern portion of the watershed (Fig. 1). The drainage area (DA) above the gully heads was determined from topographic analysis in a geographical information system (GIS) using digital elevation models with 30 m resolution.
Soil bulk density (BD) samples from different soil layers along the bank profile of each gully (the number of layers varies from 3 to 5 depending on the gully depth) were collected with a 98-cm 3 (5 cm high) cylindrical core sampler. Soil samples were dried for 24 h at 105 °C, and bulk density was calculated by dividing the mass of the oven-dried soil by the volume of the core.
Samples were also collected from each soil layer of the gully head and sidewalls for textural analysis using the hydrometer method (Day, 1965).

Statistical analysis
The statistical measures used to evaluate the goodness of fit of the empirical relationships were the coefficient of determination (R 2 , Eq. 5), the Nash-Sutcliffe efficiency (NSE, Eq. 6) and percent bias (PBIAS, Eq. 7): where xi and yi are the predicted and observed values, respectively, and the overbar indicates their arithmetic mean value. R 2 (ranges from 0~1) describes the degree of collinearity between predicted and measured data and is sensitive to extreme values and insensitive to proportional differences. NSE is a normalized statistic that determines the relative magnitude of the residual variance (ranges -∞~1). In general, NSE>0 is considered acceptable. PBIAS is the average tendency of predicted values with respect to their observed counterparts (ranges between -100 and+100). The optimal value of PBIAS is zero, with values close to zero indicating accurate model simulation.

Gully expansion rates at the watershed scale (2005 -2013)
The expansion of the gully drainage network in the Debre Mawi watershed significantly impacted the landscape (Fig. 3, Table 1).
The total length of the gully network roughly tripled in 8 years; it increased from 8.7 km in 2005 to 26 km in 2013 (Table 1)

Expansion rates of the 13 study gullies (2005-2014)
Table 2 summarizes the increase in gully surface area of the 13 study gullies between 2005 and 2014. The combined land loss area (LLA) of the 13 gullies between 2005 and 2014 was about 3 ha (increased from 0.7 to 3.8 ha). The corresponding soil loss from these gullies was estimated at 0.18 million tons. This is equivalent to 90 t ha -1 yr -1 (ranging from 9 to 400 t ha -1 yr -1 for the individual gullies). During the last 2 years of the study (2013)(2014), the land loss area produced by the 13 gullies increased by 1710 m 2 , and about 11,200 tons of soil was eroded ranging from 18 tons for gully G12 and G13 to 6074 tons for gully G6 ( Table   2).
The linear expansion by headcut retreat of the 13 gullies during the 2013 and 2014 rainy seasons is presented in Fig. 4 and Table   3. The recorded precipitation during the 2013 (44 days of rainfall) and 2014 (31 days of rainfall) monitoring periods was 917 and 1107 mm, respectively (Fig. 4c). The headward retreat of the gully heads in 2013 ranged from 0.04 to 36 m, with a combined total of 103 m increase in gully length; whereas the total headward retreat in 2014 ranged from 0 to 7 m, with a total of 19 m (Table 3). Over these two monsoon seasons (2013)(2014), about 750 m 2 of cultivated land was consumed by the headcut retreat of the 13 gullies. This is equivalent to 44% of the total surface area damaged by the 13 gullies during 2013-2014, and about 4.5% of the total surface area since their formation to 2014. During this period, the soil loss solely due to headcut migration equaled 2885 tons, which represented 26% of the total soil loss from the 13 gullies in the same period, and about 6% starting from their formation to 2014 (Tables 2 and 3). The total 2-year soil loss due to headcut retreat of the individual gullies varied from 0 (G12) to 1260 (G5) tons (Table 3, Fig. 4). During 2013, only five of the 13 gullies (G3, G5, G6, G8 and G11) actively expanded with lateral retreat (widening) varying from 3 to 11 m and headward retreat varying from 6 to 36 m, while the other gullies remained fairly stable. The headcuts of gullies G3 and G8 migrated the most, 36 m and 24 m respectively, but only during the 2013 rainy season. However, because of the relatively shallow depth (1.4 m) and narrow width (2.6 m) of gully G3, its headcut migration contributed little (only 6.7%) to the total soil loss of the 13 study gullies (Fig. 4b). As shown in Fig. 4b, the four largest gullies (G5, G6, G8, and G11) were responsible for about 94% of the total soil loss from the 13 gullies. The relationships between the lateral and longitudinal retreat and the associated volumetric soil loss are discussed in Sect.4.

Factors controlling gully headcut retreat and their relationships with gully dimensions (2013)
The linear headcut retreat of the gullies varied by more than an order of magnitude and was not explicitly related to geographic location. This variation should therefore be explained by the differences in, among others: groundwater elevation, soil physical properties (texture, bulk density, and porosity), gully head height, and drainage area (Table 3) -which will be discussed in this section.
The correlation (RL) between the observed change in linear gully headcut retreat and the precipitation recorded during the day of the gully head retreat occurrence varied between -0.23 and 0.88. Some of the big gullies such as G5, G6 and G11 showed strong correlation (RL5 = 0.88, P = 0.009 and RL6 = 0.84, P = 0.017), whereas gullies with the greatest linear retreat (LG3 = 36 m, P = 0.55 and LG8 = 24 m, P = 0.37; Fig. 4a) showed weak relationships with the daily precipitation during the retreat event (R L3 = 0.27, P = 0.55 and RG8 = 0.34, P = 0.37). The probable reason for these fairly low correlation coefficients is that daily rainfall only slowly saturates the surrounding soil, which is partly responsible for destabilizing the gully head. Due to such slow saturation processes, the daily precipitation and gully head retreat may not correlate well. However, the largest retreat rates were observed on 13 Aug 2013 ( Fig. 4a) when the watershed was already saturated, which coincided with the maximum recorded daily rainfall (94 mm). The low correlations indicate that headcut retreat may not always occur during or immediately following precipitation except for very large storm events such as on 13 Aug 2013.
The combined linear retreat (daily and cumulative) of the 13 gully heads in 2013 was compared with three different rainfall amounts: daily rainfall recorded during the retreat events, cumulative rainfall between gully head retreat events, and the cumulative rainfall since the beginning of the rainy season. The combined linear headcut retreat showed a moderate relationship with daily rainfall (RL = 0.76, p = 0.13), but fairly strong relationship with cumulative rainfall between retreat events (R L = 0.91, p = 0.01). Note that the relationship with daily rainfall was relatively high due to the retreat that occurred on 13 Aug 2013 during the largest rainfall event. When this rainfall was excluded from the analysis, the correlation was reduced to RL = 0.035 (p = 0.95).
The combined cumulative linear retreat was highly correlated with the cumulative rainfall since the beginning of the rainy season (RL = 0.99, p = 9.4E-05). This clearly indicates that cumulative rainfall, and thus gradual wetting and saturation of the soil, is more important to headcut retreat than the wetting and surface runoff from daily rainfall or individual storms.
The drainage area for the studied gullies varied from 0.7 to 68 ha with an average value of 15.4 ha and standard deviation of 18.9 ha. Simple linear regression models between cumulative headcut retreat length (LT, in m) in 2013 and drainage area (DA, in ha), and between increase in gully volume (VT, in m 3 ) in the same year and DA were: Fitting a power law relationship between both the linear (L) and volumetric (V) gully retreat and drainage area yielded Using the above equations, the relationships between the measured L and V and predicted L and V were explained as: R 2 = 0.27, NSE = 0.11, and PBIAS = 52% for L; and R 2 = 0.69, NSE = 0.47, and PBIAS = 49% for V. Figure 5 shows the water table rose above the gully bottom for all 13 gullies during the rainy season, which indicates mostly saturated gully head and bank soils. Soil saturation by a rising water table decreases the soil shear strength (Poesen, 1993;Langendoen and Simon, 2008), and therefore destabilizes banks (Simon et al., 2000;Langendoen et al., 2013). In this study, the water table measurement was carried out twice a day (i.e., in the morning and evening). The groundwater table fluctuated between these readings, but the variation was not significant (p = 0.98). The water table decreased between morning and evening readings on average by 0.7 cm with a standard deviation of 4.0 cm, and the daily water table variation between 25 July 2013 and 18 Sep 2013 ranged from -2.5 to 13.5 cm. The power-law regression model between the minimum water table depth during the rainy season (ranging from 0.02 m at G3 to 1.5 m at G1) and the linear retreat and volumetric expansion of the 13 gullies had fairly high coefficients of determination (R 2 L= 0.65 and R 2 v = 0.93). The power-law equations are presented in Table 4. By fitting a simple linear regression, the volumetric gully expansion was well explained by the height of the gully headcut (R V 2 = 0.34, p = 0.03). In contrast, the linear retreat of the gully was not well explained by the headcut height (R 2 L = 0.06, p = 0.4). This may be because of gullies G3 and G8, which had large linear retreats but small headcut heights and therefore influenced the analysis. When both gullies are excluded from the analysis, the R 2 L for the linear relationship between the gully retreat and gully head height increased from 0.06 to 0.46 (p = 0.03) and R 2 V increased from 0.34 to 0.57 (p = 0.01). Lower gully head heights reduce the mass of possible gully head failure blocks and therefore increase the stability of the gully head. An equivalent increase in gully head stability can be obtained by regrading the gully head to a lower slope.
The soil texture for all gully banks is dominated by clay-sized particles (53 to 67% with standard deviation of 4.5%), with an average bulk density of 1.19 ± 0.04 g cm -3 ( Table 3). The gully head retreat rates were only weakly correlated to these investigated soil properties (Table 4), probably because of the measured range in soil properties. Using linear regression models, the clay content did not explain the volumetric and linear headcut retreat (R 2 V = 0.07, p = 0.35; and R 2 L= 0.05, p = 0.46), whereas the coefficient of determination between the linear and volumetric headcut retreat and bulk density were R 2 L= 0.09 (p = 0.28) and R 2 v = 0.006 (p = 0.81), respectively. The power law relations and coefficients of determination are presented in Table 4.

Effects of gully erosion on agricultural lands
It is important to know the rate of gully advancement and sediment production in order to determine the economic feasibility of From our measurements, deep and wide gullies are common in the study site destroying many crop lands (Fig. 7a, b and d). In Europe, Poesen et al. (1996) discussed that the width-depth ratio (w/d) ≤ 1 represents narrow and deeper gullies which cause relatively little crop damage and the percentage of total soil loss consisting of fertile topsoil is significantly less compared to fertile topsoil losses for wide gullies. In most cases, this may not apply in wet humid areas where ground water reaches the surface. In the Debre Mawi watershed, narrow and deeper gullies expand rapidly due to the elevated water table (Tebebu et al., 2010) that saturates gully bank soils causing bank instability; hence significant portions of the watershed area have been lost by these actively expanding gullies. The rate of soil erosion due to gully expansion in the sub-humid Debre Mawi watershed (155 t ha -1 yr -1 ) is more than four times as much as the upland erosion in the watershed reported by Tebebu et al. (2010) and Zegeye et al. (2010) . The soil loss relative to the change in gully surface area is 5860 t ha -1 yr -1 or 586 kg m -2 yr -1 , which is more than 3-fold the rate reported by Daba et al. (2003) for semi-arid eastern Ethiopia over a 30-year period. The difference can likely be explained by the fact that the gullies in the sub-humid Debre Mawi watershed are much deeper than in the semi-arid area studied by Daba et al. (2003). A parallel study by Zegeye et al. (2015) in one of the studied gullies (G6) showed that more than 90% of the sediment measured at the outlet of the gully is contributed by the gully itself and the sediment concentration of the runoff increased greatly, effectively negating any positive effect of upstream erosion control practices. Table 5 lists the goodness-of-fit parameters (Eqs. 5-7) of the power-law and linear regression relations between the change in gully dimensions, such as the volumetric headcut erosion (V), the top width retreat or lateral expansion (W) and the linear headward migration of the headcut (L). The volumetric gully expansion (V) was strongly related to top width retreat (p < 0.0001), whereas no significant relation was found between V and L (p = 0.36). Similar results were obtained using a power law relationship (Table 5). Additionally, as shown in Table 1, the relative change between 2005 and 2013 in net gully area (350%) is more than 3-fold the relative change in length (109%). This indicates that sideways or lateral gully retreat is a more important mechanism of soil loss and gully expansion in the study area than linear migration of gully headcuts. Note that the R 2 value can sometimes be misleading, as indicated by the other goodness-of-fit parameters assessed. For example, the V-L power-law relationship for all 13 gullies has an R 2 = 0.83, which indicates a good fit between gully volume and linear gully extension.

The relationship between gully headcut dimensions and their controlling factors
However, based on NSE and PBIAS (Table 5) This indicates that assessing the quality of fits between gully expansion parameters cannot solely be done based on R 2 , and that good fits also require other measures like NSE and PBIAS to be in the acceptable range.
The relations between gully dimensions were generally influenced by the lateral erosion of the large gullies (G5, G6, G8 and G11). For this study, gullies whose volumetric expansion in 2013 was greater than 90 m 3 were classified as large gullies, otherwise they were classified as small gullies. To clearly identify the impact of gully dimensions (W and L) on the direction of gully expansion, two additional analyses were conducted based on the change in gully volume (see Table 5). The volumetric erosion for small gullies is significantly related to headward retreat of gully headcuts (p = 2.1E-05), whereas lateral expansion is more important (p = 0.09) than linear headcut retreat (p = 0.56) for the volumetric expansion of large and deep gullies. For large gullies the lateral erosion is mainly influenced by gravity-induced sidewall failure as the bank soil shear strength is significantly reduced by increased pore-water pressure due to elevated groundwater table.
Many studies have indicated that drainage area is a major controlling factor of gully head retreat Frankl et al., 2012). A similar result was obtained in this study. Both the linear (Eqs. 8-9) and power-law (Eqs. 10-11) regression relationships indicated that drainage area predicted the volumetric gully erosion (VT) better than the linear headward migration (LT) of the gully headcut. This suggests that the bigger the drainage area, the bigger the gully and the greater the lateral expansion by gully sidewall erosion will be, and hence the greater the sediment production will be. were also thinner. In addition, Frankl et al. (2013) stated that the establishment of relationships like Eqs. 10-11, are necessarily region-specific and only representative for similar environmental settings (climate, topography, lithology, soil, vegetation).

Viable gully erosion control measures for the sub-humid Ethiopian highlands
Gully erosion can rapidly change landscapes, see for instance Fig. 6. Gully G6 has expanded laterally into cultivated land through erosion of the right bank (west bank), whereas lateral erosion of the left bank, located on grass land, was rather limited.
This decreased gully expansion on the grassed bank may be due to the effect of the grasses either in terms of increasing the topsoil shear strength (De Baets et al, 2008) or drying out the soil through evapotranspiration (Pollen and Simon, 2005) and thereby reducing soil saturation. Also, grasses could modify overland flow and infiltration patterns, and therefore affect subsurface drainage. Gully G11 was surrounded by cultivated land on both sides, and hence expanded laterally through erosion of both left (south) and right (north) banks. In 2012 the land adjacent to the left bank was planted with eucalyptus trees to halt erosion. In 2013, erosion of this left (south) bank was significantly reduced, and gully development then occurred through extension in northeastern direction and lateral expansion in northern direction instead (see Fig. 6). The lesson learned from these two gullies (G6 and G11) is that vegetation may reduce gully expansion by increasing soil shear strength through their roots, slowing down the storm runoff and trapping sediments which was also observed by for example Gyssels and Poesen (2003) and De Baets et al. (2006). Therefore, planting suitable species on the gully face and around the boundary may reduce or slow down bank failure and water-induced erosion especially for fairly deep gullies.
Our monitoring data also contained valuable information regarding the effectiveness of soil and water conservation measures such as soil bunds that were extensively installed across the upper portion of the catchment since 2012. The measures, aimed to reduce the development of rills and gullies in the area, were implemented on saturated Vertisol areas, but have rather led to gully initiation and development (see Fig. 7f; Steenhuis et al., 2014. These soil and water conservation measures appear to be ineffective on these locations as they cannot reduce or stop upward headcut migration of gullies downslope, which requires alternative, structural measures. Similarly, diversion waterways have been tested in the watershed to arrest gully heads, but have produced new gully branches . Our data therefore supports the findings of Dagnew et al. (2014Dagnew et al. ( , 2015, which indicate that the extensive implementation of soil and water conservation measures on periodically saturated Vertisols areas may have exacerbated, rather than mitigated gully formation and expansion. In the bottomlands of the watershed where Vertisols dominate, gully formation was severe due to alternate swelling and shrinking of expanding clays resulting in deep cracks in the dry season (Fig. 7c). As was previously observed by Frankl et al. (2012), the shrink-swell behaviour of Vertisols eventually developed into pipes (Fig. 7d) and contributed to gully development.
Though pipes contribute to gully formation, we observed in this study that they are also important to drain excess subsurface water near the gully banks, thereby potentially mitigating gully expansion. For example, soil pipes in the heads of gullies G7 and G13 drained the elevated ground water table resulting in only minor headcut retreat (Table 3). This implies that gully expansion rates could be reduced by controlling the water table and therefore the pore-water pressures in the gully head. For example, drop pipes are a common practice in the United States (Field Office Technical Guide standard 587; NRCS, 2015) to control groundwater and surface water level to halt erosion of gully heads up to 15 m in height, but can be very costly (>$50,000 each).
Moreover, most of the Debre Mawi watershed gullies are deep gullies (up to 7 m) that are therefore susceptible to gravityinduced bank collapse. Regrading the gully head and bank slopes decreases their weight and reduces the probability of bank failure .
In the semi-arid Tigray region of northern Ethiopia, Frankl et al. (2012) recommended the application of a subsurface geomembrane (vertical dam) at the gully head to increase groundwater levels and subsequently decrease soil cracking and soil piping. However, this may not be effective in the (sub) humid region of Ethiopian highlands as we have shown that elevated groundwater tables increases the rate of gully expansion (Table 4, Fig. 7f). Therefore, gully mitigation measures in the subhumid Ethiopian highlands and similar climate types should target to reduce soil water content.  (Table 1). The headcut migration of the 13 gullies during the 2013 rainy season varied significantly from 0.04 to 36 m yr -1 (Table 3 and Fig.4).

Conclusions
Understanding the controlling factors of gully head migration and lateral expansion of gullies is crucial to design appropriate gully control measures. Retreat rates depended most strongly on groundwater table elevation (Table 4). The elevated water tables saturate the soils surrounding the gullies thereby reducing the soil erosion resistance. Elevated groundwater elevations may also lead to seepage-induced erosion. Additionally, the gully head depth and the drainage area, which is representative of surface runoff magnitude, were other factors controlling gully erosion in the Debre Mawi watershed (Table 4). Therefore, conservation practices that address these parameters may be most effective.  (Table 5).
Therefore, regrading the gully head and bank slopes could reduce the occurrence of gravity-induced bank collapse for deep gullies. Studies need to be designed to evaluate the effects of controlling groundwater movement, for example by subsurface drainage, on the stability of Vertisols. Vegetation may play a vital role in reducing soil water and increasing soil shear strength.
The planting of an assemblage of suitable, native plant species (both herbaceous and woody) is being tested in the watershed.

Author contributions
A.