Gully Erosion Development in Drainage Basins: A New Morphometric Approach

Abstract

The formation and evolution of management gullies is a highly intense process of soil erosion often overlooked in policies and river basin strategies. Despite the worldwide spread of the phenomenon, our ability to assess and simulate gullying and its impacts remains limited; therefore, predicting the development and evolution of these river reaches represents a significant challenge, especially in areas where the loss of productive soil or the hazards linked to landslides or floods represent critical factors. Our study demonstrates how an exclusively morphometric approach, based on the construction of the hypsometric curve and applied to small hydrographic basins that are lithologically homogeneous and hierarchized according to the Strahler classification method, is able to predict the triggering height of the gullies; this height corresponds to the mean elevation of the basin and the inflection point of the hypsometric curve itself, confirming the hypothesis that this point coincides with the point at which a sudden change in surface runoff energy occurs, The study also shows that the portion of the basin necessary to trigger these intense erosive processes is always within a small range, regardless of the size and morphology (slope) of the basin itself. Such an approach, which is quick and relatively easy to apply, could help develop hydrogeological hazard mitigation practices in land planning projects.

1. Introduction

Linear erosion caused by concentrated overland flow has been recognized as the primary cause of land degradation worldwide. Gully erosion, in particular, is an intense process of soil erosion that accounts for more than 90% of the watershed sediment load in running waters [1,2,3,4,5,6,7,8,9]. The related geomorphological features are either rills or gullies, with rills corresponding to the smaller features with incision depths below 0.3 m and gullies corresponding to features between 0.3 and 2 m (shallow gullies) or higher than 2 m (deep gullies) [2,10,11,12].

It is also known that even more intensive agricultural practices in the upper parts of river catchments lead to increased sediment input (and thus to increased nutrient loading in watercourses, which reduces the water quality of rivers and streams). The effects of sediment accumulation associated with anthropogenic structures, such as artificial reservoirs or check dams, cannot be overlooked, as they have serious consequences for water resources used for drinking or hydropower purposes or for flood hazards. Despite their importance, these processes have been little studied for a long time, as it is difficult to develop reliable models for their formation and evolution; the latter is mainly due to the variety of geomorphological features and processes involved [13,14,15,16,17].

Morphometric analysis represents one of the first quantitative approaches to understanding a river basin’s morphological and hydrological characteristics. Drainage morphometry studies, defined as the measurement and mathematical analysis of linear, areal, and relief characteristics of any drainage basin [18], were first initiated by Horton [19] and then followed by other studies, which are still the basis of modern research in this field today [20,21,22,23,24,25]. Subsequently, most studies have used simple empirical equations to analyze the morphometric and hydraulic characteristics of gullies, although many of these empirical equations are derived from typical relationships of river dynamics [26,27,28,29,30,31,32].

Many recent applications of morphometry in the study of slope erosion have mainly addressed identifying critical thresholds in terms of rainfall, topography, land use, and hydrology for the initiation, deepening, and evolution of gullies in different environments. Concerning topography, interesting results have been obtained from several studies on the relationship between the local slope gradient (S) and drainage area (A), which are considered fundamental factors in establishing the position of gully heads in river basins [8,10,33,34,35,36].

The present study, with an original approach, allows us to mark the transition between sectors of the basin characterized by the different erosive capacity of running water, confirming the intuition by Strahler [23,37], according to which the inflection point of the hypsometric curve would have a significant geomorphological value; in contrast to the above studies, the only relevant factor for predicting the relative height of activation of the main erosional processes becomes the value of the hypsometric integral and, consequently, the degree of geomorphologic evolution of the basin itself. Besides, this approach allowed us an almost automatic estimate of the upstream areas necessary to collect surface runoff capable of activating intense gullying processes, regardless of slope gradient or land use.

Using simple tools implemented in a GIS environment and ordering the stream reaches according to Strahler [38], we aim to demonstrate that the maximum concentration of gully heads, in the absence of significant tectonic–structural or anthropic conditioning, can be related to the mean elevation of the drainage basin, where the largest amount of stream reaches of a given order (depending on the resolution of the DEM used for the morphometric analysis) take their origin.

2. Materials and Methods

