Measuring rill erosion using structure from motion: A plot experiment

Catena 156 (2017) 383–392

Contents lists available at ScienceDirect

Catena journal homepage: www.elsevier.com/locate/catena

Measuring rill erosion using structure from motion: A plot experiment a

b

a

C. Di Stefano , V. Ferro , V. Palmeri , V. Pampalone a b

a,⁎

MARK

University of Palermo, Department of Agricultural and Forestry Sciences, Viale delle Scienze, 90128 Palermo, Italy University of Palermo, Department of Earth and Marine Sciences, Via Archirafi, 90123 Palermo, Italy

A R T I C L E I N F O

A B S T R A C T

Keywords: Rill erosion Field experiment Rill characteristics Image-based ground measurement technique 3D Photo-reconstruction DEM

In this paper the results of a plot experiment on rill erosion are reported. The rill network, manually incised on the soil and further shaped by a clear inflow discharge, was surveyed using the three-dimensional photoreconstruction (3D-PR) technique which allows to obtain a digital terrain model (DTM) by a large series of oblique images of the channel from consumer un-calibrated and non-metric cameras. The three-dimensional (3D) DTM and the quasi-tridimensional (2.5D) model were generated by Agisoft Photoscan software. For a single rill channel, the reliability of the 3D image-based ground measurements of morphological and hydraulic variables was positively tested by the corresponding measurements carried out on the rill gypsum cast. Moreover, the morphological and hydraulic variables were measured at 11 transects of the plot by a profilometer P, the 2.5D and 3D models. Using the 3D model as reference, the analysis showed that the reliability of the measurements by P and 2.5D methods was comparable and the error was generally lower than 15%. The rill measurements carried out by the three methods agreed with the available literature measurements and supported the applicability of both the empirical relationship between rill length and its eroded volume and the theoretical dimensionless relationship among the morphological variables describing the channelized erosion process. The analysis also showed that the effect of the distance between two consecutive cross sections on the rill volume measurement by 3D model is negligible. Finally, the shaping effect of different flow discharges on three straight rills incised into the plot was tested. In particular, three increasing flow discharges (0.15, 0.35 and 0.5 l s− 1) were used sequentially in two rills while the maximum discharge of the sequence (0.5 l s− 1) was used into the third rill. Both the grain-size distribution of the eroded sediment and the total amount of soil loss were affected by the discharge sequence.

1. Introduction Rill erosion is particularly spread on agriculturally disturbed dry areas and it has been demonstrated to be the main origin of soil loss at many sites (Bruno et al., 2008; Cerdan et al., 2002; Di Stefano and Ferro, 2011; Di Stefano et al., 2013; Hancock et al., 2008; Polyakov and Nearing, 2003; Wirtz et al., 2012). Furthermore, the transport of eroded particles from hillslopes to river channels is more efficient if a rill network is developed (Aksoy et al., 2013; Berger et al., 2010; Brunton and Bryan, 2000). Rill erosion is considered one of the most important processes affecting soil because of both the large amount of soil loss due to this phenomenon and the generation of a rill network which can produce stable and persisting landforms such as gullies (Wirtz et al., 2010). Rills are small concentrated flow paths in which water depth and flow velocity are higher than within an overland flow (Govers et al., 2007). As a consequence, rill flow is characterized by more energy for

⁎

entrainment and transport of eroded soil particles leading to an increase of erosion rates. In such shallow flow, where water depth is typically of the order of millimeters to several centimeters, bed topography significantly affects flow hydraulics (Abrahams et al., 1996; Foster et al., 1984; Gilley et al., 1990; Govers, 1992; Nearing et al., 1997; Peng et al., 2015). Bruno et al. (2008) carried out field investigations under natural rainfalls and established the proportion of both rill erosion and interrill erosion on total soil loss. The measurements demonstrated that rill erosion increases the total sediment transport efficiency because rill flow is able to transport both interrill eroded particles and sediments which are eventually detached from the rill wetted perimeter. On natural hillslopes rill and gullies can show positive and negative changes of ground level, as a consequence of a combination of strong channel incision and sediment accumulation due to bank collapsing (Vergari et al., 2013). At present the experimental work carried out in laboratory (Aksoy

Corresponding author. E-mail address: [email protected] (V. Pampalone).

http://dx.doi.org/10.1016/j.catena.2017.04.023 Received 21 November 2016; Received in revised form 11 April 2017; Accepted 21 April 2017 0341-8162/ © 2017 Elsevier B.V. All rights reserved.

Catena 156 (2017) 383–392

C. Di Stefano et al.

et al., 2013; Berger et al., 2010; Bryan and Poesen, 1989; Brunton and Bryan, 2000; Mancilla et al., 2005;) and in field condition (Bruno et al., 2008; Cerdan et al., 2002; De Santisteban et al., 2005; Wirtz et al., 2010, 2012) used both natural and simulated rainfalls applied to soils with different texture. These researches aimed to study rill network formation, to define the initial condition for rilling (Bruno et al., 2008; Mancilla et al., 2005; Torri et al., 1987), to recognize the development of rill head morphology (Bruno et al., 2008; Brunton and Bryan, 2000), to measure or estimate the main rill morphologic and hydraulic variables (cross-section area, wetted perimeter, hydraulic radius, mean velocity, shear stress) (Abrahams et al., 1996; Bruno et al., 2008; Carollo et al., 2015; Di Stefano et al., 2013; Nearing et al., 1997). Notwithstanding the efforts made for understanding rill erosion, there is still much work to be done for measuring and modelling this erosion process and to upscale all processes studied at plot scale to natural hillslopes. Accurate and repeatable measurements of rill erosion processes are required both for understanding and realizing a correct rill erosion modelling. Rill erosion measurements at plot scale by the traditional direct method (i.e. metric ruler, rillmeter, topographic instruments) modify the surface of the rilled area and are time consuming and invasive (Casalí et al., 2006), while ground measurement of rill features using a drone-based technology (Carollo et al., 2015) is considered a noninvasive method allowing a fast field survey. According to Govers et al. (2007) rill experiments “provide an opportunity to investigate to what extent the concepts used in models are a truly valid description of the erosion processes occurring”. The uncertainty of process-oriented models, due to the difficulty to estimate some parameters (such as shape and size of the cross-sections) and the unavailability of some variables (such as the soil erodibility), stimulated the proposal of empirical relationships, using simple input variables such as the channel length, for estimating the eroded volume (Bruno et al., 2008; Di Stefano et al., 2013; Toy et al., 2002). Di Stefano and Ferro (2011), using the rill and gully measurements carried out at Sparacia experimental area and available literature data, established the following equation for estimating the total eroded volume V expressed in m3:

V = ao L bo

Cumulative Frequency

Fig. 1. View of the experimental area.

Lr3, s

⎛ w H ⎞ nr = ar ⎜⎜ 2 ⎟⎟ ⎝ Lr , s ⎠

0.9 0.8 Sample_1 Sample_2 Sample_3 Sample_4

0.7 0.6 0.5 0.0001

0.001

0.01

0.1

1

10

Particle diameter (mm) Fig. 2. Grain-size distribution of samples collected in the experimental area.

based modelling creates a digital terrain model (DTM) using a set of oblique photographs taken from the same rilled surface (Frankl et al., 2015) and without using control points during the composition of the model (Gómez-Gutiérrez et al., 2014). Camera model parameters and scene geometry are simultaneously solved (James and Robson, 2012), using redundant information coming from oblique photographs. The generated point cloud has to be scaled and georeferenced using some ground control points located within the monitored rilled area. Main advantages of 3D-PR method are both that it requires little expertise, because the image processing is almost automatic (James and Robson, 2012), and that it has a high accuracy, as already demonstrated for ephemeral gullies channels by using Terrestrial Laser Scanner (TLS) method as benchmark (Castillo et al., 2012; Gómez-Gutiérrez et al., 2014; Kaiser et al., 2014) and for closed earth channels by direct measurement of their volume (Di Stefano et al., 2017). Semi-automated 3D-PR techniques, such as the coupled use of Structure from Motion (SfM) and MultiView-Stereo (MVS) workflows (Javernick et al., 2014; Seiz et al., 2006), are integrated in software packages such as Photoscan (Agisoft) (Frankl et al., 2015). Measuring rill erosion by an image-based technique requires the use of a 3D model which is able to represent the three-dimensional topography also including areas into which different z-elevations can be associated to the same couple of plane coordinates x-y. On the contrary, the 2.5D model represents the surface of the rilled area by a vertical perspective; in other words, 2.5D model cannot represent crosssections having undercut walls (Frankl et al., 2015) in which the coordinates x–y of a ground measurement point can correspond to different z elevations. In this case, the 2.5D model underestimates the area of the cross-section and, as a consequence, the eroded rill volume (Di Stefano et al., 2017). However, the 2.5D model is often used for hydrological analyses and it can be easily managed being a raster model.

(1)

in which L is the total length (m) of the channel and a0 and b0 are two numerical coefficients. Di Stefano and Ferro (2011) established that Eq. (1) can be applied for rills, ephemeral gullies (EGs) and gullies using the same exponent b0 = 1.1 and a different scale factor a0, equal to 0.0036 for rills. Bruno et al. (2008), using the dimensional analysis and the incomplete self-similarity theory, recognized that the erosive process producing a rill segment having a length Lr,s, a volume Vr,s, a width measured at the top of the cross section w, and a maximum scour depth H, can be expressed by the following dimensionless relationship

Vr , s

1.0

(2)

in which ar and nr are two constants. Di Stefano and Ferro (2011) calibrated Eq. (2) using 1929 rill, EGs and gully segment measurements and demonstrated that Eq. (2) with ar = 0.5341 and nr = 0.9379 can be indifferently applied for rills, EGs and gullies. Close-range photogrammetry can be effectively used to measure soil erosion processes using recent advancements in automatic 3D-photo reconstruction (3D-PR) techniques. These techniques use oblique images acquired from uncalibrated and non metric cameras coupled with photogrammetric software packages (Castillo et al., 2015). Terrestrial digital photogrammetry has a little impact on ground measurements, is a technique characterized by high spatial resolution (centimeters to millimeters) and it requires a time-consuming work for post-processing the digital photos acquired in the field activity. Image384

Catena 156 (2017) 383–392

C. Di Stefano et al.

Rill A

a)

Rill B

b)

Fig. 3. Plane scheme of the rill network with the numbered cross-sections (a) and rill shaping by an inflow discharge (b).

a). Comparing the 3D model of a rill channel obtained by 3D-PR method with the rill gypsum cast used for directly measuring morphological and hydraulic variables of the cross sections and the rill volume; b). Comparing the accuracy of the 2.5D model and the profilometer method, P, with respect to the 3D model for calculating morphological and hydraulic variables; c). Comparing the ground measurements obtained by the three different methods to the rill length-volume relationship and the dimensionless morphological relationship; d). Evaluating the effect of the distance interval i on rill volume measurement by 3D model; e). Studying the shaping process of three rills by different inflow discharges in terms of morphological evolution and grain size distribution of the transported sediment.

For measuring channelized (rill, ephemeral gully and gully) erosion, the channel is usually divided into a number of segments bounded by two consecutive cross sections and its eroded volume is calculated by adding the segment volumes. Therefore, it is necessary to determine the suitable distance between two consecutive measurement cross sections, i.e. the distance interval i, to properly measure channelized erosion. Recent studies (Casalí et al., 2006; Castillo et al., 2012; Di Stefano et al., 2017) assessed the effect of the distance interval i on the channel volume error. Castillo et al. (2012) introduced the variable i/L, which is the number of measured cross sections per unit channel length and proposed a model for establishing the i/L value to be used in field survey in order to obtain a given volume error magnitude with a given probability level. Di Stefano et al. (2017), analyzing the results of their investigation and those by Casalí et al. (2006) and Castillo et al. (2012), suggested that the influence of the ratio i/L on the gully volume error is limited for a large, narrow and deep gully and i/L ≤ 0.49. Soil erosion causes loss of soil fertility by a direct loss of the topsoil and of plant nutrients. In soils most chemicals are adsorbed by clay and organic matter because of their high surface area. When rainfall reaches the surface horizon of the soil, some chemicals are desorbed and go into solution while others remain adsorbed and move with the soil particles (Di Stefano et al., 2005). These eroded sediments may be transported by runoff until discharged into a water body such as a reservoir. Once agricultural pollutants enter a water system, they lower water quality and impose economic losses on water users (Di Stefano and Ferro, 2016). Therefore, attention has to be directed to the grain size distribution of the eroded sediments, especially for the usually high soil losses associated with rill erosion. The clay fraction of the eroded sediment is the site for most nutrients and chemical adsorption. The capacity of transported particles to carry chemicals is closely related to their specific surface area, which is in turn largely dependent on the clay content of the transported particles and aggregates. The specific objectives of this paper concern the measurement of rill morphology by different techniques, the analysis of the amount of eroded particles and the study of the texture of the transported sediments. In particular:

2. Materials and methods 2.1. Experimental plots Two experimental plots, located at the experimental area of the Department of Agriculture and Forestry Sciences (Fig.1) of Palermo University, 2 m wide and 7 m long were used for carrying out the investigation on incised rills. The plots, having a slope equal to 14%, are bounded by concrete walls and have a base realized by gabions. Each plot was filled with a clay soil (73% clay, 22.5% silt and 4.5% sand), collected at Sparacia experimental area (Bagarello et al., 2008, 2013). Before starting of the test, a shallow tillage was carried out on the clay soil. Fig. 2 shows the grain-size distribution of 4 samples collected in the experimental area. To avoid the detachment and transport of soil particles by inflow discharge at the rill head, an inflow device was used for breaking the water jet and allowing the water flow to enter into each incised rill as a light stream. At the plot downstream end is located a storage tank, having the same width of the plot and a cross-section 0.7 m wide and 0.4 m high. For testing the accuracy of the ground measurements carried out by 385

Catena 156 (2017) 383–392

C. Di Stefano et al.

2.2. Using of the rill gypsum cast The high accuracy of the image-based ground measurements was previously demonstrated for ephemeral gullies channels by using TLS method as benchmark (Castillo et al., 2012; Gómez-Gutiérrez et al., 2014; Kaiser et al., 2014) and for closed earth channels by direct measurement of their volume (Di Stefano et al., 2017). In this investigation, the accuracy of the image-based ground measurements for a rill channel was checked by realizing a gypsum cast (Fig.4a) of the Rill A. The gypsum was used to fill the rill volume because it is easy to be shaped when it is wet, and easy to take out of the rill and to cut when it dries. For allowing the extraction of the dried gypsum cast without destroying the rill, the rill boundary was covered by a plastic film. The dried gypsum cast was cut at 11 cross sections having a longitudinal distance interval i equal to 0.624 m (Fig.4b). The gypsum cast cross-sections were photographed and then scaled and digitized (Fig.4c) by CAD software for measuring the maximum depth H, defined as the maximum distance from bottom and surface of crosssection, the surface width w, the cross-section area A, the cross-section perimeter WP and the hydraulic radius R. The volume of each rill segment, Vr,s, was calculated by the following relationship:

Vr , s = 0.5(Ai + Ai +1 ) i

(3)

