WEPP MODEL FOR RILL EROSION ESTIMATION IN A BRAZILIAN SEMIARID WATERSHED

Soil erosion is a major environmental problem in many parts of the world and represents a serious problem for sustainable agriculture and the environment, with direct and indirect impacts on soil quality and fertility. This study aimed to use the Water Erosion Prediction Project (WEPP) model to estimate rill erosion and determine soil physical and hydraulic properties, which are essential to investigate its performance. To this end, an experiment was carried out in the Exu Creek watershed, in Serra Talhada, semi-arid region of Pernambuco State (Brazil), under increasing flow rates: T1: 5.87 L min-1; T2: 12.10 L min-1; T3: 20.33 L min1; and T4: 27.57 L min-1. Liquid and solid discharges were sampled for determination and characterization of hydraulic parameters in preformed rill flows. Reynolds numbers between 2,019 and 6,929 and Froude numbers below 1 found in this study attest to occurrence of erosive rills. Soil losses due to rill erosion increased as flow rates increased. Rill erodibility was 0.0011 kg N-1 s-1, and critical shear stress (τc) was 1.91 Pa, collapsing rill sidewalls and increasing local uplift, wet perimeter, and rill hydraulic radius. keywords: Hydraulic relations. Flow rates. Erodibility. Shear stress. ABORDAGEM WEPP EM EROSÃO DO SOLO PARA BACIAS HIDROGRÁFICAS DO SEMIÁRIDO BRASILEIRO RESUMO A erosão do solo é um dos maiores problemas em muitas partes do mundo com impactos diretos e indiretos representando um problema muito sério na agricultura e sustentabilidade ambiental, afetando a fertilidade e qualidade do solo. Elevada erosão do solo e altas cargas de sedimentos são problemas sérios em várias regiões semiáridas no Brasil. No intuito de aplicar procedimentos de mitigação, são necessários uma estimação precisa da erosão do solo e produção de sedimentos. Este trabalho visa a abordagem do Water Erosion Prediction Project (WEPP) na quantificação da erosão em sulcos e avaliação das relações físicas e hidráulicas, essenciais para investigar o desempenho do modelo. Com este objetivo, o experimento foi desenvolvido na região semiárida da bacia hidrográfica do Exu, Serra Talhada Pernambuco (Brasil), através de aplicações crescentes de diferentes vazões: T1: 5,87 L min; T2: 12,10 L min; T3: 20,33 L min e T4: 27,57 L min. Amostras de descargas líquidas e sólidas foram coletadas para determinação e caracterização dos parâmetros hidráulicos da vazão nos sulcos pré-estabelecidos. Os números de Reynolds entre 2.019 e 6.929 e de Froude abaixo de 1, atestaram a ocorrência de sulcos erosivos. As perdas de solo devido à erosão nos sulcos foram elevadas após os crescentes fluxos aplicados e o valor da erodibilidade do sulco obtido foi de 0,0011 kg N s e a tensão crítica de cisalhamento (τc) de 1,91 Pa causando o colapso das paredes laterais, elevação da área, perímetro molhado e o raio hidráulico dos sulcos experimentais. Palavras-chave: Relações hidráulicas. Taxas de fluxo. Erodibilidade. Tensão de cisalhamento. _______________________________ Corresponding author Received for publication in 06/03/2019; accepted in 04/13/2020. Department of Rural Technology, Universidade Federal Rural de Pernambuco, Recife, PE, Brazil; victorcasimiropiscoya@gmail.com ORCID: 0000-0003-1875-9771, sousa.waldemir@gmail.com ORCID: 0000-0002-3817-8766, robsonrcpm@gmail.com ORCID: 00000003-4594-5748. Department of Agronomy, Universidade Federal Rural de Pernambuco, Recife, PE, Brazil; cantalice21@hotmail.com ORCID: 00000002-8209-341X. Department of Statistics and Informatics, Universidade Federal Rural de Pernambuco, Recife, PE, Brazil; moacyr2006@gmail.br ORCID: 0000-0002-3466-8143. Department of Forestry Engineering, Universidade Federal do Tocantins, Gurupi, TO, Brazil; renisson@uft.edu.br ORCID: 0000-00029747-1276. WEPP MODEL FOR RILL EROSION ESTIMATION IN A BRAZILIAN SEMIARID WATERSHED V. C. PISCOYA et al. Rev. Caatinga, Mossoró, v. 33, n. 3, p. 835 – 843, jul. – set., 2020 836 INTRODUCTION Soil erosion is a major environmental problem in many parts of the world (EEKHOUT; VENTE, 2019). In the Brazilian semi-arid, it is one of the most serious problems threatening ecosystems, mainly in human-altered areas for their greater susceptibility to detachment and sediment transport (SILVA et al., 2019). The main impacts on these degraded areas are loss of vegetation cover and decomposed organic matter on soil surface and soil compaction due to farming practices such as overgrazing, vegetation burning, and no-till use. Moreover, the absence of vegetation cover accelerates soil erosion through water runoff (MOGHADAM et al., 2015). Rill erosion is often observed on slopes under cropped areas and contributes significantly to sediment deposition into rivers in a watershed. This type of erosion constitutes a series of small channels with a depth of around 0.30 m, being a transition between interrill and gully erosions. Several are the factors influencing rill formation such as rain and its runoff intensity, soil type, local topography, tillage system, and crop species (SUN et al., 2013). According to Heras et al. (2011), drainage by set of rills develops from flow concentration and connectivity of flows from these channels. In this context, mathematicians have been arousing great interest in rill erosion network modelling (YAN et al., 2015). Among the used models, the Water Erosion Prediction Project (WEPP) is a physically-based and most used. It was developed by the United States Department of Agriculture (USDA) and was based on the theoretical foundations of infiltration, soil physics, plant science, hydraulics, and erosion mechanics (GHOLAMI et al., 2018). This has numerous advantages over existing models, as it incorporates effects from land uses such as farming, cattle ranching, and forestry. It models spatial and temporal variability of factors that affect hydrological regime and erosion of hillsides. It has been very useful in cropping, timber harvesting, roads, and natural fire areas with predominance of Hortonian flows. WEPP is a well-established tool to simulate water erosion and sediment yield. It has already been applied in several geographical locations such as the United States (AMPOMAH et al. 2020), Australia (YANG et al. 2018), and Europe (NICOSIA et al., 2019). Bezerra et al. (2010) also used it to study hydraulic characteristics in preformed erosion rills in a Cambisol, in the semiarid region of Brazil, with different flow rates applied to erosion rills. WEPP considers flow generation as the result of excessive rainfall, which is described by the Green-Ampt infiltration model (DUN et. Al., 2009). However, studies on combinations and interactions of different factors are lacking, e.g. intrinsic mechanisms of erosion of rills. These factors are not yet clear due to their complexity, especially when under different physical conditions (SUN et al., 2013). Given the above, this study aimed to evaluate the use of WEPP model to quantify erosion of rills and to evaluate physical and hydraulic relationships, which are essential to verify the performance of this model. MATERIALS AND METHODS Experimental area location and characterization The study area is the Exu Creek watershed, which has a total area of 579,40 km and lies between the geographic coordinates of 8o 8' 2.4" and 8 9’ 7.2’’ south latitude, and 38o 24' 21.6" and 38 23’ 24.0’’west longitude (Figure 1). It is within the Pajeú River watershed, on the BR-232 federal highway km 448, and on the vicinity of the city of Serra Talhada PE (Brazil). Part of the Exu Creek watershed coincides with the water dividers in the Pajeú River watershed. According to Köppen’s classification, the climate is ‘BSwh’ type, which stands for semiarid, hot and dry (MELO et al., 2008), with summerautumn rains, and average annual rainfall of 484.06 mm year from 1992 to 2007 (PERNAMBUCO, 2015). The vegetation is basically composed of the caatinga biome, with great variety of landscapes, biological richness, and endemism. The occurrence of seasonal droughts sets intermittent river regimes and defoliates vegetation. According to the classification of Santos et al. (2013), the most common soil type found in the studied watershed includes Entisol Fluvent. The soil in the study area was physically characterized using samples collected from the 0 – 10 cm depth layer. These samples were air-dried and sieved through 2mm sieves, and the apparent results are shown in (Table 1). The morphometric parameters of the Exu Creek watershed are shown in (Table 2). WEPP MODEL FOR RILL EROSION ESTIMATION IN A BRAZILIAN SEMIARID WATERSHED V. C. PISCOYA et al. Rev. Caatinga, Mossoró, v. 33, n. 3, p. 835 – 843, jul. – set., 2020 837 Figure 1. Location of the Exu Creek watershed in the state of Pernambuco and in Brazil. Table 1. Soil physical characteristics from 0-10 cm depth layer in the Exu Creek watershed, Serra Talhada, Pernambuco State (Brazil). Grain size analysis Classification Sand Silt Clay Texture PD BD % (g cm) 59 21 20 Sandy loam 2.51 1.41 1 PD – Particles density; BD – Bulk density. Table 2. Physical, hydrographic, and morphometric parameters of the Exu Creek watershed, in the semiarid of Pernambuco State, obtained using the SRTM images (Shuttle Radar Topography Mission).