We tested our approach on 25 sample basins in two test areas: one in central Italy and the other in the central-western United States (Wyoming), using high-resolution satellite images to ensure matching with field evidence. Different climatic, morphometric, geomorphological, and land-use-related features characterize these catchment areas, as does the varying degree of anthropization; the only similarities are the less permeable bedrock, which allows the formation of a well-developed hydrographic network and the low degree of tectonic conditioning.

2.1. Study Areas

2.1.1. Adriatic Side of Central Italy

The study area includes a sector of approximately 1200 km² located on the Adriatic side of central Italy (Figure 1a). This sector is characterized by a predominantly high-hilly and hilly landscape with altitudes (generally lower than 500–600 m a.s.l.), which progressively slope towards the east and the coastal belt. Land use on the slopes is mainly agricultural, while small urban settlements and inhabited centers are primarily located on the watersheds. The bedrock, substantially homogeneous, is characterized by the presence of formations of Plio-Pleistocene age made up of alternating shales and sandy-conglomeratic levels [39,40]; these formations are locally covered by medium-fine continental deposits (colluvial and/or fluvial). The structural setting of the Plio-Pleistocene formations is characterized, on the surface, by a regular, east-dipping monocline (between 18° and 20° in the extreme western area and 5° in the coastal area), resulting from the intense tectonic uplift started in the early Pleistocene. As a whole, the monoclinal structure is displaced by rare dip-slip faults, mainly oriented NNW–SSE and WSW–ENE, the displacement of which rarely exceeds 10 m [41]. Micro- and meso-structural analyses on middle Pliocene and upper Pleistocene formations highlighted intense joint systems compatible with the abovementioned faults. These tectonic elements, however, do not produce significant morphological evidence on the surface or a particular conditioning of the hydrographic network.

From a climatic point of view, this area is a transition point between the Mediterranean climate, typical of the southern part of Italy, and the humid subtropical, or temperate, oceanic climate, widespread in most parts of the territory [42,43]. The study area shows average temperatures of almost 14 °C, with the highest values during July and August and the lowest in January. A significant increase of 0.5 °C has been observed in the last 30 years due to climate change; this has allowed the coastal areas and nearby territories to reach, and in some cases exceed, 16 °C of average annual temperature. The precipitation height, on the other hand, is between 600 mm along the coast and 1000 mm in the high-hilly sector; the relationship between temperature and rainfall values shows, from July to half of September, slight summer aridity in the low-hilly area and along the coast. The Fournier Erosion Index evidences a low erosive potential for rainfalls, except during extreme meteoric events. In these conditions, gully erosion processes are pretty common [44,45,46,47], although dense perifluvial vegetation (Figure 1b) can sometimes hide how active they are. Conversely, this vegetation makes it easy to see the transition between these processes and the parts of the slope where widespread rill erosion occurs.

2.1.2. Gillette Area, Wyoming (USA)

The second study area includes a sector of approximately 2100 km² within the Gillette coalfield in NE Wyoming (Figure 2a). The landscape of the area is characterized by the presence of a vast highland with smooth ridges and limited height differences, with altitudes ranging between 1400 and 1500 m a.s.l. The degree of anthropization in the area is extremely low and limited to single buildings located mainly along the plains; the slopes are mostly uncultivated, with little vegetation, mainly shrubby, situated within the river incisions. Also, the bedrock in this area is substantially homogeneous and mainly made up of Paleocene and Eocene siltstones and shales, locally covered by altered silty-clay bedrock (mostly on tops of hills and ridges) or fine colluvium [48,49]. From a stratigraphic point of view, the formations above show a very regular structure with almost continuous and horizontal layers and practically no tectonic conditioning.

Because of its elevation, this sector of Wyoming has a relatively cool and arid climate, with annual high temperatures of around 15 °C (and a maximum during July and August) and low yearly temperatures of around 1 °C (with a minimum in December and January). On the other hand, mountain ridges that block the humid air masses coming from the Pacific Ocean generally determine low average annual rainfalls (between 400 and 450 mm), mostly concentrated in spring and early summer (US Climate Data at usclimatedata.com). In this sector, the gully processes are very evident due to the poor protection of the soil and the almost total absence of tall vegetation (Figure 2b).