in which Ai is the initial cross-section area of the rill segment, and Αi + 1 is the final cross-section area. Eq. (3) was applied since a direct measurement of the cast volume requires the measurement of both the cast weight and specific weight. In particular, the gypsum cast presented some internal voids which determined a high variability of its specific weight. Moreover, the top surface of the gypsum cast was irregular and it could determine measurement errors comparable with the errors due to the applied method (Eq. (3)). In other words, the direct measurement of the cast volume presented uncertainties due to both the top surface irregularities and how the cast was internally made. 2.3. Comparing ground measurement techniques The incised rill network constituted by Rill A, Rill B and five tributary reaches was used for analyzing the accuracy of the morphological measurements by ground techniques. The ground measurements were carried out using 11 plot transects having a longitudinal distance interval i, measured along the plot length (6.24 m), equal to 0.624 m (Fig.3a). The transects intersects Rill A in correspondence of the 11 cross sections measured by the gypsum cast. The profilometer P was constructed by 25 pins, each having a diameter of 4 mm, installed on an aluminum bar and spaced 2 cm apart (Bruno et al., 2008). This spacing is adequate to measured rill width (7.5–20 cm) and to the roughness of the investigated rills. The pin configuration and consequently rill section geometry was photographed and pin heights were digitized from these pictures for obtaining the values of the morphological variables of the cross section. The 3D-PR technique was used for creating a DTM using a set of 70 photographs taken from the plot area by a digital camera (Samsung PL210–14.2 Mpx) enabling automatic focusing (Gómez-Gutiérrez et al., 2014; Westoby et al., 2012). Images were captured assuring that any part of the measured plot area was represented in at least 3 photographs for taking into account that the 3D algorithm is designed to work with convergent images. The image-processing software Photoscan Professional (by Agisoft, St. Petersburg, Russia) was used for obtaining the 3D model by an automatic process which couples SfM and MVS techniques (Di Stefano et al., 2017; Kaiser et al., 2014). Then, using the point cloud of the 3D model, a 2.5D surface (i.e. Digital Elevation Model; DEM) having a pixel size equal to 2 mm was generated by Photoscan. The extraction of the cross section profile was carried out by Cloud Compare software (http://cloudcompare.org) for the point cloud and by ArcMap for DEM while the geometric variables were measured by CAD.

Fig. 4. View of the gypsum cast of the main Rill A (a), two cutting lines for determining cross sections (b) and the 11 photographed and digitized cross-sections of the gypsum cast (c).

the image-based technique and the wide-spread profilometer (P) method, a rill network composed by two main rills (named Rill A and Rill B) and five tributary reaches were firstly manually cut on the plot (Fig.3a) using a small hoe. Then, the incised rills were further shaped by a clear inflow discharge (Fig.3b).

386

Catena 156 (2017) 383–392

C. Di Stefano et al.

30

8 6

20

H 3D (cm)

w 3D (cm)

+10%

-10% 10

0

4 2 0

0

10

20

30

0

2

w (Gypsum cast) (cm)

6

8

40

100

WP 3D (cm)

80

A 3D (cm2)

4

H (Gypsum cast) (cm)

60 40 20

30 20 10 0

0 0

50

100

0

10

A (Gypsum cast) (cm2)

20

30

40

WP (Gypsum cast) (cm)

5

6000

Vr,s 3D (cm3)

R 3D (cm)

4 3 2 1 0 0

1

2

3

4

4000

2000

0

5

0

R (Gypsum cast) (cm)

2000

4000

6000

Vr,s (Gypsum cast) (cm3)

Fig. 5. Comparison between the values of the cross section geometric variables and segment volume measured by 3D model and gypsum cast.

incised on the plot using a small hoe for creating a channel 5 cm wide and 2–3 cm deep. The channel was further shaped by a clear inflow discharge. For Rills 1 and 2 the experimental runs were carried out using a sequence of three inflow discharges Q equal to 0.15, 0.35 and 0.5 l s− 1 while for the Rill 3 a single discharge, equal to 0.5 l s− 1, was used. The inflow discharge was measured by the volumetric method. The runs were stopped when the same volume of the suspension (liquid + solid fractions) (308 l) was stored into the plot storage tank. Four DTMs were created using a set of oblique photographs taken after the manual cutting of the rill (R0 ground measurement) and at the end of each experimental run for a given inflow discharge (R0.15, R0.35 and R0.50 ground measurements). Using cross-sections with i equal to 0.0624 m, the generated 3D models allowed to calculate the volume of soil eroded by a known inflow discharge. This volume was converted into weight using the bulk density value equal to 1.13 kg m− 3 which was the mean of the measured values for five undisturbed samples collected into the plot before starting the experimental run. The sampling was carried out using cylindric samples having a diameter of 5 cm and a height equal to 5 cm, which were ovendried at 105 °C for 48 h. Finally, this sediment weight measured by the 3D model was compared with that stored into the tank at the end of each run. This last sediment weight was determined by oven-drying at 105 °C the 308 l of the stored suspension (liquid + solid fractions) into the tank.

For each of the 11 investigated rill cross-sections, the surface width w, the maximum depth H, the cross-section area A and perimeter WP, and the hydraulic radius R were determined using the three ground measurement methods (3D, 2.5D and P). For 3D and 2.5D methods, the surface width was determined establishing the distance between the points of the two banks where an abrupt change in slope occurred (Casalí et al., 2015). When more than one point of slope inflection was detected analyzing the cross section profile, this criterion was applied neglecting the local ground irregularities. The rill segment volume Vr,s was calculated by Eq. (3). The total volume V of a rill was calculated adding the volumes Vr,s of all segments into which the rill was divided. The estimate error E(Vr,s), expressed as percentage, of each segment volume with respect to a reference value Vr, was calculated by the following equation:

⎛ V − Vr ⎞ E (Vr , s ) = 100 ⎜ m ⎟ ⎝ Vr ⎠

(4)

in which Vm is the measured segment volume and the reference value Vr is equal to the volume measured by the 3D model. 2.4. Field measurements of rill erosion For testing the shaping effect of the flow discharge on preformed rill channels, three straight rills, named Rill1, Rill2 and Rill3, having geometrical characteristics similar to Rill A and Rill B shown in Fig.3b were incised. Each rill, having a length equal to 6.24 m, was manually 387

Catena 156 (2017) 383–392

C. Di Stefano et al.

0.3

0.1 +15%

H(m)

w (m)

0.2 -15% 2.5D

0.1

0.05

P 0

0 0

0.1

0.2

0.3

0

0.05

w 3D (m)

0.1

H 3D (m)

0.01

0.4

WP (m)

A (m2)

0.3 0.005

0.2 0.1 0

-5.2E-18 0

0.005

0

0.01

0.2

0.3

0.4

0.006

Vr,s (m3)

0.05

R (m)

0.1

WP 3D (m)

A 3D (m2)

0.025

0

0.004

0.002

0 0

0.025

0.05

0

0.002

0.004

0.006

Vr,s 3D (m3)

R 3D (m)

Fig. 6. Comparison, for the investigated variables, between the 3D reference values, and those determined by 2.5D and P methods.

depth H, cross-section area A and perimeter WP, hydraulic radius R) and the volume Vr,s of each rill segment obtained by 3D model were compared with those measured by the rill gypsum cast (Fig. 5). This test showed that 3D model is very accurate because the estimate error of each variable is generally lower than 10%. Therefore, in the following analysis the measurements obtained by the 3D model will be used as “reference” values.

Cumulative Frequency

1 0.8 0.6 0.4 0.2

2.5D

3.2. Field measurements by 2.5D and P methods

P

For each investigated variable (w, H, Α, WP, R and Vr,s) the comparison between the 3D reference values and the 2.5D and P measured values is plotted in Fig. 6. For each variable both ground measurement methods (2.5D and P) are almost always characterized by measurement errors lower than 15% and the 2.5D method is more accurate than the widely applied P method. The differences between cross section widths measured by 2.5D and the reference values could be due to the carelessness in determining abrupt slope changes implicit into DEM accuracy (Di Stefano et al., 2017). For the 2.5D method the frequency distribution of the estimate error E(Vr,s), plotted in Fig. 7, shows that the median value is quasi equal to zero, the errors vary from − 0.2 to 0.1 and 80% of them are negative. For the P method, the median error is negative (− 0.1), the errors vary from − 0.32 to 0.15 and the method tends to systematically underestimate the volume (only 4 out of 41 segment volumes are overestimated) (Fig.7). The relationship length L–volume V corresponding

0 -0.4

-0.3

-0.2

-0.1

0

0.1

0.2

E (Vr,s) Fig. 7. Cumulative frequency distribution of the estimate error E (Vr,s) for 2.5D and P methods.

3. Results 3.1. Tests by rill gypsum cast In this investigation, for testing the accuracy of the image-based ground measurements for a rill channel, the values of both the morphological and hydraulic variables (surface width w, maximum 388

Catena 156 (2017) 383–392

C. Di Stefano et al.

0.05

0.05 Rill A

Rill B

0.04

V (m3)

V (m3)

0.04 0.03

3D 0.02

2.5D P

0.01

0.03 0.02 0.01

b)

