minerals Article Soil Erosion Monitoring in Quarry Restoration Using Drones Vicenç Carabassa 1,2 , Pau Montero 1 , Josep Maria Alcañiz 1,2 and Joan-Cristian Padró 2, * 1 Ecological and Forestry Applications Research Centre (CREAF), E08193 Bellaterra, Spain; v.carabassa@creaf.uab.cat (V.C.); p.montero@creaf.uab.cat (P.M.); josemaria.alcaniz@uab.cat (J.M.A.) 2 Departament de Geografia, Edifici B, Universitat Autònoma de Barcelona, E08193 Bellaterra, Spain * Correspondence: joancristian.padro@uab.cat; Tel.: +34-935-814-343 Abstract: Mining is an essential activity that supports the provision of raw materials. However, the extraction process of mining has deep environmental impacts. For this reason, restoration actions are mandatory, and monitoring is a key step in ensuring the renaturalization of affected areas. Erosion processes are one of the main problems that affect restored areas in extractive activities due to the frequently steep slopes and the difficulty of revegetating the technosols constructed using mining debris. This research aims to develop a method for determining soil losses due to water erosion in mine-restored areas by using Geographic Information Systems (GIS) and Remote Sensing (RS) tools. For the study, images obtained using Unmanned Aerial Systems (UAS) in an open pit mine in the process of restoration are used, from which the Digital Elevation Model (DEM) of the current state of the slopes is obtained (0.10 m spatial resolution). With GIS techniques, ridges of the rills and gullies generated in the slopes are detected, whereby an estimation of a first DEM before the erosive process and a second DEM after the erosive process can be constructed. Each of these DEMs are evaluated individually in order to determine the height differences and estimate the volumetric loss. At the   same time, the results are validated with the DEM derived from official mapping agencies’ airborne Citation: Carabassa, V.; Montero, P.; Lidar data (1.00 m spatial resolution), which yield consistent data in the volumetric quantification of Alcañiz, J.M.; Padró, J.-C. Soil Erosion the erosion despite the difference in spatial resolution. In conclusion, the high spatial resolution of Monitoring in Quarry Restoration drone images facilitated a detailed monitoring of erosive processes, obtaining data from vast and Using Drones. Minerals 2021, 11, 949. inaccessible slopes that are usually immeasurable with traditional field techniques, and altogether https://doi.org/10.3390/ improving the monitoring process of mine restoration. min11090949 Keywords: mine restoration; soil degradation; DEM; drones; GIS tools; Remote Sensing Academic Editor: Javier Fernández Lozano Received: 24 July 2021 1. Introduction Accepted: 27 August 2021 Published: 30 August 2021 Mineral resources are tremendously important to economic development in terms of raw materials and energy [1,2]. However, mining can cause severe damage to the Publisher’s Note: MDPI stays neutral landscape, soil and water, especially in extraction areas. Sustainable mining is a complete with regard to jurisdictional claims in process that occurs during the lifespan of a mine, including the restoration of the affected published maps and institutional affil- land. In the restoration process, it is particularly crucial to monitor geotechnical risks, slope iations. instabilities and soil losses, where action occurs in order to control key aspects, such as slope erosion [3]. Among the various types of soil erosion, gully erosion represents the severest and has significant on-site and off-site effects [4]. In a broad sense, gully erosion can be defined as an erosion process in which deep channels are generated by runoff water that remove topsoil to a certain depth [5]. In a more detailed sense, gully erosion is a Copyright: © 2021 by the authors. threshold-dependent process that is controlled by a set of geo-environmental factors [6]. Licensee MDPI, Basel, Switzerland. This article is an open access article The eroded terrain caused by water erosion is characterized by a convergent topogra- distributed under the terms and phy that is associated with micro valleys (gullies) and separated by a series of divergent conditions of the Creative Commons ridges and drainage divisions (elevated surfaces). Therefore, the basic geomorphological Attribution (CC BY) license (https:// structure of any river landscape can be described by a net of bottom valleys with stream creativecommons.org/licenses/by/ channels and a net of ridges with drainage divisions [7]. In eroded slopes where a network 4.0/). of rills, gullies and ditches was generated, a very detailed spatial sampling is necessary to Minerals 2021, 11, 949. https://doi.org/10.3390/min11090949 https://www.mdpi.com/journal/minerals Minerals 2021, 11, 949 2 of 11 quantify and model the erosive process. It is at this scale of study where Remote Sensing (RS)-based on images captured by Unmanned Aerial System (UAS) takes on full meaning. UAS-based RS has become more modular, miniaturized and intelligent in recent years. With its unique advantages, such as high spatial resolution, low cost, short revisit cycle, efficient acquisition and easy operation, it is widely employed in various fields [8–10]. In recent years, the development of scientific research based on UAS data has enjoyed great progress and success in its application, even in combination with other Remote Sensing platforms, such as satellites or manned aircrafts [11,12]. However, the use of UAS applied on mining environments remains in constant development [13,14]. With their ease of use as well as time and cost efficiency in acquiring information, UAS equipped with cameras and sensors are increasingly replacing traditional methods. Thus, as unmanned aerial vehicle technology has advanced by leaps and bounds in the past decade, today, its use is indispensable in mining work and research. In this context, the studies focus on generating precision topography, 3D reconstruction, mining prospecting, controlling the opening of excavation fronts, monitoring geological risks and environmental control. The recent widespread availability of Digital Elevation Models (DEMs) with high spatial resolution, such as those derived from digital photogrammetry and laser altimetry (Lidar), has facilitated the mapping of morphometric changes in the terrain. The joining of classical photogrammetric methods and novel Structure from Motion (SfM) techniques from UAS imagery, now makes possible DEM generation in an affordable way, planning imagery acquisition with enough overlap to reconstruct 3D point structures and further interpolate them into grid DEMs [11–15]. Generating DEMs with photogrammetric tech- niques requires a large number of control points to ensure its quality, either on the ground or photointerpreted from a previously acquired image [16]. Surveying DEMs with centime- ter resolution derived from drone images allows the detection of elevation changes due to high precision, facilitating the monitoring of soil erosion processes [17]. The objective of the work is to generate detailed cartography of topography variations from images obtained by UAS as well as determine the loss of soil on the slopes of the mining area due to erosion by using Geographic Information Systems (GIS) and Remote Sensing tools. For this purpose, a methodology was developed starting from a DEM of high spatial resolution, which is based on the detection of watershed lines that are formed in the erosive channels on the slopes. The novelty of the method resides in the estimation of a DEM before the erosive process based on the interpolation of the points located on the ridges, taking advantage of the ultra-high spatial resolution of drone data. With these altitudinal references, the DEM prior to erosion is modeled to estimate the volume of eroded soil. 2. Materials and Methods 2.1. Study Area The area of interest in which the drone flight was carried out is the extractive activity on the Lázaro quarry, which is located in the municipality of El Vendrell (Tarragona, Spain). The exploitation has an area of 60 ha and a perimeter of 3.70 km (Figure 1). The quarry can be divided into two types of lithologies. The northern zone is formed by a mass of limestone, dolomitic and marly limestone, from the Lower Cretaceous and the southern zone by an outcrop of calcarenites belonging to the Miocene. Usually, tertiary materials do not present erosion problems. However, in cases where there is anthropic activity and clearing techniques, they can suffer erosion and caving phenomena on the slopes due to the water flow acting on it. The potency of material with commercial value is approximately 106 m of hard limestone rocks [18]. The study is focused on a subarea that is located northeast of the quarry and has an area of 23,189 m2 . This subarea, which was restored in the early 2010s, is characterized by the formation of steep slopes (>30◦ ) with a predominance of the slope-berm model. The lower slopes, which were restored Minerals 2021, 11, x FOR PEER REVIEW 3 of 11 Minerals 2021, 11, 949 3 of 11 predominance of the slope-berm model. The lower slopes, which were restored using using waste waste material material from from the the quarry quarry itself, itself, show show no change no change and and are are practically practically devoiddevoid of of vege- vegetation tation and and strongly strongly affected affected by erosive by erosive processes. processes. Figure 1.1. Location Figure Location map map ofof the the limestone limestone quarry quarry studied. studied. (a) (a) General General location location of of the the study study area area in in Catalonia, Catalonia, NE NE Iberian Iberian Peninsula. (b) Aerial space constrictions around the study area. (c) Exploitation limits and detail of the restored area. Peninsula. (b) Aerial space constrictions around the study area. (c) Exploitation limits and detail of the restored area. (d) (d) Limits of the current extractive activity and limits of the restored areas. (e) Detail of the eroded slopes focused on the Limits of the current extractive activity and limits of the restored areas. (e) Detail of the eroded slopes focused on the study. study. (f) Detail of the slope-berm restoration and erosion. (f) Detail of the slope-berm restoration and erosion. Minerals 2021, 11, x FOR PEER REVIEW 4 of 11 Minerals 2021, 11, 949 4 of 11 2.2. 2.2.Unmanned UnmannedAerial AerialSystems Systems(UAS) (UAS) The Theunmanned unmannedaircraft aircraftused used waswas an an Atmos-7 Atmos-7 [19] [19] with with flexible flexible fixed fixed wings. Having Having an anautonomy autonomy ofofupupto to70 70min, min,thetheAtmos-7 Atmos-7allows allowsaaflight flightceiling ceiling of of4000 4000m mand andaatelemetry telemetry range range of of 5000 5000 m. The sensor used used for fordata datacollection collectionwas wasthetheMicaSense MicaSense RedEdge-M RedEdge-M cam- camera, era, withwith a professional a professional multispectral multispectral sensor sensor capable capable of simultaneously of simultaneously capturing capturing fivefive dis- discrete crete spectral spectral bandsbands to generate to generate precise precise and quantitative and quantitative information information on theon the vegetation vegetation cover cover on theon the slopes slopes (Figure(Figure 2). The2).flight The flight was carried was carried out onout on 8 January 8 January 2020, 2020, at at a height a height of 120 of m 120 abovem above ground,ground, taking taking off andoff and landing landing on an esplanade on an esplanade located alocated a few kilometers few kilometers from the quarry. from theWith no clouds quarry. With noorclouds wind, or thewind, goodtheweather conditions good weather allowedallowed conditions the flight thetoflight take place between 11:15 and 12:00 without repetitions (Figure to take place between 11:15 and 12:00 without repetitions (Figure 2). 2). Figure2. Figure (a)Fixed-wing 2.(a) Fixed-wingdrone drone platform platform used used for for the the mission. mission. (b) (b) Sensor Sensor used used to to capture capture aerial aerial data. data. From the overflight of the Lázaro quarry, a multispectral image and the DEM with a From the overflight of the Lázaro quarry, a multispectral image and the DEM with a resolution of 0.10 m were obtained. Since the main purpose of the flight carried out was to resolution of 0.10 m were obtained. Since the main purpose of the flight carried out was co-register the resulting cartography with the official cartography, georeferencing was car- to co-register the resulting cartography with the official cartography, georeferencing was ried out by using 10 control points obtained by photointerpretation. The current orthophoto carried out by using 10 control points obtained by photointerpretation. The current ortho- (0.25 m) of the Cartographic and Geologic Institute of Catalonia (ICGC) [20,21] for the plani- photo (0.25 m) of the Cartographic and Geologic Institute of Catalonia (ICGC) [20,21] for metric coordinates (XY), and the DEM of 2.00 m × 2.00 m was also taken as a reference the planimetric coordinates (XY), and the DEM of 2.00 m × 2.00 m was also taken as a from the ICGC for the altimetric coordinate (Z). Three test points were used for achiev- reference from the ICGC for the altimetric coordinate (Z). Three test points were used for ing geometric accuracy in line with the reference pixel measurement (RMSEx = 0.11 m; achieving geometric accuracy in line with the reference pixel measurement (RMSEx = 0.11 RMSEy = 0.19 m; RMSEz = 0.45 m). m; RMSEy = 0.19 m; RMSEz = 0.45 m). Therefore, with these materials, it is possible to obtain detailed enough imagery and Therefore, digital models towith these materials, accurately it is possible detect ridges to obtain and water detailed troughs enoughthem and manage imagery with and GIS digital models to accurately detect ridges and water troughs tools to obtain volumetric information of the eroded terrain. and manage them with GIS tools to obtain volumetric information of the eroded terrain. 2.3. Data Processing 2.3. Data TheProcessing photogrammetric processing of the mosaic of the drone images, the georeferencing, The photogrammetric the orthophoto generation, processing the 3D pointofcloud the mosaic and theofDEM the drone images,out was carried thewith georeferenc- Structure ing, fromthe orthophoto Motion generation,using (SfM) techniques the 3D point Phostoscan Agisoft cloud and the DEM was carried photogrammetric out with software [11]. Structure It is worth from Motion noting that(SfM) techniques the DEM using Agisoft was obtained usingPhostoscan only the 3D photogrammetric points classified soft- as ware ground[11]. inItAgisoft is worth noting that Photoscan, a the DEM was obtained categorization performedusing onlySfM in the theprocedure 3D points classified thorough as ground in differential geometrical Agisoft Photoscan, slopes anda categorization performed pixel radiometric in the features. SfMthese From procedure data, thor- three ough geometrical differential slopes and pixel radiometric features. From these different ways of calculating the eroded area were applied: (1) relief inversion and drainage data, three different network ways of calculating detection, the eroded i.e., detection areainwere of ridges applied: elevation the original (1) relief inversion model; (2)andthedrain- Flow Accumulation age Model (FAM) network detection, [22]; and i.e., detection of (3) the Relative ridges Slopeelevation in the original Position model model;(RSP) (2) the[23]. FlowA Minerals 2021, 11, x FOR PEER REVIEW 5 of 11 Minerals 2021, 11, 949 5 of 11 Accumulation Model (FAM) [22]; and (3) the Relative Slope Position model (RSP) [23]. A workflow consistingof workflow consisting oftwelve twelvechained chainedprocesses processeswas was applied applied (Figure (Figure 3).3). ForFor thethe GISGIS andand Remote Sensing Sensingprocessing, processing,MiraMon MiraMon Version 8.2 [24] and QGIS (SAGA GIS) [25] Version 8.2 [24] and QGIS (SAGA GIS) [25] software soft- ware werefor were used used for generating generating RSP and RSP FAMand FAM models, models, which can which can be implemented be implemented using modelusing builderbuilder model tools (Figure 3). tools (Figure 3). Figure Figure 3.3. Comprehensive workflowofofthe Comprehensive workflow thedescribed describedmethodology, methodology, starting starting from from thethe drone drone imagery imagery and its Structure from Motion (SfM) processing, following with the obtention of the and its Structure from Motion (SfM) processing, following with the obtention of the intermediate intermediate products (DEMof products (DEM ofthe thecurrent currentrelief, relief,3D 3Dpoint pointcloud cloudand and orthomosaic), orthomosaic), thethe GIS GIS treatment treatment of the of the DEMDEM to obtain the ridge network, the estimation of the DEM before the erosive processes and the DEM to obtain the ridge network, the estimation of the DEM before the erosive processes and the DEM of of Differences (DoD) between the current relief DEM and the estimated original DEM. Differences (DoD) between the current relief DEM and the estimated original DEM. The The domains domains of of the the different different drainage drainagethresholds thresholdsin inthe theridge ridgepattern patternmake makeititpossible pos- to delimit sible the accumulated to delimit the accumulated runoff flowflow runoff according according to atowell-defined a well-defined criterion. criterion. With Withlow thresholds, low thresholds, parallel ridges parallel ridges without withoutintermediate intermediatecells cells occur withwithgreater greaterperiodicity, periodicity, whereas with with high high thresholds, thresholds,they theyomit omitdrainage drainagedividers, dividers, and the threshold and the threshold selection selection was such such that thatititdetected detectedthe thelargest largest number number of ridges. of ridges. Here, Here, it should it should be noted be noted thattodue that due thethe to complex complex nature of the nature surface, of the the the surface, process is only process an approximation is only an approximation of theofposition of the position thethe of ridges. In this ridges. complex In this complexand very detailed and very drainage detailed network, drainage orthophotography network, orthophotography is used is to check used which scene-dependent to check which scene-dependent thresholdthreshold value is best to represent value is best tothe specific landscape represent the specific features. landscape features. To obtain To obtainthetheridge ridgenetwork network inin vector vectorformat, format,thethe drainage drainageorder was was order defined, resulting defined, result- ing in a raster that assigns a value of first order to the primary sections and second the in a raster that assigns a value of first order to the primary sections and second order to order secondary to sections. the secondary The ridge sections. Themesh ridgewasmesh thenwas vectorized. To subsequently then vectorized. build a DEM To subsequently build a without ruts or gullies, a buffer was generated around the ridges, and subsequently, a DEM without ruts or gullies, a buffer was generated around the ridges, and subsequently, selection was made of the 3D point cloud obtained from the flight with respect to these a selection was made of the 3D point cloud obtained from the flight with respect to these polygons, thereby obtaining a new point cloud limited only to the points near the crest polygons, thereby obtaining a new point cloud limited only to the points near the crest (without the points in the troughs of the slopes). To quantify the volume of rills and (without the points in the troughs of the slopes). To quantify the volume of rills and gul- gullies, two DEMs estimated before erosion were generated, which were defined as the lies, two DEMs Potential Elevation estimated before erosion Model (DEMP) and thewere generated, Filtered ElevationwhichModelwere defined (DEMF). TheasDEMPthe Po- tential Elevation Model (DEMP) and the Filtered Elevation was modelled from the 3D points adjacent to the ridges using a Triangular Irregular Model (DEMF). The DEMP was modelled from the 3D points adjacent to the ridges using Network (TIN) and Multilevel B-Spline [26], which subsequently transformed into a rastera Triangular Irregular Net- work (TIN) and Multilevel B-Spline [26], which subsequently format resulting in a DEMP with the approximate height values as those before the erosion transformed into a raster format incisions. resulting Apart fromin a this, DEMP a DEMwithfiltering the approximate of the current height statevalues as those of the slopes was before the ero- performed sion with incisions. a median Apart from this, morphological a DEM filtering operator. The median of the currentisstate operator of the slopes a non-linear was is filter that per- formed with a median used to eliminate morphological noise from images, wherebyoperator. The median optimal use of thisoperator is a non-linear filter eliminates positivefilter and negative that is used tovalued spikes, eliminate noisebeing from also an effective images, whereby depression optimalfillusetechnique. of this filter Each pixel eliminates of the original positive raster isvalued and negative replaced in the spikes, destination being raster bydepression also an effective the outcome fill of the central technique. Each tendency pixel value, of the thus raster original resulting in the DEMF. is replaced in theFinally, a Lidarraster destination pointby cloud the was outcomefiltered, from of the cen- tral tendency value, thus resulting in the DEMF. Finally, a Lidar point cloud was filtered, Minerals 2021, 11, 949 6 of 11 which the ground points were extracted and duplicate points were eliminated. Once the classified point cloud was obtained, the Elevation Model was created from the derived Lidar data (DEML). The estimation of the height difference is a quantification of the value of the elevation changes within the slopes caused by erosion and deposition. It is calculated from the difference between the DEM of the current state of the ground and the DEM prior to erosion (Equation (1)). DoD = DEM − DEM (p) (1) To calculate the volumes of gullies, the method applied was the subtraction of the DEM from the current state and the DEM prior to the erosive process. The volume calculation was carried out considering a global detection of the changes in the slopes, an uncertainty equal to the spatial resolution of the data (0.10 m) and a minimum detection level, which is obtained from the estimation of the error of the DEMs [27]. There are two different directions of volume change present in the study area: on the one hand, the continuous erosion leading to a negative volume change (material losses), and on the other, the positive values which represent interpolation artefacts and do not necessary correspond to a gain of material. This is because the methodology proposed in this work interpolates areas with concave terrain morphology where erosion has occurred. Therefore, convex areas of the relief (i.e., ridges between gullies, colluvial fan or other landforms) are detected as part of the relief unaffected by erosion processes. The validation of such detailed data is difficult to carry out with field measurements. Therefore, we take as reference values for the erosion the DoD between airborne Lidar data obtained in 2010 and 2017 by the local official mapping agency (Institut Cartogràfic i Geològic de Catalunya, ICGC, Barcelona, Spain) [20,21]. Lidar data allowed the generation of two DEM of 1 m spatial resolution to quantify the volumetric differences in the 2010–2017 period. 3. Results 3.1. Gully Detection Methods The automated detection of ridges using the proposed methodology for relief inversion is consistent with the visualization of the DEM and its derived products. As can be seen in Figure 4, buffer polygons derived from the drainage network of the inverted relief model are related to the areas of higher elevation along streams. These polygons represent the estimation of the topography before the erosion of the relief (Figure 4a). In this study, the flow accumulation model has values ranging from 0.01 for the lowest accumulation pixels to 185 for highest accumulation pixels. A threshold of up to 5 was used to select DEM pixels that were not affected by the effects of erosion from surface runoff (Figure 4b). Choosing these parameters and limits is variable, depending on the morphological characteristics of the relief. The threshold determination is scene dependent and needs the human criteria to check with the ad-hoc orthophotography the ridges the algorithm performance. In the RSP model, values range from 0 to 1. The value 0 represents the lowest relative height, while 1 represents the part with the highest relative topography. Pixels greater than 0.7 were selected, which represent the area of the relief that was not modified by erosion. Figure 4c shows the relative position of the relief. Minerals2021, Minerals 2021,11, 11,949 x FOR PEER REVIEW 77 of of 11 11 Figure4.4.Detail Figure Detailof ofthe thecartographic cartographicrepresentation representationof ofthe thedifferent differentmethods methodsfor fordetecting detectinggully gullyzones. zones.(a) (a)Inverse InverseMDE MDE method, method, in ochre, the buffer polygon. (b) Flow accumulation method, values range from 0 to 5 is represented intransparent in ochre, the buffer polygon. (b) Flow accumulation method, values range from 0 to 5 is represented in transparent whichcorresponds which correspondstotothe thezones zoneswith withaalow lowflow flowaccumulation accumulationor orequivalent equivalenttotothe thezone zonewithout withouterosion. erosion.(c) (c)Relative Relativeslope slope positionmethod, position method,where wherevalues valuesrange rangefrom from0.7 0.7toto11represents representsareas areasnot notaffected affectedbybyerosion. erosion. 3.2. DEM 3.2. DEM of of Differences Differences(DoD) (DoD) Theresult The resultof ofthe thedifference differencebetween betweenthe themodels models(DoD) (DoD)isisshown shownin inFigure Figure5.5.The Thefirst first case shows case shows thethe interpolation interpolation of of the the pre-erosion pre-erosion DEM DEM calculated calculated by bythe therelief reliefinversion inversion method(Figure method (Figure5a,d). 5a,d). The The second second case caseshows showsthetheinterpolation interpolationofofthe themodel modelthrough throughthe the FAM (Figure 5b,e). Finally, the result of the model interpolated with the RSP FAM (Figure 5b,e). Finally, the result of the model interpolated with the RSP is shown in is shown in Figures 5c,f. In these height difference models, the areas of greatest erosion are Figure 5c,f. In these height difference models, the areas of greatest erosion are highlighted highlighted incolor, in color,whereby wherebyred redindicates indicatesaacorrelation correlationwith withthe thedeepest deepestgullies gullies(Figure (Figure5).5).The Theflow flow accumulation method accumulation method shows shows the thelowest lowesterosion erosionvolume, volume,while whilethe themodel modelmademadewithwiththe the reliefinversion relief inversiontogether togetherwith withthe theRSPRSPmodel model interpolated interpolated byby M-SPLINES M-SPLINES shows shows thethe high- highest est amount amount of eroded of eroded volume volume as well as well as eroded as eroded areaarea (Table (Table 1). 1). Table 1. Statistical data for the different gully detection methods used. The reference data for validation are the DEM of Differences (DoD) at 1.00 m spatial resolution obtained from airborne Lidar data sensed in 2010 and 2017 (local official mapping agency [20,21]). Raster Interpolation Volume (m3 ) Area (m2 ) Min. Value (m) Mean Value (m) Stand. Dev. (m) M-Spline Relative Slope Position −5679.42 20,125.88 −5.44 −0.23 0.38 TIN Relief Inversion −5343.70 15,829.37 −4.41 −0.22 0.38 M-Spline Relief Inversion −5194.99 16,302.30 −4.41 −0.21 0.39 TIN Relative Slope Position −3342.82 13,594.35 −5.71 −0.13 0.29 M-Spline Flow Accumulation Model −1155.10 15,649.33 −2.05 −0.05 0.08 TIN Flow Accumulatio Model −541.96 7987.75 −3.19 −0.02 0.06 Reference LIDAR 2010–2017 DoD −2689.09 10,696.00 −3.87 −0.11 0.27 Minerals2021, Minerals 2021,11, 11,949 x FOR PEER REVIEW 8 8ofof11 11 Figure5.5.DEM Figure DEMofof Differences Differences(DoD) using (DoD) different using combination different of relief combination modelling of relief and interpolation modelling methods: and interpolation (a) through methods: (a) relief inversion through and Splines relief inversion andinterpolation; (b) flow accumulation Splines interpolation; and Splines (b) flow accumulation andinterpolation; (c) relative Splines interpolation; (c)slope position relative slope and po- Splines interpolation; sition and (d) relief inversion Splines interpolation; (d) reliefand TIN interpolation; inversion (e) flow accumulation and TIN interpolation; and TIN interpolation; (e) flow accumulation (f) relative and TIN interpolation; (f) relative slope slope position andposition and TIN interpolation. TIN interpolation. The colors The colors the represent represent depth the depth of the of the gully. gully. Values Values below 0.2 below 0.2shown to 0 are to 0 are in shown in transparent. transparent. 4. Discussion The results obtained in this study confirm the possibility of mapping erosion processes on slopes using different techniques for automatic gully detection. However, the quality of the results is strongly linked to the spatial resolution of the DEM used to perform the Minerals 2021, 11, 949 9 of 11 work as well as the terrain control points used to adjust the DEM itself. Precise erosion estimation depends on the relative erosion rates on the slopes, and the mapping precision depends on the spatial resolution of the data as well as the vertical accuracy. So far, other studies have determined that a resolution of 0.05 m is sufficient for obtaining more than 85% of the information. However, this will ultimately depend on the study site and the geomorphology [2]. Among the models studied in this work, the flow accumulation model sufficiently detects the small gullies of the slopes, but it does not accurately detect the sectors where intense erosion occurred. On the other hand, by using both the RSP and the detection of crests by relief inversion, it is possible to detect the areas where greater erosion occurred, but the areas with more incipient erosion are not well detected. This information could be complemented with the photo interpretation of the images for delimiting the gullies [28]. However, photointerpretation can be time-consuming and reduces the efficiency of the protocol, especially for monitoring large areas, such as those existing in extractive activities. Moreover, the use of Lidar data, if available, from rutinary airborne flights could also improve the proposed methodology [15]. However, the low point density of some rutinary flights (<2 point m−2 ) make difficult the use this methodology, at least for detecting incipient erosion processes. Drone data give much more 3D point density, allowing us to work up to scales of 1:500 [12], although severely affected areas could be detected using Lidar data [21]. Vegetation creates problems when calculating the volume eroded on the land, as it generates artefacts that do not correspond to eroded areas. This problem could be partially solved with the exclusion of vegetation patches through vegetation classification using certain classification protocols [3,12]. It should also be noted that in the areas where the slope changes, i.e., from slope to plain, an artefact is produced. In models using RSP and inverted relief, such artefacts can be detected and measured. However, in flow accumulation models, such artefacts are not detectable (Figure 5). Interpolation by M-Splines tends to generate a greater difference in heights and occupy more area affected by erosive processes. However, with TIN interpolation, the surfaces are more conservative, and therefore, the difference between DEMs is usually smaller. The comparison of the drone-acquired photogrammetric results with the airborne Lidar reference shows that the volumetric result of Relative Slope Position method is the most similar to airborne Lidar reference among the methods tested, but regarding the pixel minimum value, the Flow Accumulation and the Relief Inversion methods are closer to the reference. The validity and uncertainties of the results generated by this methodology should be evaluated in depth with field data surveys of the gullies on the slopes, so as to determine the success of the results in more detail. However, fieldwork can be difficult or even impossible in areas having steep and unstable slopes. For future research on debris erosion detection in slopes, one can envisage using higher spatial resolution together with control points on the ground in order to minimize the uncertainty in volume calculation. 5. Conclusions In this study, a practical method to estimate the volume of eroded material in restored areas of an open pit mine is presented, using UAS imagery, photogrammetric and Remote Sensing techniques. Geographical Information Systems are the basis of the presented methodology, and we provide a workflow easily implementable in most GIS software thorough model builders. From the results, it is possible to obtain maximum and minimum estimations of the total volume on restored but eroded slopes and sedimentation, being especially relevant to identify where the main erosion processes take place. This approach facilitates the quantitative assessment of erosion processes, becoming a potential reliable tool for restoration programs, mining landscapes and many more slope erosion scenarios. Minerals 2021, 11, 949 10 of 11 Author Contributions: Conceptualization, J.-C.P., J.M.A. and V.C.; methodology, J.-C.P. and P.M.; software, J.-C.P. and P.M.; validation, V.C., J.-C.P. and P.M.; formal analysis, V.C., J.-C.P. and P.M.; investigation, V.C., J.-C.P. and P.M.; resources, J.-C.P., J.M.A. and V.C.; data curation, J.-C.P. and P.M.; writing—original draft preparation, V.C., J.-C.P. and P.M.; writing—review and editing, V.C., J.-C.P., J.M.A. and P.M.; visualization, J.-C.P. and P.M.; supervision, J.-C.P., J.M.A. and V.C.; project administration, V.C. and J.M.A.; funding acquisition, V.C. and J.M.A. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the Government of Catalonia through the project “Research and Innovation in the process and restoration of extractive activities”, by the European Union through the LIFE project “New Life for Drylands” (LIFE20 PRE/IT/000007), by Exodronics SL and by the Spanish Government through the Centro para el Desarrollo Tecnológico Industrial (CDTI) under the StratoTrans Project [IDI-20191270] (Proyecto cofinanciado por el Centro para el Desarrollo Tecnológico e Industrial (CDTI) en el marco de la convocatoria de Proyectos de I + D de Transferencia Tecnológica CERVERA). Acknowledgments: The authors acknowledge the support of COMERCIAL LÁZARO SA. Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. References 1. Ren, H.; Zhao, Y.; Xiao, W.; Hu, Z. A review of UAS monitoring in mining areas: Current status and future perspectives. Int. J. Coal Sci. Technol. 2019, 6, 320–333. [CrossRef] 2. Xiao, W.; Chen, J.; Zhao, J.; Hu, Z.; Zhang, S. Identify Maize chlorophyll impacted by coal mining subsidence in high groundwater table area based on UAS remote sensing. J. China Coal Soc. 2019, 44, 302–313. [CrossRef] 3. Carabassa, V.; Montero, P.; Crespo, M.; Padró, J.C.; Pons, X.; Balagué, J.; Brotons, L.; Alcañiz, J.M. Unmanned aerial system protocol for quarry restoration and mineral extraction monitoring. JEMA 2020, 270, 110717. [CrossRef] 4. Vanmaercke, M.; Poesen, J.; Van Mele, B.; Demuzere, M.; Bruynseels, A.; Golosov, V.; Bezerra, J.F.R.; Bolysov, S.; Dvinskih, A.; Frankl, A.; et al. How fast do gully headcuts retreat? Earth-Sci. Rev. 2016, 154, 336–355. [CrossRef] 5. Zabihi, M.; Mirchooli, F.; Motevalli, A.; Khaledi Darvishan, A.; Pourghasemi, H.R.; Zakeri, M.A.; Sadighi, F. Spatial modelling of gully erosion in Mazandaran Province, northern Iran. Catena 2018, 161, 1–13. [CrossRef] 6. Nicu, I.C. Is digital shoreline analysis system “fit” for gully erosion assessment? Catena 2021, 203, 105307. [CrossRef] 7. Lindsay, J.B.; Seibert, J. Measuring the significance of a divide to local drainage patterns. IJGIS 2013, 27, 1453–1468. [CrossRef] 8. Esposito, G.; Mastrorocco, G.; Salvini, R.; Oliveti, M.; Starita, P. Application of UAS photogrammetry for the multi-temporal estimation of surface extent and volumetric excavation in the Sa Pigada Bianca open-pit mine, Sardinia, Italy. Environ. Earth Sci. 2017, 76, 103–118. [CrossRef] 9. Guirado, R.; Padró, J.-C.; Zoroa, A.; Olivert, J.; Bukva, A.; Cavestany, P. StratoTrans: Unmanned Aerial System (UAS) 4G Communication Framework Applied on the Monitoring of Road Traffic and Linear Infrastructure. Drones 2021, 5, 10. [CrossRef] 10. Jaramillo, R.; Padró, J.-C. Generación de cartografía a partir de imágenes captadas con dron de ala fija, asociada a proyectos hidráulicos fluviales. GeoFocus 2020, 26, 93–117. [CrossRef] 11. Padró, J.C.; Muñoz, F.J.; Avila, L.A.; Pesquer, L.; Pons, X. Radiometric Correction of Landsat-8 and Sentinel-2A Scenes Using Drone Imagery in Synergy with Field Spectroradiometry. Remote Sens. 2018, 10, 1687. [CrossRef] 12. Padró, J.C.; Muñoz, F.J.; Planas, J.; Pons, X. Comparison of four UAV georeferencing methods for environmental monitoring purposes focusing on the combined use with airborne and satellite remote sensing platforms. IJAEO 2019, 79, 130–140. [CrossRef] 13. Honarmand, M.; Shahriari, H. Geological Mapping Using Drone-Based Photogrammetry: An Application for Exploration of Vein-Type Cu Mineralization. Minerals 2021, 11, 585. [CrossRef] 14. Jackisch, R.; Lorenz, S.; Zimmermann, R.; Möckel, R.; Gloaguen, R. Drone-Borne hyperspectral monitoring of acid mine drainage: An example from the Sokolov Lignite District. Remote Sens. 2018, 10, 385. [CrossRef] 15. Niculita, M.; Mărgărint, M.C.; Tarolli, P. Using UAV and LiDAR data for gully geomorphic changes monitoring. Dev. Earth Surf. Process. 2020, 23, 271–315. [CrossRef] 16. Padró, J.C.; Carabassa, V.; Balagué, J.; Brotons, L.; Alcañiz, J.M.; Pons, X. Monitoring opencast mine restorations using Unmanned Aerial System (UAS) imagery. Sci. Total Environ. 2019, 657, 1602–1614. [CrossRef] 17. Eltner, A.; Mulsow, C.; Maas, H.G. Quantitative Measurement of Soil Erosion from TLS and UAS Data. ISRPS 2013, 40, 119–124. [CrossRef] 18. INECO. Ministerio de Transporte, Movilidad y Agenda Urbana. Available online: https://www.mitma.gob.es/recursos_mfom/ paginabasica/recursos/12t3750_03-geologia-y-proc-de-materiales.pdf (accessed on 21 June 2021). 19. CATUAV. Available online: https://www.catuav.com/drones/atmos7-3/ (accessed on 21 June 2021). Minerals 2021, 11, 949 11 of 11 20. ICGC (Institut Cartogràfic i Geològic de Catalunya). Full 273–135 de l’Ortofoto de Catalunya 1:2500 (OF–25C) (1:2500 Catalonia Orthophoto Tile 273–135). Available online: http://www.icc.cat/vissir3/ (accessed on 21 June 2021). 21. ICGC (Institut Cartogràfic i Geològic de Catalunya). Lidar Data. Available online: https://www.icgc.cat/en/Downloads/ Elevations/Lidar-data (accessed on 11 July 2021). 22. Jenson, S.K.; Domingue, J.O. Extracting Topographic Structure from Digital Elevation Data for Geographic Information System Analysis. ASPRS 1988, 54, 1593–1600. 23. Boehner, J.; Selige, T. Spatial prediction of soil attributes using terrain analysis and climate regionalisation. In SAGA–Analysis and Modelling Applications; Boehner, J., McCloy, K.R., Strobl, J., Eds.; Goettinger Geographische Abhandlungen: Göttingen, Germany, 2006; pp. 13–28. 24. Pons, X. MiraMon. Sistema d’Informació Geogràfica i Software de Teledetecció. Versió 8.2d. Centre de Recerca Ecològica i Aplicacions Forestals, CREAF, Bellaterra 84-931323-4-9. 2019. Available online: http://www.creaf.uab.cat/miramon/ (accessed on 21 June 2021). 25. Conrad, O.; Bechtel, B.; Bock, M.; Dietrich, H.; Fischer, E.; Gerlitz, L.; Wehberg, J.; Wichmann, V.; Boehner, J. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. Geosci. Model Dev. 2015, 8, 1991–2007. [CrossRef] 26. Lee, S.; Wolberg, G.; Shin, S.Y. Scattered Data Interpolation with Multilevel B-Splines. IEEE Trans. Vis. Comput. Graph. 1998, 3, 228–244. [CrossRef] 27. Williams, R. DEMs of difference. Geomorphol. Technol. 2012, 2, 1–17. 28. Fiorucci, F.; Ardizzone, F.; Rossi, M.; Torri, D. The use of stereoscopic satellite images to map rills and ephemeral gullies. Remote Sens. 2015, 7, 14151–14178. [CrossRef]