2.2. Morphometric Analysis

We carried out the morphometric analysis of the individual basins in a GIS environment using the ArcGIS Pro software (ESRI, v.3.3.1); in contrast, for the graphic representation of the results, we employed a simple spreadsheet (Microsoft Excel v. 365). More specifically, we structured the analysis into four different steps:

• Step 1: Tracing of the drainage network, definition of the watercourse orders, and calculation of the basic morphometric parameters by GIS procedures.
• Step 2: Construction of the frequency distribution histogram (FDH).
• Step 3: Construction of the hypsometric curve and the flow contribution curve (FCC).
• Step 4: Detection of gullying initiation areas by satellite image interpretation.

Step 1: Tracing of the drainage network, definition of the watercourse orders, and calculation of basic morphometric parameters by GIS procedures

We performed the morphometric analysis of the Italian and US basins using 10 m-resolution DTMs, freely downloadable from the links https://doi.org/10.13127/tinitaly/1.1 [50] and https://portal.opentopography.org/, respectively. For analyzing the selected basins, whose size varied between 1.15 km² and 37.5 km², we used 10 m resolution DTMs.

Choosing a DTM with a suitable resolution is a key issue. Nowadays, global databases with a relatively high resolution are available (e.g., ~10–30 m), and the local/regional sources can even surpass them (e.g., LiDAR). Although our basic goal was to use DTMs with the best possible resolution, we also wanted to develop a processing method that could be used for as many different databases as possible. Many studies in the literature have shown that a DTM with a resolution of 10 m is an excellent choice for morphometric analyses [51,52,53]. Furthermore, in this work, we avoided using LiDAR data with a resolution of 1 m because in partially anthropized areas, “false” elevation differences can often occur, related, for example, to the presence of roads or paths; such “anomalies” can then lead to errors in the tracking of the “flow direction” by GIS methods.

After we had obtained the DTMs and “improved” them by clipping the basin outlines, we removed all sinks with the Fill tool (Spatial Analyst) to obtain “hydrologically” correct DTMs. Following, we calculated some basic morphometric parameters, such as area, perimeter, average slope, maximum and minimum elevation, and circularity ratio, which is the ratio between the area of the single basin and that of the circle with the same perimeter [54]. To define the hydrographic network of each basin, we chose the Strahler [38] classification method, a “top-down” system where the first orders are the outermost tributaries. Following this method, if two reaches of the same order merge, the resulting reach will have a number that is one higher, while if two reaches with different stream orders merge, the resulting stream will have the higher of the two numbers. Subsequently, starting from the reconditioned DTM, we created, in a GIS environment, a Flow direction raster and a following Flow accumulation raster, representing the weight of all cells that flow into each downslope cell. Using a conditional (if/else) evaluation and a threshold of 10 cells (compatible with the DTM resolution and, as will be shown later, irrelevant to the final result in any case), we made a new raster with a specific stream definition. Next, we used the relative ArcGIS Stream Order tool to create a stream order raster. Finally, we used the Stream to Feature tool to make a polyline feature showing all the river reaches grouped by order (Figure 3a,b). As previously mentioned, we used the same procedure on a DTM with 1 m resolution (LiDAR, Figure 3c); however, we could not correctly show the hierarchization of the hydrographic network with such fine details. Therefore, we chose to use a DTM with a lower resolution that was still compatible with the average size of the basins we looked at.

Step 2: Construction of the Frequency Distribution Histogram (FDH)

After designing the hydrographic network, we obtained, for each hierarchical order, the frequency distribution histogram (FDH), which describes the number of river segments arranged along different elevations. This number, coded in this study with the symbol Nuz, expresses the number of Nu sections intercepted at a predefined height, z (Figure 4a). We obtained this parameter by counting the points of intersection between the vector file of the hydrographic network created in the previous step and the vector file of the contour lines extracted from the DTM (Figure 4b); for consistency with the quantities determined previously, we chose the equidistance between the contour lines to be equal to 10 m.

Similarly, we could also represent the number of reaches of each hierarchical order (Nuz_I for the first-order reaches, Nuz_II for the second-order, and so on). These representations allowed us to verify the arrangement and, specifically, the number of reaches for each hierarchical order as a function of their elevation (represented along the x-axis) (Figure 4a).