a)

0.00

0.00 0

2

4

6

8

0

2

L (m)

4

6

8

L (m)

Fig. 8. Relationship L-V for rill A (a) and rill B (b), determined by 3D, 2.5D and P methods.

Fig. 11 demonstrates, for the 3D model and for the investigated cases (0.0624 ≤ i ≤ 1.248 m, L = 6.24 m) for which the i/L ratio ranges from 0.01 to 0.20, that the distance between two subsequent crosssections affects the calculated volumes. In particular, using as reference condition the minimum investigated value i = 6.24 cm, when the distance interval increases cases of volume overestimation (Rill A) or underestimation (Rill B) can be detected. In any case the distance interval has a moderate effect on the L-V relationship. Respect to the reference condition (i = 6.24 cm), the mean of absolute measurement errors of the two total volumes was equal to 10.5% for i/L = 0.1 and 15.6% for i/L = 0.2. Casalí et al. (2006), investigating six natural rills by a profilometer and tape, found either positive and negative volume errors and established that using i values < 4 m guarantees that the errors remain below 10%. Taking into account that the length of the investigated rills ranged between 20.4 m and 29.4 m, the corresponding i/L values were equal to 0.2 and 0.14, respectively. Rill A and Rill B were characterized by values of the cross section area coefficient of variation equal to 0.33 and 0.35, respectively, which are practically coincident with those determined by Casalí et al. (2006). Therefore, the differences between the volume errors determined by Casalí et al. (2006) and in the present investigation cannot be explained by cross section variability. In any case, according to both investigations, for i/ L ≤ 0.2 the measurement of rill volume can be considered slightly dependent on the ratio i/L.

1

V(m3)

0.1 0.01 0.001 Di Stefano et al. (2013) Rill A Rill B eq. 1

0.0001 0.00001 0.1

1

10

100

L(m) Fig. 9. Comparison between Eq. (1) and the L-V pairs measured by 3D method and by Di Stefano et al.. (2013).

to the three different methods (P, 2.5D and 3D) is shown in Fig. 8. For Rill A, this figure points out a light underestimation of the rill volume calculated by profilometer measurements. Fig. 9 shows that the pairs (L, V) measured by 3D model are near the curve of Eq. (1) calibrated by Di Stefano and Ferro (2011) (a0 = 0.0036 and b0 = 1.1), and fall within the data scattering corresponding to the calibration data set of Eq. (1). This expected result supports the idea that the severity of the rill erosion process can be estimated by the channel length even if the cross-section area variability can affect this estimate. The rill measurements carried out by 3D, 2.5D and P methods were used to test the applicability of Eq. (2) with ar = 0.5341 and nr = 0.9379. Fig. 10 confirms that Eq. (2), as calibrated by Di Stefano and Ferro (2011) is applicable and that this equation is independent of the used measurement method (3D, 2.5D, P). The 3D-PR technique allowed to measure rill cross sections choosing different distance interval values (i = 6.24, 62.4 and 124.8 cm) and to evaluate the effect on the volume measurement of the distance interval.

3.3. Tests with different rill flow discharges The general aim of the last part of this investigation was to test the effect on rill channel evolution of different rill flow discharges by using 3D-PR technique. The corresponding soil losses measured in the runs are indicated as SL1, SL2 and SL3 in Table 1. The results listed in Table 1 show that the soil loss measurements by 3D model are affected by errors ranging from − 6.4% to 16.8% as compared to the actual weight of the stored sediment in the tank. Moreover, using of 3D model yielded to a low estimate error of the total soil loss SLT (= SL1 + SL2 + SL3), ranging from −6.4% to 6.7%. The measurement errors by 3D model

0.1

0.1 Rill B

0.01

Vr,s/Lr,s3

Vr,s/Lr,s3

Rill A

3D 2.5D P eq.2

a) 0.001 0.005

b) 0.05

wH/Lr,s

0.01

2

0.001 0.005

0.05

wH/Lr,s

2

Fig. 10. Comparison, for Rill A (a) and Rill B (b), between Eq. (2) and the pairs (wH/Lr,s2, Vr,s/Lr,s3) measured by 3D, 2.5D and P methods.

389

Catena 156 (2017) 383–392

C. Di Stefano et al.

0.05

0.05

Rill A

Rill B

0.04

0.03

V(m3)

V(m3)

0.04

0.02 i = 6.24 cm i= 62.4 cm i= 124.8 cm

0.01

0.03 0.02 0.01

a)

b)

0

0 0

2

4

6

0

8

2

4

6

8

L (m)

L (m)

Fig. 11. Effect of the distance between two subsequent cross-sections on the L-V relationship for Rill A (a) and Rill B (b).

Cumulative Frequency

Table 1 Values of the soil losses, SLs, measured in the Rills 1, 2 and 3 by weighing and 3D model and measurement errors, E, of the 3D model with respect to weighing assumed as reference. Rill1

Weighing 3D model E (%)

Soil loss SL1 (g) 1180.2 1166.5 −1.2

Soil loss SL2(g) 1136.4 1317.4 15.9

Soil loss SL3 (g) 519.0 541.3 4.3

Total soil loss SLT (g) 2835.6 3025.2 6.7

Soil loss SL1 (g) 400.8 407.4 1.6

Soil loss SL2(g) 761.4 735.2 −3.4

Soil loss SL3 (g) 252.7 295.3 16.8

Total soil loss SLT (g) 1414.9 1437.8 1.6

–

–

– – –

–

Soil loss SL3 (g) 4136.2 3870.0 − 6.4

Rill2

Rill3

Weighing 3D model E (%)

0.9

a)

0.8 0.7

Q=0.15 l/s, rill1 Q=0.35 l/s, rill1 Q=0.5 l/s, rill1 Q=0.5 l/s, rill3

0.6 0.5 0.4 0.0001

0.001

0.01

0.1

1

10

Particle diameter (mm) Cumulative Frequency

Weighing 3D model E (%)

1.0

Total soil loss SLT (g) 4136.2 3870.0 − 6.4

