Groundwater-induced seasonal slumps in gullies of the Bação Complex, Southeastern Brazil

Areas with a high density of large-scale gullies are frequent in regions of the crystalline basement of Southeastern Brazil, such as the Bação Complex. These gullies pose a high hazard for people and properties, and cause loss of agricultural land, among other impacts. Gullies originating as erosion by channelized runoff can be controlled in their first stages by ordinary containment practices. However, when erosive channels reach the groundwater, erosive processes conditioned by subsurface flows start acting, causing mass movements and their control becomes more difficult. Field monitoring shows that these mass movements occur not only in the rainy season, as expected, but also in the dry season. To understand the dynamics of mass movement in the evolution of these features, a representative unstable gully in this Complex was selected. As it is common in this region, this gully presents unstable slopes, especially due to slumps. Numerical simulations of saturated and unsaturated flow have shown that, in this region of high seasonality, the aquifer is anomalously recharged in two stages. Safety factor analysis by limit-equilibrium method indicates a condition near to instability, with slumps occurring during the rainy season, when the aquifer at the toe of the slope is recharged, and in the dry season, when the upper slope is recharged after a few months' lags due to thicker unsaturated zone. As these gullies are very unstable and difficult to access, a viable means of stabilization could be the sediment retention with gabion walls, backfilled with soil or inert materials.


Introduction
Erosion in channels may have serious impacts, particularly in tropical regions with deep weathering profiles, where they can attain large dimensions. In Brazil, these kinds of erosion are referred to as ravines i.e., those formed only by surface erosion mechanisms, and gullies (or "voçorocas") that also involve subsurface erosion (Guerra 1995). This classification has a practical basis, because ravines can be easily controlled by surface water 1 3 containment, while gullies require controlling mass movements and subsurface erosion caused by the exfiltration of groundwater, which is much more complex.
Gullies are quite common in regions that have deep and erodible soil profiles, such as those on crystalline basement or psammitic sedimentary rocks (Ruxton and Berry 1957;Okagbue and Ezechi 1988;Lal 1990;Scholten et al. 1995). Gullies usually develop rapidly, constituting a geodynamic process that may cause immense economic, social, and environmental damage. This includes loss of agricultural soil, destruction of buildings and habitat of fauna and flora, change in hydrological regime, siltation of channels and hydroelectric reservoirs, in addition to lowering of the phreatic surface, with baseflow reduction and drying of springs (Lal 1990;Bacellar 2000;Poesen et al. 2003;Ionita et al. 2015;Lima et al. 2018).
The Bação Complex, an area with crystalline basement rocks, is one of the most intensely affected areas by gullying in the world, with 225 large-sized gullies mapped in an area of 375 km 2 (Bacellar et al. 2005;Lana et al 2022), causing the elimination of approximately 1.35 × 10 7 m 2 of farmland (Lima et al. 2018). Erosion in this region is facilitated by the thick soil profile, with a very erodible C pedological horizon (henceforth called saprolite) (Bacellar 2000;Bacellar et al. 2005). The erosion in the Bação Complex is initiated by surface processes and when the water table is reached, the erosion processes induced by groundwater can be initiated. As the soil profile has a great water storage capacity, its removal by erosion reduces the aquifer recharge by about 7.5 × 10 7 m 3 , decreasing water availability in the dry season in this region (Lima et al. 2018). The hundreds of gullies in the region constitute a serious hazard for the surrounding residents of this densely populated region of high economic and historical relevance for Brazil, in addition to major environmental impacts, such as reduction in base flow and silting of channels and reservoirs.
To prevent the progression of gullies, the common practice in the region is to restore through revegetation and divert the surface waters with terracing on contour lines in the contour of the gullies. However, these practices are usually ineffective, since gullies grow mainly due to mechanisms controlled by the exfiltration of the water table at the base of slopes, with emphasis on mass movements, especially slumps (Bacellar 2000;Lima et al. 2018). Futai (2002) inferred that these slumps develop because of the loss of apparent cohesion (Abramson et al. 2001) of the exposed soils on gully slopes due to an increase in the saturation during the rainy season. Therefore, Futai (2002) suggested that these slumps could be stabilized through measures of surface protection on unstable gully slopes, using revegetation with geotextile blankets on the face of the gullies.
However, it has been found that such slumps are also frequent in the dry period, maybe due to the late recharge of the unconfined aquifer, as suggested by Drumond (2006). The groundwater recharge time lag of several months, measured in piezometers by this author, is due to the high thickness of unsaturated soil on the subvertical slopes of the gullies. Hence, the above-mentioned gully-stabilization measures that target only surface saturation of soils and runoff control are ineffective. Gullies evolve by backward erosion and tend to stabilize as the upstream contribution area (Hayas et al. 2019) is progressively reduced. One evidence of the stabilization of these gullies is the progressive colonization of pioneer vegetal species in the bottom and toe slope of gullies (Farias 1992). This author demonstrated that the vegetal succession in these gullies only develops when the groundwater is at a certain depth to allow the root development of the plants. The extreme runoff in the rainy season and the aquifer rise after the recharge can explain why several attempts to stabilize gullies with revegetation have failed in the region. Therefore, the main objective of this study was to investigate the influence of the seasonal water flow regime in the saturated and unsaturated zones of the soil profile on 1 3 the degree of stability of the slopes of these gullies. A representative gully in the Bação Complex was selected and detailed data were collected to carry out this analysis, which included numerical simulation of flow by finite element method, coupled with the analysis of the degree of stability to slumping by limit-equilibrium method. By identifying the seasonal mechanisms of the evolution of gullies, this study can assist projects of stabilization of these features and in the control of sediment transport.