For each basin, only one function describes the structure of the relative hydrographic network along the vertical. In the simplest case of first-order basins, only one absolute minimum falls on the watershed (Nuz_I = 0).

Furthermore, the resulting graphic solution highlights that in correspondence with the peak values, the elevation belts within which the erosive action of water flow was particularly effective.

Step 3: Construction of the Hypsometric Curve and the Flow Contribution Curve (FCC)

To estimate the mean elevation of the individual basins and evaluate the mean height of the starting points of the reaches for each hierarchical order, we constructed the hypsometric curve [23], which represents the fraction of the area of a river basin above a given elevation. These curves are presented in non-dimensional form (i.e., plots of the relative height (h/H) versus the relative area (a/A) of the region under investigation), which enables the comparison between different hydrographic basins. For the specific objectives of this work and to make the correlations existing between the relative elevations of the basins analyzed more evident, we created the hypsometric curve in a semi-dimensional form with the “normalized” basin area values (a/A between 0 and 1, Figure 5) using a specific ArcGIS tool called “Calhypso” [55]. Also, the hypsometric curve equation, expressed with a 3rd-degree polynomial function, shows two particular points (oblique inflection points—OIPs), one corresponding to the change of the curve slope at the upper parts of the basin and the other placed outside of the diagram area (Figure 5).

The inflection point is an essential concept in differential calculus and represents the point at which the curve of a function changes its concavity; mathematically, it corresponds to the point at which the second derivative of this function changes sign, with a positive second derivative indicating a concave upward curve and the negative second derivative a downward concave curve. In real-world applications, where a system is modeled using a curve, finding the inflection point is critical to anticipating the system’s behavior. In the specific case of a hydrographic basin, Strahler [23] highlighted that this change in concavity marks the level at which the rate of decrease of mass upwards changes from an increasingly rapid rate of decline to a diminishing rate of decline; in other words, this point represents the transition from an area of greater to one of lower energy on the slope. The second OIP, on the other hand, represents only one of the mathematical solutions of the derivative function but does not correspond to any real location within the river basin; in altimetric terms, it would be positioned at a lower altitude than the outlet of the basin itself (Figure 5).

Subsequently, we compared the semi-dimensional hypsometric curve with the FDH obtained in the previous step; for this purpose, we normalized the Nuz values of each hierarchical order to their maximum value, obtaining the flow contribution curve (FCC) (Figure 6).

The curve in the graph represents the number of reaches of first, second, and third order as a function of the elevation. We normalized this number to 1 to make it comparable with the hypsographic curve (light green in the graph). Moreover, it was possible to identify the highest altimetric belt, where the greatest number of hierarchically classified reaches of a given order (in our case, the third order reaches) drain a greater basin area.

Step 4: Detection of gully erosion by satellite image interpretation

Finally, we carried out a detailed analysis of the selected basins to estimate the average altitude of gullying initiation areas using high-resolution (0.3–0.6 m) satellite images (source: ESRI, Maxar, Earthstar Geographics, and the GIS User Community) (Figure 7a,b).

In the case of the US basins, the recognition of the trigger elevation of the gullies was rather simple, taking into account the complete absence of vegetation. In the case of the Italian basins, since they are areas characterized by a slight degree of anthropization, the recognition of the trigger elevation, when not clearly visible, has been hypothesized to coincide with the presence of perifluvial vegetation, which tends to colonize incisions along the slope which cannot be erased by typical agricultural practices.

3. Results