(Table 1) are generally consistent with those obtained for Rill A and the rill gypsum cast as reference, as shown in the first part of present investigation. This result further confirms the reliability of the 3D model for measuring rill erosion. In the comparison between Rill 1 and Rill 2, the ratio SL2/SL1 varies from 0.96 to 1.9 while the ratio SL3/SL2 is similar and equal to 0.46 and 0.33, respectively. In both cases, about 80% of total soil loss SLT occurred in the first two runs (SL1 + SL2). This result can be explained by the “bulldozer effect” (Regűés et al., 2000; Seeger et al., 2004) which is “a mobilization of the loose material available in the rill before an experiment starts” (Wirtz et al., 2012). In particular, in the first run the flow discharge (0.15 l s− 1) and the corresponding transport capacity Tc,1 were low, as a consequence the rill flow was saturated (the actual sediment load was equal to the transport capacity Tc,1) by the sediment particles available for the “bulldozer effect” and the soil detachment due to rill flow was negligible. In the subsequent run (Q = 0.35 l s− 1) less loose material was available (the bulldozer effect was quasi finished) and the erosion process of the rill wetted perimeter was very active producing soil losses SL2 which ranges from SL1 to 2SL1. In the third sequential discharge (0.50 l s− 1) the bulldozer effect was finished, the transport capacity increased with the discharge and the soil loss decreased (Table 1). Therefore, since the transport capacity of the third run, Tc,3, was greater than the second one, Tc,2, and the corresponding soil loss SL3 was lower than SL2, the erosion process became detachment-limited. Increasing the discharge (from 0.15 to 0.50 l s− 1) the transport capacity increased too and the flow was able to transport particles having a larger size: the grain-size distribution of the transported sediments (sampled in the plot storage tank), shown in

1.0 0.9

b)

0.8 0.7 sediment

0.6

rill bed

0.5 0.4 0.0001

0.001

0.01

0.1

1

10

Particle diameter (mm) Fig. 12. Comparison among the grain-size distributions of sediment particles transported in the three sequential runs (Q = 0.15, 0.35 and 0.50 l s− 1) in the Rill 1 and that of the sediments transported in the Rill 3 (Q = 0.50 l s− 1) (a), and comparison between the grain-size distributions of the transported sediments and the rill bed after the run in the Rill 3 (b).

Fig. 12a as an example for the Rill 1, moves towards large-size particles when the flow discharge increases. The soil loss measurement carried out in Rill 3 using a single discharge equal to 0.50 l s− 1 (Table 1) was compared with those carried out in Rill 1 and Rill 2 for the third run (having the same values of discharge and runoff volume) of the sequence 0.15–0.350.50 l s− 1. This soil loss was always largely greater (from 8 to 16 times) than the corresponding value in the sequential run and this result can be justified taking into account that in the Rill 3 the flow was always characterized by the highest transport capacity Tc,3, while in the Rill 1 and 2 the transport capacity increased during the run. The flow in Rill 3 was able to simultaneously transport the loose material available in the rill before the experiment started and the detached particles from the rill wetted perimeter. The total soil loss for the Rill 3 was also greater than the total soil losses of the Rill 1 and 2, although its associated total runoff was 1/3 of those relative to Rill 1 and 2. Some process oriented models, such as WEPP, use the peak discharge calculated by the hydrologic component as steady discharge in the erosion component and the event duration is such that the total runoff is maintained. 390

Catena 156 (2017) 383–392

C. Di Stefano et al.

Cumulative Frequency

1

Therefore, according to the above results establishing that the soil loss corresponding to an increasing inflow discharge in sequence is less than that of a single discharge equal to the maximum value of the sequence, WEPP assumption is highly conservative for soil conservation purposes. From the soil loss point of view, the steady state condition is more burdensome than an increasing discharge sequence. Fig. 12a shows that the size of the particles transported by the flow in the Rill 3 increased respect to the size of the sediments transported in the three sequential runs of the Rill 1; in other words, the grain-size effect due to the single discharge 0.50 l s− 1 was greater than the sequence of three discharges (0.15–0.35-0.50 l s− 1). The flow moving in Rill 3 was able to transport larger particles originating from the bulldozer effect and the scour phenomena along the wetted perimeter. Furthermore, Fig. 12a shows that in the three sequential runs of Rill 1 the flow transport was selective and for increasing discharge values (i.e. for increasing transport capacity values from Tc,1 to Tc,3) more and more large particles were entrained into the flow. This mechanism is similar to an “armouring” effect produced by the small discharges that prevent the erosion of the greater ones. For the 0.50 l s − 1 discharge flowing in the Rill 3 the grain-size distribution of the transported sediment was coincident with that of the particles sampled into the rill channel at the end of the run (Fig. 12b). Therefore the transport capacity in the rill 3, Tc,3, was so high that the flow was able to transport all available particles independently of their size. The flow determined an equal-mobility condition of the particles, i.e. all particles, independently of their size, had the same probability to be entrained and transported by the flow. For assessing the morphological evolution of the rill channel, the initial bed configuration, i.e. before starting the run, was compared with the final one (at the end of the single run for Rill 3 or the sequence of runs for Rill 1 and 2). The empirical frequency distributions of the ratio between the cross-section area in the final condition, Α0.5, and in the initial condition, Α0, are compared in Fig. 13a. For Rill 3 this figure shows that 80% of the Α0.5/Α0 values is > 1 while this percentage is equal to 76% and 62% for Rill 1 and 2, respectively. The empirical frequency distributions of the ratio between the maximum depth in the final condition H0.5 and that corresponding to the initial condition H0 are compared in Fig. 13b. For Rill 3 this figure shows that 97% of the H0.5/H0 values is > 1 while this percentage is equal to 95% for Rill 1 and 81% for Rill 2. In other words, rill shaping by a single discharge (0.50 l s− 1) produced a more frequent deepening of the cross-section and increase of the cross-section area than the sequence of three discharges (0.15–0.35–0.50 l s− 1). Fig. 14 shows, as an example for a cross-section of Rill 1 (a), Rill 2 (b) and Rill 3 (c), the comparison between the initial condition R0 (Q = 0) and the final one R0.5 (Q = 0.50 l s − 1). The sequence of discharges (0.15–0.35–0.50 l s− 1) flowing into Rill 1 and 2 produced bed scouring and the consequent deepening of the cross-section (Fig. 14a, b); in other words the morphological evolution was characterized by the increase of the area by bed deepening and determined a deep cross-section. On the contrary, the discharge Q = 0.50 l s − 1 flowing into Rill 3 produced the increase of the cross-section area by widening and the morphological evolution determined a wide crosssection (Fig. 14c). In this case, the channel shaping processes acted prevalently on the banks and produced a cross-section enlargement by bank collapse.