INTRODUCTION
Soil erosion is a major environmental problem in many parts of the world (EEKHOUT; VENTE, 2019). In the Brazilian semi-arid, it is one of the most serious problems threatening ecosystems, mainly in human-altered areas for their greater susceptibility to detachment and sediment transport (SILVA et al., 2019). The main impacts on these degraded areas are loss of vegetation cover and decomposed organic matter on soil surface and soil compaction due to farming practices such as overgrazing, vegetation burning, and no-till use. Moreover, the absence of vegetation cover accelerates soil erosion through water runoff (MOGHADAM et al., 2015).
Rill erosion is often observed on slopes under cropped areas and contributes significantly to sediment deposition into rivers in a watershed. This type of erosion constitutes a series of small channels with a depth of around 0.30 m, being a transition between interrill and gully erosions. Several are the factors influencing rill formation such as rain and its runoff intensity, soil type, local topography, tillage system, and crop species (SUN et al., 2013). According to Heras et al. (2011), drainage by set of rills develops from flow concentration and connectivity of flows from these channels. In this context, mathematicians have been arousing great interest in rill erosion network modelling (YAN et al., 2015).
Among the used models, the Water Erosion Prediction Project (WEPP) is a physically-based and most used. It was developed by the United States Department of Agriculture (USDA) and was based on the theoretical foundations of infiltration, soil physics, plant science, hydraulics, and erosion mechanics (GHOLAMI et al., 2018). This has numerous advantages over existing models, as it incorporates effects from land uses such as farming, cattle ranching, and forestry. It models spatial and temporal variability of factors that affect hydrological regime and erosion of hillsides. It has been very useful in cropping, timber harvesting, roads, and natural fire areas with predominance of Hortonian flows.
WEPP is a well-established tool to simulate water erosion and sediment yield. It has already been applied in several geographical locations such as the United States (AMPOMAH et al. 2020), Australia (YANG et al. 2018), andEurope (NICOSIA et al., 2019). Bezerra et al. (2010) also used it to study hydraulic characteristics in preformed erosion rills in a Cambisol, in the semiarid region of Brazil, with different flow rates applied to erosion rills. WEPP considers flow generation as the result of excessive rainfall, which is described by the Green-Ampt infiltration model (DUN et. Al., 2009). However, studies on combinations and interactions of different factors are lacking, e.g. intrinsic mechanisms of erosion of rills. These factors are not yet clear due to their complexity, especially when under different physical conditions (SUN et al., 2013).
Given the above, this study aimed to evaluate the use of WEPP model to quantify erosion of rills and to evaluate physical and hydraulic relationships, which are essential to verify the performance of this model.