We reported all the results from the processing above in a single summary table (

Table 1. A summary table with all the processing results described in the text.

BASIN	Area (km2)	Perimeter (km)	Average Slope	Hmax (m)	Hmin (m)	Δh (m)	Circularity Ratio	Hypsographic Curve Function	OIP (m)	Hypso_Int	Hmean_gullie Field (m)	Hmean_III Start (m)	Hmean_III FCC (m)	Hmean (m)	Std. Dev. σ	Amean_III Contr (m2)
Basin 1 ITALY	6.18	14.46	0.30	377.00	114.00	263.00	0.372	y = −474.13x^3 + 696.78x^2 − 435.58x + 332.08	230.16	0.431	230.50	223.25	250.00	227.35	9.25	13,406.00
Basin 2 ITALY	2.65	8.65	0.30	304.00	78.00	226.00	0.445	y = −252.79x^3 + 396.49x^2 − 345.44x + 289.78	181.40	0.473	196.50	176.07	190.00	184.90	7.03	15,008.00
Basin 3 ITALY	2.91	11.62	0.29	368.00	47.00	321.00	0.271	y = −236.06x^3 + 378.18x^2 − 436.5x + 353.27	192.07	0.483	195.91	171.00	170.00	202.04	13.22	10,656.00
Basin 4 ITALY	2.68	9.25	0.28	283.00	68.00	215.00	0.394	y = −290.14x^3 + 441x^2 − 332.68x + 258.58	165.45	0.459	177.20	156.67	170.00	166.69	6.66	12,666.00
Basin 5 ITALY	14.63	22.56	0.23	331.00	49.00	282.00	0.361	y = −419.28x^3 + 689.84x^2 − 492.66x + 283.28	151.46	0.401	169.12	157.27	140.00	162.08	9.87	11,575.00
Basin 6 ITALY	6.56	14.60	0.29	432.00	151.00	281.00	0.387	y = −277.06x^3 + 371.65x^2 − 328.99x + 397.73	300.18	0.487	305.15	276.10	250.00	287.85	19.71	13,369.00
Basin 7 ITALY	16.48	24.42	0.27	414.00	118.00	296.00	0.347	y = −387.87x^3 + 628.63x^2 − 487.95x + 368.95	227.69	0.402	239.11	224.45	230.00	236.99	5.56	14,587.00
Basin 8 ITALY	33.83	37.50	0.25	453.00	58.00	395.00	0.303	y = −713.95x^3 + 1207.5x^2 − 793.84x + 372.03	180.31	0.355	196.15	187.58	150.00	198.23	17.44	11,974.00
Basin 9 ITALY	6.09	13.42	0.34	496.00	68.00	428.00	0.426	y = −569.25x^3 + 993.29x^2 − 786.07x + 434.41	201.14	0.378	245.30	205.84	190.00	229.78	20.18	14,443.00
Basin 10 ITALY	6.88	14.71	0.31	333.00	67.00	266.00	0.400	y = −455.64x^3 + 713.34x^2 − 471.06x + 289.79	173.46	0.417	185.40	157.90	210.00	177.92	17.09	15,261.00
Basin 11 ITALY	1.36	6.85	0.23	435.00	171.00	264.00	0.364	y = −292x^3 + 464.27x^2 − 407.38x + 417.52	288.55	0.472	284.64	298.74	320.00	295.61	12.30	12,201.00
Basin 12 ITALY	2.61	9.30	0.28	340.00	143.00	197.00	0.380	y = −300.77x^3 + 486.83x^2 − 359.21x + 321.95	222.57	0.439	193.81	215.24	180.00	229.48	18.50	13,952.00
Basin 13 ITALY	2.94	12.22	0.21	364.00	119.00	245.00	0.247	y = −407.48x^3 + 620.58x^2 - 436.85x + 345.53	230.34	0.461	244.94	225.28	240.00	231.95	7.04	15,123.00
Basin 14 ITALY	5.17	13.14	0.21	312.00	105.00	207.00	0.376	y = −288.65x^3 + 476.31x^2 − 352.24x + 272.39	174.71	0.376	171.35	172.42	190.00	182.83	7.12	14,883.00
Basin 1 USA	2.49	7.80	0.14	1527.00	1416.00	111.00	0.514	y = −160.06x^3 + 200.51x^2 − 133.42x + 1515	1482.57	0.529	1485.21	1470.92	1470.00	1474.72	6.15	12,807.00
Basin 2 USA	8.79	13.87	0.14	1551.00	1440.00	111.00	0.574	y = −166.77x^3 + 235.34x^2 − 165.71x + 1548.2	1504.98	0.559	1512.44	1497.12	1490.00	1502.05	7.53	18,468.00
Basin 3 USA	3.13	7.89	0.14	1494.00	1400.00	94.00	0.631	y = −111.84x^3 + 160.36x^2 − 126.23x + 1483.4	1447.48	0.481	1453.69	1441.32	1440.00	1445.21	4.87	16,462.00
Basin 4 USA	6.88	11.80	0.14	1587.00	1460.00	127.00	0.622	y = −162.46x^3 + 250.72x^2 − 185.27x + 1566	1514.94	0.443	1523.71	1513.24	1510.00	1516.26	4.55	15,345.00
Basin 5 USA	1.15	4.38	0.14	1586.00	1492.00	94.00	0.753	y = −188.45x^3 + 302.21x^2 − 198.28x + 1581.2	1532.76	0.465	1544.26	1530.62	1550.00	1535.71	7.32	16,146.00
Basin 6 USA	6.86	12.04	0.14	1502.00	1395.00	107.00	0.595	y = −167.37x^3 + 237.75x^2 − 160.21x + 1494.4	1454.10	0.530	1449.98	1442.76	1450.00	1451.71	3.79	16,798.00
Basin 7 USA	5.56	10.54	0.11	1483.00	1390.00	93.00	0.630	y = −122.32x^3 + 160.9x^2 − 111.8x + 1478.4	1450.01	0.598	1452.55	1441.60	1450.00	1445.61	3.88	15,529.00
Basin 8 USA	3.86	10.51	0.15	1459.00	1338.00	121.00	0.439	y = −97.791x^3 + 126.93x^2 − 127.11x + 1440.9	1401.72	0.472	1400.84	1391.17	1400.00	1395.11	4.02	16,187.00
Basin 9 USA	1.27	5.96	0.22	1454.00	1347.00	107.00	0.448	y = −120.92x^3 + 122.01x^2 − 104.26x + 1453	1427.15	0.601	1419.11	1401.67	1400.00	1411.31	10.31	14,992.00
Basin 10 USA	4.18	9.21	0.07	1531.00	1471.00	60.00	0.619	y = −70.996x^3 + 107.49x^2 − 88.141x + 1527.4	1501.15	0.506	1511.69	1497.61	1500.00	1501.36	4.85	15,596.00
Basin 11 USA	2.46	6.93	0.09	1529.00	1465.00	64.00	0.644	y = −50.376x^3 + 64.438x^2 − 67.169x + 1523.3	1502.48	0.526	1504.36	1496.82	1500.00	1498.66	2.68	16,104.00

) and verified the existence of possible correlations between the elevations obtained with different methods. The following paragraph will show the results with significant examples for synthesis. The complete list of the documents is provided in the attachment to this work (Supplementary Material).

