Drone-Based Identiﬁcation of Erosive Processes in Open-Pit Mining Restored Areas

: 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 identifying 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 identiﬁed. Subsequently, an estimated DEM before the erosive processes was reconstructed by interpolating the gully ridges; this DEM serves as a reference for the relief before the erosion. To calculate the volume of eroded material, the DEM of Differences was calculated, which estimates the volume difference between the previously estimated DEM and the current DEM. Additionally, we calculated the material necessary for the geomorphological adaptation of the quarry and the slope map, which are two valuable factors 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-ﬁlled quarries, it would be important to characterize the subsurface relief. Essentially, the presented method can be applied with affordable and non-invasive materials to create digital grid maps at 10 cm resolution, obtaining data ready for 3D metrics, being a very practical landscape modelling tool for characterizing the restoration evolution of open-pit mining spaces.


Introduction
Extracting raw materials such as sand, gravel, crushed stone, brick clays, gypsum, and ornamental stones requires intensive and extensive efforts as all of these products are obtained from open-pit mining operations [1]. The open-pit mining method is used where large volumes of materials are mobilized, which gives rise to the formation of quarries [2]. Open-pit mining is one of the most land-degrading and landscape-modifing activities. Public awareness of sustainable land management has been growing, and since the 1990s, environmental policies and restoration actions have been increasingly focused on addressing environmental impacts [3]. In this type of mining, exploitation occurs in phases 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 highquality 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 (DEM T0 ), from which the current DEM (DEM T1 ) 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.

Study Area
The study area corresponded to the "La Constancia" quarry (Lat. 41 • 33 8.2547 Lon. 1 • 59 31.7548 ), 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 buildings, 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.

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 georeferencing 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 m 2 . 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 Drone imagery acquired during this study was processed by the authors.

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 georeferencing 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 m 2 . 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 rotation with ± 0.01° control driven by a dedicated processor. Vibration is reduced by using three shock-absorbing balls.

Methodology
Widely used software was employed for the photogrammetric treatment of drone images (Agisoft Metashape PhotoScan 1.5.2 [36]) and GIS processing (ArcGIS 10.8 [37], MiraMon 8.2e [38], QGIS 3.10 [39], and SAGA 7.9.0 [40]). The main novelty of this workflow was the method of estimating the DEM in 2019, before the erosion process (DEM T0 ), from the DEM after the erosion process (DEM T1 ), in 2021. Subtracting the current DEM (DEM T1 ), which is the result of erosion processes, from the estimated DEM T0 makes it possible to obtain a DoD and quantify the volume of eroded material. A Digital Surface Model (DSM) was obtained with the photogrammetric processing of the drone-acquired images. Thanks to the overlapping of the images, a 3D point cloud (XYZ coordinates) with a very dense sampling can be constructed. The point cloud was classified to keep only the terrain points, thus avoiding vegetation and building surfaces ( Figure 2). Therefore, only ground points' Z coordinate (altitude) was interpolated to model the terrain surface and obtain a grid map containing the terrain altitude (DEM T1 ). This DEM was used to obtain the current drainage network by applying GIS tools based on the aspect and slope models of DEM T1 , which defines the water flow direction, water flow accumulation, the channel flow network, and the watershed ridge divisions. The estimated DEM before the erosion process (DEM T0 ) was obtained from the interpolation of the Z coordinate (altitude) of the 3D point cloud points located on the ridges of the DEM T1 . By selecting the points close to the watershed divisions and interpolating them, the terrain surface was modelled and rasterized to obtain a primary grid map with the terrain altitude (DEM T0 ), being a reasonable approximation of the original relief before the erosion processes considering that there was no previous direct source of data. If the original 3D point cloud was not available, a variant could be reached if the DEM T1 cells intersecting the ridge vector were converted to points, for cloud comparison, interpolation, and rasterization.

Methodology
Widely used software was employed for the photogrammetric treatment of drone images (Agisoft Metashape PhotoScan 1.5.2 [36]) and GIS processing (ArcGIS 10.8 [37], MiraMon 8.2e [38], QGIS 3.10 [39], and SAGA 7.9.0 [40]). The main novelty of this workflow was the method of estimating the DEM in 2019, before the erosion process (DEMT0), from the DEM after the erosion process (DEMT1), in 2021. Subtracting the current DEM (DEMT1), which is the result of erosion processes, from the estimated DEMT0 makes it possible to obtain a DoD and quantify the volume of eroded material. A Digital Surface Model (DSM) was obtained with the photogrammetric processing of the drone-acquired images. Thanks to the overlapping of the images, a 3D point cloud (XYZ coordinates) with a very dense sampling can be constructed. The point cloud was classified to keep only the terrain points, thus avoiding vegetation and building surfaces ( Figure 2). Therefore, only ground points' Z coordinate (altitude) was interpolated to model the terrain surface and obtain a grid map containing the terrain altitude (DEMT1). This DEM was used to obtain the current drainage network by applying GIS tools based on the aspect and slope models of DEMT1, which defines the water flow direction, water flow accumulation, the channel flow network, and the watershed ridge divisions. The estimated DEM before the erosion process (DEMT0) was obtained from the interpolation of the Z coordinate (altitude) of the 3D point cloud points located on the ridges of the DEMT1. By selecting the points close to the watershed divisions and interpolating them, the terrain surface was modelled and rasterized to obtain a primary grid map with the terrain altitude (DEMT0), being a reasonable approximation of the original relief before the erosion processes considering that there was no previous direct source of data. If the original 3D point cloud was not available, a variant could be reached if the DEMT1 cells intersecting the ridge vector were converted to points, for cloud comparison, interpolation, and rasterization.
In this processing, 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 In this processing, 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 the subsequent comparison of Digital Elevation Models. Finally, an arithmetic subtraction operation was executed (DEM T0 -DEM T1 ) to obtain a difference raster that corresponds to the volume of material eroded.
The drone imagery acquisition was fully planned, automatized, and replicable by using mission planner software. The photogrammetric and SfM processing was automatized using batch processes. The GIS process to extract hydrologic and morphometric information from a DEM was fully automatized with the graphical model builder tools of ArcGIS 10.8 or QGIS 3.10 ( Figure 3).