Experimental area location and characterization
The study area is the Exu Creek watershed, which has a total area of 579,40 km 2 and lies between the geographic coordinates of 8º 8' 2.4" and 8 o 9' 7.2'' south latitude, and 38º 24' 21.6" and 38 o 23' 24.0''west longitude ( Figure 1). It is within the Pajeú River watershed, on the BR-232 federal highway -km 448, and on the vicinity of the city of Serra Talhada -PE (Brazil). Part of the Exu Creek watershed coincides with the water dividers in the Pajeú River watershed.
According to Köppen's classification, the climate is 'BSwh' type, which stands for semiarid, hot and dry (MELO et al., 2008), with summerautumn rains, and average annual rainfall of 484.06 mm year -1 from 1992 to 2007 (PERNAMBUCO, 2015). The vegetation is basically composed of the caatinga biome, with great variety of landscapes, biological richness, and endemism. The occurrence of seasonal droughts sets intermittent river regimes and defoliates vegetation.
According to the classification of Santos et al. (2013), the most common soil type found in the studied watershed includes Entisol Fluvent. The soil in the study area was physically characterized using samples collected from the 0 -10 cm depth layer. These samples were air-dried and sieved through 2mm sieves, and the apparent results are shown in (Table 1). The morphometric parameters of the Exu Creek watershed are shown in (Table 2)  1 PD -Particles density; BD -Bulk density.

Determination of erosion rates in preformed rills
Experimental plots consisted of preformed rills constructed using the cutting edge of a hoe, following the natural slope of the land. Sixteen preformed rills were constructed and divided into four blocks of four. They were immediately subjected to erosion tests, applying different flow levels at random in each block, one test at a time.
Treatments consisted of different flow levels, as follows: T1) 5.87 L min -1 ; T2) 12.10 L min -1 ; T3) 20.33 L min -1 ; and T4) 27.57 L min -1 . At the upper end of the rill, energy sinks were buried in form of circular plastic containers. Then, water-conducting hoses were into these containers, reaching the rills after overflowing the containers. At the lower end of the rills, a metal collecting trough was installed to assist in sampling of liquid and solid discharges. At the end of each gutter, flow samples were collected every 5 minutes in 1000 mL beakers and then transferred to properly identified 1-L plastic pots, sealed, and conditioned in racks for subsequent transport and evaluation in the laboratory (Figure 2).