As mentioned in the previous chapter, the fourth step of this work (Step 4) made it possible to obtain a complete mapping of the streams classifiable as “gullies” in all 25 basins analyzed. Subsequently, we calculated the average maximum elevations (trigger heights) of gully heads for each basin (Hmean_gullies field in Table 1).

The results were initially compared with those derived from the basic morphometric analysis conducted in a GIS environment (Step 1), thus allowing us to observe in the analyzed basins a fair amount of variability in terms of area size (between 4.38 km² and 37.50 km²), average slope (between 7% and 34%), maximum and minimum elevation and height difference (this last was between 60 and 428 m), and circularity ratio (between 0.247 and 0.753). This confirms the basins’ representativeness and the results’ non-dependence on specific morphometric conditions. Subsequently (Step 2), through geoprocessing operations, we reconstructed the hydrographic network of each basin, and the hierarchical order was defined (generally between the fourth and fifth orders) through the methodology proposed by Strahler [38].

By comparison with the mapping of gullies obtained in Step 4, we could observe a clear correspondence between the latter and the third-order reaches in all the basins considered; consequently, we focused the morphometric analysis on these stream reaches. First of all, similarly to what we did in the interpretation of the satellite images, we calculated the average of the maximum elevations of the third-order reaches (Hmean_III start) and the average of the feeding areas of every single reach (Amean_III contr) the latter, understood as the sub-basin with outlet coinciding with the head of the reach itself (Table 1). Once we set up the hydrographic network, we constructed a semi-dimensional hypsometric curve (with the elevation along the y-axis and the areas normalized to 1 along the x-axis). Then, we compared this curve with the FCCs relating to the third-order reaches nNuz_III. Solely for comparison, the FCCs relating to the first- and second-order reaches are shown in the graphs (i.e., Figure 6 and Supplementary Material).