Results
From the photogrammetric processing, an orthophotomosaic of 2.4 cm/pixel native resolution was obtained, which was resampled to 2.5 cm/pixel to facilitate raster processing and conform with the Catalan official mapping agency's cartography grid (multiple of 2.5 cm). A point cloud with a density of 122 points/m 2 was generated, which allowed us to elaborate a DEM with a resolution of 10 cm/pixel using only the points classified as ground. In the AOI the vegetation coverage is negligible, as can be observed in the orthoimage, thus avoiding possible errors related to 3D point cloud classification errors. At this spatial resolution, the photointerpretation of the relief at very high detail is possible, particularly when combining 2D and 3D models (Figure 4).

Results
From the photogrammetric processing, an orthophotomosaic of 2.4 cm/pixel native resolution was obtained, which was resampled to 2.5 cm/pixel to facilitate raster processing and conform with the Catalan official mapping agency's cartography grid (multiple of 2.5 cm). A point cloud with a density of 122 points/m 2 was generated, which allowed us to elaborate a DEM with a resolution of 10 cm/pixel using only the points classified as ground.
In the AOI the vegetation coverage is negligible, as can be observed in the orthoimage, thus avoiding possible errors related to 3D point cloud classification errors. At this spatial resolution, the photointerpretation of the relief at very high detail is possible, particularly when combining 2D and 3D models (Figure 4).

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

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 DEM T1 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 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). 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).  Once detected, the channels, together with the DEM T1 , 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 DEM T0 (Figure 6).

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 material of the slopes is not compacted and the topographic slope in some areas is very steep. The calculation of eroded material was obtained by means of an arithmetic subtraction operation (DEMT0 -DEMT1) (Figure 7).

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 material of the slopes is not compacted and the topographic slope in some areas is very steep. The calculation of eroded material was obtained by means of an arithmetic subtraction operation (DEM T0 -DEM 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 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 15,414.9 m 3 was estimated, corresponding to an erosion rate of 364 Mg/ha/year, considering a substrate density of 1.7 Mg/m 3 [46].

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 material of the slopes is not compacted and the topographic slope in some areas is very steep. The calculation of eroded material was obtained by means of an arithmetic subtraction operation (DEMT0 -DEMT1) (Figure 7).

Slopes and Erosion Identification
Slope maps are useful in hydrological studies and erosion models. To produce slope maps, it is necessary to calculate the slope of the land with respect to a horizontal plane. In this case study, the slope calculation was made based on the following categories [41]: planar (0-3 • ), slightly inclined (3-12 • ), inclined (12-30 • ), sliding (30-45 • ), and free fall (>45 • ). In terms of area, half of the land (49%) presents slopes greater than 30 • , which in theory would facilitate the drag phenomenon mainly when there is drizzle, while only 3% of the area is part of the land categorized as plain (Figure 8). 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 15,414.9 m 3 was estimated, corresponding to an erosion rate of 364 Mg/ha/year, considering a substrate density of 1.7 Mg/m 3 [46].

Slopes and Erosion Identification
Slope maps are useful in hydrological studies and erosion models. To produce slope maps, it is necessary to calculate the slope of the land with respect to a horizontal plane. In this case study, the slope calculation was made based on the following categories [41]: planar (0°-3°), slightly inclined (3°-12°), inclined (12°-30°), sliding (30°-45°), and free fall (>45°). In terms of area, half of the land (49%) presents slopes greater than 30°, which in theory would facilitate the drag phenomenon mainly when there is drizzle, while only 3% of the area is part of the land categorized as plain ( Figure 8). To analyze the relationship between the eroded area and the slope of the land, the eroded area calculation was performed for each slope category at pixel scale (Table 1).  To analyze the relationship between the eroded area and the slope of the land, the eroded area calculation was performed for each slope category at pixel scale (Table 1).

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 studies [24,25,27]. The photogrammetric and SfM processing allowed the generation of a 3D dense point cloud, orthoimage, and DEM T1 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 estimated DEM (DEM T0 ) using the points of the dense point cloud located over the watersheds of DEM T1 . 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 according 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.
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.

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 (DEM T0 ) must therefore be estimated from the current DEM (DEM T1 ). 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 DEM T0 , 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. 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).