land Article Drone-Based Identification of Erosive Processes in Open-Pit Mining Restored Areas Joan-Cristian Padró 1, * , Johnsson Cardozo 1 , Pau Montero 2 , Roger Ruiz-Carulla 3 , Josep Maria Alcañiz 1,2 , Dèlia Serra 4 and Vicenç Carabassa 2 1 Grumets Research Group, Departament de Geografia, Edifici B, Universitat Autònoma de Barcelona, E08193 Bellaterra (Cerdanyola del Vallès), Catalonia, Spain; johnsson.cardozo@e-campus.uab.cat (J.C.); josemaria.alcaniz@uab.cat (J.M.A.) 2 Centre de Recerca Ecològica i Aplicacions Forestals (CREAF), E08193 Bellaterra (Cerdanyola del Vallès), Catalonia, Spain; p.montero@creaf.uab.cat (P.M.); v.carabassa@creaf.uab.cat (V.C.) 3 Division of Geotechnical Engineering and Geosciences, Department of Civil and Environmental Engineering, Universitat Politècnica de Catalunya-BarcelonaTech, E08034 Barcelona, Catalonia, Spain; roger.ruiz@upc.edu 4 RPAS Unit, Aerial Support Service, Rural Agents Corp Government of Catalonia, E08130 Santa Perpètua de Mogoda, Catalonia, Spain; delia.serra@gencat.cat * Correspondence: joancristian.padro@uab.cat; Tel.: +34-935-814343 Abstract: Unmanned Aerial Systems, or drones, are very helpful tools for managing open-pit mining operations and developing ecological restoration activities. This article presents a method for identi- fying water erosion processes in active quarries by means of drone imagery remote sensing, in the absence of pre-existing imagery or mapping for comparison. A Digital Elevation Model (DEM) with a spatial resolution (SR) >10 cm and an orthophoto with an SR >2.5 cm were generated from images captured with a drone and their subsequent photogrammetric processing. By using Geographical  Information Systems tools to process the DEM, a detailed drainage network was obtained, the areas  of detected water erosion were separated, and the watersheds in the gullies identified. Subsequently, Citation: Padró, J.-C.; Cardozo, J.; an estimated DEM before the erosive processes was reconstructed by interpolating the gully ridges; Montero, P.; Ruiz-Carulla, R.; Alcañiz, this DEM serves as a reference for the relief before the erosion. To calculate the volume of eroded J.M.; Serra, D.; Carabassa, V. material, the DEM of Differences was calculated, which estimates the volume difference between the Drone-Based Identification of Erosive Processes in Open-Pit Mining previously estimated DEM and the current DEM. Additionally, we calculated the material necessary Restored Areas. Land 2022, 11, 212. for the geomorphological adaptation of the quarry and the slope map, which are two valuable factors https://doi.org/10.3390/land11020212 closely related to the monitoring of erosive processes. The results obtained allowed us to identify the erosion factors quickly and accurately in this type of mining. In the case of water-filled quarries, it Academic Editors: Gergely Tóth, would be important to characterize the subsurface relief. Essentially, the presented method can be Eva Ivits and Bozena Smreczak applied with affordable and non-invasive materials to create digital grid maps at 10 cm resolution, Received: 31 December 2021 obtaining data ready for 3D metrics, being a very practical landscape modelling tool for characterizing Accepted: 26 January 2022 the restoration evolution of open-pit mining spaces. Published: 29 January 2022 Publisher’s Note: MDPI stays neutral Keywords: UAS; Unmanned Aerial Systems; mining; restoration; erosion; digital terrain models; with regard to jurisdictional claims in GIS tools published maps and institutional affil- iations. 1. Introduction Extracting raw materials such as sand, gravel, crushed stone, brick clays, gypsum, Copyright: © 2022 by the authors. and ornamental stones requires intensive and extensive efforts as all of these products Licensee MDPI, Basel, Switzerland. are obtained from open-pit mining operations [1]. The open-pit mining method is used This article is an open access article where large volumes of materials are mobilized, which gives rise to the formation of distributed under the terms and quarries [2]. Open-pit mining is one of the most land-degrading and landscape-modifing conditions of the Creative Commons activities. Public awareness of sustainable land management has been growing, and since Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ the 1990s, environmental policies and restoration actions have been increasingly focused on 4.0/). addressing environmental impacts [3]. In this type of mining, exploitation occurs in phases Land 2022, 11, 212. https://doi.org/10.3390/land11020212 https://www.mdpi.com/journal/land Land 2022, 11, 212 2 of 13 where the extraction of materials is linked with the ecological restoration of previously exploited areas, which generally includes geomorphological adaptation [4]. Open-pit mining can trigger erosive processes and landslides depending mainly on four principal factors (soil or substrate characteristics, vegetative cover, topography, and climate) [5]. Thus, to protect the environment, international organizations such as the United Nations (UN) and the Organization for Economic Co-operation and Development (OECD) have established guidelines [6,7]. At the same time, national policies (e.g., China [8], Canada [9], South Africa [10], Colombia [11], Australia [12], and Spain [13]) are becoming more and more committed to the responsible management of waste from extractive industries, and to protecting and rehabilitating the area affected by mining activities. These rehabilitating frameworks aim to prevent, or at least mitigate and rectify, land erosion as much as possible. Open-pit mining causes the destruction of soil vegetation, the formation of steep slopes, and the accumulation of waste and production materials. Consequently, these detrimental factors catalyze erosive processes that are further aggravated by meteorological conditions. Given the urgency of addressing and finding solutions to environmental impacts, decision makers are becoming increasingly interested in monitoring and modelling [14]. So far, soil erosion characterization [15] has been traditionally modeled with the Universal Soil Loss Equation (USLE) [16] and the Revised USLE (R-USLE) [17], by using topographic information, land cover information, Digital Elevation Models (DEMs), field observations, and environmental measurements applied by technicians [18,19]. These specialists rely on official cartography, satellite and/or airborne imagery, lidar data, and many atmospheric sensors to characterize the erosion evolution of a given quarry. Remote sensing is one of the most important data sources for land management, land use and land cover mapping, and change detection. Since the 1970s, airborne and satellite remote sensing has provided large amounts of spatial and thematic data for environmental monitoring. Since the mid- 2010s, the sensor miniaturization and commercialization of affordable Unmanned Aerial Systems (UASs) or drones have provided new possibilities of remote sensingdata capture, especially very detailed aerial images, for relief modelling and waterflow [20,21]. By using photogrammetric techniques, such as Structure from Motion (SfM), and Geographical Information System (GIS) tools, it is possible to generate topographic maps and three- dimensional models for monitoring mining production in stages [22], the calculation of stockpiles and dumps [23], as well as the analysis of the fragmentation of extracted material. Therefore, they are instrumental in planning and administrating future work in mines, both in terms of exploitation [24] and restoration [25]. A UAS includes an Unmanned Aerial Vehicle (UAV), a ground station, a communi- cations link, and often a sensor to capture data. In relatively small areas (around 10 km2 ), UASs are the best tools for obtaining information by remote sensing since they are inex- pensive, deploy quickly, and can fly at low altitudes, thereby limiting atmospheric effects such as clouds and fog. UAVs can carry different types of sensors depending on their mis- sion, including acoustic, visual, chemical, and biological sensors for environmental remote sensing purposes [26]. Compared to previous methods, the planned acquisition of optical imagery from a UAS joining the image processing with SfM and GIS results in various monitoring and volume calculations, which are essential for time and cost reduction. The use of UASs facilitates the obtainment of Digital Surface Models (DSMs), orthoimagery, land cover mapping, and the generation of three-dimensional models with centimetric resolution [27]. In addition to supporting the monitoring of mining production, risk maps can be obtained by processing and analyzing the images captured with a UAS. With these maps of areas that are difficult to access due to steep slopes and unstable ground, one can examine the rock geometry and analyze the topography in areas where explosives were used, among other extractive activities [28]. The monitoring of soil erosion in quarries with UASs has been evolving since 2015. In 2018 Padró et al. [29] presented a whole series of high-quality products derived from drone flights, which were intended to support ecological restoration plans in the 483 open-pit mining operations in Catalonia (Spain, Europe). With regional government assistance, a Land 2022, 11, 212 3 of 13 detailed workflow was proposed to support both the mining process and post-mine-closure activities using UASs [30,31]. The aim was to offer a set of techniques for generating spatial information, including vegetation indices, soil erosion, drainage networks, and digital models, among other measures supporting restoration policies [25]. The main objective of this work is to define in detail and demonstrate, with high- quality and high-resolution images, the applicability of a methodology for identifying areas eroded by surface runoff, focusing on the case where there is no reference DEM of before the erosion. The novelty of this work lies in the strategy of estimating the DEM before the erosion process (DEMT0 ), from which the current DEM (DEMT1 ) is subtracted to obtain a DEM of Differences (DoD) and the eroded volume of soil between time date 0 (T0) and date 1 (T1) is quantified. Apart from the method itself, interesting results could be derived from the slope map, whereby erosion rates can be estimated and erosion ranges can be determined according to slope categories. We focused on open-pit mine scenarios located in Catalonia, Spain. Based on the results obtained using this monitoring methodology, private companies and public administrations could be significantly supported in their decision-making concerning quarry operations and further restorations with the help of renewable and accurate mapping products. 2. Materials and Methods 2.1. Study Area The study area corresponded to the “La Constancia” quarry (Lat. 41◦ 330 8.254700 Lon. 1◦ 590 31.754800 ), which is part of the municipality of Terrassa (Catalonia, Spain, Europe) (Figure 1). The open-pit mine extracts inert aggregate materials, such as sand and gravel of different sizes, which are currently used for construction work. According to the geological map of the area, Miocene non-consolidated conglomerates with a silty-sandy matrix are the geological substrate which mainly constitutes those granular materials [32]. The study area has a Mediterranean sub-humid climate, where intense rainfalls during summer and autumn are frequent. According to the data from the closest meteorological station (Sabadell-Parc Agrari), the area was affected by two intense rainfall episodes (>50 mm/h) during October and December 2019 [33]. The drone-acquired data coversd a wide and heterogeneous area that includes build- ings, roads, and areas surrounding the front side, where the aggregates are extracted. The area which corresponds to the bottom of the mining basin has water deposits and was delimited and omitted to avoid errors in processing. Thus, we had a study area of approximately 36 ha, with a polygonal shape of 725 m × 550 m, wherein the land height varied between 288 m and 367 m above sea level. Specifically, there was a difference of near 80 m between the quarry bottom and its highest part (Figure 1). The most eroded part, which is the focus of the study and where more significant processes were undertaken, had a polygonal shape of 200 m × 150 m. 2.2. Unmanned Aerial Vehicle (UAV), Image Sensor and Flight Planning A photogrammetric survey was carried out on 5 May 2021, between 11:55 a.m. CET and 12:45 p.m. CET under excellent weather conditions. Georeferencing was carried out by locating 24 Ground Control Points uniformly distributed throughout the flight area and subsequently corroborated by photointerpretation. The Ground Control Point georefer- encing method has been tested and compared with other methods (e.g., Post-Processed Kinematic) and offers optimal results [34]. To cover the study area, 753 photographs were taken with a transverse and longitudinal overlap of 80%, at an average flight altitude of 115 m above the ground, covering an area of 353,000 m2 . The Unmanned Aerial Vehicle used was a DJI Inspire 2 [35], which has a magnesium-aluminum alloy shell and carbon fiber arms, can fly at speeds of up to 94 km/h with a thrust of 2 kg per rotor, 15” propellers, and a maximum flight time of 27 min. The sensor used was DJI Zenmuse X5S (FC6520) RGB of 20.8 MP (5280 × 3956 pixels). The sensor is stabilized by three axes of precise motor Land 2022, 11, 212 4 of 13 Land 2022, 11, x FOR PEER REVIEW 4 of 14 rotation with ± 0.01◦ control driven by a dedicated processor. Vibration is reduced by using three shock-absorbing balls. Figure 1. (a) Location of the study area. Source: Modified from Google Earth. (b) Regional location Figure 1. (a) Location of the study area. Source: Modified from Google Earth. (b) Regional location of “La Constancia” quarry near Barcelona, Spain, Europe. Source: Modified from Google Earth. (c) of “La Constancia” quarry near Barcelona, Spain, Europe. Source: Modified from Google Earth. Local location of “La Constancia” quarry in Terrassa, Catalonia. Source: Modified from Google (c) Local Earth. (d)location of “La Constancia” Orthoimagery quarry in quarry. of the “La Constancia” Terrassa, Catalonia. Source: Drone Source: Modified imagery from acquired Google during this Earth. processed study (d) Orthoimagery of the “La by the authors. (e)Constancia” quarry. Digital Elevation Source: Model Drone of the “Laimagery acquired Constancia” during quarry. this Source: study processed Drone by the authors. imagery acquired (e) Digital during this Elevation study was Model processed byof thethe “La Constancia” quarry. Source: authors. Drone imagery acquired during this study was processed by the authors. 2.2. Unmanned Aerial Vehicle (UAV), Image Sensor and Flight Planning 2.3. Methodology A photogrammetric survey was carried out on 5 May 2021, between 11:55 a.m. CET Widely used software was employed for the photogrammetric treatment of drone and 12:45 p.m. CET under excellent weather conditions. Georeferencing was carried out images (Agisoft Metashape PhotoScan 1.5.2 [36]) and GIS processing (ArcGIS 10.8 [37], by locating 24 Ground Control Points uniformly distributed throughout the flight area MiraMon 8.2e [38], QGIS 3.10 [39], and SAGA 7.9.0 [40]). The main novelty of this workflow and subsequently corroborated by photointerpretation. The Ground Control Point was the method of estimating the DEM in 2019, before the erosion process (DEMT0 ), from georeferencing method has been tested and compared with other methods (e.g., Post- the DEM after the erosion process (DEMT1 ), in 2021. Subtracting the current DEM (DEMT1 ), Processed Kinematic) and offers optimal results [34]. To cover the study area, 753 which is the result of erosion processes, from the estimated DEMT0 makes it possible to photographs were taken with a transverse and longitudinal overlap of 80%, at an average obtain a DoD and quantify the volume of eroded material. A Digital Surface Model (DSM) flight altitude with was obtained of 115 mphotogrammetric the above the ground,processing covering an areadrone-acquired of the of 353,000 m2. The Unmanned images. Thanks Aerial Vehicle used was a DJI Inspire 2 [35], which has a magnesium-aluminum to the overlapping of the images, a 3D point cloud (XYZ coordinates) with a very alloydense shell and carbon fiber arms, can fly at speeds of up to 94 km/h with a thrust of 2 kg sampling can be constructed. The point cloud was classified to keep only the terrain points, per rotor, 15” thuspropellers, and a maximum avoiding vegetation flight time and building of 27(Figure surfaces min. The 2). sensor used was DJI Zenmuse X5S (FC6520) RGB of 20.8 MP (5280 × 3956 pixels). Therefore, only ground points’ Z coordinate (altitude) The sensorwas is stabilized by three interpolated axes to model of theprecise terrainmotor surfacerotation with ±a 0.01° and obtain control grid map driven by containing thea dedicated processor. terrain altitude (DEMVibration T1 ). This is reduced by using three shock-absorbing balls. Land 2022, 11, x FOR PEER REVIEW 5 of 14 Land 2022, 11, 212 5 of 13 2.3. Methodology DEM was used to obtain the current drainage network by applying GIS tools based on Widely used software was employed for the photogrammetric treatment of drone the aspect and slope models of DEMT1 , which defines the water flow direction, water images (Agisoft Metashape PhotoScan 1.5.2 [36]) and GIS processing (ArcGIS 10.8 [37], flow accumulation, the channel flow network, and the watershed ridge divisions. The MiraMon 8.2e [38], QGIS 3.10 [39], and SAGA 7.9.0 [40]). The main novelty of this estimated DEM before the erosion process (DEMT0 ) was obtained from the interpolation workflow was the method of estimating the DEM in 2019, before the erosion process of the Z coordinate (altitude) of the 3D point cloud points located on the ridges of the (DEMT0), from the DEM after the erosion process (DEMT1), in 2021. Subtracting the current DEMT1 . By selecting the points close to the watershed divisions and interpolating them, DEM (DEMT1), which is the result of erosion processes, from the estimated DEM T0 makes the terrain surface was modelled and rasterized to obtain a primary grid map with the it possible to obtain a DoD and quantify the volume of eroded material. A Digital Surface terrain altitude (DEMT0 ), being a reasonable approximation of the original relief before Model (DSM) was obtained with the photogrammetric processing of the drone-acquired the erosion processes considering that there was no previous direct source of data. If the images. Thanks to the overlapping of the images, a 3D point cloud (XYZ coordinates) with original 3D point cloud was not available, a variant could be reached if the DEMT1 cells a very dense sampling can be constructed. The point cloud was classified to keep only the intersecting the ridge vector were converted to points, for cloud comparison, interpolation, terrain and points, thus avoiding vegetation and building surfaces (Figure 2). rasterization. Figure 2. Figure 2. (a) (a) Camera coverage and Camera coverage and Ground Ground Control Control Point Point (GCP) (GCP) locations. locations. (b) (b) 2D 2D nadiral nadiral perspective perspective of RGB colored point cloud. (c) 2D nadiral perspective of point cloud classification (ground points of RGB colored point cloud. (c) 2D nadiral perspective of point cloud classification (ground points were used for the DEM construction). (d) 3D perspective of RGB colored point cloud. (e) 3D were used for the DEM construction). (d) 3D perspective of RGB colored point cloud. (e) 3D perspective of point cloud classification. Source: Drone imagery acquired during this study perspective of point cloud classification. Source: Drone imagery acquired during this study processed processed with Metashape Photoscan by the authors. with Metashape Photoscan by the authors. Therefore, only ground points’ Z coordinate (altitude) was interpolated to model the In this processing, special care was given to ensure that the cell-size value of the DEM terrain surface and obtain a grid map containing the terrain altitude (DEM T1). This DEM before erosion was equal to the cell-size value of the current DEM, so as to avoid errors in was used to obtain the current drainage network by applying GIS tools based on the the subsequent comparison of Digital Elevation Models. Finally, an arithmetic subtraction aspect andwas operation slope models executed of DEM (DEM T1, which defines the water flow direction, water flow T0 –DEMT1 ) to obtain a difference raster that corresponds to accumulation, the the channel volume of material eroded.flow network, and the watershed ridge divisions. The estimated Taking into accountthe DEM before thaterosion process (DEM slope inclination ) was obtained andT0length are factors from the interpolation affecting the erosion of the Z coordinate (altitude) of the 3D point cloud points located process, the slope calculation was carried out using the classification criteria on the ridges of of the Lugo DEM T1. By selecting the points close to the watershed divisions and interpolating them, Hubp [41], which classifies the slopes in five categories given the well-known correlation the terrain between surface slope was modelled and erosion: and◦ ),rasterized planar (0–3 to obtain slightly inclined a primary (3–12 grid ◦ ), inclined map◦ ),with (12–30 the sliding terrain◦ ), (30–45 altitude (DEM and free T0), being fall (>45 ◦ ). a reasonable approximation of the original relief before the erosion processes considering The drone imagery acquisition thatwas there fullywas no previous planned, directand automatized, source of data. replicable If the by using original 3D point cloud was not available, a variant could be reached mission planner software. The photogrammetric and SfM processing was automatized if the DEM T1 cells intersecting the ridge vector were converted to points, for using batch processes. The GIS process to extract hydrologic and morphometric informationcloud comparison, interpolation, from a DEM was and rasterization. fully automatized with the graphical model builder tools of ArcGIS 10.8 or QGISIn this(Figure 3.10 processing, 3). special care was given to ensure that the cell-size value of the DEM before erosion was equal to the cell-size value of the current DEM, so as to avoid errors in Hubp [41], which classifies the slopes in five categories given the well-known correlation between slope and erosion: planar (0°−3°), slightly inclined (3°–12°), inclined (12°–30°), sliding (30°–45°), and free fall (>45°). The drone imagery acquisition was fully planned, automatized, and replicable by using mission planner software. The photogrammetric and SfM processing was Land 2022, 11, 212 automatized using batch processes. The GIS process to extract hydrologic and 6 of 13 morphometric information from a DEM was fully automatized with the graphical model builder tools of ArcGIS 10.8 or QGIS 3.10 (Figure 3). Figure 3. Complete workflow from the data caption to the eroded volume estimation. Figure 3. Complete workflow from the data caption to the eroded volume estimation. 3. Results 3. Results From From thethe photogrammetric processing, photogrammetric processing,ananorthophotomosaic orthophotomosaic of 2.4 of cm/pixel nativenative 2.4 cm/pixel resolution was obtained, which was resampled to 2.5 cm/pixel to facilitate rasterraster resolution was obtained, which was resampled to 2.5 cm/pixel to facilitate processing processing and conform with the Catalan official mapping agency’s cartography grid and conform with the Catalan official mapping agency’s cartography grid (multiple of (multiple of 2.5 cm). A point cloud with a density of 122 points/m2 was generated, which 2.5 cm). A point cloud with a density of 122 points/m2 was generated, which allowed us to allowed us to elaborate a DEM with a resolution of 10 cm/pixel using only the points elaborate a DEM classified with aInresolution as ground. the AOI theofvegetation 10 cm/pixel using is coverage only the points negligible, classified as can as ground. be observed In the AOI the vegetation coverage is negligible, as can be observed in the in the orthoimage, thus avoiding possible errors related to 3D point cloud classification orthoimage, thus avoiding errors. At thispossible errors related spatial resolution, to 3D point cloudofclassification the photointerpretation errors. the relief at very highAt thisisspatial detail Land 2022, 11, x FOR PEER REVIEW resolution, 7 of 14 possible, the photointerpretation particularly when combining of the relief 2D and 3Datmodels very high detail (Figure 4). is possible, particularly when combining 2D and 3D models (Figure 4). Figure 4. (a) Orthoimage resulting from the UAS flight and delimitation of the Area of Interest 4. (a)(b) Figure(AOI). Orthoimage resulting Digital Elevation Modelfrom the resulting (DEM) UAS flight fromand thedelimitation UAS flight. (c)ofDetail the Area of 3D of the Interest (AOI). (b) Digital Elevation representation Model of the (DEM) orthoimage in resulting the AOI. (d)from the UAS flight. 3D representation of the(c) DEMDetail of the and the 3D representation orthoimage of theinorthoimage the AOI, where in the theeroded AOI. terrain (d) 3Dcan be clearly seen. Source: representation of the Drone DEM imagery and theacquired during in the AOI, orthoimage this study processed with QGIS by the authors. where the eroded terrain can be clearly seen. Source: Drone imagery acquired during this study processed with QGIS by 3.1. Morphometric the authors. Analysis and Drainage Network To carry out this analysis, the 10 cm resolution DEM T1 was used as the main input. The slope and aspect models were calculated for initial morphometric analysis. Firstly, the sinks were corrected by using an algorithm [42] to modify the elevation values in those places of the DEM where water could accumulate. In this correction, the drainage network flowed without interruption through the model. Secondly, the drainage network derived Land 2022, 11, 212 7 of 13 3.1. Morphometric Analysis and Drainage Network To carry out this analysis, the 10 cm resolution DEM T1 was used as the main input. The slope and aspect models were calculated for initial morphometric analysis. Firstly, the sinks were corrected by using an algorithm [42] to modify the elevation values in those places of the DEM where water could accumulate. In this correction, the drainage network flowed without interruption through the model. Secondly, the drainage network derived from the DEMT1 was characterized by the natural course of water and sediment transport by gravity on the ground, using the catchment area algorithm [43]. Subsequently, the flow direction was calculated to determine the direction of the water flow from one pixel to another, considering each pixel as eight adjacent cells (D8), following the direction where there was a more pronounced slope from the cell being analyzed. Cells with concentrated flow accumulation were used to identify the channels of the drainage network [44] in order Land 2022, 11, x FOR PEER REVIEW 8 of 14 to determine the amount of water that flowed through each cell from the adjacent cells and allow the extraction of the sectors of the land where water accumulates (Figure 5). Figure 5. (a) Flow water catchment (Log10 for better visualization) of the Area of Interest (AOI) Figure 5. (a) Flow water catchment (Log10 for better visualization) of the Area of Interest (AOI) based based on DEMT1. (b) Zoom of (a), with the channels and the XYZ points located over the watershed, on DEMwere which T1 . (b) Zoom used of (a),DEM to build withT0the . (c)channels and the XYZ 3D representation of points DEMT1located over with the the watershed, channel network. which were used to build DEMT0 . (c) 3D representation of DEMT1 with the channel network. Once detected, the channels, together with the DEMT1, were used to delimit the watershed network. The selection of the XYZ points of the dense point cloud, located over the watershed lines, was the key to obtaining the reference altitudes on the ridges, being the T0 estimation. Interpolating the Z coordinate at the ridges, it was possible to approximate a DEMT0 (Figure 6). Land 2022, 11, 212 8 of 13 Once detected, the channels, together with the DEMT1 , were used to delimit the watershed network. The selection of the XYZ points of the dense point cloud, located over the watershed lines, was the key to obtaining the reference altitudes on the ridges, being the Land 2022, 11, x FOR PEER REVIEW 9 of 14 T0 estimation. Interpolating the Z coordinate at the ridges, it was possible to approximate a DEMT0 (Figure 6). Figure6.6. (a) Figure (a)DEM DEMT0 of the T0 of the Area Area ofof Interest Interest (AOI), (AOI),estimating estimating the theterrain terrainmorphology morphology before beforethe the erosionprocess. erosion process.(b) (b)DEM DEMT1 of of the the AOI AOI based basedon ondrone-acquired drone-acquiredimagery. imagery.(c) (c)3D 3Drepresentation representationofof T1 DEMT0. (d) 3D representation of DEMT1. DEMT0 . (d) 3D representation of DEMT1 . 3.2.Erosion 3.2. ErosionEstimation Estimation TheDEM The DEMof of Differences Differences between between two two DEMsDEMs is a technique is a technique used to geomor- used to estimate estimate geomorphological phological changes changes and quantifyand quantify volumesvolumes of material, of material, mainly mainly in land in land exposed exposed to water to erosion [45]. Erosion water erosion in the area [45]. Erosion occurs in the areamainly occursin the form mainly of gullies in the form ofbecause gulliesthe material because the of the slopes material is slopes of the not compacted and the and is not compacted topographic slope inslope the topographic someinareas someisareas very issteep. very The calculation steep. of erodedofmaterial The calculation eroded was obtained material was by means of obtained byanmeans arithmetic of ansubtraction arithmetic operation (DEM subtraction T0 –DEM operation T1 ) (Figure (DEM T0 – DEM7).T1) (Figure 7). Difference raster values of 0 referred to areas with no change in volume, the positive values corresponded to land-eroded areas, and the negative values corresponded to land accumulation areas. Likewise, the slope raster was reclassified and clipped to the same scope of the difference raster to obtain adequate results by conducting the zonal statistical analysis, which generated tables of slope-category data sets by the area and volume of the eroded material. In Figure 4b, the browned areas are those without evidence of soil loss due to erosion, followed by areas with orange and yellow colors. On the other side, magenta areas are those with greater land accumulation, followed by areas with blue and green col- ors. As can be seen in the previous image, the effects caused by erosion and the subsequent material accumulation at lower slopes were observed and could be quantified. Laminar erosion processes, mainly the formation of gullies of different depths, are represented as variations in coloration. Through the gullies the material was dragged in granulated form towards the bottom of the quarry. In total, a removal amount of approximately 15,414.9 m3 was estimated, corresponding to an erosion rate of 364 Mg/ha/year, considering a substrate density of 1.7 Mg/m3 [46]. 3.2. Erosion Estimation The DEM of Differences between two DEMs is a technique used to estimate geomorphological changes and quantify volumes of material, mainly in land exposed to water erosion [45]. Erosion in the area occurs mainly in the form of gullies because the Land 2022, 11, 212 material of the slopes is not compacted and the topographic slope in some areas is9very of 13 steep. The calculation of eroded material was obtained by means of an arithmetic subtraction operation (DEMT0 – DEMT1) (Figure 7). Land 2022, 11, x FOR PEER REVIEW 10 of 14 Difference raster values of 0 referred to areas with no change in volume, the positive values corresponded to land-eroded areas, and the negative values corresponded to land accumulation areas. Likewise, the slope raster was reclassified and clipped to the same scope of the difference raster to obtain adequate results by conducting the zonal statistical analysis, which generated tables of slope-category data sets by the area and volume of the eroded material. In Figure 4b, the browned areas are those without evidence of soil loss due to erosion, followed by areas with orange and yellow colors. On the other side, magenta areas are those with greater land accumulation, followed by areas with blue and green colors. As can be seen in the previous image, the effects caused by erosion and the subsequent material accumulation at lower slopes were observed and could be quantified. Laminar erosion processes, mainly the formation of gullies of different depths, are represented as variations in coloration. Through the gullies the material was dragged in granulated form towards the bottom of the quarry. In total, a removal amount of approximately Figure 7. (a) 15,414.9 (a) Erosion Erosion m3 wasfrom map resulting resulting estimated, corresponding thedifference difference betweenDEMDEMtoT0 an anderosion rate of DEMT1 (DoD) (DoD) 364 at the the Figure 7. map from the between 3 [46]. T0 and DEM T1 at Mg/ha/year, Area of Interestconsidering a substrate (AOI). (b) Detail of (a). (c)density of 1.7 Mg/m 3D representation of the erosion map at the AOI. Area of Interest (AOI). (b) Detail of (a). (c) 3D representation of the erosion map at the AOI. 3.3.Slopes 3.3. Slopesand andErosion ErosionIdentification Identification Slopemaps Slope mapsare areuseful usefulininhydrological hydrologicalstudies studiesandanderosion erosionmodels. models.To Toproduce produceslope slope maps, it is necessary to calculate the slope of the land with respect to maps, it is necessary to calculate the slope of the land with respect to a horizontal plane. a horizontal plane. InInthis thiscase casestudy, study,thetheslope slopecalculation calculationwas wasmade madebased basedon onthe thefollowing followingcategories categories[41]: [41]: planar(0–3 planar ◦ (0°–3°), slightlyinclined ), slightly inclined(3–12 ◦ (3°–12°), inclined(12–30 ), inclined ◦ (12°–30°), sliding(30–45 ), sliding ◦ (30°–45°), andfree ), and freefall fall (>45 ◦ ). In (>45°). Interms termsof ofarea, area,half halfofofthe theland land(49%) (49%)presents presentsslopes slopesgreater greaterthan ◦ , which than3030°, whichinin theorywould theory wouldfacilitate facilitatethe thedrag dragphenomenon phenomenonmainlymainlywhenwhenthere thereisisdrizzle, drizzle,while whileonly only3% 3% ofofthe thearea areaisispart partofofthe theland landcategorized categorizedasasplain plain(Figure (Figure8).8). Figure8.8.DEM Figure DEMT1slopes slopesanalysis analysisatatthe theArea AreaofofInterest Interest(AOI) (AOI)(water (waterisisnot notcomputed). computed). T1 Toanalyze To analyzethe therelationship relationshipbetween betweenthetheeroded erodedarea areaand andthe theslope slopeofofthe theland, land,thethe erodedarea eroded areacalculation calculationwas wasperformed performedfor foreach eachslope slopecategory categoryatatpixel pixelscale scale(Table (Table1).1). Table 1. Area eroded by range of slopes. Range (degrees) Area (m2) Volume (m3) Area of Slope Range (%) 0–3 191.1 357.7 3 Land 2022, 11, 212 10 of 13 Table 1. Area eroded by range of slopes. Range (Degrees) Area (m2 ) Volume (m3 ) Area of Slope Range (%) 0–3 191.1 357.7 3 3–12 1117.6 2462.3 19 12–30 1569.1 3923.0 29 30–45 1903.1 6495.5 40 >45 640.6 2176.2 9 Total 18,224 15,414.7 100 4. Discussion Water erosion occurs when it rains and the infiltration capacity of the land is saturated, or when the intensity of rain exceeds the rate of soil permeability, and subsequently the water runs superficially (runoff) and drags soil particles, which are mostly removed by the direct impact of raindrops [47]. Erosion in mining activities can cause structural risks such as landslides, and therefore it is important to understand erosion behavior. The analysis carried out accurately identifies the erosion of the slopes and calculates the volume of displaced material [48]. The use of UASs for aerial image acquisition represents an affordable and accurate source of data for mining environments, as is well documented in recent research stud- ies [24,25,27]. The photogrammetric and SfM processing allowed the generation of a 3D dense point cloud, orthoimage, and DEMT1 at a very high sampling distance for the entire area of interest, which is consistent with recent studies [22,28]. The methodology used for estimating the pixels corresponding to channels and watersheds, according to a method previously tested and implemented in morphometric [43] and hydrologic [44] analysis, was applied in this research with success in drone imagery, and enabled us to construct an esti- mated DEM (DEMT0 ) using the points of the dense point cloud located over the watersheds of DEMT1 . The estimated erosion was calculated by applying DoD technique [45]. The raster product of this processing was obtained with GIS software, which is widely used for morphometric and water balance analysis, as the analytical capabilities of the algorithms implemented on Digital Elevation Models and imagery can yield results in very high detail. To improve the process of detecting the ridges or bottoms of gullies, it was very important to take into account the meteorological conditions when carrying out flights using UASs. It was necessary to pay attention to these factors in order to minimize the presence of shadows that could cause errors; capturing images during cloudy days reduces the shadow presence, when diffuse solar irradiance is dominant over direct sun irradiance [29]. Shadow reduction techniques can also be applied [49]. Additionally, it is important to determine the geomorphology under water bodies to avoid underestimating erosion records because it is the lowest area of the land where the entrained material ends up accumulating. The method implemented for the quantitative estimation of erosion may vary accord- ing to some aspects at the researcher’s discretion, such as the choice of pixels considered to represent different levels of erosion in the gullies or when the type of interpolation is chosen (different methods generate different results). Another aspect to consider would be how to determine the submerged terrain, where the dragged material ends up accumulating. In calculating volumes and areas, results may vary depending on whether only the DEM is used, or if the vector file is included to identify the area volume for each range of slopes. In the latter case, the calculation algorithms are different, since statistical subsets are used and this causes differences in results. It is also important to consider the difference in the geometry of the pixels (square shape), and the polygons of the vector file (mainly triangular irregular shapes). The calculation of the erosion rate was carried out considering events of great activity on bare soil and with little resistance to weathering, which would explain the high values obtained. Land 2022, 11, 212 11 of 13 In general terms, we consider that the implemented methodology is robust and easily replicable in other scenarios with similar geomorphological conditions, either in natural areas or in the framework of infrastructure works where water erosion can be a factor for evaluation. Additionally, the methodology would facilitate the risk analysis of structural problems caused by erosion. 5. Conclusions As demonstrated in this study as well as previous literature, remote sensing using drone-acquired imagery is a very useful tool for the morphometric characterization of terrains in open-pit mining activities, especially for erosion processes related to restored areas. However, drone data for the initial moment of the erosion process are not often available, and the DEM before the erosion process (DEMT0 ) must therefore be estimated from the current DEM (DEMT1 ). Here we described a new strategy for estimating the DEM before the erosion process by using well-known GIS techniques based on hydrologic analysis (channel network and watershed detection). From the estimated DEMT0 , the current DEM T1 was subtracted to obtain a DEM of Differences (DoD) and identify the estimated eroded area and volume of soil between the date 0 (T0) and date 1 (T1). In this case study, an erosion detection and calculation workflow was applied, providing valuable information for quarry restoration management, particularly regarding those areas with higher erosion rates and subsequent land accumulations. Controlling the slope and calculating erosion rates enable the monitoring of areas with active erosion processes and facilitate decisionmaking that is crucial for implementing stabilization measures in areas of high landslide susceptibility. Therefore, the presented method is proven to be valid and highly suitable for mining operations and applicable to other geomorphologically similar scenarios. Author Contributions: Conceptualization, V.C., P.M. and J.-C.P.; data acquisition, R.R.-C. and D.S.; methodology, P.M. and J.-C.P.; software, J.C., P.M. and J.-C.P.; validation, V.C. and J.C.; formal analysis, J.C.; investigation, V.C.; resources, V.C. and J.M.A.; data curation, V.C., P.M., R.R.-C. and J.M.A.; writing—original draft preparation, J.C. and J.-C.P.; writing—review and editing, R.R.-C., J.C., J.-C.P., V.C. and J.M.A.; visualization, P.M., J.C. and J.-C.P.; supervision, V.C.; project administration, V.C.; funding acquisition, V.C. 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” and by the European Union through the LIFE project “New Life for Drylands” (LIFE20 PRE/IT/000007). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: The data that support the findings of this study are available from the corresponding author upon reasonable request. Conflicts of Interest: The authors declare no conflict of interest. References 1. European Commission. Internal Market, Industry, Entrepreneurship and SMEs. Available online: https://ec.europa.eu/growth/ sectors/raw-materials/eip/strategic-implementation-plan_en (accessed on 19 January 2022). 2. Altiti, A.H.; Alrawashdeh, R.O.; Alnawafleh, H.M. Open Pit Mining, Mining Techniques—Past, Present and Future, 1st ed.; Abhay, S., Ed.; IntechOpen: London, UK, 2021; Available online: https://www.intechopen.com/chapters/71931 (accessed on 19 January 2022). 3. United Nations, and Deutsche Stiftung für Internationale Entwicklung. Mining and the Environment: The Berlin Guidelines: A Study Based on an International Round Table in June 1991, Organized by the Department of Technical Co-Operation for Development, United Nations, and the Development Policy Forum of the German Foundation for International Development; Mining Journal Books: London, UK, 1992. 4. Narendra, B.H.; Siregar, C.A.; Turjaman, M.; Hidayat, A.; Rachmat, H.H.; Mulyanto, B.; Maharani, R. Managing and Reforesting Degraded Post-Mining Landscape in Indonesia: A Review. Land 2021, 10, 658. [CrossRef] Land 2022, 11, 212 12 of 13 5. Segura-Salazar, J.; Tavares, L.M. Sustainability in the Minerals Industry: Seeking a Consensus on Its Meaning. Sustainability 2018, 10, 1429. [CrossRef] 6. United Nations Environment Programme. Sustainability Reporting in the Mining Sector. Current Status and Future Trends. 2020. Available online: https://www.unep.org/resources/report/sustainability-reporting-mining-sector (accessed on 19 January 2022). 7. Organization for Economic Co-operation and Development. OECD Brochure: A Global Standard: Towards Responsible Mineral Supply Chains. 2020. Available online: https://mneguidelines.oecd.org/Brochure_OECD-Responsible-Mineral-Supply-Chains. pdf (accessed on 19 January 2022). 8. Government of China. The 13th Five Year Plan for Economic and Social Development of the People’s Republic of China (2016–2020). 2016. Available online: https://en.ndrc.gov.cn/policies/202105/P020210527785800103339.pdf (accessed on 19 January 2022). 9. Canadian Government. The Canadian Minerals and Metals Plan. 2020. Available online: https://www.minescanada.ca/sites/ default/files/cmmp-actionplan2020_rev52_feb_29_2020-a_en.pdf (accessed on 19 January 2022). 10. Minerals Council South Africa. Modernization: Towards the Mine of Tomorrow. 2020. Available online: https://www.mine2030. co.za/ (accessed on 19 January 2022). 11. Government of Colombia. Decreto 1076 de 2015 Sector Ambiente y Desarrollo Sostenible. 2015. Available online: https://www. funcionpublica.gov.co/eva/gestornormativo/norma.php?i=78153 (accessed on 19 January 2022). 12. Australian Government. Mine Rehabilitation. Leading Practice Sustainable Development Program for the Mining Industry. 2016. Available online: https://www.industry.gov.au/sites/default/files/2019-04/lpsdp-mine-rehabilitation-handbook-english.pdf (accessed on 19 January 2022). 13. Government of Spain. Real Decreto-ley 975/2009, del 12 de Junio, Sobre Gestión de los Residuos de las Industrias Extractivas y de Protección y Rehabilitación del Espacio Afectado por Actividades Mineras. 2009. Available online: https://www.boe.es/ buscar/act.php?id=BOE-A-2009-9841 (accessed on 19 January 2022). 14. Huang, Y.; Espley, S. A 3D mine simulation model for decision-making in mine design and production. International Journal of Surface Mining. Reclam. Environ. 2005, 19, 251–259. [CrossRef] 15. Ramli, M.; Thamrin, M.; Asrafil, M. Analysis of Soil Erosion on Mine Area. IOP Conf. Ser. Mater. Sci. Eng. 2020, 875, 012052. [CrossRef] 16. Wischmeier, W.H.; Smith, D.D. Predicting Rainfall Erosion Losses: A Guide to Conservation Planning with Universal Soil Loss Equation (USLE); Agriculture Handbook; Department of Agriculture: Washington, DC, USA, 1978. 17. Renard, K.G.; Foster, G.R.; Weesies, G.A.; Porter, J.P. RUSLE: Revised universal soil loss equation. J. Soil Water Conserv. 1991, 46, 30–33. 18. Duarte, L.; Cunha, M.; Teodoro, A.C. Comparing Hydric Erosion Soil Loss Models in Rainy Mountainous and Dry Flat Regions in Portugal. Land 2021, 10, 554. [CrossRef] 19. Barrena-González, J.; Rodrigo-Comino, J.; Gyasi-Agyei, Y.; Pulido Fernández, M.; Cerdà, A. Applying the RUSLE and ISUM in the Tierra de Barros Vineyards (Extremadura, Spain) to Estimate Soil Mobilisation Rates. Land 2020, 9, 93. [CrossRef] 20. 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] 21. 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] 22. Yaqiu, Y.; Cunhao, J.; Jing, L.; Jie, W.; Xing, J.; Hao, W.; Yu, X.; Lushou, Z. Mining Ground Surface Information Extraction and Topographic Analysis Using UAV Video Data. In Proceedings of the E3S Web of Conferences, Shanghai, China, 18–20 September 2020; Volume 194, p. 05030. [CrossRef] 23. Andras, M.; Domozi, Z. Tracking Production Volumes of Open-Pit Mines with Photogrammetry. In Proceedings of the 2017 European Modelling Symposium (EMS), Manchester, UK, 20–21 November 2017; pp. 100–106. [CrossRef] 24. Honarmand, M.; Shahriari, H. Geological Mapping Using Drone-Based Photogrammetry: An Application for Exploration of Vein-Type Cu Mineralization. Minerals 2021, 11, 585. [CrossRef] 25. Carabassa, V.; Montero, P.; Alcañiz, J.M.; Padró, J.-C. Soil Erosion Monitoring in Quarry Restoration Using Drones. Minerals 2021, 11, 949. [CrossRef] 26. Manfreda, S.; McCabe, M.; Miller, P.; Lucas, R.; Pajuelo Madrigal, V.; Mallinis, G.; Ben-Dor, E.; Helman, D.; Estes, L.; Ciraolo, G.; et al. On the Use of Unmanned Aerial Systems for Environmental Monitoring. Remote Sens. 2018, 10, 641. [CrossRef] 27. Shahmoradi, J.; Talebi, E.; Roghanchi, P.; Hassanalian, M. A Comprehensive Review of Applications of Drone Technology in the Mining Industry. Drones 2020, 4, 34. [CrossRef] 28. Ren, H.; Zhao, Y.; Xiao, W.; Hu, Z. A review of UAV monitoring in mining areas: Current status and future perspectives. Int. J. Coal Sci. Technol. 2019, 6, 320–333. [CrossRef] 29. Padró, J.C.; Carabassa, V.; Balagué, J.; Brotons, L.; Alcañiz, J.; Pons, X. Monitoring opencast mine restorations using Unmanned Aerial System (UAS) imagery. Sci. Total Environ. 2018, 657, 1602–1614. [CrossRef] 30. Carabassa, V.; Montero, P.; Crespo, M.; Padró, J.C.; Balagué, J.; Alcañiz, J.M.; Brotons, L.; Pons, X. UAS Remote Sensing Products for Supporting Extraction Management and Restoration Monitoring in Open-Pit Mines. Proceedings 2019, 30, 4. [CrossRef] Land 2022, 11, 212 13 of 13 31. Carabassa, V.; Montero, P.; Crespo, M.; Padró, J.C.; Pons, X.; Balagué, J.; Brotons, L.; Alcañiz, J. Unmanned aerial system protocol for quarry restoration and mineral extraction monitoring. J. Environ. Manag. 2020, 270, 110717. [CrossRef] 32. ICGC. Mapa Geològic de Catalunya. Institut Cartogràfic i Geològic de Catalunya. 2021. Available online: http://betaportal.icgc. cat/visor/client_utfgrid_geo.html (accessed on 19 January 2022). 33. SMC. Xarxa Agrometeorològica de Catalunya (Estació Sabadell-Parc Agrari). Servei Meteorològic de Catalunya. 2022. Available online: https://ruralcat.gencat.cat/agrometeo.estacions (accessed on 19 January 2022). 34. 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. Int. J. Appl. Earth Obs. Geoinf. 2019, 79, 130–140. [CrossRef] 35. DJI. Inspire 2 User Manual v1.0 2016.12. 2016. Available online: https://dl.djicdn.com/downloads/inspire_2/INSPIRE+2+User+ Manual+.pdf (accessed on 19 January 2022). 36. Agisoft LLC. Agisoft Metashape User Manual: Professional Edition, Version 1.5. 2019. Available online: https://www.agisoft. com/pdf/metashape-pro_1_5_en.pdf (accessed on 19 January 2022). 37. ESRI. ArcGIS Desktop: Release 10.8; Environmental Systems Research Institute: Redlands, CA, USA, 2016. 38. Pons, X. MiraMon. Sistema d’Informació Geogràfica i Software de Teledetecció. Versió 8.1f [MiraMon. Geographical Information System and Remote Sensing Software. Version 8.2e]. Centre de Recerca Ecològica i Aplicacions Forestals, CREAF. Bellaterra. 2021. Available online: http://www.creaf.uab.cat/miramon/Index_usa.htm (accessed on 19 January 2022). 39. QGIS.org. QGIS Geographic Information System. QGIS Association. 2021. Available online: http://www.qgis.org (accessed on 19 January 2022). 40. Conrad, O.; Bechtel, B.; Bock, M.; Dietrich, H.; Fischer, E.; Gerlitz, L.; Wehberg, J.; Wichmann, V.; Böhner, J. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. Geosci. Model Dev. 2015, 8, 1991–2007. [CrossRef] 41. Lugo-Hubp, J.I. Elementos de Geomorfología Aplicada (Métodos Cartográficos): Universidad Nacional Autónoma de México; Instituto de Geografía: México City, Mexico, 1988; 128p. 42. Wang, L.; Liu, H. An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. Int. J. Geogr. Inf. Sci. 2006, 20, 193–213. [CrossRef] 43. Wood, J. Geomorphometry in LandSerf. In Geomorphometry: Concepts, Software, Applications. Developments in Soil Science; Hengl, T., Reuter, H.I., Eds.; Elsevier: Amsterdam, The Netherlands, 2009; Volume 33, pp. 333–349. 44. Fairfield, J.; Leymarie, P. Drainage Networks from Grid Digital Elevation Models. Water Resour. Res. 1991, 27, 709–717. [CrossRef] 45. Williams, R.D. DEMs of Difference. In Geomorphological Techniques (Online Edition); Cook, S.J., Clarke, L.E., Nield, J.M., Eds.; British Society for Geomorphology: London, UK, 2012; pp. 1–17. 46. Chesworth, W. (Ed.) Encyclopedia of Soil Science; Springer: Berlin/Heidelberg, Germany, 2008; ISBN 978-1-4020-3994-2. 47. Farek, V.; Unucka, J.; Ponfzilova, I.; Gergerova, M.; Židek, D.; Pallos, R. Assessment of the runoff conditions of small ungauged catchments using GIS and fully distributed hydrologic models. Acta Montan. Slovaca 2014, 19, 22–30. 48. Krusekopf, H.H. The effect of slope on soil erosion. In University of Missouri Agricultural Experimental Station Research Bulletin; University Archives of the University of Missouri-Columbia: Columbia, MI, USA, 1943; p. 363. 49. Pons, X.; Padró, J.-C. An Operational Radiometric Correction Technique for Shadow Reduction in Multispectral UAV Imagery. Remote Sens. 2021, 13, 3808. [CrossRef]