From each graph, we obtained the altimetric belt with the greatest number of hierarchically ordered reaches (Hmean_III FCC).

Finally, Step 3 allowed us to figure out the trend of the hypsometric curve for each basin and calculate the relative hypsometric integral (Hypso_Int, ranging between 0.355 and 0.601), the mean elevation of the basin (Hmean) and, deriving the mathematical function that describes the trend of each single curve (Hypso_curve function in Table 1), the value of the oblique inflection point (OIP) located at highest elevation in the curve.

4. Discussion

The analyses above allowed us to obtain significant indications of the evolution of gullies and identify, through simple morphometric analysis, the areas of their first activation. This approach can be used on basins that are very different in terms of size and shape (mean elevation, difference in height, slope, etc.) but are similar in terms of hydrogeomorphology (well-developed hydrographic network), low tectonic-structural conditioning, and reduced environmental changes by anthropic activity); as far as land use is concerned, the method used shows that it can be effectively applied in areas with low soil resistance and low hydraulic roughness (e.g., the K coefficient described in [34,56]), as characterized by the presence of bare soils, grasslands or arable lands.

The most evident aspect is the high correlation between the “elevations” calculated with different methods; this correlation, verified through Pearson’s correlation coefficient “r”, is always higher than 0.99 (Figure 8a) as well as the probability value p which always shows extremely small values (Figure 8b). Even the standard deviation σ calculated for the “heights” of each basin is consistently below or around 10m (acceptable error for the resolution of the DTM chosen). Only five cases among the basins selected in the Italian territory show a standard deviation between 15 and 20 (Table 1 and Figure 8c,d).

Concerning the Italian basins, the causes of this “anomaly” can be related in three cases with the presence of badlands (basins 6, 8, and 9) and two cases (basins 10 and 12) with a slightly higher degree of anthropization; in all cases, as observed by various authors [12,34,57], these factors may have influenced the development and evolution of the initial portions of the hydrographic network and, consequently, the elevation of the starting points of the gullies (Figure 9a–e).

A practically zero degree of anthropization, among other things, is undoubtedly the basis of the systematically lower standard deviation values found in the US basins, where, in only one case, the value of 10 is reached (Table 1); the presence of this “external” factor (compared to morphological and climatic ones) can influence the natural evolution of a hydrographic network. The second aspect to consider is the height of activation of the gully erosion processes, verified, as mentioned, by analyzing satellite images. The fact that this elevation corresponds precisely to the height of activation of third-order reaches represents a secondary factor linked exclusively to the resolution of the DTM chosen; this hierarchical order, however, is constant in the basins of the two test areas and appears independent from the general morphometric parameters of the basins (area, perimeter, shape, slope, etc.). Concerning this elevation, it is interesting to note that it can also be extrapolated from the FCCs of each basin and generally corresponds to the highest elevation with the greatest number of reaches (in our case, those of third order; Figure 6). Another interesting aspect is the correspondence between these relative elevations and the mean elevation of the basin (calculated with GIS tools).

The explanation for these analogies is that gullies undoubtedly represent the dominant morphodynamic process in the highest portions of the river basins. It is no coincidence that the inflection points of the hypsometric curve of the basins, which marks the transition from a concave to a convex shape of the curve itself, are also placed at the same elevation. These points represent the threshold between low-impact erosive processes (sheet and rill erosion) and deep incision phenomena (gullies) and, probably, the point at which the maximum decrease in potential energy is linked to the maximum increase in kinetic energy. This aspect, as mentioned previously, had already been highlighted by Strahler [23], although he merely evidenced the transition point between an area with a higher change in mass rate and one with a lower one.