rill 1 rill 2 0.1

0.01

rill 3

a) 0.5

1

A0.5/A0

1.5

2

Cumulative Frequency

1

0.1

0.01

b) 0.5

1

H 0.5/H0

1.5

2

Fig. 13. Empirical frequency distribution of the ratio Α0.5/Α0 (a) and the ratio H0.5/H0 (b).

R0 R0.5

Rill 1

a)

Rill 2

b)

4. Conclusions

Rill 3

The present investigation showed the results of a plot experiment aimed at measuring rill erosion using the 3D-PR technique. In particular, morphological and hydraulic variables of rill cross sections of a rill network and rill segment volumes were measured. With reference to a single main rill of the network, the measurements by the 3D model were practically coincident with the ones obtained by the cross section profiles of its gypsum cast. Therefore, the 3D method can be assumed as

c)

Fig. 14. Comparison, as an example for three cross-sections of Rill 1 (a), Rill 2 (b) and Rill 3 (c), between the initial condition (R0) and the final one (R0.5).

391

Catena 156 (2017) 383–392

C. Di Stefano et al.

583–594. Cerdan, O., Le Bissonnais, Y., Couturier, A., Bourennane, H., Souchère, V., 2002. Rill erosion on cultivated hillslopes during two extreme rainfall events in Normandy France. Soil Tillage Res. 67, 99–108. De Santisteban, L.M., Casalí, J., Lòpez, J.J., Giràldez, J.V., Poesen, J., Nachtergaele, J., 2005. Exploring the role of topography in small channel erosion. Earth Surf. Process. Landf. 30, 591–599. Di Stefano, C., Ferro, V., 2011. Measurements of rill and gully erosion in Sicily. Hydrol. Process. 25, 2221–2227. Di Stefano, C., Ferro, V., 2016. Establishing soil loss tolerance: an overview. J. Agric. Eng. XLVII 560, 127–133. Di Stefano, C., Ferro, V., Palazzolo, E., Panno, M., 2005. Sediment delivery processes and chemical transport in a small forested basin. Hydrol. Sci. J. 50, 697–712. Di Stefano, C., Ferro, V., Pampalone, V., Sanzone, F., 2013. Field investigation of rill and ephemeral gully erosion in the Sparacia experimental area, South Italy. Catena 101, 226–234. Di Stefano, C., Ferro, V., Palmeri, V., Pampalone, V., Agnello, F., 2017. Testing the use of an image-based technique to measure gully erosion at Sparacia experimental area. Hydrol. Process. 31, 573–585. Foster, G.R., Huggins, L.F., Meyer, L.D., 1984. A laboratory study of rill hydraulics: I. Velocity relationships. Trans. ASAE 27, 790–796. Frankl, A., Stal, C., Abraha, A., Nyssen, J., Rieke-Zapp, D., De Wulf, A., Poesen, J., 2015. Detailed recording of gully morphology in 3D through image-based modelling. Catena 127, 92–101. Gilley, J.E., Kottwitz, E.R., Simanton, J.R., 1990. Hydraulics characteristics of rills. Trans. ASAE 27, 797–804. Gómez-Gutiérrez, A., Schnabel, S., Berenguer-Sempere, F., Lavado-Contador, F., RubioDelgado, J., 2014. Using 3D photo-reconstruction methods to estimate gully headcut erosion. Catena 120, 91–101. Govers, G., 1992. Relationship between discharge, velocity and flow area for rills eroding loose, non layered materials. Earth Surf. Process. Landf. 17, 515–528. Govers, G., Giménez, R., Van Oost, K., 2007. Rill erosion: exploring the relationship between experiments, modelling and field observations. Earth-Sci. Rev. 84, 87–102. Hancock, G.R., Crawter, D., Fityus, S.G., Chandler, J., Wells, T., 2008. The measurement and modelling of rill erosion at angle of repose slopes in mine spoil. Earth Surf. Process. Landf. 33, 1006–1020. James, M.R., Robson, S., 2012. Straightforward reconstruction of 3D surfaces and topography with a camera: accuracy and geosciences applications. J. Geophys. Res. 117, 1–17. Javernick, L., Brasington, J., Caruso, B., 2014. Modeling the topography of shallow braided rivers using structure-from-motion photogrammetry. Geomorphology 213, 166–182. Kaiser, A., Neugirg, F., Rock, G., Müller, C., Haas, F., Ries, J., Schmidt, J., 2014. Smallscale surface reconstruction and volume calculation of soil erosion in complex Moroccan gully morphology using structure from motion. Remote Sens. 6, 7050–7080. Mancilla, G.A., Chen, S., McCool, D.K., 2005. Rill density prediction and flow velocity distributions on agricultural areas in the Pacific Northwest. Soil Tillage Res. 84, 54–66. Nearing, M.A., Norton, L.D., Bulgakov, D.A., Larionov, G.A., West, L.T., Dontsova, K.M., 1997. Hydraulics and erosion in eroding rills. Water Resour. Res. 33, 865–876. Peng, W., Zhang, Z., Zhang, K., 2015. Hydrodynamic characteristics of rill flow on steep slopes. Hydrol. Process. 29, 3677–3686. Polyakov, V.O., Nearing, M.A., 2003. Sediment transport in rill flow under deposition and detachment conditions. Catena 51, 33–43. Regüés, D., Balasch, J.C., Castelltort, X., Soler, M., Gallart, F., 2000. Relación entre las tendencias temporales de producción y transporte de sedimentos y las condiciones climáticas en una pequeña cuenca de montaña mediterránea (Vallcebre, Pirineos Orientales). Cuadernos de Investigación Geográfica 26, 41–65. Seeger, M., Errea, M.P., Beguería, S., Arnáez, J., Martí, C., García-Ruiz, J.M., 2004. Catchment soil moisture and rainfall characteristics as determinant factors for discharge/suspended sediment hysteretic loops in a small headwater catchment in the Spanish pyrenees. J. Hydrol. 288, 299–311. Seiz, S.M., Curless, B., Diebel, J., Scharstein, D., Szeliski, R., 2006. A Comparison and Evaluation of Multi-View Stereo Reconstruction Algorithms. IEEE Conference on Computer Vision ad Pattern Recognition. IEEE Computer Society, New York. Torri, D., Sfalanga, M., Chisci, G.M., 1987. Threshold conditions for incipient rilling. In: Bryan, R.B. (Ed.), Rill Erosion. Catena Suppl Vol. 8. pp. 97–105. Toy, T.J., Foster, G.R., Renard, K.G., 2002. Soil Erosion: Processes, Prediction, Measurement and Control. John Wiley & Sons Inc., New York. Vergari, F., Della Seta, M., Del Monte, M., Fredi, P., Lupia Palmieri, E., 2013. Long- and short-term evolution of several Mediterranean denudation hot spots: the role of rainfall variation and human impact. Geomorphology 183, 14–27. Westoby, M.J., Brasington, J., Glasser, N.F., Hambrey, M.J., Reynolds, J.M., 2012. ‘Structure-from-motion’ photogrammetry: a low-cost, effective tool for geoscience applications. Geomorphology 179, 300–314. Wirtz, S., Seeger, M., Ries, J.B., 2010. The rill experiment as a method to approach a quantification of rill erosion process activity. Z. Geomorphol. 54, 47–64. Wirtz, S., Seeger, M., Ries, J.B., 2012. Field experiment for understanding and quantification of rill erosion processes. Catena 91, 21–34.

