The Use of Stereoscopic Satellite Images to Map Rills and Ephemeral Gullies

Abstract

Accurate mapping and measurement of erosion channels is necessary to accurately estimate the impact of channeled erosion in an area. Field surveys can provide optimal quantitative results, but they are only applicable to small areas. Recently, photogrammetric techniques have been applied to small format aerial photographs that were taken by UAVs. Few studies have applied photogrammetry for mapping and measuring single permanent gullies using very high resolution stereoscopic satellite images. We explore the use of such images to map rills and ephemeral gullies and to measure the length, width and depth of individual erosion channels to estimate the eroded volumes. The proposed methodology was applied to the Collazzone area of Central Italy. All of the channel characteristics were determined using GeoEye-1^® panchromatic stereoscopic satellite images of the 48-km² study area and a 3D floating cursor. We identified, mapped, and measured the lengths of 555 channel segments. The top width and depth could be measured in only a subset of the channel segments (the SMC subset). The SMC data were used to determine the coefficients of the power law relationship between the rill/gully volume and length (V = aL^b) and the uncertainties due to the channel depth measurements and the cross-sectional shape. The field data of the rill and gully volumes were within the estimated uncertainty. We defined a decision rule to distinguish rills from gullies on the basis of the segment length and applied the corresponding power law relationship that was derived from the SMC subset to estimate the eroded volume of the entire dataset. The erosion values that were calculated at different scales (0.680 Mg∙ha⁻¹ at the catchment scale, 28.4 Mg∙ha⁻¹ on the parcels affected by erosion) are consistent with values found in the literature. Our results indicate that erosion at the catchment scale can be considered moderate, whereas the erosion at the field scale exceeds the tolerance limit, which is consistent with data that have been summarized and/or discussed by several authors.

1. Introduction

During intense rainfall events, soil erosion often excavates two types of channels: rills and ephemeral gullies. Rills are usually small, with typical widths of approximately 10–20 cm and depths of 10–15 cm. Ephemeral gullies are larger channels that are approximately 50–100 cm wide and 40–100 cm deep. Rills and ephemeral gullies are typical erosion products that are associated with cropland and are obliterated by normal tillage practices as well as by land levelling when the excavated channels become too large. Permanent gullies are characterized by larger channels that can be as much as 500 m wide and up to 300 m deep. Detailed discussions on the definition of and the differences between rills and gullies can be found in Poesen et al. [1] and Torri and Borselli [2].

The importance of rills and ephemeral gullies is not limited to erosion but is also related to sediment and water connectivity in the landscape. Rills and ephemeral gullies increase the connectivity of the slope to the base slope and to the drainage network, which affects the chances of local floods and contributes to damage to infrastructure and private property [3,4]. For these reasons, assessing their presence and monitoring their changes may contribute to a better understanding of floods and post-disaster management. Hence, increasing the number of techniques for detecting and evaluating rill and ephemeral gully erosion is beneficial to disaster management and prevention. Because permanent gullies are permanent features of the landscape, they will not be discussed in this paper, and the term gully will refer only to ephemeral gullies.

The mapping of rills and ephemeral gullies and the measuring of soil loss when rills and ephemeral gullies are generated or further excavated are not simple procedures. Ephemeral gullies are usually eradicated soon after the triggering rainfall [5] and can even go undetected when gauges in the instrumented basin become clogged by excessive sediment [6]. As a result, most intense events are measured based on the dimensions of rills and gullies that are carved into the slopes [7]. Collecting data on rill and ephemeral gully erosion is made even more complex by the temporary nature of the phenomenon, which is easily removed by tillage equipment or bulldozers or masked by the growth of vegetation, making prompt and fast mapping surveys of rills and gullies extremely important [5]. Different techniques have been used to compile rill and gully maps. Traditionally, reconnaissance field surveys [8,9,10,11,12] or interpretations of topographic maps and aerial photographs have been used to obtain information on the number, density, size, type and location of head basins [8,13]. Over the last 20 years, various remote sensing data and techniques have been used for rill and gully detection and mapping, including stereoscopic aerial photographs that were taken shortly after an event, sets of multi-temporal aerial photographs [14,15,16] and the visual or digital analysis of high-resolution DEMs that were obtained from kite, airborne Lidar sensors and Unmanned Aerial Vehicles (UAVs) [17,18,19,20] or from terrestrial photographs [21,22,23].

Examples of rill and gully erosion measurements with photogrammetric techniques abound in the literature. Nachtergaele and Poesen [5] assessed the use of stereoscopic aerial photographs to determine the erosion rates of ephemeral gullies in Belgium. The authors suggested a correction factor that correlates the measured length on aerial photographs to the actual length of the channels. Watson and Evans [24] compared the measurements of eroded volumes from a field survey with the estimated volumes that were obtained through perspective photographs for an area in northeastern Scotland. Despite the limitations, these techniques offer a quick and easy method to estimate soil loss in areas that are dominated by channeled erosion processes. Keech [25] used multi-temporal stereoscopic aerial photographs at a scale of 1:25,000 that were taken in 1946, 1956, 1968, 1971, 1976, and 1984 to map gullies in the Mhondoro Communal Area of Zimbabwe. The depth of the gullies and eroded volumes were estimated using photo-analytical plotters.

Satellite images can be used to map erosion features. In addition to filters and classification techniques [26,27,28], low- and medium-resolution images were used in the past to map large- and medium-sized gullies [29,30] or to estimate the total area that is affected by gully erosion.

Over the past decade, the technological development of satellite sensors, the availability of images with ever higher spatial resolution and the improvement of digital visualization and analysis techniques have encouraged investigators to utilize satellite images to detect and map slope processes. Previously, automatic and semiautomatic techniques were used to map channeled erosion features. Bouaziz et al. [31] used a pixel-based maximum likelihood classification, whereas Vrieling et al. [32] used MODIS and Aster NDVI to evaluate erosion risk. New high-resolution satellite images and improved digital imaging technology also open up new possibilities for object-based image analysis [33] in the mapping of rills and gullies. Shrutri et al. [34] developed a knowledge-based generic method for gully assessment by coupling object oriented analysis and spatially dynamic erosion modeling in the sub-humid to semi-arid region of Sehoul, Morocco.

An inexpensive and rapid method for mapping permanent gullies based on the exploitation of very high-resolution (VHR) satellite images (GeoEye-1^® and Cnes SPOT^® images) that are accessed via Google Earth^® has been tested in Ethiopia [35,36]. The availability of VHR satellite stereoscopic images (Cartosat-1, 2.5 m resolution) has recently been used as an alternative to stereoscopic aerial photographs to quantify the erosion of permanent gullies in northern Australia [37].

In the abovementioned studies, the gullies that were examined using satellite images were permanent medium- to large-size features of the landscape. Ephemeral gullies and rills were ignored because of both their small dimensions and low permanence in the landscape. Most research focused on mapping erosion features using monoscopic satellite images rather than measuring the morphological characteristics (i.e., width, depth, length) of erosion channels. Of the previous studies, only Sattar et al. [37] used 2.5 m stereoscopic satellite images to measure the depth of large permanent gully channels.

In this work, we present the results of an experiment that tests the efficiency of using digital stereoscopic satellite images to rapidly quantify ephemeral gully erosion over a large area. The main objectives of this paper are to (1) evaluate the capability of mapping and measuring ephemeral gullies and rills using stereoscopic satellite images and (2) assess the reliability and associated error of satellite-based rill and gully measurements to estimate eroded volumes over large areas.

2. Study Area

The study area is located in Umbria, Central Italy (Figure 1A,B), at the coordinates 42°49′N and 12°27′E, and covers an area of approximately 48 km², with an elevation that ranges from 146 m to 634 m. The area is characterized by outcrops of sedimentary rocks (Figure 2A), including (i) alluvial deposits (5%), chiefly along the main valley bottoms; (ii) Plio-Pleistocene continental deposits (52%), sand (10%) and clay (2%); (iii) Pleistocene travertine deposits (6%); (iv) Miocene layered sandstone and marl (12%); and (v) Lias to Oligocene thinly layered limestone (13%) [38,39,40,41].