Last but not least, the indication that emerges from calculating the mean contribution areas of the reaches mentioned above (Amean_III contr); all these areas, which belong to basins of very different sizes and morphometric characteristics, are located in a relatively limited range, with values between 10,000 and 20,000 m² (between 1 and 2 hectares, Table 1). This fact, considering the different rainfall regimes that characterize the two test areas, represents a crucial point in understanding the mechanisms that regulate the erosion processes at the slope scale, since it allows us to consider the “position” factor as the only actual “driver” of these processes. The assessment of topographic thresholds for gully head development is certainly not a new aspect, and many scientific articles have dealt with this issue. However, it is interesting to note that although the topographic threshold conditions in such articles are usually given as double logarithmic plots of slope area (A) and slope gradient (s), the method used in this paper yielded similar values [58,59,60]. On the other hand, it differs greatly from other methods that use various topographic indices to describe the hydro-geomorphological characteristics of river catchments; in the latter, in particular, they emphasize specific aspects, such as the sediment transport capacity (RUSLE LS factor [61]), the potential channel erosion (Stream Power Index—SPI [62]), the transport capacity (Terrain Characterization Index—TCI [63]) or even the topographic control over hydrological processes (Topographic Wetness Index—TWI [62,64]). The above indices identify, among other things, areas prone to triggering different types of erosion and deposition processes that do not always correspond to the triggering height of the gullies’ heads.

The advantages of this approach are, as already mentioned, indisputable. The possibility to apply an exclusively morphometric method to DTMs with a medium depth of detail (e.g., 10 m), which can be easily acquired over large areas, makes the method itself extremely versatile and useful in an initial phase of investigation or territorial planning; compared to other approaches, no high-resolution images, field measurements or UAV surveys are required [36,57]. Finally, it allows comparisons between environments and geographical areas, even if they are very different and distant from each other.

An obvious limitation is that it deals with essentially homogeneous basins (e.g., with extensive debris cover or with outcropping clayey bedrock) that have little or no tectonic conditioning and a low degree of anthropization. Furthermore, as already mentioned, it is not readily applicable to DTMs with a low level of detail (20–30 m) or those with a high level of detail (e.g., lidar); in the latter case, automatic reconstruction in a GIS environment often results in an unrealistic and overly hierarchical hydrographic network.

5. Conclusions

The present work aimed to demonstrate how it is possible, through an exclusively morphometric analysis and the use of a digital terrain model with adequate resolution, to hypothesize the elevation and the average contribution area to activate accelerated erosion processes classified as gullies in relatively homogeneous basins (from a lithological point of view) and with little or no geological–structural conditioning. This approach is much more versatile and easier to handle than other methods that require a higher number of parameters and/or specific field investigations, and also with regard to the use of increasingly powerful IT tools (GIS) and detailed digital terrain models.

Although verified on a limited but still significant number of drainage basins, the results from the present study demonstrate that gully erosion processes require a minimum threshold of “energy” to be activated and that this threshold can be associated with a specific area or altitude. Consequently, morphometric analysis can be considered a valuable and versatile tool for rapidly assessing the location of erosive processes and joined to other techniques, for their quantification and the design of soil erosion remedial works.

More specifically, it highlighted that:

• Using an open-source digital terrain model with 10 m resolution, easily downloadable from public sites, it is possible to obtain the hydrographic network hierarchized according to the Strahler method of each river basin through simple geoprocessing operations in a GIS environment.
• The third-order reaches, classified with the method described above, correspond to segments of the hydrographic network characterized by intense linear erosion processes (gullies); this correspondence is easily verifiable through aerial photo interpretation.
• The mean starting elevation of these features coincides with the mean elevation of the basin (calculated from the hypsometric curve with GIS procedures) and with the highest elevation among those where the reaches described above show a peak in the frequency distribution: a similar coincidence exists between these relative elevations and the more elevated inflection point of the hypsometric curve of the basin.
• The elevation of the inflection point obtained from the hypsographic curve, which coincides with the mean elevation of the drainage basin, is confirmed to be a transition point between low-impact erosive processes (sheet and rill erosion) and deep incision phenomena (gullies).
• The portion of the basin necessary to activate these intense erosive processes is always included in a range of 1–2 hectares or, in any case, in a relatively narrow range regardless of the size and the morphometric characteristics of a specific basin.
• The almost completely automated method described in this study represents both an advantage and a limitation for its application since, as mentioned, it can only be used in homogeneous and practically unconditioned drainage basins.

Software: All data processing necessary to produce the maps was performed using ESRI-ArcGIS Pro software (ESRI, v.3.3.1) and Microsoft Excel v.365. The Figures in the document were created using CorelDraw Essential 2021.