Figure 2.
Preformed rill following the natural slope of the terrain, energy sinks buried at the upper end, and liquid and solid discharge collector installed at the lower end.

Determination of hydraulic parameters
Runoff velocity (Rv) was calculated by measuring the time spent for a methylene blue dye to reach 3 m away from the rill, with values expressed in ms -1 . These measures were taken every 3 minutes since formation of sheet flow. The same method was used to determine liquid and solid discharges. Average runoff velocity (m s -1 ) was obtained by multiplying the records during surface velocity tests by α = 0.6, which is a correction factor for runoff velocity variations with flow depth, due to soil friction (LAFAYETTE; CANTALICE; COUTINHO, 2011).

Geometrical parameters of hydraulic characterization of surface runoff in erosion rills
Wet area, wet perimeter, and hydraulic radius were estimated using a profilometer, through cross sectional measurements of the rill flow. Two measurements were taken, the first after the first 4 minutes of the test, and the second 4 minutes before the end of the test. Rill areas were measured using a mechanical planimeter. Channel perimeter was measured using an analogue curvimeter. If combined to the cross-sectional area of rills, the section perimeter provided the hydraulic radius used to calculate shear stress flow, as described below (Equation 1): Wherein: Rh is the hydraulic radius, A is the crosssectional area of the rill (m 2 ), and Pm is the wetted perimeter of the cross-section (m).
Hydraulic roughness in the rills was obtained by determining the Darcy-Weisbach Friction coefficient (ƒ), as in the (Equation 2) below: Wherein: g corresponds to the acceleration of gravity (m s -2 ), R h is the hydraulic radius (m), S is the slope of the terrain (m m -1 ), and V is the average velocity of runoff (m s -1 ), which was estimated by Manning's equation (n), as follows (Equation 3): Wherein: V is the average velocity of surface flow (m s -1 ), N is the Manning roughness coefficient (m s -1 ), R h : is the hydraulic radius (m), and S is the rill slope (m m -1 ).
The Reynolds number (Re) was determined according to the equation described by Simons and Senturk (1992), as follows (Equation 4) Wherein: R e is the Reynolds number (dimensionless), V m is the average runoff velocity (m s -1 ), R h is the hydraulic radius of the rill cross-section (m), and is the kinematic viscosity of water (m 2 s -1 ).
The kinematic viscosity of water was determined according to Bezerra et al. (2010), as in (Equation 5); and temperature (ºC) was measured using a thermometer in each experiment.

(5)
Wherein: ʊ is the kinematic viscosity of water (m 2 s -1 ) and T is the water temperature (°C).
The number of Froude (F r ) was obtained through the (Equation 6), according to Simons and Senturk (1992): Wherein: F r is the number of Froude (dimensionless), V m is the mean runoff velocity (m s -1 ), g is the gravity acceleration (m s -2 ), R h is the radius of the rill cross-section (m).
Flow shear stress can be determined after recognizing as true that, in rill erosion, as flow rates increases, sediment load increases to a level much larger than its transport capacity. Based on this, (Equation 7) of Elliot et al. (1989) was used to determine flow detachment capacity.
Wherein: D c is the flow detachment capacity in rills (kg m -2 s -1 ), K r is the soil erodibility in rills (Kg N -1 s -1 or s m -1 ), τ is the flow shear stress ʊ ʊ= 1.14-0.031 T-15 +0.00068 T-15 2 x10 -6 m 2 s -1 F r = V m g R h D c =K r τ − τc (N m -2 or Pa), and τc is the soil critical shear stress (N m -2 or Pa), which was expressed according to (Equation 8), as follows: Wherein: γ is the specific weight of water (N m -3 ), R h is the hydraulic rill radius (m), S is the rill slope (mm 1 ). (Equation 9) was used to estimate total soil losses from instantaneous sediment concentration data of runoff and liquid discharge rate (CANTALICE et al., 2013).
Wherein: P S is the total soil loss (kg m -2 ), Q in is the liquid discharge rate (L m -1 ), C in is the sediment concentration (kg L -1 ), t is the interval between samples (min), and A is the rill area (m 2 ).