Figure 1. Location map. (A) The Umbria Region is in pink (inset), and the shaded relief map shows the location of the study area (DEM 20 meter resolution obtained from interpolation of counter lines). The blue dot represents the Todi rain gauge. (B) GeoEye-1^® image taken on 27 May 2010.

Figure 1. Location map. (A) The Umbria Region is in pink (inset), and the shaded relief map shows the location of the study area (DEM 20 meter resolution obtained from interpolation of counter lines). The blue dot represents the Todi rain gauge. (B) GeoEye-1^® image taken on 27 May 2010.

The climate of the study area is Mediterranean, with most of the precipitation falling between October and December. The historical rainfall time series (1951–2013) for the Todi rain gauge (42°46′36″N and 12°24′18″E; elevation: 283 m asl), (Figure 1A) indicates that the average annual rainfall is 841.1 mm and that the maximum and minimum average monthly precipitation levels occur in November (116 mm) and July (38 mm), respectively. The maximum observed monthly rainfall was 456 mm, recorded in October, and the maximum rainfall amount in 24 h was 124 mm, recorded in September.

The soil moisture regime in the study area is a typical xeric moisture regime of areas with cold and damp winters and dry summers. The soil texture is medium-fine, and the soil thickness ranges from 0.25 m to 1.5 m. Five profiles are available in the study area, and these profiles exhibit a well-developed A horizon (45–50 cm) with deterioration due to ploughing and a B horizon up to 90 cm deep (Soil map of the Regione Umbria, 1:250,000 (http://www.agricoltura.regione.umbria.it)).

The land cover (Figure 2B) is characterized by croplands (60.22%), forests (23.08%), urban areas (8.99%), pastures (6.03%), vineyards and orchards (1.46%) and water (0.23%). Farming in the area favors the development of slope failures (landslides, channeled and sheet erosion). The slope dynamics are characterized by landslides, mainly soil slides and earth flows [42,43], which are triggered by prolonged rainfall and rapid snow melting [44]. No studies have focused on gully erosion, despite the fact that ephemeral gullies in croplands are a common feature in this area.

Figure 2. (A) Simplified lithological map and pie chart for the study area, modified after Servizio Geologico d’Italia [39]. The figures represent the percentages of the lithologies that are present in the study area. (B) Land use map and pie chart for the study area. The figures represent the percentages of the soil use in the study area.

Figure 2. (A) Simplified lithological map and pie chart for the study area, modified after Servizio Geologico d’Italia [39]. The figures represent the percentages of the lithologies that are present in the study area. (B) Land use map and pie chart for the study area. The figures represent the percentages of the soil use in the study area.

3. Materials

Two pairs of stereoscopic panchromatic satellite images are available for the test area: one taken by the WorldView-1^® satellite on 8 March 2010 and the other from the GeoEye-1^® satellite on 27 May 2010 (Table 1). The GeoEye-1^® and WorldView-1^® images are in the panchromatic band with a 0.5 × 0.5 m ground sample distance. All the images feature Rational Polynomial Coefficients (RPCs), which indicate the image’s position relative to the ground.

Table 1. Characteristics of the stereoscopic satellite images that were acquired by the GeoEye-1^® and WorldView-1^® satellites.

	Date	Overlap (%)	Azimuth Angle (°)	Elevation Angle (°)
WorldView-1®	8 March 2010	70	77.30	57.90
WorldView-1®	8 March 2010	70	141.1	52.70
GeoEye-1®	27 May 2010	95	8.41	72.26
GeoEye-1®	27 May 2010	95	199.51	71.99

Table 1. Characteristics of the stereoscopic satellite images that were acquired by the GeoEye-1^® and WorldView-1^® satellites.

	Date	Overlap (%)	Azimuth Angle (°)	Elevation Angle (°)
WorldView-1®	8 March 2010	70	77.30	57.90
WorldView-1®	8 March 2010	70	141.1	52.70
GeoEye-1®	27 May 2010	95	8.41	72.26
GeoEye-1®	27 May 2010	95	199.51	71.99

The ERDAS IMAGINE and Leica Photogrammetry Suite (LPS) software were used to prepare the oriented stereo-pair. The WorldView-1^® stereo-pair was oriented using the model sensor for the internal orientation and the RPC files for the external orientation. The stereo-pair GeoEye-1^® was oriented using the sensor model and the associated RPC. The oriented model was improved in accuracy by adding 26 ground control points (GCPs). The estimated error of the oriented model given by the root mean square error (RMSE) is 0.11 pixels. In detail, the GCPs’ RMSEs are as follows: X is 1.75 m, Y is 0.99 m, and Z is 2.01 m. Files of this block type are compatible with the Stereo Analyst ArcGIS^® extension.

4. Methods

4.1. Rainfall Data and Event Reconstruction

To investigate the triggering of rills and gullies, we analyzed the rainfall that was recorded by the Todi rain gauge station (Figure 1). Two periods were analyzed: (i) one preceding the date of the images acquisition (27 May 2004), and (ii) one preceding the date of the field survey activities (30 April 2013). This is needed to understand whether the rainfall events in the two periods are similar.

Because abundant rills and gullies are visible in the 27 May 2004, GeoEye-1^® stereo image but not in the 8 March 2004, WordView-1^® stereo image, we assume that they developed between 8 March and 27 May.

An analysis of the 30-min cumulative rainfall data series allowed us to identify multiple rainfall events within the selected time period. This analysis was performed using the automated classification procedure for rainfall events that was proposed by Rossi et al. [45,46,47], which classifies a rainfall event as a period of consecutive hours with a cumulative rainfall of more than 1 mm. The total cumulative rainfall (C), the duration (D) and the average intensity (I) were calculated for each event in the selected period.

4.2. Use of Planar’s SD StereoMirror^TM Technology

A 3D view of the GeoEye-1^® pan-sharpened images was generated in the Stereo Analyst module [48], and the Planar StereoMirror^™ [49] technology was used to visualise the 3D view.

Planar’s SD StereoMirror™ technology works with two active matrix liquid crystal display (AMLCD) monitors that are oriented at a 110° angular distance. A passive beam splitter mirror bisects the angle between the two monitors. One side of the glass mirror has a reflective coating, and the other side has an anti-reflective coating, thereby allowing the user to see the two stereo-images of the monitors at the same time [50]. When stereo-pair images from the two monitors are viewed through crossed polarizing glasses, the user only sees the left-eye image with one eyepiece and the right-eye image with the other eyepiece due to different polarization angles (Figure 3). The result is a single, fused stereoscopic image. The Stereo Analyst extension works with a floating cursor; in this case, the cursor resides on the topographic surface and is used to draw vector features. A 3D floating cursor consists of an independent cursor that is displayed on the left image and an independent cursor that is displayed on the right image of the stereo-pair. When images are not viewed in stereo, the 3D floating cursor appears to be two separate cursors that may or may not rest on the same feature. However, when viewed in stereo, the two cursors fuse to create the perception of a 3D floating cursor. In this type of workstation, the scroll wheel moves up or down to either increase or decrease the elevation (Z) of the 3D floating cursor. With this digitizing device, the technician has full control over the 3D floating cursor and can position it accurately and reliably within the stereo-pair. Every time the 3D floating cursor is adjusted, new 3D coordinates (E and N coordinates and elevation) are computed using the sensor model information that is associated with each oriented image in the pair. Importantly, an elevation model is not required to collect reliable 3D GIS data when oriented images are used. However, the manual adjustment of the position of the 3D floating cursor requires the continuous attention of the operator.

Figure 3. Schematic of the StereoMirror™ technology (http://www.planar.com/media/211327/mn-planar-sd2220w.pdf).

Figure 3. Schematic of the StereoMirror™ technology (http://www.planar.com/media/211327/mn-planar-sd2220w.pdf).

Using the 3D floating cursor it is possible to measure relative heights of objects from the difference of elevation. Stereoscopic satellite images provide low-accuracy estimates of the absolute elevation (z) (1–2 m) but high-precision measurements of the relative heights of objects [51].

4.3. Mapping and Measuring Rills and Gullies on Satellite Stereo Images

Rills and gullies can be recognized in stereoscopic panchromatic satellite images through the heuristic interpretation of photographic characteristics, such as shape, color (grey shade), tone and texture. In particular, linear erosion features, such as rills and gullies, are characterized by long, narrow shapes and light tones. The light tones that highlight recently eroded channels are caused by increased reflectance, which results from several differences between the channel surface and the inter-channel degraded seed-bed surface (e.g., surface roughness, organic matter content, soil moisture, and vegetation cover).

The visibility of the ephemeral gullies and rills within the images depends on morphometric factors, such as the width of the channel, the illumination conditions and the characteristics of the image (contrast and resolution). The contrast is the most important characteristic for image interpretation; features can only be distinguished if they have noticeably different levels of brightness. The advantage of using satellite images (and digital images in general) is the ability to dynamically vary the contrast (Figure 4).

Figure 4. (A) Example of a ploughed field on an unstretched image and (B) on the enhanced image that was obtained by stretching the histogram values of the scene.

Figure 4. (A) Example of a ploughed field on an unstretched image and (B) on the enhanced image that was obtained by stretching the histogram values of the scene.

All of the visible channels were mapped and edited in a 3D GIS environment through the photointerpretation of stereoscopic satellite images. The central lines of the channels were digitized in the flow direction. Each channel segment is delimited by two consecutive nodes, where a node is either an intersection point between channels or a channel head/end (Figure 5A). The length of every channel segment was obtained as an attribute of the GIS database. The top widths and depths of the channels could be measured only in a subset of the channels (SMC) that was chosen based on their visibility on the images. The bottom widths of the channels could not be measured without ambiguity.

Figure 5. (A) Sketch of the channel mapping on satellite images. The yellow dots represent the nodes. Each channel is delimited by two consecutive nodes. (B) Sketch of the channel depth measurements. The crosses represent the positions of the floating cursor.

Figure 5. (A) Sketch of the channel mapping on satellite images. The yellow dots represent the nodes. Each channel is delimited by two consecutive nodes. (B) Sketch of the channel depth measurements. The crosses represent the positions of the floating cursor.

All of the attributes of the SMC channel segments were measured using the 3D floating cursor [50] directly on the stereo-pair in a 3D environment. For each segment, we measured a series of top widths (TW) and depths (h). The top width of the channel was measured as the length of the segment that is transverse (orthogonal) to the channel and could be identified from the differences in the shades of grey inside and outside the channel. The channel depth is the difference between the elevations of the side (top) of the channel segment and the channel bottom (Figure 5B).

4.3.1. Rill and Gully Volume Estimation from Satellite Measurements

The volume of each of the SMC channel segments was calculated. From the series of TW and h values that were measured along the channel segment, we calculated the average values of $\overline{TW}$ and $\overline{h}$. Hence, each channel segment was characterized by its length (L), $\overline{TW}$ and $\overline{h}$. We assumed a trapezoidal cross-sectional shape of the channel with a bottom width BW that is estimated to be a fraction of $\overline{TW}$. This proportion was estimated based on the average BW/TW ratio that was measured in the field (Section 5.3).

The CS of each channel segment was then estimated using the following equation:

$$\text{CS}=\left[\frac{\left(1+k\right)\times \overline{TW}}{2}\right]\times \overline{h}$$

(1)

where k =BW/TW, to obtain the volume,

$$\text{V}=\text{CS}\times \text{L}$$

(2)

4.4. Rill and Gully Field Measurements

In the spring of 2013, a field survey was conducted in the study area to measure the characteristics of the rill and gully channels and to calculate their volumes. The gullies and rills were measured using the procedure that was described in Casalí et al. [10] for the “detailed measurement of section characteristic lengths with a tape.” The sections were chosen and measured every meter or every five meters along the channel depending on the channel length (Figure 6A). In each section, the depth of the channel was measured with a variable pitch according to the width of the channel (from 5 cm to 50 cm for wider channels) (Figure 6B). The sections were then drawn in a CAD system to automatically calculate the area of each channel’s cross section (CS). Finally, the volume of each channel segment (V_c) was calculated using the equation

$${\text{V}}_{\text{c}}={{\displaystyle \sum }}^{\text{}}\left(\frac{{\text{CS}}_{\text{i}}+{\text{CS}}_{\text{i}+1}}{2}\right)\times {\text{L}}_{\text{i}}$$

(3)

where CS_i is the area of a section, CS_{i + 1} is the area of the next section, and l_i is the length of the channel between CS_i and CS_{i + 1} (Figure 6C).

Figure 6. Sketch of the field mapping procedure. (A) Measurement points in the field. (B) Cross section measurements. (C) Measurements of the channel volumes.

Figure 6. Sketch of the field mapping procedure. (A) Measurement points in the field. (B) Cross section measurements. (C) Measurements of the channel volumes.

4.5. Channel Classification

We adopted the cross section dimension threshold [52,53] as the criterion to distinguish between the two erosion phenomena (rills and ephemeral gullies), which states that the limit between rills and gully is a cross-sectional area of 0.0929 m² (1 square foot). This criterion was applicable to the field dataset and to the SMC subset.

4.6. Uncertainty of the Measurements

4.6.1. Accuracy of Satellite Height Measurements

We selected specific locations in the field as reference points (hereafter denoted as the height measurement survey) to assess the accuracy (i.e., the errors) of the height measurements that were derived from the satellite images using the floating cursor (see Section 4.2). The survey was conducted in October 2010 when we measured the heights of 24 stable objects (e.g., pavement, walls, and buildings) because they are constant over time, unlike the rills and gully channels. The field measurements were compared with the heights of the corresponding objects that were calculated from the satellite images (Section 5.4.1).

Figure 7. The pink dots represent the rainfall intensity in mm∙h⁻¹, the blue bars indicate the daily rainfall, and the violet line represents the cumulative rainfall. Only events that exceed 4 mm∙h⁻¹ are shown.

Figure 7. The pink dots represent the rainfall intensity in mm∙h⁻¹, the blue bars indicate the daily rainfall, and the violet line represents the cumulative rainfall. Only events that exceed 4 mm∙h⁻¹ are shown.

4.6.2. Volume-Length Relationship from Satellite Data and Associated Uncertainty

Several authors [54,55,56,57,58] have proposed that the length, L, of a gully or rill channel can serve as an estimator of the eroded volume according to the following power law equation:

$$V=a{L}^b$$

(4)

where $a$ is a coefficient and $b$ is a power law exponent. In a bi-logarithmic diagram, the power law equation is linear, with $b$ controlling the angular coefficient and expressing the rate of variation of V depending on L, and a is the scale factor.

During the analysis of the SMC subset that were measured from the satellite images, we first calculated the best-fitting values of a and b in the V-L relationship and then estimated the uncertainty that was associated with the fit using a bootstrap procedure [59,60] and by considering the error that was associated with the satellite depth measurements (see Section 4.6.1). The analysis was performed using R, a free software environment for statistical computing and graphics [61].

The best-fitting procedure for the calculation of the power law parameters was performed using a classical least-squares approach.

We executed the following steps to estimate the uncertainty that was associated with the fit:

• The uncertainty that was associated with the channel depth satellite measurements (Section 5.4.1) was modelled by assuming a normal distribution for the residuals in the linear regression model.
• We used the Kolmogorov-Smirnov test (KS) and the quantile-quantile (QQ) plot of the modelled and measured errors to determine whether the normal distribution was appropriate for modeling the uncertainty that was associated with the channel depth satellite measurements.
• The normal distribution was used to generate n synthetic series of channel heights that were associated with each channel segment. Then, knowing the width and the length of the channel segments that were measured from the satellite stereo-images, we derived a synthetic series of channel volumes.
• The n synthetic series of channel volumes that were associated with the satellite channel length measurements were fitted using a least-squares approach to estimate the n power law parameters (a and b).
• The 35th and the 65th percentile values of the n power law parameters (a and b) were assumed to be representative of the uncertainty that was associated with the satellite V-L relationship estimates.

The procedure was applied to the rills and gullies separately, considering n = 10,000 synthetic series to evaluate the V-L relationship parameters and the associated uncertainty.

4.6.3. TW/BW Ratio and Associated Uncertainty.

It is only possible to measure the TW dimension on the satellite images. The cross sections are evaluated based on an average value of the ratio BW/TW, which is calculated on the rills and gullies that were surveyed in the field; this process introduces a simplification and is a source of error. To assess its importance, we can assume that the shape of the cross section is approximated by different geometrical shapes. The CS formula can be varied continuously to change the section shape from rectangular to triangular through all possible trapezoids. CS can be calculated as follows:

$$\text{CS}=\left[\frac{\left(1+k\right)\times TW}{2}\right]\times h$$

(5)

where k = BW/TW.

This is a general form that results in a rectangle when BW/TW = 1 and a triangle when BW/TW = 0. Thus, we calculated the effect of changing the value of the ratio from 0 to 1 on the calculated volume (Section 5.4.2).

4.7. Eroded Volumes and Denudation

The V-L power law relationships that were obtained from the analysis of the SMC segments were used to estimate the eroded volumes for all 555 channels that were mapped on the satellite images. For the 500 channels for which the widths, depths and associated cross sections were not measured, we could not define whether the channels were rills or gullies and thus could not choose which equation to use. An alternative solution consists of evaluating the volumes for the extreme cases in which the channels were all rills or all gullies. The two components of the total soil erosion (rill and gully) can be added to obtain the total erosion if the weights are known. Because we want to evaluate the erosion, the weighing factor should be the channel volume. Unfortunately, we cannot estimate the ratio (component volume)/(total volume) from the SMC or field survey data because in both cases, the rill and gully samples are not representative of the population; the selection was aimed at determining the volume/length relationship. Because the segment length is the only parameter that is available for all of the segments, we must define a criterion that is based on the segment length to be applicable to the entire satellite dataset. Hence, we propose a mixed approach that is based on segment length: we use the rill volume equation if the segment is shorter than the minimum gully length and the gully volume equation if it is longer the maximum rill length. If the segment is between minimum gully length and maximum rill length, we arbitrarily use the average of the rill and gully volumes.

In both extreme cases (all rills and all gullies), we used the V-L power law relationships that were derived using the 35th, 50th and 65th percentile values of the a and b parameters to simulate the total eroded volumes and their uncertainties. The values that correspond to the mixed approach were calculated only for the 50th percentiles. We estimated the erosion based on an approximate density (γ) of the eroded material of γ = 1.6 Mg∙m⁻³ [62].

5. Results

5.1. Rainfall Events

The total recorded rainfall from 8 March to 27 May 2010, was 217 mm. Intense rainfall events occurred between 3 May and 8 May with 32 mm of cumulative rainfall in 138 h and between 11 May and 17 May with 66 mm of cumulative rainfall in 160 h (violet line in Figure 7). During the 3–8 May event, the most intense rainfall was recorded on 6 May and amounted to 20 mm (blue bars in Figure 7), which corresponds to an intensity of 10 mm∙h⁻¹ (pink dots in Figure 7). During the 11–17 May event, the most intense rainfall was recorded on 13 May and totaled 21 mm, whit an intensity of 14 mm∙h⁻¹ (pink dots in Figure 7).

The rainfall conditions during the period that preceded the field survey (1 January to 30 April 2013) were analyzed (Figure 8). The total recorded rainfall in this period was 338 mm. The most intense rainfall events occurred on 12 March, with 45 mm of rainfall in 24 hours (blue bars in Figure 8), whit an intensity of 14 mm∙h⁻¹. Another intense event was recorded on 28 April with a maximum intensity of 12 mm⁻¹.

Figure 8. Rainfall measurements from the Todi rain gauge from 1 January to 30 April 2013. The purple dots represent the rainfall intensity in mm∙h⁻¹, the blue bars represent the daily rainfall, and the green line represents the cumulative rainfall. Only events that exceed 4 mm∙h⁻¹ are shown.

Figure 8. Rainfall measurements from the Todi rain gauge from 1 January to 30 April 2013. The purple dots represent the rainfall intensity in mm∙h⁻¹, the blue bars represent the daily rainfall, and the green line represents the cumulative rainfall. Only events that exceed 4 mm∙h⁻¹ are shown.

5.2. Mapping and Measuring Rills and Ephemeral Gullies on Satellite Images

We mapped 555 channels using the method that was described in Section 4.3. The erosion channels were distributed over 19 cropland under seed-bed conditions (Figure 9 and Figure 10).

Figure 9. Map that shows the locations of the field measurements of the rills and gullies (blue dots), the cropland under seed-bed conditions that were affected by channeled erosion (red polygons), and the cropland under seed-bed conditions not affected by rill and gully erosion (green polygons).

Figure 9. Map that shows the locations of the field measurements of the rills and gullies (blue dots), the cropland under seed-bed conditions that were affected by channeled erosion (red polygons), and the cropland under seed-bed conditions not affected by rill and gully erosion (green polygons).

Figure 10. Examples of rills and ephemeral gullies that were mapped on the satellite images.

Figure 10. Examples of rills and ephemeral gullies that were mapped on the satellite images.

The average channel length is 26.5 m (σ = 26.2 m), and the total length of the channels is 1.4 × 10⁴ m. We measured the top widths and depths of a subset of 55 channels segments (see Table 2) of the SMC subset (14 rills and 41 gullies). Rills and gullies are classified based on the cross-sectional dimensions [53]. The rill lengths range from 7.33 to 43.4 m, the depths range from 0.08 to 0.23 m, top widths range from 0.50 to 1.3 m, and the cross-sectional areas range from 0.052 to 0.092 m². The gully lengths vary between 10.61 and 113.44 m, the depths range between 0.08 and 0.33 m, top widths range from 0.52 to 2.6 m, and the cross-sectional areas range from 0.11 to 0.36 m².

For each of the 55 channels, the volume was calculated as described in Section 4.3.1. The volumes of the rills range between 0.39 m³ and 3.61 m³, and the volumes of the gullies range between 1.84 m³ and 32.80 m³.

Table 2. Measurements of the 55 channels that were obtained from the stereo satellite images. L: channel segment length; CS: cross-sectional area; $\overline{TW}$: major base (top width); $\overline{h}$ channel depth, BW: minor base (bottom width; equal to 0.6 TW); R: rill; EG: ephemeral gully.

Channel #	L m	$\overline{TW}$ m	BW m	$\overline{h}$ m	CS m2	Volume m3	Type
1	7.34	0.82	0.49	0.08	0.05	0.39	R
2	9.40	0.82	0.49	0.08	0.05	0.49	R
3	15.56	0.8	0.48	0.10	0.06	1.00	R
4	28.03	0.5	0.30	0.16	0.06	1.79	R
5	9.70	1.05	0.63	0.08	0.07	0.65	R
6	14.17	1.05	0.63	0.08	0.07	0.95	R
7	18.57	1.05	0.63	0.08	0.07	1.25	R
8	25.17	0.5	0.30	0.17	0.07	1.71	R
9	29.18	0.5	0.30	0.17	0.07	1.98	R
10	13.76	1.05	0.63	0.09	0.08	1.04	R
11	11.23	1.20	0.72	0.08	0.08	0.86	R
12	11.79	1.30	0.78	0.08	0.08	0.98	R
13	43.43	1.30	0.78	0.08	0.08	3.61	R
14	32.15	0.50	0.30	0.23	0.09	2.96	R
15	36.42	1.46	0.87	0.08	0.09	3.40	EG
16	22.63	1.20	0.72	0.1	0.10	2.17	EG
17	38.02	1.50	0.90	0.08	0.10	3.65	EG
18	39.12	1.50	0.90	0.08	0.10	3.76	EG
19	20.81	0.64	0.38	0.19	0.10	2.02	EG
20	27.16	1.50	0.90	0.09	0.11	2.93	EG
21	29.70	1.73	1.04	0.08	0.11	3.29	EG
22	55.96	0.93	0.56	0.16	0.12	6.66	EG
23	44.99	0.52	0.31	0.3	0.12	5.61	EG
24	17.86	2.00	1.20	0.08	0.13	2.29	EG
25	16.05	0.8	0.48	0.22	0.14	2.26	EG
26	26.48	0.9	0.54	0.2	0.14	3.81	EG
27	12.18	1.05	0.63	0.18	0.15	1.84	EG
28	12.98	0.87	0.522	0.22	0.15	1.99	EG
29	37.00	1.2	0.72	0.16	0.15	5.68	EG
30	34.85	1.00	0.60	0.2	0.16	5.58	EG
31	18.34	2.6	1.56	0.08	0.17	3.05	EG
32	30.96	0.88	0.52	0.25	0.18	5.45	EG
33	11.00	0.90	0.54	0.25	0.18	1.98	EG
34	44.81	1.35	0.81	0.17	0.18	8.23	EG
35	49.58	0.80	0.48	0.29	0.19	9.20	EG
36	25.93	0.90	0.54	0.26	0.19	4.85	EG
37	31.34	1.24	0.74	0.2	0.20	6.22	EG
38	40.70	1.10	0.66	0.23	0.20	8.24	EG
39	28.35	1.50	0.90	0.17	0.20	5.78	EG
40	20.92	1.20	0.72	0.23	0.22	4.62	EG
41	62.63	1.40	0.84	0.2	0.22	14.03	EG
42	70.42	1.80	1.08	0.16	0.23	16.23	EG
43	10.61	1.00	0.6	0.29	0.23	2.46	EG
44	113.43	1.58	0.94	0.2	0.25	28.67	EG
45	80.74	1.54	0.92	0.21	0.26	20.89	EG
46	20.87	0.98	0.58	0.33	0.26	5.40	EG
47	16.90	1.7	1.02	0.2	0.27	4.60	EG
48	26.02	1.8	1.08	0.19	0.27	7.12	EG
49	18.25	1.56	0.93	0.22	0.27	5.01	EG
50	25.01	1.46	0.87	0.26	0.30	7.60	EG
51	105.80	1.25	0.75	0.31	0.31	32.80	EG
52	43.60	1.48	0.88	0.27	0.32	13.94	EG
53	54.77	1.50	0.90	0.27	0.32	17.74	EG
54	57.06	2.15	1.29	0.2	0.34	19.63	EG
55	53.16	1.56	0.93	0.29	0.36	19.24	EG

Table 2. Measurements of the 55 channels that were obtained from the stereo satellite images. L: channel segment length; CS: cross-sectional area; $\overline{TW}$: major base (top width); $\overline{h}$ channel depth, BW: minor base (bottom width; equal to 0.6 TW); R: rill; EG: ephemeral gully.

Channel #	L m	$\overline{TW}$ m	BW m	$\overline{h}$ m	CS m2	Volume m3	Type
1	7.34	0.82	0.49	0.08	0.05	0.39	R
2	9.40	0.82	0.49	0.08	0.05	0.49	R
3	15.56	0.8	0.48	0.10	0.06	1.00	R
4	28.03	0.5	0.30	0.16	0.06	1.79	R
5	9.70	1.05	0.63	0.08	0.07	0.65	R
6	14.17	1.05	0.63	0.08	0.07	0.95	R
7	18.57	1.05	0.63	0.08	0.07	1.25	R
8	25.17	0.5	0.30	0.17	0.07	1.71	R
9	29.18	0.5	0.30	0.17	0.07	1.98	R
10	13.76	1.05	0.63	0.09	0.08	1.04	R
11	11.23	1.20	0.72	0.08	0.08	0.86	R
12	11.79	1.30	0.78	0.08	0.08	0.98	R
13	43.43	1.30	0.78	0.08	0.08	3.61	R
14	32.15	0.50	0.30	0.23	0.09	2.96	R
15	36.42	1.46	0.87	0.08	0.09	3.40	EG
16	22.63	1.20	0.72	0.1	0.10	2.17	EG
17	38.02	1.50	0.90	0.08	0.10	3.65	EG
18	39.12	1.50	0.90	0.08	0.10	3.76	EG
19	20.81	0.64	0.38	0.19	0.10	2.02	EG
20	27.16	1.50	0.90	0.09	0.11	2.93	EG
21	29.70	1.73	1.04	0.08	0.11	3.29	EG
22	55.96	0.93	0.56	0.16	0.12	6.66	EG
23	44.99	0.52	0.31	0.3	0.12	5.61	EG
24	17.86	2.00	1.20	0.08	0.13	2.29	EG
25	16.05	0.8	0.48	0.22	0.14	2.26	EG
26	26.48	0.9	0.54	0.2	0.14	3.81	EG
27	12.18	1.05	0.63	0.18	0.15	1.84	EG
28	12.98	0.87	0.522	0.22	0.15	1.99	EG
29	37.00	1.2	0.72	0.16	0.15	5.68	EG
30	34.85	1.00	0.60	0.2	0.16	5.58	EG
31	18.34	2.6	1.56	0.08	0.17	3.05	EG
32	30.96	0.88	0.52	0.25	0.18	5.45	EG
33	11.00	0.90	0.54	0.25	0.18	1.98	EG
34	44.81	1.35	0.81	0.17	0.18	8.23	EG
35	49.58	0.80	0.48	0.29	0.19	9.20	EG
36	25.93	0.90	0.54	0.26	0.19	4.85	EG
37	31.34	1.24	0.74	0.2	0.20	6.22	EG
38	40.70	1.10	0.66	0.23	0.20	8.24	EG
39	28.35	1.50	0.90	0.17	0.20	5.78	EG
40	20.92	1.20	0.72	0.23	0.22	4.62	EG
41	62.63	1.40	0.84	0.2	0.22	14.03	EG
42	70.42	1.80	1.08	0.16	0.23	16.23	EG
43	10.61	1.00	0.6	0.29	0.23	2.46	EG
44	113.43	1.58	0.94	0.2	0.25	28.67	EG
45	80.74	1.54	0.92	0.21	0.26	20.89	EG
46	20.87	0.98	0.58	0.33	0.26	5.40	EG
47	16.90	1.7	1.02	0.2	0.27	4.60	EG
48	26.02	1.8	1.08	0.19	0.27	7.12	EG
49	18.25	1.56	0.93	0.22	0.27	5.01	EG
50	25.01	1.46	0.87	0.26	0.30	7.60	EG
51	105.80	1.25	0.75	0.31	0.31	32.80	EG
52	43.60	1.48	0.88	0.27	0.32	13.94	EG
53	54.77	1.50	0.90	0.27	0.32	17.74	EG
54	57.06	2.15	1.29	0.2	0.34	19.63	EG
55	53.16	1.56	0.93	0.29	0.36	19.24	EG

5.3. Rill and Gully Field Measurements

The field measurements and relative volume estimates of the rills and gullies are reported in Table 3.

Table 3. Measurements of the rills and ephemeral gullies that were surveyed in the field. L: channel segment length; CS: cross-sectional area; BW: minor base (bottom width); TW: major base (top width); EG: gully; R: rill.

Id	Type	L m	# CS	CS Average m2	CS Max m2	CS Min m2	TW Max m	TW Min m	Volume m3	Average BW/TW ratio
1	R	47	20	0.069	0.11	0.02	0.90	0.34	3.32	0.7
2	R	8.8	7	0.050	0.08	0.03	1.20	0.40	0.55	0.7
3	R	19	7	0.069	0.04	0.01	0.50	0.34	0.44	0.5
4	R	30	8	0.061	0.10	0.03	1.00	0.40	1.50	0.7
5	R	3.2	2	0.085	0.09	0.08	0.50	0.50	0.30	0.5
6	R	32.5	7	0.062	0.14	0.08	1.30	0.51	1.88	0.7
7	EG	34	7	0.103	0.25	0.07	1.70	1.00	4.14	0.6
8	R	30	7	0.082	0.10	0.04	1.20	0.60	2.13	0.6
9	EG	22.4	7	0.144	0.20	0.09	1.20	0.48	2.68	0.7
10	EG	13	4	0.205	0.27	0.12	1.08	0.50	2.23	0.5
11	R	22.5	8	0.091	0.11	0.20	1.17	0.30	1.15	0.5
12	EG	18.8	6	0.250	0.36	0.16	1.19	0.77	3.80	0.3
13	EG	10	5	0.198	0.24	0.14	1.00	0.70	1.85	0.3
				Average BW/TW ratio						0.6

Table 3. Measurements of the rills and ephemeral gullies that were surveyed in the field. L: channel segment length; CS: cross-sectional area; BW: minor base (bottom width); TW: major base (top width); EG: gully; R: rill.

Id	Type	L m	# CS	CS Average m2	CS Max m2	CS Min m2	TW Max m	TW Min m	Volume m3	Average BW/TW ratio
1	R	47	20	0.069	0.11	0.02	0.90	0.34	3.32	0.7
2	R	8.8	7	0.050	0.08	0.03	1.20	0.40	0.55	0.7
3	R	19	7	0.069	0.04	0.01	0.50	0.34	0.44	0.5
4	R	30	8	0.061	0.10	0.03	1.00	0.40	1.50	0.7
5	R	3.2	2	0.085	0.09	0.08	0.50	0.50	0.30	0.5
6	R	32.5	7	0.062	0.14	0.08	1.30	0.51	1.88	0.7
7	EG	34	7	0.103	0.25	0.07	1.70	1.00	4.14	0.6
8	R	30	7	0.082	0.10	0.04	1.20	0.60	2.13	0.6
9	EG	22.4	7	0.144	0.20	0.09	1.20	0.48	2.68	0.7
10	EG	13	4	0.205	0.27	0.12	1.08	0.50	2.23	0.5
11	R	22.5	8	0.091	0.11	0.20	1.17	0.30	1.15	0.5
12	EG	18.8	6	0.250	0.36	0.16	1.19	0.77	3.80	0.3
13	EG	10	5	0.198	0.24	0.14	1.00	0.70	1.85	0.3
				Average BW/TW ratio						0.6

The 13 rills and ephemeral gullies that were mapped in the field range between 3.2 and 47 m in length. A total of 95 cross sections were measured; the areas ranged from 0.01 to 0.36 m², and the average area was 0.11 m². To evaluate the mean value of the BW/TW (bottom width/top width) ratio, each of the cross sections that were measured in the field was approximated by a trapezoid.

5.4. Measurement Uncertainty

5.4.1. Accuracy of Satellite Height Measurements

Figure 11C shows the correlation between the heights that were measured in the field and those that were measured from the satellite stereo images using the 3D floating cursor. The inset in Figure 11C shows the uncertainty that is associated with the channel depth satellite measurements that are modelled by assuming a normal distribution (mean value of 0.002 and standard deviation of 0.059).

5.4.2. Volume Estimation and Associated Error

The volumes and lengths of the SMC subset were used to estimate the V-L relationships for the rills and the ephemeral gullies. Figure 12 shows the results of the regressions. Both regressions are characterized by high coefficients of determination (R² = 0.5 for rills and R² = 0.47 for gullies). Additionally, the probability of rejection (p value), which was calculated using the one-tail Pearson’s correlation coefficient test, is less than 0.0005.

Figure 11. (A) Example of an object that was selected as a reference during the field survey (the points that were used to calculate and measure the height difference are shown in green). (B) Measurement of the height difference of a stable object in the field. (C) Correlation between the heights that were measured in the field and the corresponding heights that were derived from the satellite stereo images. The inset shows the error, which is modelled by a normal distribution.

Figure 11. (A) Example of an object that was selected as a reference during the field survey (the points that were used to calculate and measure the height difference are shown in green). (B) Measurement of the height difference of a stable object in the field. (C) Correlation between the heights that were measured in the field and the corresponding heights that were derived from the satellite stereo images. The inset shows the error, which is modelled by a normal distribution.

Figure 12. Power law relationship ($V=a{L}^b$) in a bi-logarithmic diagram for the rills and ephemeral gullies and the associated uncertainties. The small blue and red dots are the ephemeral gully volumes and rill volumes, respectively, which were measured using the stereo satellite images, and the large blue and red dots are the volumes of the rills and gullies, respectively, which were measured in the field. The uncertainties that are associated with the power law fits for the rills (red shaded area) and ephemeral gullies (blue shaded area) are defined based on the 35th percentile (m − 15) and the 65th percentile (m + 15), respectively, of the power law parameters a and b.

Figure 12. Power law relationship ($V=a{L}^b$) in a bi-logarithmic diagram for the rills and ephemeral gullies and the associated uncertainties. The small blue and red dots are the ephemeral gully volumes and rill volumes, respectively, which were measured using the stereo satellite images, and the large blue and red dots are the volumes of the rills and gullies, respectively, which were measured in the field. The uncertainties that are associated with the power law fits for the rills (red shaded area) and ephemeral gullies (blue shaded area) are defined based on the 35th percentile (m − 15) and the 65th percentile (m + 15), respectively, of the power law parameters a and b.

To consider the uncertainty of the channel depth estimates (which contributes to the uncertainty of the volume estimates), we compared the data that were collected during the height measurement survey (the height calibration curve in Section 5.4.1). Despite the high values of the R² coefficients, which ensure a high significance of the relationships with p < 0.0005, normally distributed errors are still associated with the measurements. By applying the bootstrap techniques that were described in Section 4.6.2, we generated error lines that correspond to the 35th and 65th percentiles of the errors that are associated with the estimated parameters a and b for the two volume-length relationships (Table 4). Figure 13 shows the quantile-quantile plot of the modelled and measured heights and the Kolmogorov-Smirnov test results. The shaded area in Figure 12 represents the uncertainties in the volumes that were estimated using the best-fitting values for the rills and ephemeral gullies. The value of b controls the length-dependent variation in the cross-sectional area. Larger values of b are associated with larger cross-sectional areas and longer channel lengths [36]. The b coefficients range from 1.04 to 1.30 (35th and 65th percentiles) for the rills and from 1.15 to 1.24 (35th and 65th percentiles) for the gullies (Table 4).

Table 4. Values of the coefficients a and b from the interpolation of the data that were collected using the satellite stereo-pairs. The values of the coefficients a and b include the 35th percentile (minimum value), the 65th percentile (maximum value) and the median value of the parameter distribution. The total eroded volume is estimated as the sum of the volumes of the 555 channels. The weight of the eroded volume is based on a medium with a density of γ = 1.6 Mg∙m⁻³ [62].

	Rill			Gully				Mixed Approach
	Min (35th perc.)	Max (65th perc.)	Average Median		Min (35th perc.)	Max (65th perc.)	Average Median	Average Median
a	0.027	0.061	0.044		0.072	0.101	0.082	--
b	1.04	1.30	1.20		1.15	1.24	1.17	--
V (m3)	4.62 × 102	2.791 × 103	1.321 × 103		1.844 × 103	3.655 × 103	2.317 × 103	2.041 × 103
Mg (γ = 1.6 Mg∙m−3)	7.39 × 102	4.465 × 103	2.113 × 103		2.950 × 103	5.848 × 103	3.707 × 103	3.260 × 103

Table 4. Values of the coefficients a and b from the interpolation of the data that were collected using the satellite stereo-pairs. The values of the coefficients a and b include the 35th percentile (minimum value), the 65th percentile (maximum value) and the median value of the parameter distribution. The total eroded volume is estimated as the sum of the volumes of the 555 channels. The weight of the eroded volume is based on a medium with a density of γ = 1.6 Mg∙m⁻³ [62].

	Rill			Gully				Mixed Approach
	Min (35th perc.)	Max (65th perc.)	Average Median		Min (35th perc.)	Max (65th perc.)	Average Median	Average Median
a	0.027	0.061	0.044		0.072	0.101	0.082	--
b	1.04	1.30	1.20		1.15	1.24	1.17	--
V (m3)	4.62 × 102	2.791 × 103	1.321 × 103		1.844 × 103	3.655 × 103	2.317 × 103	2.041 × 103
Mg (γ = 1.6 Mg∙m−3)	7.39 × 102	4.465 × 103	2.113 × 103		2.950 × 103	5.848 × 103	3.707 × 103	3.260 × 103

Figure 13. Kolmogorov-Smirnov test (KS) and the quantile-quantile (QQ) plot of the modelled and measured errors.

Figure 13. Kolmogorov-Smirnov test (KS) and the quantile-quantile (QQ) plot of the modelled and measured errors.

The BW/TW ratio influences the dimensions of the CS and the estimated volume. The error that is associated with the selection of the BW/TW ratio was evaluated. The best fit that was obtained for the 50th percentile for several values of BW/TW is shown in Figure 14. Because the field observations have BW/TW ratios between 0.3 and 0.7, the error due to the assumption of a constant BW/TW ratio is approximately half of the error due to the error that is associated with the depth (shaded area).

The volumes of the rill and gully channels that were measured in the field were compared with those that were estimated from the stereo images. In Figure 12 and Figure 14, the large dark blue and red dots represent the volumes of the rills and gullies in the field, whereas the smaller dots represent the volumes that were estimated from the satellite data. The field values refer to new channels, which were measured in a different year but under similar land-use and rainfall conditions (Section 5.1). The field values can serve as a general reference for evaluating the reliability of satellite measurements. In particular, the field-based volume measurements fall well inside the shaded area that represents the uncertainty that is associated with the volume-length relationship from the satellite stereo images.

5.5. Eroded Volumes and Denudation

The total eroded volumes for the study area were simulated (i) for the extreme cases in which the channels were all rills or all gullies and (ii) for the mixed approach (in which a fraction of the channel segments are considered as rills and the remaining as gullies, based on the length of the channel segment). Table 2 and Table 3 show that rill segments are never longer than 43.4 m (SMC data) or 47 m (field survey). The gully segments are never less than 10 m (field survey) or 10.6 m (SMC data) long. Hence, in the mixed approach we use the rill volume equation if the segment is shorter than 10 m and the gully volume equation if it is longer than 47 m. If the segment is between 10 m and 47 m long, we arbitrarily use the average of the rill and gully volumes.

Table 5. Erosion in Mg∙ha⁻¹ and denudation in mm for the entire catchment (total surface area: 48 km²), cropland/seed-bed conditions (5.92 km²) and the parcels that were affected by rill and gully erosion (1.149 km²) based on the values of the coefficients a and b for the rills and gullies for the 35th percentile (minimum value), 65th percentile (maximum value) and median value.

		Catchment A = 48 km2		Cropland Seed-Bed Condition A = 5.92 km2		Affected Parcels A = 1.149 km2
		Erosion (Mg∙ha−1)	Denudation (mm)	Erosion (Mg∙ha−1)	Denudation (mm)	Erosion (Mg∙ha−1)	Denudation (mm)
Rill	Min (35th perc.)	0.153	0.010	1.248	0.078	6.427	0.402
	Max (65th perc.)	0.930	0.0581	7.544	0.471	38.833	2.429
	Average (median)	0.440	0.027	3.571	0.223	18.382	1.150
Gully	Min (35th perc.)	0.614	0.038	4.985	0.311	25.661	1.605
	Max (65th perc.)	1.218	0.076	9.880	0.617	50.860	3.181
	Average (median)	0.772	0.048	6.263	0.391	32.242	2.017
Mixed Approach	Average (median)	0.680	0.043	5.517	0.344	28.401	1.776

Table 5. Erosion in Mg∙ha⁻¹ and denudation in mm for the entire catchment (total surface area: 48 km²), cropland/seed-bed conditions (5.92 km²) and the parcels that were affected by rill and gully erosion (1.149 km²) based on the values of the coefficients a and b for the rills and gullies for the 35th percentile (minimum value), 65th percentile (maximum value) and median value.

		Catchment A = 48 km2		Cropland Seed-Bed Condition A = 5.92 km2		Affected Parcels A = 1.149 km2
		Erosion (Mg∙ha−1)	Denudation (mm)	Erosion (Mg∙ha−1)	Denudation (mm)	Erosion (Mg∙ha−1)	Denudation (mm)
Rill	Min (35th perc.)	0.153	0.010	1.248	0.078	6.427	0.402
	Max (65th perc.)	0.930	0.0581	7.544	0.471	38.833	2.429
	Average (median)	0.440	0.027	3.571	0.223	18.382	1.150
Gully	Min (35th perc.)	0.614	0.038	4.985	0.311	25.661	1.605
	Max (65th perc.)	1.218	0.076	9.880	0.617	50.860	3.181
	Average (median)	0.772	0.048	6.263	0.391	32.242	2.017
Mixed Approach	Average (median)	0.680	0.043	5.517	0.344	28.401	1.776

The results are given in Table 5, which summarizes the eroded soil and denudation respectively in terms of tonnes per hectare (Mg∙ha⁻¹) and volume per unit area (mm). The total erosion and denudation refer to three areas: (i) the entire catchment (48 km²), (ii) only to the land that is affected by rills and ephemeral gullies (cropland under seed-bed conditions; 5.92 km²), and (iii) only the fields where channel erosion took place (1.149 km²)

6. Discussion

Previous attempts to use stereoscopic satellite images (e.g., [36,37]) have consisted of the recognition, measurement and mapping of medium to large gullies. In our experiment, VHR (i.e., 0.5 m spatial resolution) stereoscopic satellite images were used to detect, map and measure ephemeral gullies and rills that are easily removed by agricultural practices or hidden by vegetation growth; hence, this study identified and measured channel features at a more detailed scale.

The lengths and widths of channels have been measured in the previous studies because they are usually derived from monoscopic satellite images [63]. This has been overcome by low-altitude UAV photogrammetry because it uses stereo pairs that allow millimetric resolutions [64]. These types of extremely detailed measurements are generally not applicable to medium size catchment (i.e., >20–25 km²) because of the extensive field surveys and complex post-processing operations that are required to orient the large number of aerial photographs. The use of VHR satellite stereo pairs allows us to study large catchment with only one pair and generally requires only a few GCPs to improve the accuracy of the orientation.

Because rill and gully erosion networks are often anastomosing, it was necessary to subdivide channels into segments between successive nodes (head/end channel points and channel branches).

In this experiment, the rill and gully depths could only be measured along some of the channel segments (SMC subset). The width, length and depth measurements allowed the calculation of the channel volumes and the values of the coefficients a and b that were the result of the interpolation of the channel volume-length (V-L) data. Because the depth measurements depend on the ability of the photo-interpreters to use floating cursor techniques, the proposed approach evaluates the depth-measurement uncertainty and its impact on the estimated eroded volumes. The shaded areas in Figure 12 represent the uncertainties in the volumes that were estimated using the best-fitting values for the rills and ephemeral gullies. The uncertainty of the rill volume is greater than the uncertainty of the gully volume, which indicates that the error in the depth measurements has a larger effect on the volume estimates of rills.

Among the other sources of error, the effect of changing the shape of the channel cross section from triangular (BW/TW = 0) through trapezoidal (0 < BW/TW < 1) to rectangular (BW/TW = 1) was also evaluated. The cross-sectional shape (i.e., BW/TW) was found to have a smaller effect on the volume estimation than the depth error. Hence, the use of the field-data derived average TW/BW ratio to estimate cross section to obtain gully and rill volume from satellite images has a negligible effect on the volume estimation.

The values of the coefficients a and b for the ephemeral gully that were obtained in this study were in the range of values that have been obtained by several authors from field data (Table 6) [54,55,56,57,58], which suggests that the satellite stereoscopic images can provide reliable measurements to characterize the sizes of gullies and rills. When these trends were compared with the rill and gully channel data that were collected during a field survey in the same area, we obtained a reasonable agreement with the calculated trends. The comparison is possible because the rainfall conditions that triggered the erosion channels in the two periods were similar (maximum rainfall intensity in the two periods is 14 mm∙h⁻¹).Therefore, satellite images can be used to acquire information in areas that are affected by channeled erosion and to characterize transient erosion features, such as ephemeral gullies and rills in agricultural areas.

Figure 14. Analysis of the influence of BW/TW on the estimation of the volumes. A sketch of cross sections with different values of BW/TW is shown in the inset, different colors of the lines correspond to different value of the BW/TW. The thin lines represent the fit for the 50th percentile of the V-L relationship that was obtained for different BW/TW values, and the thick lines are the fits that were obtained for BW/TW = 0.6. The shaded areas represent the uncertainty of the volume estimates of the power law fits that is related to the error in the depth evaluation (Figure 12). The dots represent the volumes and lengths of the rills and ephemeral gullies that were measured in the field.

Figure 14. Analysis of the influence of BW/TW on the estimation of the volumes. A sketch of cross sections with different values of BW/TW is shown in the inset, different colors of the lines correspond to different value of the BW/TW. The thin lines represent the fit for the 50th percentile of the V-L relationship that was obtained for different BW/TW values, and the thick lines are the fits that were obtained for BW/TW = 0.6. The shaded areas represent the uncertainty of the volume estimates of the power law fits that is related to the error in the depth evaluation (Figure 12). The dots represent the volumes and lengths of the rills and ephemeral gullies that were measured in the field.

Table 6. Values of the coefficients a and b from the interpolation of the data that were collected from satellite stereo-pairs and the coefficients for ephemeral gullies that were obtained by several authors.

Authors	[54,56]	[55]	[56]	[57]	[58]	This Experiment
Study area	Spain	Belgium	Sicily(Italy)	Sicily(Italy)	Sicily(Italy)	Umbria(Italy)
a	0.054	0.048	0.082	0.034	0.812	0.082
b	1.0.37	1.29	1.416	1.42	1..062	1.174

Table 6. Values of the coefficients a and b from the interpolation of the data that were collected from satellite stereo-pairs and the coefficients for ephemeral gullies that were obtained by several authors.

Authors	[54,56]	[55]	[56]	[57]	[58]	This Experiment
Study area	Spain	Belgium	Sicily(Italy)	Sicily(Italy)	Sicily(Italy)	Umbria(Italy)
a	0.054	0.048	0.082	0.034	0.812	0.082
b	1.0.37	1.29	1.416	1.42	1..062	1.174

Distinguishing rills from gullies is a key step for correctly estimating the eroded volume. Table 4 shows that the rill volume that was estimated using the a and b parameters of the V-L relationship that corresponds to the 65th percentile is 6 times the volume that was estimated for the 35th percentile, while for the gullies, the volume that corresponds to the 65th percentile is 1.9 times the volume that was estimated for the 35th percentile. This result demonstrates the importance of the error in the depth estimation, especially for small depths (in the SMC the average rill depth is 11 cm, while the average gully depth is 20 cm).

The differences between considering the all-rill case versus the all-gully case (a and b parameters of the 50th percentile) result in a gully/rill ratio of 1.8; hence, the two cases add a further non-negligible uncertainty to the volume estimate. This difference decreases if we use the mixed approach, which recognized rills and gullies based on the segment length. The mixed approach suggests eroded values that are 1.5 times greater than the all-rill case and 0.9 times smaller than the all-gully case.

Assuming that the mixed approach provides the most realistic estimate, we can compare the estimate from the mixed approach to values of erosion from the literature. In particular, a value of approximately 0.5 Mg∙ha⁻¹ can be deduced from the data in Torri et al. ([65]) over the same period (i.e., 2.5 months), which is comparable to the erosion value at the catchment scale (0.680 Mg∙ha⁻¹). At the field scale (i.e., the affected parcels), the mixed approach suggests erosion of 28.4 Mg∙ha⁻¹, which falls within the values of gully erosion from the literature (i.e., 9–90 Mg∙ha⁻¹ over 2.5 months [7,66]). The fields that are affected by linear erosion are subject to intense erosion with a total denudation that is 1.77 mm, greater than the expected pedogenization for a substantially mature soil (e.g., Alexandrovskiy [67], who reported erosion of approximately 0.1–0.07 mm in 2.5 months for a 1–2 ka luvisol).

Our results indicate that erosion at the catchment scale is moderate [68,69], whereas the erosion at the field scale exceeds the tolerance limit, which is consistent with data that have been summarized and/or discussed by several authors [70,71].

Despite these advantages, the methodology that we described for the measurement of rills and gully channels has certain limitations. This analysis requires that operators be skilled at using the 3D floating cursor mapping approach. In particular, the operator must position the floating cursor correctly on the ground to obtain realistic measurements of an object’s height. Second, the analysis must be performed during fixed illumination conditions of the day of the satellite acquisition. As suggested by Giménez [72], certain illumination conditions can result in an inaccurate depth estimates of narrow erosional channels (i.e., when the TW/D ratio is less than 2.5). These problems were limited in our analysis because the satellite images were taken under adequate illumination conditions (27 May 2010, at 10:06 am GMT, 11:06 am local time in Italy; the local sun azimuth was 121.06° with an elevation of 56.63°), and the images were used to map shallow rills and gully channels with an average TW/D ratio of 8.9. However, these conditions limit the time window during which images should be acquired.

The use of very high-resolution stereoscopic satellite images allows the mapping of erosional phenomena over a short period of time compared to the time that is usually required for field survey mapping. Nachtergaele and Poesen [5] estimated that a 2 km² field survey requires two people and a full day to obtain basic data. Conversely, one person can map and measure the erosional features in an area of approximately 20 km² in one day using a satellite stereo-pair. Moreover, VHR stereoscopic satellite images can cover a much larger area in a single frame than traditional aerial photographs. In fact, the minimum area of a stereoscopic commercial satellite image is 100 km². When studying these surface phenomena at the basin scale or at least across very large areas, the wider overview that is provided by satellite images can limit the omission of features that are typically missed during field investigations [5]. Moreover, this method is more useful than field surveys for estimating ephemeral erosion channels in remote areas or in places that are hard to reach; these estimates can even be made months after the occurrence of ephemeral gullies and rills and when the features have been destroyed by farming activities.

The proposed method involves many aspects, including the aforementioned skill of the operator, the spatial resolutions and the quality of the stereoscopic satellite images, and the time span during which the correct illumination is achieved. The spatial resolution has a strong influence on the size of the linear erosion features that can be detected [63]. The quality of the image is fundamental for detecting small features, such as rills or ephemeral gullies, and includes photographic aspects such as the contrast, brightness and haze. In particular, the ability to detect and map rills and gullies is related to the contrast between the channel and the surrounding objects and can be emphasized dynamically changing the contrast. Furthermore, the stereovision allows data to be obtained at an apparent resolution that is better than the nominal resolution (i.e., pixel size; 50 cm in our case). Despratts [63] showed values that were smaller than the nominal resolution. This may be due to the simple geometrical effect of a 3D view with respect to a 2D view (single image). During the internal process of reconstructing an image, the brain can actually increase the resolution. Unfortunately, we do not have sufficient data to control this aspect with proper ground truth. Hence, we decided to only explore channels larger than 50 cm.

7. Conclusions

In this study, we proposed a new methodology for mapping and measuring rills and ephemeral gullies using very high-resolution GeoEye-1^® panchromatic stereo images. This new methodology was successfully applied and tested to recognize, map and measure relatively small channels, such as rills and ephemeral gullies. The proposed method is faster than classical field surveys, improves the ability to map these features over large areas, and complements UAV-based surveys, which are applicable to detailed scales and analyses, and other more traditional methods. Moreover, the availability of archive images allows this method to be used to map ephemeral channels that were removed by farming activities, which improves the analysis of erosion over time.

The results that were obtained in this study were consistent with a local field survey and, more generally, with results from the literature, which confirms that satellite stereo pairs can provide reliable estimates of ephemeral gully and rill erosion. The field-surveyed rill and gully data plot within the error band around the satellite-based volume-length relationship. Furthermore, the values of the coefficients of the volume-length relationship (a = 0.082; b = 1.174) are within the range of ephemeral gully values that were observed by other authors. The method was also applied to estimate the eroded volume from a large area (48 km²) by incorporating the uncertainty into the V-L relationship based on the uncertainty that is associated with the channel depth measurements. The calculated erosion value (mixed approach) at the scale of fields that were affected by channeled erosion (1.14 km²) was 28.4 Mg∙ha⁻¹ over 2.5 months. At the catchment scale (48 km²), the erosion value is 0.68 Mg∙ha⁻¹. In addition to the methodological results, this study provides the first estimate of gully and rill erosion in a medium-size mountain catchment in Central Italy and complements numerical studies that have been conducted in other areas, such as Sicily [56].

The proposed methodology can easily be applied wherever stereoscopic satellite images are available provided that the illumination conditions and channel width/depth constraints are satisfied. We consider this methodology of using very high-resolution stereoscopic satellite images to be particularly useful for building extensive and reliable databases of ephemeral erosional features over large areas.