Topographic thresholds for gully head formation in badlands

Gully erosion is a particularly damaging process which is not yet sufficiently understood and parameterized. Gully head topographic threshold relative to Hortonian runoff have been studied in cropland, rangeland and forest. This study extends such modelling approach to badlands. Different badlands (eight sites) have been studied in the Mediterranean environment in Italy and Spain, characterized by diversified climatic, lithological, and geological settings under different anthropogenic conditioning. Many badlands have been characterized by their specific human history in addition to their geomorphological properties. Land use, as part of the human history, strongly affected many badland formation and development, through extremely impacting land exploitation (usually overgrazing). The effect of geological and geomorphological processes are usually particularly well visible. While the weakening effect of joints is confirmed, the different geological layer bedding orientation with respect to the slope aspect generates a different development of badland morphologies and different values of gully head thresholds values (as shown in two badlands sites on the same geological material and climate.


| INTRODUCTION AND THEORETICAL BACKGROUND
Gully development is a particularly damaging erosion process which is not yet sufficiently well understood and parameterized . It can lead to badland development requiring intense topographic reshaping to remove gully channels and prevent the formation of new channels. Evidently, land use, slope gradient, soil and the underneath geological layers with their characteristics and those inherited by the local orogenetic processes (e.g., joints and faults) contribute to the resistance offered by the land to the erosional potential of concentrated overland flow.
Gully channels extend over lengths that can be measured in metres up to kilometres, with widths spanning over tens of metres to less than one and depths ranging from a few centimetres (e.g., ephemeral summer gullies, Nachtergaele et al., 2002) to tens of metres (e.g., Radoane et al., 1995;Poesen et al., 2003). Gullies occur in soils with a large textural range. Possible effects of soil texture on gully channel initiation and development have been observed and discussed in the literature (Radoane et al., 1995;Shahrivar & Christopher, 2012;Valentin et al., 2005). Evidently, a deep gully can erode by subsurface flow, seepage, piping and tunnel erosion, which can undercut more resistant top layers such as those colonized by plant roots (De Baets et al., 2008). Furthermore, tension cracks develop easily on vertical walls (particularly those of recently formed gully channels) which facilitate mass movements, leading to a change in slope gradient of the new channel walls towards gradients evolving towards the typical repose angle of the material. Lithology and its features may also strongly influence gully locations and patterns, particularly lithological layer bedding and orientation whereas density of joints may control drainage network, preferential flows and in turn surface and subsurface erosion (Beavis, 2000;Colica & Guasparri, 1990;Parkner et al., 2007). Soil material with its organic matter content (especially glomalin; Rillig & Mummey, 2006) makes soil aggregates more resistant to degradation and the pore system better structured for diffusing water into the soil profile. Also plant roots add to the soil resistance (Gyssels et al., 2005;De Baets et al., 2007). Subsurface flow can complement soil erosion by overland flow, resulting in a further upslope retreat of the gully head (GH). The more resistant top layer can reduce GH retreat rate. Despite important advances in understanding gully erosion processes over the last decades, predicting gully erosion still remains a major challenge (Poesen, 2017;Vanmaercke et al., 2021).
Despite large efforts made to map gullies and gully susceptibility, dealing with all the gully subprocesses remains difficult because of their interactions and because it is difficult to obtain adequate instrumentation to cover all possible field conditions .
Nevertheless, several models have been developed during the past 50 years. The empirical line that spans from the pioneering works of Beer and Johnson (1965) and Seginer (1966) to the recent review by Vanmaercke et al. (2016) suggests that rainfall intensity (average, maximum) and catchment area are two important factors for estimating GH retreat rates. More sophisticated models, such as those presented by Sidorchuk (1998Sidorchuk ( , 1999Sidorchuk ( , 2005Sidorchuk ( , 2006 with the dynamic, static, deterministic and stochastic modelling of gullies, or by De Ploey (1992), Willgoose (2005), Tucker et al. (2001), Campo-Bescos et al. (2013) and Torri and Poesen (2014), have been developed, based on a strong or a semi-empirical physical base. These models introduce parameters which are difficult to extract from public databases or time-consuming and costly to determine for each case study. The GH threshold model, introduced by Patton and Schumm (1975), refined by Montgomery and Dietrich (1994) and further expanded by Torri and Poesen (2014) to include land use, falls midway between a totally empirical and a physically-based approach. Such model predicts where GHs will form and then subsequently may retreat further upslope. Such sites in the landscape are described as a function of the local soil slope gradient at the GH, the surface area draining into it and the land use expressed through the maximum potential runoff production (which is predicted by the runoff Curve Number (CN); Hawkins et al., 2009). Next, rock fragment content in the topsoil, presence of joints (faults) and piping susceptibility are other (secondary) factors which help improving the prediction of the landscape position where a gully channel may form and retreat. This approach allows land-use types, such as rangeland, pasture, cropland and forest to be considered in the evaluation of the probability that a GH will be observed in a given landscape location.
Given the relatively simple data required for running the model proposed by Torri and Poesen (2014), it is worthwhile to further improve and expand this model to other land-use types. The dataset on which this model is based, contains field observations on GHs made in every continent and for several land-use types, including forest, rangeland, and cropland. Badlands were not included yet, despite the fact that they typically have a high gully density (e.g., Joshi, 2014;Joshi & Nagare, 2013;Torri, Poesen, et al., 2018). The Torri, Poesen, et al. (2018) study covers a Plio-Pleistocene marine sediment, silty clay, not cemented and overconsolidated, revealing two topographic threshold values, due to small differences in management for the same land use (overgrazed rangeland).
Given the scarce data on GHs in these environments, this study aims to better understand GH formation in badland-type environments by expanding the GH database for badlands that developed in other lithological and geological settings. Study sites have been selected in such a way that they are all located in a similar climatic zone, defined following the Köppen-Geiger classification as Mediterranean Cs type (see Table 1 for details on symbols and acronyms), in order to allow for the study of meteorological differences, without mixing their effects with climatic effects. Study sites were therefore selected in Italian and Spanish Mediterranean regions where badlands are common landscape features.
The study additionally explores the possible differences in gully occurrence between different badlands.
The study focus is on GH thresholds following the standardized methodology proposed by Torri and Poesen (2014) and revised by Torri, Poesen, et al. (2018) in order to avoid the risk of producing biased trends for topographic thresholds . The study will also investigate the contribution of soil rock fragment content to the processes linked to the resistance to GH retreat (Torri & Poesen, 2014).
To address all these items in a coherent way, this article (i) presents the theoretical framework behind the topographic threshold concept, (ii) compares this with field observations collected during this study, (iii) identifies the shortcomings of such a theoretical framework, and finally (iv) proposes a generalization of the theory to account also for gullies in badlands.

| Gully head threshold theoretical background
This paragraph summarizes the physics behind the gully-head threshold model, its empirical components, weaknesses and aspects to be T A B L E 1 Symbols and acronyms used in the text, with descriptions, units of measure and references for a more detailed information Symbol further developed. Montgomery and Dietrich (1994) proposed the basic inequality (Equation 1). Table 1 summarizes all symbols and   acronyms used in the text later, providing a proper description and   specifying the units of measure and references for more detailed information. Rossi et al. (2015) and Torri and Poesen (2014) Torri, Poesen, et al. (2018) expanded the observed range of land-use types with GHs to include one badland site. Furthermore, they proposed a modification of the basic equations in order to avoid problems for areas that have low CN values.
The threshold concept for GH formation which then retreats further upslope derives from the fact that concentrated overland flow should exert flow shear stresses in excess of a critical value to erode the soil. This critical flow shear stress is the value at which the most resistant fraction of the soil material is entrained by the flow and either wall or bed armouring is impeded. This condition should last long enough to allow the flow at the GH to cut a cross-sectional area exceeding a critical value (Hauge, 1977;Poesen et al. 2003). Montgomery and Dietrich (1994) noted that the flow shear stress (τ f ) produced by runoff must exceed a critical threshold value (τ cr expressed in pascals) to erode a gully channel. Hence, it follows: where ρ f is the fluid density measured in kg/m 3 , g is acceleration due to gravity in m/s 2 , h is flow depth (or hydraulic radius) measured in metres, and γ is slope angle of the soil surface.
Considering that: where Q is flow discharge rate measured in m 3 /s, u is mean flow velocity measured in m/s, w is the width of the channel head measured in metres, n is Manning's hydraulic roughness coefficient expressed in s/m 1/3 .
According to Montgomery and Dietrich (1994), the discharge rate Q, which assumes that the entire upslope catchment area A contributes to the overland flow (i.e., when depression storage is filled, hence R > I), is given by: where R is rainfall intensity and I is infiltration capacity, both expressed in m/s, and A is measured in m 2 .
Introducing Equation (8) and sinγA 0:4 > k with k ¼ k 0:6 w τ cr gρ f n 0:6 R À I ð Þ 0:28 with k (reflecting the resistance of a site to gully development) is a coefficient accounting for several factors, some of which relate to the soil resistance (e.g., τ cr ) and others relate to flow erosivity (e.g., R À I).
Despite this relation is supported by the literature (this is substantially derived from equation (7) in Montgomery and Dietrich, 1994), there are still parameter interdependencies to be explored. For instance, k decreases with increasing hydraulic roughness (n), which is counterintuitive. This apparent inconsistency is due to the fact that in Equation (10), the physical threshold is τ cr while k is a dummy constant used to simplify a multi-dimensional threshold into a two dimensional one. Other approaches removed this apparent inconsistency, but introduced less well-established empirical equations.
An alternative k formulation can be based on the stream power per unit of volume W in W/m 3 (SI system; Torri et al., 2012): where u is mean flow velocity, which in sediment laden flows depends only on discharge: where k u and β are an empirical constant with β 0:3;0:45 ½ (respectively for rills and gullies, see Govers et al., 2000;Nearing et al., 1999).
Using Equations (3), (12) and (11) we obtain Equation (10) with Torri and Poesen (2014) developed an equation for predicting k values for GHs caused by rainfall excess runoff as a function of S 0:05 (Hawkins et al., 2009). This parameter is a quantity derived from the S 0:20 is derived directly from CN and it is a fundamental quantity in the NRCS method which appears in its basic equation for runoff (Q CN , in millimetres) evaluation given the daily precipitation P (in millimetres) and the land use (S 0:20 ): The equation for the threshold k, which was proposed by Torri and Poesen (2014), is the following: where ψ (originally named c in Torri and Poesen [2014]) represents all factors that may influence k but could only be identified for too a few cases to develop a possible correction factor (e.g., for GHs forming on joints or fault lines, or for the presence of piping at some depth). RFC is rock fragment content (estimated visually, hence it is a surface cover).
Unfortunately, Equation (17)  1.2 | Tan γ ð Þ or sin γ ð Þ: Additional theoretical basis for a new insight into the threshold model Following the approach for the estimation of k, supported by Montgomery and Dietrich (1994) and expressed by Equations (9) and (10), the GH threshold depends on τ cr which is the threshold stress exerted by the flow on the soil particles and aggregates to entrain these. This critical stress is assumed to be proportional to the soil strength. A first model can be constructed assuming the forces expressed by vectors instead of tensors. The latter would require more basic data, which are usually not available, especially for such small soil depths (i.e., at the very surface) as discussed by Brunori et al. (1989). The resistance vector (R ! expressed in pascals) of a soil grain per unit of soil surface may be expressed as follows. First, there is a component parallel to the slope (R x ) due to the gravity component normal to the soil surface (pressure exerted over the surface of contact between the soil mass and the soil grain). This component must be decreased by the gravity component parallel to the soil surface (which is not included in the potentially detaching fluid stress): where h f (dimensionless) is the fraction of the grain (or rock fragment) that is actually submerged in the flow (h f = 1 represents total submersion) and ζ (measured in metres) represents the diameter of the grain normal to the soil surface (usually the shortest grain or rock fragment axis), φ is the friction angle, ρ s is the saturated soil grain density and ρ f is the fluid density both expressed in kg/m 3 .
Note that with increasing slope gradient the subtractive gravity component grows while the frictional part decreases, hence R x necessarily decreases with slope.
The R ! component (R z ) normal to the soil surface is the component that is opposed to the runoff lift forces and equals: The incipient motion of soil grains will occur in a certain direction following which of the two components of the drag/lift forces, produced by overland flow at the GH, dominates. Along this direction and opposite to the detaching force, another resistance force, due to various bonds between the grain and the surrounding soil grains, adds to R ! and this is what is usually called cohesion (c).
The only assumption made here for estimating the intensity of the resisting pressure, and then of the critical flow detaching force, is that of using the following estimation of the total resistance intensity, supposing that cohesion and R ! have the same direction and sense: where (20) is substantially a Ranking Coulomb approximation with a frictional and a cohesive component, the main difference being that it also contains the gravity detaching component.
As stated earlier, the critical fluid pressure (drag and lift) can be assumed to be proportional to the soil resisting forces per soil surface unit, that is τ cr $ R tot . Hence, the threshold equation can be rewritten to account for the additional dependence on γ. Let us now introduce Equation (20) into Equation (10): The k value of Equation (21) corresponds to the k 1st part, a constant: we have substituted a fluid shear stress with soil resistance but, as Nearing (1991) noted, the fluid shear strength is much smaller than the soil resistance probably because only the most intense eddies could actually entrain the soil grains (hence, turbulent flow conditions).
Equation (21) introduces four new parameters to be estimated, that is cohesion, friction angle, average grain or rock fragment size, and submerged grain or rock fragment density. Some of these can be estimated from the literature, as shown later in the section. Now let us take a look at the theoretical threshold trends. The resistance term R tot is made of two parts, only one changing with the slope angle γ.
Hence, the threshold line will vary in shape between two extremes following whether the cohesion term is larger than the friction term or F I G U R E 1 Relation between hillslope gradient and drainage area at gully heads. Four curves are shown for four c/R ratios (0, 1, 3300, 3). The three ratios correspond to three total resistance values (c + R = 50, 100, 200 Pa). The third value (200 Pa) is obtained for two situations, one valid for a fine soil and one for a coarse soil but still rich in clay, hence with two different ratios of cohesive to frictional forces. The two corresponding curves differ only for steep gradient values. This bending is present in all the curves with a large frictional component (i.e., frictional component > 0.3 cohesion) [Color figure can be viewed at wileyonlinelibrary.com] vice versa ( Figure 1). If the cohesion term dominates over the frictional one (i.e., large values of c/R in Figure 1), the threshold line will be a straight line in the log-log area-slope plot (the hillslope gradient is expressed as the sine of the slope angle). Conversely if the frictional resistance becomes more and more relevant (i.e., low values of c/R in Figure 1), the threshold line will bend for low drainage area values. In such a case, the threshold straight line is obtained expressing the hillslope gradient as the tangent of the slope angle.
We can run some simulations varying the ratio between the vector intensities c/R from zero (only friction forces) to infinity (only cohesive forces) and see how the threshold curve varies.
The trends shown in Figure 1 indicates that one must expect a convex threshold line when the frictional component is sufficiently large, which is the case when R > 0.3c.
Let us now describe a way for estimating the frictional parameters. In order to estimate the friction angle, we can make use of some relationships reported by Poesen & Torri (1989) and Torri et al. (1990), in agreement with Nearing (1991). They found that for a grain or rock fragment of a given size and a given shape, its friction angle over a fixed bed is a function of the relative size of the grain with respect to the grain size of the bed. Re-examining their data and those reported in other studies (Biancalani, 1988(Biancalani, , 1989Dewilde, 1986;Torri & Poesen, 1988) it was possible to draw the trends shown in Figure , where the angular and the round rock fragment data where presented in Biancalani (1988), while the intermediately shaped rock fragment data are shown in Poesen & Torri (1989) and Torri et al. (1990).
These friction angles, as reported by Poesen and Torri (1989) and Torri et al. (1990), are distributed roughly following a gaussian distribution. As the friction angle does vary between 0 and 90 the equation of the gaussian curve should be modified for transforming the finite interval into an infinite one. As we are not looking for situations very close to 0 or to 90 , we will consider the usual form of the normal Gauss distribution. Figure 3 represents the variation of the friction angle standard deviation for the same relative grain size as shown in Figure 1: Combining the results shown in Figures 2 and 3, we can now calculate the mean value of the friction angles for soils, using the following formula: where φ 50 and φ 80 are respectively the 50th and the 80th percentiles; S d,φ is the standard deviation shown in Figure 3.

| The ratio between the fluid stresses and the soil resistance strengths
Usually one uses soil resistance characteristics which apparently do not match the stresses and pressures exerted by the fluid on the soil grains and aggregates to detach and entrain these. Actually, when one measures the soil resistance, whatever measuring device is used, the resulting resistance values are always larger than the fluid stresses measured using equations such as the one for the bottom shear stress. This is a problem that was addressed by Nearing (1991). Here, we follow this approach and add some means to evaluate the critical fluid stresses in the absence of any measurements of peak overland flow.
The following three statements summarize this and allow to further develop a model. Nearing (1991) suggested that the distribution of fluid stresses for a given value of average fluid stress has a tail toward high stress values, that overlaps the tail of the distribution of soil grain resistances towards low soil strength values; this causes the apparent contradiction that lower fluid stresses erode soils with higher strengths (i.e., a low average fluid stress, F stress , will erode a soil with a much F I G U R E 2 The 80th percentile of the friction angle distribution (φ 80 ) as a function of the relative rock fragment size (D p , rock fragments) with respect to the median grain or rock fragment size (D 50 ) of the channel bed. The friction angle φ 80 expresses the probability of motion equal to 80%. The three curves show the effect of the rock fragment shape and roundness (angular means tending towards prismatic, with sharp corners, intermediates are more ellipsoidal in shape, while the round ones have bulging faces and rounded corners). The angular fragments are usually characterized by a ratio between the maximum (Dmax) and the minimum (Dmin) diameter larger than 2.5, the intermediate fragments have Dmax/ Dmin ratios between than 2.5 and 1, while for the round fragments this ratio varies between 1 and 1.3. [Color figure can be viewed at wileyonlinelibrary.com] F I G U R E 3 Relation between the standard deviation of the soil grain friction angle (φ) calculated from two data sets reported by Poesen and Torri (1989) and Torri et al. (1990) and the relative rock fragment size (D p , rock fragments) with respect to the median grain or rock fragment size (D 50 ) of the channel bed. [Color figure can be viewed at wileyonlinelibrary.com] higher average strength, R soil ); this is often represented through a ratio, FS ratio , between F stress and R soil ; i. FS ratio resulted to be much lower than unity, with an order of magnitude given for rills by Torri et al. (1987) (corrected according to Brunori et al., 1989) and by Abdel-Rahman (1963) reported in Torri et al. (1987), having values of 0.0002 and 0.0004, respectively; those values can be considered as an approximate lower boundary value for gullies; ii. Beasley (1959, 1961) proposed an equation for evaluating the average critical fluid stress using the clay fraction (C) (see also Nachtergaele et al., 2002): Note that when C ¼ 0 then F stress ¼ 0:5 Pa, independently from the shape and other variables driving the friction angle, as discussed pre-

viously (Figures 2 and 3).
Hence, we expect that the ratio FS ratio between the fluid forces acting at the GHs and the soil resistance forces varies between 0:0004 < FS ratio < 1.

| Characteristics, origin and evolution of badlands
Verghereto and Montione badlands (Figure 4) share the same landuse history, given their proximity and the status of dependence of Montione with respect to Verghereto. The short description that follows is based on books such as Amati (1868), Rampoldi (1832), and Repetti (1841) Loamy-skeletal variations are also common.
Obviously, these soils characterize the area not affected by gully erosion. Where gully erosion dominates, the soils substantially consist of rock fragments, with a scarcely-developed soil horizon structure (A/C over C or directly over R with A/C horizon very shallow), particularly in Montione (Figure 4).
The Leonina badlands are located inside the protected Natura2000 area 'Crete di Camposodo e Crete di Leonina' (IT5190004), surrounded by croplands (Figure 4). Land-use history is presented and discussed in Amici et al. (2017) and Torri, Poesen, et al. (2018), . tion rates determined by Chiaverini et al. (1999), most probably just after the first black death epidemic, which was reinforced by a series of plague upsurges during the following decades. This situation, which persisted for 50-70 years, kept the countryside poorly managed by lack of peasants and allowed soil erosion to expand and to transform a fertile territory into badlands. These biancana badlands are presently just a small remnant of what they were at the beginning of the 20th century. They were mostly transformed to cropland through a series of techniques of which the most effective was land levelling using bulldozers, which ended in the 1990s. The present land use (protection of the biancana geoform and the biodiversity it hosts) is such that vegetation has encroached the eroded slopes and the bare pediments for the last 20-30 years, substantially transforming most gullies into vegetated dormant gullies (Figure 4) Hence, the GHs are currently mostly inactive. Torri et al. (2013) showed that vegetation covers was much lower in the 1950s, when the bare biancana slopes were common (LB) and the interbiancana gullies (LH) had almost bare GH. Biancana badland remains are now located in the middle of croplands and invaded by grass and bushes, such as the ruderal species among which is the Avena fatua (Chiarucci et al., 1995;Maccheriniet al., 2000). Plants such as  (Calzolari, Ristori, Busoni, et al., 1993;Calzolari, Ristori, Sparvoli, et al., 1993).
According to the CORINE Land Cover classification of 2012 the land cover in the study area consists mainly of forest and semi-natural areas with at some places coniferous or mixed forest and at other places open or sparsely vegetated areas with shrubs ( Figure 4). There is also some arable land, but it is scattered and mixed with natural vegetation (matorral). However, this has not always been the case.  (Suárez, 2008).  Martín-Moreno et al., 2014;Ternan et al., 1996). However, after reforestation, as Ballesteros Cánovas et al. (2017) and Ternan et al. (1996) noticed, a phase of renewed channel incision occurred due to a concentration of runoff and thus erosion in the main gully channel because of a declined sediment connectivity between slope and gully channel (clear water effect). This indicates that the recent land-use changes are still affecting positively badland formation.
The resulting landscape is characterized by Eutric Cambisol soils (WRB, 2015). However, due to the strongly weathered character of the soils with disintegrated pebbles and clayey material, well developed soils can also be found and classified as Ultisols or Alfisols following Soil Taxonomy (Espejo Serrano, 1985;Gutierrez-Elorza et al., 2002). Badlands usually occur in the poorest soil situation.

| Gully head selection
Within each badland site, the gully sampling locations were selected to guarantee the correct characterization and measuring of GH features. According to the definition of Hauge (1977), GHs with a crosssection perpendicular to flow direction larger or equal to about 1 ft 2 (i.e., approximately equal to 0.093 m 2 ) were selected directly in the field. This criterion was adopted because it has been largely used in many studies for standardization purposes despite the fact that there are no physical reasons for this criterion (Poesen et al., 2003). Single and still active GHs selected were (i) the identification and measuring of the GH cross-sections was straightforward, (ii) the contributing area was easy to determine, (iii) the GH vegetation cover and land use was almost uniform in the contributing area; (iv) the human disturbances (e.g., roads, drainage systems, etc.) were absent in the contributing area and in the GH proximity. At each site, the collection procedure aimed at maximizing the range of GH slope gradients and contributing areas selecting among those where the drainage area and slope gradient around the GH were more clearly detectable.  Table 3. The comparison RFC by mass (after sieving) with the visually estimated RFC in the field allowed to evaluate the reliability of field estimations.

| Gully head threshold value and Curve Number parameterization criteria
In this section, the criteria applied for the determination of both k and S 0:05 are described. The former will be always calculated as the first percentile (k 1st ) of the observed distribution of the k-values calculated as k ¼ sin γ ð ÞA 0:4 , with γ and A being measured for a set of GHs in the field. The latter, that is S 0:05 , will be selected considering that the resistance of the land-use soil system can change during the year, following seasons or generally following time and the many causes of modifications such as plants growth, symbiotic behaviour (e.g., mycorrhiza), soil biomass modification, anthropic direct interventions (e.g., tillage), grazing, changes in grazing management, to mention just a few of what is an almost never ending list of sources. The situations corresponding to the conditions for maximum hazard of GH retreat will be selected, for example, usually keeping the possible choices within the range of vegetation cover below 0.3 as this value brings soil resistance already pretty high thanks to the herbaceous root effect (De Baets et al., 2008;Thornes, 1985;. Presence of cropland will bring in the seed bed situation.
Usually, potential crop field that are difficult to protect from erosion damages are abandoned and land use shift to poor pasture or to rangeland, with frequent episodes of overgrazing, hence rangeland in poor to medium soil hydraulic conditions.
The S 0:05 is then calculated through Equations (14) and (15) Still from Branson et al. (1981), we have CN ¼ 91 À 13 V c for pasture, rangelands and annual grass ðfA1fÞ From Romero et al. (2007), Taguas et al. (2015) and further crop cover descriptions by Meerkerk et al. (2008): As an example, for a vegetation cover V C This method reduces arbitrariness once the land-use type and cover values are identified and quantified. Having the researchers aware that being wrong or scarcely defined in the evaluation of the earlier listed quantitative characteristics does reduce the certainty of correctly evaluating thresholds, then they will be more careful in their observations (lack of information on vegetation cover or on rock fragment cover was the base for the reducing to a quarter the total number of GH field data in the review paper by Torri and Poesen (2014)).

| RESULTS AND DISCUSSION
This section summarizes and discusses the results of the analysis per-    (Poesen et al., 2003). This is related to the difficulty to select a proper gully cross-section in the field, which may vary significantly along the gully channel even over a short distance. However, these small variations corresponding in general to shifts of the gully cross-section measurement locations along the channel of few centimetres, do not correspond to significant changes in GH catchment area and soil surface slope values. The slightly larger variation of the GH cross-section perimeter and average depth (Figure 7f,g) reflects the complexity and the roughness of the gully channels, which is particularly larger in badland environments.
The RFC by mass is highly variable in the different sites The correspondence between RFC visual estimated in the field and the RFC by mass values determined by sieving, were evaluated through linear fitting and the relative results are summarized in Table 5. For each subset, Table 5 reports the sample size, the fitting coefficients and their significance as well as the R 2 and its significance.

Differences between the RFC by mass values and the values visu-
ally estimated in the field are expected, since the latter is necessarily biased when the distance between the observer and the soil surface is larger than c. 50 cm and additionally conditioned by the observer eyesight quality and by visual defects. In addition, another possible explanation of such differences is related to the fact that the RFC by mass determinations, in some sampling conditions, may account not only the RFC at the soil surface but also incorporates subsurface fragments. Indeed, in general the number of RFC at the soil surface is larger than below the surface, because of erosion pavement formation. Obviously, the rock fragment cover is estimated with a surface as unit reference while the by mass content is automatically referring to a volume (or better to a mass of soil, always three dimensional).
However, in the study areas the sampling for RFC by mass determination was done collecting material from the first 10 cm of topsoil, hence being representative for the surface and topsoil condition. The R 2 results in Table 4 reveal that visual estimation better relates to RFC diameters larger than 4 mm. The correlation improves when the visual estimations are averaged over all those corresponding to the same by mass sample. In particular, for the average values the angular coefficient is the closest to one and the intercept to zero indicating a better degree of correlation. We acknowledge that the RFC determinations either measured by sieving or visually estimated in the field, are rarely described in detail in the literature. This may potentially introduce biases in the analyses. We therefore suggest to use the more objective RFC by mass values specifying the sampling depth and the sampling and measuring conditions.
The soil surface slope and catchment area values measured in the field of the different badland study areas are shown in Figure 9(a). In F I G U R E 6 Comparison of minimum (a) and maximum (b) air temperature and precipitation (c) monthly average statistics calculated for each gully head site (see Table 2). [Color figure can be viewed at wileyonlinelibrary.com] T A B L E 4 Values of the climatic indexes calculated for each gully head site in the study areas (see Table 2) Site code  With respect to the other main land-use classes, the badland data plot lower, corresponding to a larger erodibility. The partial overlap of the badlands data with the cropland confirms that tillage is a powerful mean to undermine the soil and land-use system resistance to gully erosion. Also rangeland is partially overlapping the badlands (the median line plots already below two badland GH points, suggesting that the lower 50% of the rangeland observations are mixed with the badland ones), which suggest that badlands can develop directly from this land use. The variability of k values (Equation 10) is shown for the different study sites (Figure 9b) and for the different study areas ( Figure 9c). The tangent of the slope is used rather than the sine value to better represent the k value variability for high slope gradients.  Figure 9b). The lower boundary k-values (i.e., the most probable estimators of the threshold) seem to vary within a small range over the different areas (Figure 9c).

| The grain size distribution and the soil resistance to detachment
In Figure 9(a) the tangent instead of the sine was used for representing the effects of slope angle γ ð Þ on GH thresholds. This is not the expected trend (i.e., sine γ versus area) which results from the dependence of the flow characteristics (shear stress) on the sine of the slope angle. Figure 9(a) shows that the trend for the Spanish data are rectilinear if the gradient is expressed as a tangent instead of a sine. According to what is presented earlier the friction angle component seems to dominate. Hence, one needs to compare the field data and trends with evaluations of the two components which depend on the effective friction angle (Equation 21).
There are different elements to be considered for such evaluations. The first is the particle diameter representative for the eroded soil. This must be calculated using the grain size distribution extended to include the RFC by mass sampled in the top 5-10 cm of soil. As representative diameter, we use D 50 , hence the ratio D p /D 50 = 1. This can be calculated for all our sites. Each site will be represented by the median grain diameter averaged over all the samples.
A second element to consider is the average bulk density, which depends on the soil characteristics: if the soil is aggregated, the bulk density (at saturation) is a reasonable value, if it is dispersed/flocculated because of the presence of high sodium content, the values of density may be much lower than for an undispersed clay. If the soil is fresh of F I G U R E 7 Variability of gully catchment (a, b) characteristics and gully head (GH) geometric and morphometric features (c, d, e, f, g) within different GH sites (see Table 2). [Color figure can be viewed at wileyonlinelibrary.com] weathering, its density will be very close to the density of the unaltered rock. The more the weathering the lower the density of the grains.
Values can be found in Poesen and Lavee (1994).
To evaluate the position of the threshold curves for the different study sites, we adopted the values of the parameters reported in Table 7, where the values of cohesion were used for obtaining a perfect match between the k 1st value and the one calculated with Equation (21).
Results are shown in Figure 10. Regarding the LH and LB sites, the friction component was absolutely negligible and the position and shape of the threshold curve was defined only by the cohesion value.
Note that Torri et al. (2013),  evaluated the cohesion value of the saturated material to be below 500 Pa due to the dispersive action of the exchangeable sodium. The material was transported as floccules due to the high salt concentration in the runoff (2550-5000 μS/cm) hence the low values for the Δρ in Table 7.
The relationship between the measured k 1st and R tot (i.e., the average soil resistance over the observed surface gradient range, calculated with Equation (20) and data in Table 7) is obtained using Equation (21) for all the eight study sites: where the constant 0.000015 is evaluated as the linear best fitting coefficient between k 1st and R tot 10=7 assuming a nil intercept. In other terms.
As discussed in Section 1.3, we expect that the ratio between critical fluid stresses and soil resistance values (FS ratio Þ should be included in the interval 0.0004 and 1. The estimated soil resistance R tot is not the average soil resistance, as it would have been measured using a torvane or other device, but this is rather an estimation of the resistance offered by the grains that can be detached and entrained by the flow, and hence it is smaller than the average soil resistance R soil . Consequently, we expect a FS ratio much larger than 0.0004. If we estimate the critical fluid stress (F stress ) values using the Smerdon and Beasley (1961)   F I G U R E 9 (a) Soil surface slopes (tan slope) at the gully heads and their corresponding catchment areas (a) measured in the different badland study areas (see Table 2); solid and dashed coloured lines show, respectively, the median and average threshold trends determined for cropland (dark red), rangeland (purple), and forest (green) (Torri & Poesen, 2014 Figure 11 with the solid symbols depicting the dispersive soils at the LH and LB sites. Additionally, in Figure 11, open red circles report the FS ratio calculated for the non-dispersive badland soils and the red dashed curve show the corresponding best fitting equation and the relative R 2 value. Such relation, which is significant at a 0.01 level, seems to suggest that we can rely on Smerdon and Beasley (1961) relation or on its improvements, at least for non-dispersive soils.

| Badland thresholds characterization
In order to investigate the effects of RFC on the resistance to the GH regression (k-value), the RFC by mass and RFC cover data for the GH sites have been classified into five equally spaced RFC classes and plotted against their corresponding k values ranges. Leonina data have been omitted in this analysis due to its extremely diversified lithological characteristics, being over consolidated marine dispersive silty-clay to clayey silt sediments without any type of cement, while the others are mostly cemented or partly cemented sediments with distinct layering. Figure 12 shows a slight increase of k in relation to the increase of RFC.
The trend is better shown if the RFC is by mass, where the effect of RFC tends to a horizontal asymptote. When RFC is given by visual estimation then a slight decrease is shown after a peak at about (60, 80].

| Curve Number value selection for badlands
Once the distribution of the k-values for each study site was established, a possible CN value and the derived S 0:05 (i.e., depending on the land use and vegetation) parameter for each site was identified in order to evaluate if the badland GH data do agree with the k 1st model first proposed by Torri and Poesen (2014) and then further developed by Torri, Poesen, et al. (2018). With this aim in mind, the CN values and possibly how these change with vegetation type and cover need to be defined.
For this purpose, a criterion was used that is explained and moti- Regarding the experimental sites presented here, the first 50 cm of the soil profile are usually interrupted by an impervious layer, hence these soils can be classified as hydrologic soil group D. As shown in To compare the k 1st values determined with data observed in the field with S 0:05 as proposed by Torri and Poesen (2014) and Torri, Poesen, et al. (2018), ψ factors also need to be considered because values are different from one for the different test sites.
3.4 | The ψ factor Table 8 shows the steps taken to obtain the S 0:05 parameter (Equation 17). Table 9 contains the S 0:05 , the ψ factor and the observed k 1st (Equation 17). The ψ factor has so far not yet been discussed. The ψ factor accounts for the local effects of processes that may be important in a site but absent in another. For instance, such effects can be related to the presence of joints or faults (ψ ¼ 0:76 T A B L E 7 Data used as input for calculating the position of the threshold curve: The values of cohesion were used for obtaining a perfect match between the k 1st value and the one calculated with Equation (21). When cohesion was nil (i.e., when the amount of clay was lower than 0.05) then the ρ-value of the representative grain was used for the fine tuning of the two k-values. The shape 'rounded' was chosen when the material was eroded mostly as soil aggregates or sand grains (smallest D 50 size) from Torri and Poesen [2014]). In addition, there is another correction factor for the bare slopes of the biancanas in the Leonina site (LB) which was derived by laboratory experiments discussed by Torri et al. (1994). This factor accounts for the over-consolidation of silty clay marine sediment biancana deposits, which are characterized by a very high resistance to erosion when dry but which drops to almost zero when thoroughly wet. The boundary line between the erodible phase (topsoil) and the still too resistant subsoil layer, moves slowly into this sediment at a very low velocity, which corresponds to the velocity of the water front penetration (u wfp ). This rate is smaller than the rate at which the erodible layer is potentially eroded, that is its potential denudation rate (d pr ). When the erosion is driven by raindrop impact and sheet flow dominates, d pr < u wfp , but whenever concentrate flow erosion occurs, it cannot erode at a rate that is larger than u wfp and the flow detachment is substantially limited.
A similar situation is observed in the two sites VR and MN, with VR characterized by a deeper surface alterite layer than the MN site due to peculiar geological conditions. These two sites are characterized by similar lithologies, predominantly constituted of regular alternations of sandstone and pelitic marls (i.e., siltstone and mudstone, see T A B L E 9 Final values of Curve Number (CN) and S 0.05 calculated for the different gully head (GH) sites obtained with row averages of values in shaded columns corresponding to the vegetation characteristics (V c ) of the corresponding GH sites (see Table 8). The ψ factor values account for various factors adding or reducing resistance to GH formation, for example in presence of joints ψ = 0.76 (Equation (17), Torri & Poesen, 2014), and other local characteristics, such as the water penetration front rate as limiting factor for erosion to proceed, that is ψ =3.5 (Torri et al., 1994(Torri et al., , 2018a Final averaged parameter values (for Equations (25) and (26)) b based on selected values [values with asterisk (*) in Table 8 3.5 | The land-use interpolating function (S 0:05 ) Torri, Poesen, et al. (2018) suggested a slightly different equation than the one suggested previously (Torri & Poesen, 2014) for including the land-use effect using S 0:05 . As the badlands k 1st data are among the smallest values recorded, these can be compared and two interpolating functions identified and evaluated. The two equations, ignoring RFC effects (RFC eff = 1), become: and k 1st ¼ 0:00113ψ 1 À exp À0:0137S 0:05 ð Þ ½ S 0:05 ð 26Þ with k 1st in ha À0.4 , S 0:05 in millimetres and ψ the corrector factors reported in Table 8.
The k 1st values calculated using Equation (26) are reported in Table 7. Given the importance of the ψ factor values, the ratio k1st ψ is plotted versus S 0:05 in Figure 14. Figure 14 shows how the data fit the two curves (Equations 25 and 26). The only datum that clashes with the Torri, Poesen, et al. (2018) curve is the one indicated with a green cross overlapping a small red circle. If we accept this curve then this F I G U R E 1 3 Badland development in relation to layer bedding orientation in the Verghereto (a,b,c,d) and in the Montione (e,f,g,h) study sites. [Color figure can be viewed at wileyonlinelibrary.com] datum should be considered as an outlier, either because there was an error in reporting it or due to some local peculiarities. Hence, the Nazari Samani et al. (2009) data (green crosses in Figure 14) were reexamined. These reported data are for cropland (a very small group) and rangeland. The latter were not used by Torri and Poesen (2014), but they now are included in this analysis. For this purpose, data reported in the original article were digitized and the k 1st values calculated for rangeland and cropland, the latter being in total agreement with Torri and Poesen (2014) data. The values of S 0:05 were calculated following the methodology proposed in the material and methods section, obtaining a different value for the S 0:05 for cropland. The difference between the two evaluations of the CN and hence S 0:05 values for cropland, was not due to the methodology. It is related mainly to another paper by Nazari Samani et al. (2016), where they conducted runoff experiments on plots for cropland, abandoned, and rangeland soils in the same locality. The photographs showed that the cropland soil was bare, crusted and probably also with a biological crust (lichens or cyanobacteria). This pushes the CN values further towards the largest values, and consequently S 0:05 towards very low values. The two recalculated data points are shown in Figure 14 with green crosses. Now also the data reported by Nazari Samani et al. (2009) along with the other data plots very close to the Torri, Poesen, et al. (2018) interpolating curve and support the finding that the Torri and Poesen (2014) Equation (25) must be abandoned.

| Insight on rock fragment content effects
The data reported in Table 9 relative to badland sites, allow to calculate the effect of rock fragments on the ratios (k ratio ) between measured and predicted k 1st values. The comparison of the predicted k 1st values without RFC effect with those including a correction for RFC are shown Figure 14. For this comparison, the k ratio values were obtained considering Equations (26) and (18). If the k ratio equals 1.0 then there is no RFC effect. If the ratio shows an increasing trend with RFC then there is an RFC effect and the function interpolating the data can be multiplied by the k calculated with Equation (26) to obtain the best estimate of the measured k 1st .
Note that the ratio is one at an RFC larger than zero because the k 1st reflects the average rock fragment cover over the 55 gully sites used in Figure 14 (average RFC between 0.1 and 0.3).
As shown in Figure 15, with increasing RFC values, predicted k 1st values tend to underestimate the measured k 1st values. Overall, the general behaviour of kratio suggests an increase of k 1st with RFC. Anyhow, the number and scatter of data do not allow to calculate a reliable model for accounting the effect of RFC.
Most of the variance regarding the threshold model (Equation 26) is explained by the land use as quantified by the CN-S 0.05 runoff model hence only a relatively small data scatter still remains to be characterized and explained. It is well known that weather and climate data as well as the RFC have an impact on the threshold conditions. This means that it would be better to examine weather/climate and RFC at the same time in order to extract these factors. To do that a much larger data base then the one relative to badlands is obviously needed.

| CONCLUSIONS
The data presented and discussed in this article extend the GH topographic and land-use thresholds to Mediterranean badlands, whose data align to those discussed by Torri and Poesen (2014).
A strong accent is here posed on the repeatability and reproducibility (RR) of the parameters needed for calculating the gully threshold curves and for characterizing the GH data bases. To ensure this issue, the procedure to assign a CN value to a given representative F I G U R E 1 4 Ratio k1st ψ versus S 0:05 for badlands (triangles). The green crosses depict the same data for cropland (lower value) and rangeland (Nazari Samani et al., 2009), while the red circles are data reported in Torri and Poesen (2014). Solid green line and dashed orange line, show respectively the interpolating fitting Equations (25) and (26). [Color figure can be viewed at wileyonlinelibrary.com] F I G U R E 1 5 The ratio (k ratio ) between the k 1st observed and the k 1st calculated as a function of S 0.05 alone is plotted versus the rock fragment content (RFC) by mass at the studied badland sites. Calculated k 1st values were obtained using Equation (26). [Color figure can be viewed at wileyonlinelibrary.com] land use was improved. It is obvious that under any possible condition, when runoff concentrates, a rill or gully will be initiated when and where the soil/land-use system is the weakest (e.g., scarce vegetation protection, presence of various disturbances, water concentration paths).
Hence a land-use type should also be characterized by its weak spots such as animal paths, machinery trails, bare soil patches and not just through its average situation. This results in a selection of CN values representing the poorest hydraulic conditions. A set of equations, linking vegetation cover to CN was collected from the literature for characterizing the conditions at the boundary between a safe, well protected soil situation and a totally bare soil surface, under different types of vegetation. This CN estimation approach reduces arbitrariness and increase the RR of the evaluation of the land-use effect on GH thresholds.
The study also re-evaluated the effects of rock fragments proposing the use of a by mass rather than a visual estimation of rock fragment cover at the soil surface by sampling the top soil over a depth of 5 to 10 cm (depth depending on the grain size of the surface: if rock fragments are less than 5 cm in size then a soil sample of a few centimetres thick is sufficient; for larger rock fragment sizes a soil sample taken over a larger depth is needed).
The local factors, inherited by the geological and the orogenetic history (here mainly joints and local structural attributes) have been shown to affects the threshold curve. These effects still need further studies and confirmation to be safely used in the threshold prediction.
Another aspect that has been addressed is the importance of the history of the area under study: what one observes today is a landscape that was formed under the current conditions, or it may be a legacy of the past centuries or even millennia. The analysis of historical documents provides a better view of the situation under which the gullies at the study sites formed. The analysis of past geoenvironmental records (e.g., paleopalynological data, dendrochronological series, radionuclide data, etc.) may complement such investigations. The conditions for the formation of the currently dormant GHs at the LH study area is the result of past local conditions similar to the ones that the LB study area was experiencing until 40-50 years ago. Similar considerations can be made for Verghereto/Montione and for the Guadalajara sites. In many of these study sites one cannot use the present-day conditions and obtain theoretical k-values in agreement with the measured k. History also shows that the cause for these badland sites to develop had mainly anthropogenic drivers while present conditions testify a today abandonment (which is not necessarily the cause of GH formation).
The effect of the rock fragments component on the GH threshold model was confirmed: that is with increasing rock fragment content by mass, the resistance to gully formation increases. Similar confirmations were proposed for the presence of joint networks. Another probable cause of differentiation in threshold values was found due to the layer bedding orientation: while Verghereto badland gullies are almost all formed into a given bedding orientation, the Montione site is characterized by gullies that developed on the opposite bedding.
The GH topographical threshold model gains its present form following studies by Patton and Schumm (1975) and the generalization proposed by Montgomery and Dietrich (1994). The Torri and Poesen (2014) literature review on the Hortonian runoff generated gullies brought the model to a better clarification and included the effects of land use to expand the model from a topography-based to a topography and land-use based one. The often ignored Rossi et al. (2015) paper warns about the bias of the classical approach to calculating threshold lines. The data collected in this study required an improvement of the model. It was necessary to express the critical flow shear stress explicitly. This was done using a classical soil strength concept with a component depending on a friction term and another on a cohesive term. When cohesion is dominant the topographic threshold line is straight in a (A, sin γ) bilog representation of the GH conditions. When friction is dominant the threshold line bends at slope gradients above 60 . Hence, the condition for GH initiation and retreat is the result of the tradeoff between the frictional and cohesive components of the resistance forces. Additional studies are needed to completely exploit this new modelling perspective.