Statistical analysis
The experiment was carried out under a randomized block design. Rill erosion data were statistically analysed using the SAS software. Means were compared by the Tukey's test at 5% significance. Microsoft Excel 2010 was used for regression analysis, refitting the equations to each regression model. Figure 3 shows soil detachment within the entire wet perimeter and rill sidewall collapse, forming an irregular rill bed. τ = γ R h S P S = Q in C in t A Figure 3. Soil detachment in the entire wet perimeter and rill sidewall collapse, forming an uneven rill bed. (Table 3) shows the hydraulic behaviour of preformed rills after application of increasing flow rates. Liquid discharges were statistically similar for flow rates of 5.87 and 12.10 L min -1 , whereas larger flow rates of 20.33 and 27.57 L min -1 showed a significant difference for this parameter.  Means followed by uppercase letters in the same line do not differ from each other at 5% probability by the Tukey's test. Q = liquid discharge; V m = mean flow velocity; S = mean rill slope; R e = Reynolds number; F r = Froude number: log ƒ = hydraulic roughness (Darcy-Weisbach coefficient).

RESULTS AND DISCUSSION
According to Reynolds numbers, a transitional flow regime was observed in the lowest liquid discharges, while the two highest flow rates had a turbulent regime, as observed by Piscoya et al. (2018), Knapen and Poesen (2010), and thus characterizing the hydraulic behaviour of runoff erosion in rills. Analysis of slope data (S) showed no significant differences, demonstrating\ that the preformed rills had the same location.
The mean velocities (V m ) of the different flow rates did not differ statistically from each other at 5% probability by the Tukey's test. It might have occurred because flow rates were tested on the surface of uncovered rills and under the same roughness conditions. This was confirmed by the statistically undifferentiated values of Darcy-Weisbach. However, the applied flow rates showed a potential increase in \velocity (Figure 4). Similar relationships were found by Piscoya et al. (2018) in a Haplic Cambisol soil, while simulating rainfall erosion in the Jacu River watershed, also in the semiarid region of northeaster Brazil. According to the Froude number, all flows were below 1 (one), and hence considered slow.
( Figure 5) shows the differentiation in hydraulic radius of rills as flow rate is raised, showing rill deepening and geometric response.
1 Figure 5. Exponential relationship between area and hydraulic radius in rills generated by increasing flow rates applied in an Entisol Fluvent soil.
Shear stresses increased with the rise in flow rates applied to preformed rills, with a significant difference between the lowest flow rate and the others, due to hydraulic radius variation. In turn, such shear stresses led to differences in soil detachment rates and losses due to the applied flow rates.
Likewise, soil losses increased from the second applied flow rate onwards, until values considered high for a semi-arid Entisol, poorly developed, and under high temperatures, and variable rainfall. Rill top erodibility (K r ) was 0.0011 kg N -1 s -1 , which was obtained by coefficient b or slope of the linear regression line ( Figure 5) relationship between shear stress and soil detachment rate for a fit with R 2 = 0.80 (KNAPEN; POESEN, 2010). The critical shear stress in the preformed rills, for τc = -a / b or the intercept at x for y = 0, was 1.91 Pa.
In the semiarid region, Bezerra et al. (2010) obtained a rill erodibility in a Cambisol of 0.0051 Kg N -1 s -1 and τc of 1.18 Pa. Also, in a Cambisol of the Pernambuco semiarid, Piscoya et al. (2018) found rill erodibility and critical shear stress of 0.00211 Kg N -1 s -1 and 2.34 Pa, respectively ( Figure  6).

Figure 6.
Rill soil erodibility (k r ) and critical shear stress (τ c ) obtained by regression of soil detachment rates (D r ) and respective shear stresses generated by the applied flow rates.
Rev. Caatinga, Mossoró, v. 33, n. 3, p. 835 -843, jul. -set., 2020 842 The rill erodibility obtained in our study can be justified because the profile of an Entisol Fluvent soil has a sandy grain-size composition. (Table 4) shows the erosion rates in the preformed rills in an Entisol Fluvent soil. In this table, one can observed significant differences between the smaller and larger applied flow rates in terms of soil detachment rates (D r ), while soil losses (SL) had greater differentiation due to an integrated effect. Means followed by uppercase letters in the same line do not differ from each other at 5% probability by the Tukey's test.

CONCLUSIONS
Shear stresses from the applied flow rates caused rill side wall collapse and local uplift, and increases in wet perimeter and hydraulic radius. Erodibility in rills was 0.0011 kg N -1 s -1 , while critical shear stress (τc) was 1.91 Pa. These results corroborate those in the literature for all the soil erosion models. The WEPP model proved to be accurate in predicting erosion, which is crucial for further development of new approaches.