reference to evaluate the reliability of 2.5D and P methods. Both ground measurement methods gave errors lower than 15%. The three measurement methods yielded to the same relationship length – volume, although with a light underestimation of the rill volume measured by P exclusively for the downstream reach of a main rill. Moreover, according to literature results, this relationship was slightly affected by the distance between two consecutive measurement cross sections and the length – volume measurements were consistent with the literature ones used for calibrating the relationship. The dimensionless morphological relationship was applicable to the measurements regardless of the used method (3D, 2.5D, P). Finally, the dynamics of three rill channels was investigated by using increasing inflow discharges in sequence or a single discharge equal to the maximum value of the sequence. The total soil loss in this last case was higher than the other two cases and the shaping channel mechanism was also different. In fact, in the run with a single discharge it acted prevalently on the banks and produced cross-section enlargement by bank collapse, while in the three sequential runs, erosion processes were mainly due to bed scouring and produced a consequent deepening of the cross-section. In the three sequential runs, the flow transport of the eroded sediments was selective with more and more large particles for increasing discharges, while an equal-mobility condition of the particles was determined in the run with a single discharge. In conclusion, the present analysis established that, with reference to both the sediment size and the total amount of soil loss, rill erosion is affected by the hydrograph. As a consequence, for improving accuracy of rill erosion prediction models, runoff temporal variability should be adequately simulated. Acknowledgements All authors contributed to outline the investigation, analyze the data, discuss the results, and write the manuscript. This research did not receive any specific grant from funding agencies in the public, commercial or not-for-profit sectors. References Abrahams, A.D., Gang, L.I., Parsons, A.J., 1996. Rill hydraulics on a semiarid hillslope, southern Arizona. Earth Surf. Process. Landf. 21, 35–47. Aksoy, H., Unal, N.E., Cokgor, S., Gedikli, A., Yoon, J., Koca, K., Inci, S.B., Eris, E., Pak, G., 2013. Laboratory experiments of sediment transport from bare soil with a rill. Hydrol. Sci. J. 58, 1505–1518. Bagarello, V., Di Piazza, G.V., Ferro, V., Giordano, G., 2008. Predicting unit plot soil loss in Sicily, south Italy. Hydrol. Process. 22, 586–595. Bagarello, V., Ferro, V., Pampalone, V., 2013. A new expression of the slope length factor to apply USLE-MM at Sparacia experimental area (Southern Italy). Catena 102, 21–26 (Special Issue “The effect of scale on soil erosion”). Berger, C., Schulze, M., Zapp, D.R., Schlunegger, F., 2010. Rill development and soil erosion: a laboratory study of slope and rainfall intensity. Earth Surf. Process. Landf. 35, 1456–1467. Bruno, C., Di Stefano, C., Ferro, V., 2008. Field investigation on rilling in the experimental Sparacia area South Italy. Earth Surf. Process. Landf. 33, 263–279. Brunton, D.A., Bryan, R.B., 2000. Rill network development and sediment budgets. Earth Surf. Process. Landf. 25, 783–800. Bryan, R.B., Poesen, J., 1989. Laboratory experiments on the influence of slope length on runoff, percolation and rill development. Earth Surf. Process. Landf. 14, 211–231. Carollo, F.G., Di Stefano, C., Ferro, V., Pampalone, V., 2015. Measuring rill erosion at plot scale by a drone-based technology. Hydrol. Process. 29, 3802–3811. Casalí, J., Loizu, J., Campo, M.A., De Santisteban, L.M., Àlvarez – Mozos, J., 2006. Accuracy of methods for field assesment of rill and ephemeral gully erosion. Catena 67, 128–138. Casalí, J., Giménez, R., Campo-Bescós, M.A., 2015. Gully geometry: what are we measuring? Soil 1, 509–513. Castillo, C., Pérez, R., James, M.R., Quinton, J.N., Taguas, E.V., Gómez, J.A., 2012. Comparing the accuracy of several field methods for measuring gully erosion. Soil Sci. Soc. Am. J. 76, 1319–1332. Castillo, C., Pérez, R., James, M.R., Gómez, J.A., 2015. SF3M software: 3-D photoreconstruction for non-expert users and its application to a gully network. Soil 1,

392