Study area
The Bação Complex is the crystalline basement of the Ferriferous Quadrangle ("Quadrilátero Ferrífero"), which is one of the most important minerals provinces worldwide, with large reserves, particularly of iron and gold, in the supracrustal rocks of the Rio das Velhas and Minas supergroups (Door 1969;Alkmim and Marshak 1998). This complex of Archean age covers approximately 385 km 2 , and outcrops in a structural window between the supracrustals (Fig. 1). It is composed predominantly of banded, variably migmatized gneisses from tonalitic-trondhjemitic-granodioritic composition, intruded by granitoids, and subordinately, by mafic and ultramafic rocks (Gomes 1986;Noce et al. 1998;Lana et al. 2013). A polyconvex (mar de morros) or slightly aligned hilly relief predominates the area, with a dendritic or angular dendritic drainage (Bacellar 2000). The weathering profile can attain a thickness of up to 50 m in gently sloping areas, where Ferralsols develop. Here, the saprolite thickness is tens of meters, the pedological Bw horizon (diagnostic horizon of Ferralsols) is 2-10 m and that of A horizon (major eroded) is only a few decimeters (Bacellar 2000). The saprolite has a silt loam texture according to the soil texture classification diagram and is characterized by low hydraulic conductivities (3.7 × 10 -6 cm/s < K < 6.6 × 10 -5 cm/s) and high erodibility (Parzanese 1991;Bacellar 2000;Morais 2003;Bacellar et al. 2005), even for piping erosion, as demonstrated by pin-hole tests (Morais et al. 2004). Although the Bw horizon has a predominantly clay texture, its hydraulic conductivity is high (1.0 × 10 -4 cm/s < K < 9.2 × 10 -4 cm/s ) and has low erodibility due to aggregation of fine particles by iron oxides and hydroxides, providing a granular structure to the soil, typical of highly evolved tropical soils (Lal 1990, Brito Galvão and Shulze 1996, Futai et al. 2005. Cambisols develop in the steepest areas and Fluvisols and Histosols in valley floors (Bacellar 2000).
The climate of the Bação Complex is classified as CWa (Koppen classification), with a rainy season between October and March, and a dry season between April and September (Sá Júnior 2009). The average annual precipitation is 1348 mm, within the range 1024-1744 mm (Bacellar 2000), which exceeds the average annual evapotranspiration of 963 mm (Amorim et al. 1999), resulting in an effluent drainage system. The vegetation is classified as semi-deciduous seasonal forest and is considered a transition zone between the "Cerrado" (Brazilian savannah) and the Atlantic forest domains (Santos et al. 2002).
The flatter areas of the Bação Complex have a high concentration of gullies of large dimensions (Fig. 1b) that are up to 400-500 m long, 150 m wide, and 50 m deep (Bacellar et al. 2005). This study was conducted on one of the largest gullies in the region, located in an area known as the Holanda monitoring station (Bacellar et al. 2005) on the southern part of the complex (Fig. 1). Like most other gullies in the complex, this gully lies in an urban expansion in the Santo Antônio do Leite County posing a potential risk to the occupation (Sobreira 2000).
These gullies begin as small rectilinear rills due to surface runoff. Although the Bw horizon is resistant to erosion, the concentration of runoff water, exacerbated by various human activities since the end of the seventeenth century (Bacellar 2000), can erode it to form ravines. When the phreatic surface normally situated at the saprolite base is reached, the gully is extended by the centripetal force of percolating groundwater due to the increase in intake area with the progression of backward erosion (Terzaghi et al, 1996) and gradually acquires a more rounded form (Sobreira, 2000). The erosion-resistant Bw horizon forms a protective cover for the saprolite, resulting in very steep gully slopes that are often subvertical in the middle and upper thirds (Bacellar 2000).
Water infiltration into the soil is high during wet months (Morais 2007) due to the region's climate. However, due to the high thickness of the unsaturated zone consisting largely of saprolite, the recharge of unconfined aquifers involves a high lag time, as can be seen by the rise (up to 3 m) of the groundwater exfiltration zone on the slopes of gullies (Drumond and Bacellar 2006). These authors used Casagrande piezometers installed in a gully nearby (Mangue Seco gully), where the water table was approximately 20 m deep (Fig. 2a) and identified a lag time of ~ 5 months between the onset of rains (October) and the beginning of the phreatic surface recharge (March). This high lag time was also apparent in the hydrograph of the drainage channel of this gully outlet (Fig. 2b) and was also recognizable in annual hydrographs of drainage channels in some catchments of this region (Costa 2005; Costa and Bacellar 2007) (Fig. 2c). In fact, in catchments with gullies, the water table tended to be lowered by these erosive features that functioned as drains, thereby delaying the recharge and consequent base flow in comparison to gully-free catchments (Fig. 2c).
Hence, gullies initially evolve by surface erosion processes, but when the phreatic surface is reached, various groundwater-induced processes become more important, including several types of mass movements, particularly slumps, which involve the entire soil profile and have the bedrock as a basis (Drumond and Bacellar 2006;Lima et al. 2018).

Fig. 2
Hydrological behavior of the region. a Variation in the water level (hydraulic head) of three piezometers located upstream of Mangue Seco gully. b Hydrograph representing the flow over a hydrological year of the drainage channel at the watershed mouth of Mangue Seco gully (scale on the right). The base flow was separated from total flow. Note that the base flow rises after the rainy season (adapted from Lima et al 2018). c Hydrographs representing the specific daily flow (flow/area-Q/A) of daily flows over a hydrological year of two catchments in the vicinity of the area. The specific base flow allows comparing catchments with different areas. Note that the specific base flow increases at the end of the rainy season and lasts until April/May in the dry season (adapted from Costa and Bacellar 2007) Hence, controlling surface flow does not prevent the progression of gullies in the region (Sobreira 2000;Bacellar et al 2005). Futai (2002) suggested that the steep gully slopes become unstable by the loss of apparent cohesion caused by saturation during rainy periods. However, field observations and laboratory test results showed that soil saturation only occurs at the toe of the slopes, with the exfiltration of groundwater (Bacellar 2000).
Drumond and Bacellar (2006) also point out that while large slumps prevail during the rainy season, small mass movements, such as creep, small slumps and flow-slides, develop at the base of gully slopes in the dry season probably due to the delayed ascent of the phreatic surface. According to these authors, such instabilities would explain the large volume of soil transported in suspension through internal drainage channels during the first low-intensity summer rains, indicating soil mobilization on gully slopes in the dry season. This means that during the dry season the slopes are progressively undercut, and the soil dislodged inside the gully is eliminated with the first rains. This is a cause of concern since the seasonal pulsing of eroded soil from gullies into the streams may have strong effect on the quality of the abstracted water as well as on the aquatic ecosystems of the region (Wantzen 2006).
Standard penetration test (SPT) carried out in boreholes in the Holanda station area shows that the NSPT resistance (Terzaghi et al 1996) is lower (10-20 blows) at the top of the saprolite and that it gradually grows (20-30 blows) and tends to stabilize below 15-20 m deep (Bacellar 2000). If the overburden correction factor (Liao & Whitman 1986) is applied to these NSPT data, the results indicate an approximately constant penetration resistance for the saprolite throughout the profile. With the exception of the first upper meters, where the saprolite has a slightly finer texture and lower mechanical resistance due to weathering (Futai 2002), the rest of the thick profile presents relative uniformity in texture and physical indices (Bacellar 2000).
The shear strength parameters of the Bw horizon and saprolite of the Holanda station were obtained from the detailed geomechanical characterization carried out by Futai (2002). As A horizon is shallow and is usually eliminated by sheet erosion (Bacellar 2000) near gullies, they were not analyzed by this author. Based on the large volume of geotechnical characterization data of the area, Futai (2002) considered that the results of the shear tests obtained with soils at a depth of 5 m would adequately represent the mechanical strength of most of the thick saprolite profile on the slopes of the gullies. On the other hand, the soils of the Bw horizon are characteristically very homogeneous and with more uniform resistance parameters. Performing triaxial compression tests with servo-controlled suction on soils in the gully areas, Futai (2002) found a cohesion value of approximately 38 kPa, with an angle of friction of 26.4° for the saprolite and cohesion of 8 kPa and friction angle of 27.7º for the Bw horizon.

Methods
Gullies of the Bação Complex have experienced similar genesis and evolutionary process (Bacellar et al 2005). For this study, it was selected the area known as the Holanda monitoring station, which has a large repository of geological, geotechnical, and hydrogeological data (Parzanese 1991;Bacellar 2000;Silva 2000;Sobreira 2000;Santos et al 2002;Futai 2002). There are two large gullies at this station (Fig. 3a), with slopes affected by mass movements, especially slumps. An aerophotogrammetric map of the area was carried showing that erosion in the last 14 years has advanced a lot upstream (north and northwest), but that it has also increased against the topographic slope (south and southwest). Note that the unstable gully slopes are those devoid of vegetation out on 05/20/2020 with a Phantom 3 advanced drone, at a height of 80 m and with a resolution of 1.8 cm/px. For this study, the focus was on the eastern gully (Gully 1, Fig. 3a).
Despite the large volume of existing soil geotechnical data, basic characterization tests were performed, including soil texture and plasticity, following Brazilian (NBR) and international standards (ASTM). For this, deformed and intact samples were collected from trenches and gully slopes from Bw and C (saprolite) horizons, since A horizon was largely eroded. Hydraulic conductivity was measured in undeformed samples using the variablehead permeability test. The variation in soil suction with soil moisture in the drying path was assessed with the filter paper technique of D5298-16 ASTM (2016) standard. Laboratory data for the soil-water characteristic curve were adjusted following Gitirana and Fredlund (2004), a method suitable for tropical soils, which often exhibit porous space with a bimodal distribution. The hydraulic conductivity function was then adjusted with the method proposed by Fredlund et al. (1994), using Geostudio (2020a) software SEEP/W module.
The shear strength parameters of the soils were obtained from Futai (2002). For the analysis of the shear strength the ϕ b linearity in the extended Mohr-Coulomb failure envelope of unsaturated soils, as defined by Fredlund et al. (1978), was assumed.
The climate data were obtained from INMET (Brazilian weather center) weather stations. Precipitation was measured with a pluviometer located proximately and potential evapotranspiration determined by the Thornthwaite (1948) method with photoperiod, average temperature, and relative air humidity data. The actual evaporation was calculated by the method of Wilson et al. (1997) available in the SEEP/W module (Geostudio 2020a).
For stability analysis, a representative section (Fig. 3b) was selected according to historical data, in addition to in situ observations, according to the advance of the slump gullying process. The NW slope of Gully 1 was highly unstable, as reported previously (Bacellar 2000).
The geological model of the section was based on the description of gullies slopes and previous data from auger and percussion (with SPT) borings and geophysical surveys of electrical resistivity and GPR (Ground Penetration Radar) (Bacellar 2000). Some vertical electrical soundings (VES) were also conducted, since this geophysical method is recommended (Telford et al, 1990) to better identify the position of the water table and the top of the bedrock. The obtained VES were inverted using the IPI2Win-1D software (2000) to obtain the geoelectrical layers that represent the soil profile.
This representative section was used for numerical simulations of flow analysis by finite elements in the SEEP/W module of the Geostudio (2020a). The analysis was initiated with a steady-state flow simulation for establishing the initial flow conditions (Anderson and Woessner 2002). Subsequently, the conditions of transient flow in the hydrological year between 10/01/18 and 09/30/19 were simulated.
The steady-state and transient flow numerical simulations were validated with the water table position assessed using electro-resistivity and by comparing with seasonal piezometric data compiled by Drumond (2006) and simulated by Lima (2016) in Mangue Seco gully, which presents remarkably similar characteristics to the studied gully.
The slope stability analysis was conducted on the SLOPE/W module on Geostudio (2020b), using the limit-equilibrium method, with Morgenstern and Price (1965) technique. The stability analyses were coupled with flow analyses to assess the seasonal variations of the safety factor. The import of flow simulation results, as proposed here, is preferable to the simple addition of the piezometric line in slope stability analyses, as usually done. This is because the numerical simulation of unsaturated and saturated flow makes it possible to better analyze the influence of hydraulic gradients, often with upward components, on slope stability (Puy and Ramesch 2006). The results of safety factors over the hydrological year were smoothed using the moving average method to better observe trends.

Results and Discussion
The climate data for the hydrological year 2018-19 illustrate the local climate seasonality, with rainfall concentrated between October and March (972 mm), and a dry season between April and September (199.1 mm). The annual precipitation in this hydrological year was 1171.1 mm, below the regional average (1348 mm), but within the measured range, between 1024 mm and 1744 mm (Bacellar 2000). This drought was a consequence of an abnormally dry January, a phenomenon regionally called "veranico" (dry spell). The potential evapotranspiration estimated using the Thornthwaite (1948) method was 1099.5 mm, with higher values in the dry season, when actual evapotranspiration is lower. The cumulative water balance, based on the actual evaporation (Wilson et al. 1997), showed the climatic seasonality, with an excess of water during rainy months and deficit in the dry months.
The two gullies in the Holanda station are on the side of a hill gently sloping to the south and southeast (Fig. 3a). On the gully slopes, it is possible to see that the thickness of the saprolite may reach tens of meters in size, while Bw horizon is 2-10 m. A Horizon is almost completely eroded.
The rills on the gully slopes occur only in the saprolite, confirming its greater erodibility. As a result of the lower erodibility of the Bw horizon, the slopes are subvertical in the upper third part. Gully 1 has a maximum height of ~ 50 m in the most critical section, an average width of 250 m, and length over 310 m. Bedrock occurs at the bottom of the gully and the phreatic surface fluctuates (~ 1-3 m) at the base of the saprolite in response to recharge events. In the satellite image of 05/31/19, in the dry season and after 15 days with no precipitation, there is a large effluent flow at the base of the slopes of the two gullies feeding the main channel (Fig. 3b). This behavior is confirmed by hydrographs of daily flows of catchments in the area (Fig. 2c).
The intact vegetation indicates that some slopes are more stable, while bare slopes are unstable, such as those located in the west-northwest section of both gullies of Holanda station (Fig. 3b). Mass movements, especially slumps, and occasional exfiltration of groundwater are frequently observed in these active slopes in the field. A comparison of two remote sensing images (August 2006 and June 2020, Fig. 3b) shows that both gullies have expanded not only in the direction of the topographic slope but also in other directions due to the influence of subsurface processes induced by the centripetal force of groundwater percolation. In other words, although the topographical slope is northward, the erosion features are expanding even to the south, mainly by the action of slumps. It is important to note this behavior is common in the region, showing the high relevance of mass movements in the development of these gullies.
The macroscopic characterization of the soils of the gullies slopes in the Holanda station area confirmed that the Bw horizon and saprolite present remarkably homogeneous texture. The saprolite is mainly classified as a silt loam according to the soil texture classification diagram that justifies its high erodibility. The reddish-brown or yellowish Bw horizon is classified as a clay but it has a granular structure due to the aggregation of the particles by clay and iron oxyhydroxides. Because of this aggregation that is typical of tropical soils (Lal 1990, Lacerda 2010) the Bw horizon presents higher porosity and hydraulic conductivity ( Table 1). The results obtained in this study (Table 1)  The soil-water characteristic curves and the hydraulic conductivity functions of both soil horizons are represented in Fig. 4 and the adjustment parameters are represented in Table 2. The obtained results were in agreement with those by Futai (2002) in the same area. The Bw horizon exhibits a bimodal behavior, typical of well-evolved pedological horizons (Carvalho and Leroueil 2000). The saprolite also presents bimodal behavior, unusual for this type of soil, but already identified in the area by Futai and Almeida (2005). If the minimum limit to drain macropores is set as 10 kPa suction (Marshall 1959), the Bw horizon, unlike saprolite, is distinguished by the large volume of macropores. This favors infiltration and decreases the surface flow, consequently reducing erosivity, in contrast to the behavior of the saprolite.  Systematic monitoring of several gullies in the region during the last 20 years (Bacellar 2000;Futai 2002;Bacellar et al. 2005;Futai et al. 2005;Drumond and Bacellar 2006;Costa and Bacellar 2007;Lima et al, 2018;Lana et al, 2022) has shown that surface erosion processes are important in the early stages of erosion in channels; however, when the water table is reached, subsurface processes become more prominent. The influence of surface flow becomes less relevant as measures to contain floods such as terraces on contour lines, also adopted upstream of the gullies of the Holanda station (Fig. 3b), are often unsuccessfully deployed in the region.
The gullies in the study area exhibit the following mechanisms of erosion (Bacellar 2000;Drumond and Bacellar 2006): rotational slides (slumps) mobilizing the whole slope (Fig. 5a); falling soil blocks from Bw horizon, which break by tension (tension joints) over the saprolite that is eroded by plunge-pool erosion; wash, rill (Fig. 5b, c) and splash erosion of the highly erodible saprolite exposed on slopes; boiling and piping erosion and moisture-driven seasonal creep or small flow slides of saprolite in the area of exfiltration of the phreatic surface at the toe of the slopes (Fig. 5b and c).
As already verified by Drumond and Bacellar (2006), the creep and flow slide processes in these gullies intensify in the dry season, when the phreatic surface rises due to delayed recharge. Consequently, the gully walls are vulnerable to the action of these two processes, facilitating subsequent development of large slumps that affects the entire slope. Therefore, the undercutting of the slope base during the dry period is fundamental to the continuous advancement of gullies.
To better understand the role of water dynamics in the seasonal behavior of slope stability, a 2D equilibrium-limit analysis was performed coupled with a numerical model of groundwater flow in a geologically representative section of the unstable slope in the northwest sector of gully 1 (Fig. 3b). Data from previous surveys and geophysical mapping, such as electrical resistivity and GPR (Bacellar 2000) were used to build the conceptual model. To resolve ambiguities, particularly regarding the actual position of the groundwater and the top of the bedrock, additional geoelectric data were acquired (vertical electrical soundings-VES).
The joint processing of these data allowed to build a conceptual model of the soil horizon distribution, the top of the bedrock, and the position of the phreatic surface, from which the analysis section was built (Fig. 6a). The natural (not eroded) slope has an average slope of 12.2%. At the gully face, the eroded slope has sub-vertical declivities towards the top (maintained by the more resistant Bw horizon) that softens towards the inner gully channel, with a permanent drainage. In the direction of the topographic divider, it is possible to notice that the overall thickness of the weathering cover decreases, from 6 to 2 m for Bw horizon and from 38 to 22 m for saprolite. The water table is on average 18 m deep in the divider and 26 m under the erosion scar ridge, with an average hydraulic gradient 0.1-0.4% under the natural slope, increasing in the eroded portion of the slope caused by the lowering of the aquifer by the gully.
The hydrodynamic parameters for simulating the saturated and unsaturated flows were obtained in hydraulic conductivity tests (Table 2), soil-water characteristic curves and hydraulic conductivity function (Fig. 4). The boundary conditions specified for the finite element simulation (Anderson and Woessner 2002) were: specified flux (Type 2), such as impermeable border at the base (contact with the gneissic rock mass) and on the right (hydraulic boundary on the drainage divider) and at the top (specified flux by weather conditions); specified hydraulic head (Type 1) at the left border, where occurs the drainage channel inside the gully.
Results of the subsurface flow analysis under transient regime along the hydrological year revealed two distinct patterns: one on the eroded slope and the other on the natural slope (Fig. 6a). The phreatic surface in the corresponding hydrological year varied more on the eroded slope (2.97 m) than on the natural slope (maximum variation 1.24 m) and the hydraulic gradients for the eroded slope have lower values during the wet season (0.66) than at the peak of the dry season (0.77). The decrease in the hydraulic gradient and the higher variation in the hydraulic head on the eroded slope are caused by the faster recharge due to the lower thickness of the unsaturated zone in this section (Fig. 6a). It is also worth mentioning the high vertical upward flow component throughout the hydrological year at the base of the eroded slope (Fig. 6a), which is manifested in the field by the presence of soil boiling points, especially in the dry season.
A water balance enabled an estimate of the water exfiltration flows on the gully face, showed the highest value in January (3.57 × 10 -6 m 3 /s), and a second high in May (3.30 × 10 -6 m 3 /s), which falls again until the end of the dry season. These results are similar to the ones presented by Drumond and Bacellar (2006) for the flow in the drainage channels of a gully in the region, with large flows in the rainy season and an increase in the dry period due to the increase in the base flow (Fig. 2b).
Slope safety factor analysis was conducted in the same hydrological year, by importing the flow simulations (head distribution) and not adding the phreatic surface, as is usual (Puy and Ramesch 2006;Ventura and Bacellar 2011). It then became possible to analyze the influence of the hydraulic head magnitude and direction in the saturated and unsaturated zone on slope stability. When the phreatic surface is just added in a limit equilibrium analysis, it is wrongly assumed that the water flows are horizontal and the head equipotentials are vertical. It was observed that the critical slipping surface throughout the hydrological year showed little variation in the safety factor (Fig. 6b,   Fig. 6 a Simulated geological cross section (location in Fig. 3), illustrating the position of the water table at the beginning of the simulation 10/01/2018 at the end of the dry season (dashed line), and on 01/01/2019 in the rainy season. The water table rises most in the eroded sector of the slope in the rainy season because of the thinner unsaturated zone. The detail shows the upward flow components that occur throughout the year at the toe of the slope. b Results of the safety factor (brown line) over the hydrological year and the corresponding smoothed curve (red line). The dry period is marked with a blue rectangle. The small increase in the safety factor in the second half of January is caused by the reduction in recharge due to dry spell (DS). Note that the average safety factor is low at the beginning of the dry season, when the water table is high, and gradually rises in subsequent months. c Slump rupture surface (in white) at the beginning of the hydrological year (10/01/18) and the potential zone of rupture (in red) c), always remaining in the potential range of rupture, as is expected for such unstable slopes. A bimodal behavior of the safety factor is observed, with two minimum values, one at the height of the rainy season in January (01/01/19), influenced by the rise of the phreatic surface in the eroded slope, and another in April (04/17/19), when this surface rises inside the slope. The lowest value persists during the first months of the dry season, influenced by delayed water table recharge. The safety factor starts to increase from May when the water table begins to recede again. The small increase in the safety factor during the second half of January (Fig. 6b) is a consequence of the dry spell (veranico).
The small variation in the safety factor may be a consequence of the lower rainfall of the year studied, which was approximately 82% of the local historical average (1348 mm y −1 ). The lower rainfall also explains why the water table does not return to the initial level at the end of the simulation (Fig. 6b). Slope undercutting during the dry season may also result in a decrease in the safety factor. However, stability analyses of the safety factor in this condition have not been done, as the forms of erosion at the slope base are very variable in space and time.
The analysis of these data and previous work (Bacellar et al. 2005;Drumond and Bacellar 2006;Costa and Bacellar 2007;Lima et al. 2018) allows the construction of an evolutionary model for these gullies. It was possible to confirm two patterns of water percolation and aquifer recharge, which reflect the characteristics of the soils and the geometry of the gully slopes. On the portion of the natural slope, unaffected by the gully, the Bw horizon is more porous and permeable and water infiltrates more easily. However, it takes longer to recharge the aquifer at a depth of over 18 m due to the thick, relatively less permeable saprolite. In this section, the phreatic surface reaches its peak between July and August in the period of drought, confirming the findings of Drumond and Bacellar (2006), who found a 5-month lag time between the peak of the rain and the recharge of groundwater located 20 m deep.
Conversely, on the gully face, the saprolite emerges and the phreatic surface is shallower (< 3 m), enabling faster recharge, especially when considering the rising capillary fringe which is high in this silty material. As a result of this distinct behavior, the safety factor showed low values throughout the year, however, is especially lower at two different times: the peak of the wet season when the water table rises at the bottom of the slope, and the middle of the dry season when the water table rises inside the slope. In the rainy season, large slumps, involving the entire slope may occur, associated with other surface erosion processes (Fig. 7a). In the dry season, slope undercut prevails because of the increase in the hydraulic gradient (Fig. 7b).
The rising hydraulic gradient at the toe of the slope is not sufficient to cause the static liquefaction of the saprolite but contributes to the decrease in the effective stress and, therefore, in the stability of the slope. However, the small slips observed on this stretch cause dynamic liquefaction, transforming them into smaller silt flow slides (Fig. 7b). The channelized flow through small discontinuities also generates small-scale piping erosion. All these mechanisms lead to the removal of slope-supporting saprolite in the dry season, facilitating the reactivation of large slumps affecting the whole slope in the subsequent wet season.
Therefore, the rise of the phreatic surface in the dry months mobilizes the sediments, but the water flow in the drainage channels of the gully does not have the capacity to transport all of them. The suspended sediment transport in these channels increases during the first summer rains (October/November), as demonstrated by Drumond and Bacellar (2006). The drawdown of the phreatic surface in the vicinity of gullies causes a two-stage recharge pattern, faster at the toe of the gully slopes and slower towards the 1 3 slope's interior (not eroded). This differs from the pristine morphology, where recharge in the non-eroded headwaters was gradual and sediment transport by the drainage system is significantly lower.
Since the instability of these slopes is influenced by the groundwater regime, providing efficient drainage at the slope base is an effective means of stabilization. However, the difficult access to the gully bottom, due to deep, long and tortuous channels, makes it impossible to use most conventional slope stabilization methods, such as the installation of toe drains in trenches or horizontal drains (Abramson et al. 2001). Any excavation at the gully bottom is dangerous since it can increase the upward hydraulic gradient and may trigger quicksand behavior of the soil or even a large rupture.
However, solutions for an unstable slope depend on a number of factors, such as slope geometry, mechanical and hydraulic behavior of soils, climate, etc. A low-cost option that has proven effective in stabilizing these kinds of gullies slopes is sediment containment by permeable walls, such as gabions. In this sense, Carvalho et al. (2018) employed gabion walls with geotextile on the back face in order to retain the sediments and allow water seepage. The gabion wall was backfilled with construction and demolition waste and Fig. 7 Erosion process operating on the slopes of the gullies. a In the rainy season, large slumps, involving the entire slope may occur, associated with others surface erosion processes such as plunge-pool erosion, sheet erosion, and shallow ravines; b in the dry season, as a result of delayed water table recharge, several mechanisms can act by undercutting the slope base, preparing for further slumps in the subsequent rainy season. c A possible solution to stabilize the slopes of gullies, with the construction of toe and blanket drains backfilled with soil from the slope benching (not shown in the picture) subsequent raises are made in the upstream direction until the original topography is partially or even completely restored. The high shear strength and hydraulic conductivity of this kind of waste ensures its stability. However, this method requires that waste is available at short range, which is not always feasible. Another option is the method proposed by Prandini et al. (1974) which involves the construction of a landfill, with material from the slope benching, thus making the execution more economic and faster. As the slope soils show low hydraulic conductivity, drains are necessary. Prandini et al. (1974) suggested an alternative drainage system, which consists of throwing bags of soil-cement at the toe of the slope to work as a toe drain and a drainage blanket, which is subsequently covered by the landfill (Fig. 7c). The throwing of bags is a good option, since the base of these gullies slopes is an area of high geological risk and of difficult to access. The stabilization is complemented by a retaining wall designed with gabions and connected to the drainage blanket. Except for the gabion screens and cement, all other materials can be obtained locally and executed by unskilled labor, reducing costs.
In both stabilization methods, the recommended time to build the contention structures is at the end of the dry season, when the slope is more drained and more stable. For complete stabilization, it is recommended to continue with the control of runoff. These passive containment methods make it possible to maintain the water table deeper, at lower levels, which enables vegetation recovery and sediment retention (Fig. 7c). It is not necessary to backfill the entire gully to ensure stability. Equilibrium-limit analyses show that at about 40% of the backfilling, stable conditions have already been attained, with a safety factor higher than 1.30.
Therefore, these can be feasible solutions to stabilize the slopes and prevent gullies from expanding. This is of great importance to the study area, which has suffered severe economic, social, and environmental impacts because of the accelerated evolution of hundreds of gullies.

Conclusions
Like many other areas with crystalline basement rocks in Southeastern Brazil, the Bação Complex is greatly affected by channelized erosion. These erosion features begin with the concentration of surface waters, forming ravines, but when the groundwater at the base of the erodible saprolite is reached, the processes of subsurface erosion begin, giving rise to hundreds of large gullies. These gullies pose a high hazard for people and properties, as they cause loss of agricultural soil, destruction of buildings and habitats, change in hydrological regime, siltation of channels and hydroelectric reservoirs, in addition to lowering of the phreatic surface, with baseflow reduction and drying of springs. For this reason, the usual practice of controlling surface runoff by terracing in contour lines and revegetation has proven ineffective in stabilizing these gullies. Among the processes of subsurface erosion, the mass movements stand out, especially the slumps.
The climate in the Bação Complex is seasonal, with 5-6 months of rainy season alternating with months of low rainfall. The monitoring of these gullies along the hydrological year indicates that slumps occur not only in the rainy season but also in the middle of the dry season. This phenomenon can be explained by the slow rise of the water table identified in the field because of the recharge, delayed to the rainfall peak. To verify this possibility, a representative gully with unstable slopes in the southeast of the Bação Complex was selected. The numerical simulations of saturated and unsaturated flow on one of the slopes of this gully for a hydrological year, showed that the water table presents seasonal fluctuations, with a maximum five to six months' delay in relation to the peak of the rainy season. This recharge delay is caused by the low percolation speed through the thick unsaturated zone, largely formed by the saprolite of silty texture. On the eroded stretch of the slope, where the unsaturated zone is shallower, the phreatic surface responds more quickly to rainy events. There is also strong upward flow at the toe of the slope throughout the year, with greater intensity in the dry season. These gradients explain the erosive processes that undercut the slope in the dry season, like boiling, piping erosion, creep, and flow slides.
Stability analyses of this slope by limit-equilibrium coupled with numeric flow simulation revealed a meta-stable situation, with low safety factor values, with little variation throughout the hydrological year. The safety factor is lower in the rainy season due to the rapid saturation of the toe of the slope and then decreases again in the middle of the dry season, when the aquifer is recharged more slowly on the rest of the slope.
These data confirm that the stabilization of these erosive features can be achieved by the control of runoff and with the improvement of the water exfiltration condition at the toe of the gully. However, conventional stabilization measures are difficult and costly, due to the unstable conditions and the difficulty of access the base of these gullies. Therefore, slope stabilization with gabion walls backfilled with local soil or with construction and demolishing wastes may be a good alternative stabilization method.
The monitoring of the evolution of hundreds of large gullies in the Bação Complex region is paramount because of the high social, economic, and environmental impacts that de destabilization of such landscape features can cause. Therefore, the proposition of a low-cost stabilization method for easy implementation is a major issue for regional sustainable development.