UAV-Based Soil Erosion Assessment in Mediterranean Agricultural Orchards

1. Introduction

The Mediterranean region stands out as one of Europe’s most sensitive areas to land degradation in [1,2], due to its distinctive arid climate, land abandonment, and unsustainable agricultural practices [3,4]. While soil erosion rates in the Mediterranean are generally lower than in other parts of Europe, this often reflects the already degraded state of soils, underlining that further degradation would have severe impacts [5,6]. Mediterranean erosion rates are generally highest in croplands, where it causes a loss of fertile soil and leads to reduced productivity in agricultural systems [1]. Vineyards and olive groves are the agricultural systems that are most prone to erosion, with vineyards having the highest erosion rates by far [5,7]. The main causes for high erosion rates in these systems are low vegetation cover, tillage, and steep slopes [8]. Soil erosion in olive groves has increased over time due to the intensification of olive production and is a major threat to the sustainability of olive groves [9,10]. Furthermore, conventional management leads to reduced sustainability and yield in vineyards [11]. Meanwhile, increased periods of prolonged drought are expected to lead to an intensification of irrigation in olive groves and vineyards, and in turn, this intensification is causing more soil erosion [12].

As one of the most important drivers behind land degradation in many ecosystems, soil erosion negatively affects soil properties, reduces agricultural productivity and available agricultural land [1,13], and leads to a decrease in organic carbon and nutrients as the soil is lost [14]. Furthermore, soil erosion can lead to negative off-site effects such as decreasing water quality [15]. On the local scale, soil erosion causes agricultural systems to become unsustainable, negatively impacting the resilience of farmer livelihoods [16]. Furthermore, erosion-induced decreases in agricultural productivity are a driver of land abandonment [17]. Land abandonment can exacerbate soil erosion and stimulate rural depopulation, negatively affecting rural economies [18]. In turn, rural depopulation eventually leads to the loss of cultural values and traditions of local communities [19]. As such, the prevention of soil erosion is important for both ecological and socio-economic reasons. These factors show that there is a need for research into measures to reduce soil erosion and increase the sustainability of agricultural orchards in the Mediterranean [20].

Remote sensing is a useful tool for the spatial assessment of soil erosion and its different components and is increasingly being used in soil erosion assessment and monitoring [21]. Satellite remote sensing is often used to estimate erosion parameters such as vegetation cover and to map erosion from catchment to global scale by using the measured parameters in erosion models [1,21,22]. On smaller scales, UAV-based (Unmanned Aerial Vehicles, e.g., drones) remote sensing is a suitable alternative [23]. Whereas satellite remote sensing can only be used to measure erosion parameters, UAV-based remote sensing can directly measure soil erosion over time [24]. Furthermore, UAV-based remote sensing can assess soil erosion at a much higher resolution and can calculate soil erosion at a higher accuracy [25] and allows for the accurate assessment of microtopology [26].

UAV-based remote sensing is increasingly used to investigate the components of soil erosion on the field and catchment scale. For example, UAV-based remote sensing has been used to study spatial patterns of erosion [27], calculate sediment budgets [28], detect rill formation in agricultural fields [29], and estimate micro-topological soil parameters such as soil roughness [30]. However, for even smaller scales, such as small erosion plots within fields, UAV-based remote sensing has only been applied in a small number of studies [25], making it relevant to study the use of UAV-based remote sensing on smaller scales.

Many studies have already used field experiments to study soil erosion with conventional methods, such as sediment traps or rainfall simulations [5,31]. However, these methods often do not identify the individual components of inter-rill, rill, and gully erosion and cannot be used to study spatial patterns of sedimentation and deposition. This causes uncertainty about the actual impacts of soil erosion and the influence of different erosion processes, such as the beforementioned rill or gully erosion [32]. Meanwhile, exact measurements of surface changes with high temporal and spatial resolution are needed to investigate the development of erosion patterns over time, which is necessary to understand rill and inter-rill erosion processes [24].

UAV remote sensing can offer new methods for studying soil erosion at the plot scale [24,29,33,34]. Previous research includes only a limited number of studies using UAV remote sensing to assess soil erosion in vineyards [35], and none have focused on olive orchards. Moreover, few studies have compared the effects of different soil conservation measures between plots using UAV remote sensing [36,37], and none have tested the application of biochar; a product that, beyond its effects to reduce erosion, can also improve soil habitat and crop growth [38,39].

To address these research gaps, the present study focusses on vineyards and olive orchards in Alentejo, Portugal, a region where climate change and agricultural intensification have been shown to increase erosion risks [40], making the region an interesting location for studies of soil erosion in Mediterranean agricultural systems. The general objective of this research is to evaluate the applicability and limitations of UAV-based photogrammetry for soil erosion assessment in Mediterranean vineyards and olive groves.

Specifically, this study aims to: (i) assess the accuracy of UAV-derived DTMs and their suitability for quantifying soil erosion through DoDs, (ii) compare erosion estimates derived from DoDs with in-field erosion measurements at the plot scale, (iii) evaluate the effects of mulch and mulch–biochar treatments on soil erosion, and (iv) identify key topographic and vegetation-related variables derived from UAV data that explain spatial variability in measured soil erosion. By addressing these objectives, this study provides insight into both the potential and the current limitations of UAV photogrammetry for erosion monitoring in Mediterranean agricultural orchards.

2. Materials and Methods

2.1. Study Area and Experimental Design

This research was conducted across four study areas in Portugal’s Alentejo region: two olive orchards and two vineyards. The olive orchards were the irrigated, integrated-management Vale Formoso (VFO) area near Granja (40 km east of Évora, GPS coordinates 38°18′13.3”N 7°16′37.6”W) and the organic, rainfed Mouchão Olive (MOO) area near Casa Branca (20 km north of Évora, 38°55′03″ N 7°47′20″ W). Approximately 600 m from MOO lies the integrated-management Mouchão Vineyard (MOV) area, with the second vineyard, the Fundação Eugenio de Almeida (FEA) area, situated 15 km east of Évora (38°31′03″ N 7°44′17″ W). In the integrated-management study areas, ground cover was permitted to grow throughout, with herbicide application focused near the tree lines [10,41]. The regional climate is characterized as hot-summer Mediterranean (Köppen Csa), with a mean annual temperature of 16 °C and annual precipitation ranging from 500 to 600 mm. Soil parent materials were slate/conglomerate at VFO, conglomerate at MOO and MOV, and schist at FEA. Soils were classified as Cambisols, except for Regosols and Luvisols identified at the MOO and MOV areas. Nine erosion plots with areas ranging between 150 and 200 m² were set up in each study area, in three blocks of three plots each. Three different treatments were randomly applied to each of the plots: untreated, a mulch layer, and a mulch layer over a biochar layer applied previously. The mulch application rates were 3 Mg/ha of olive leaves in the olive groves and 2 Mg/ha of straw in the vineyards. Biochar was applied at a rate of 10 Mg/ha in both land uses. The treatments were applied in a completely randomized block design, with three blocks of three plots per location, resulting in 36 plots. Figure 1 shows a map with the locations of the study sites, and Table 1 shows a summary of the plots and treatments per study area.

2.2. In-Field Erosion and Bulk Density Measurements

Sediment fences [42] were installed at the downslope end of every erosion plot. Sediment was removed from the fences and weighed during monthly field visits. The dry weight of the sediments was determined by drying sediment samples in the lab.

Soil bulk density of each plot was determined to calculate the weight of eroded sediment estimated with the DoDs. Per plot, nine soil samples were taken: in the bottom, middle, and top of the plot, and for three micro-environments: below the tree or vine, in areas compacted by wheel passages, and in vegetated areas. All samples were weighed and then dried for 24 hours at 105 degrees Celsius. After drying, the samples were once again weighed. The bulk density was calculated as the mass of dry soil divided by the total volume of the sample.

2.3. Deriviation of Digital Elevation Models

Two sets of drone flights were conducted, in early November 2024 and in early April 2025. A Mavic 2 Pro Camera was used, with a resolution of 0.8 cm per pixel. DTMs were created from the photographs of each drone flight through a Structure from Motion photogrammetry workflow, using the Agisoft Metashape Pro software [43]. Ground Control Points (GCPs) were used to geo-reference the DTMs. Geo-referencing was necessary to correct camera calibration parameters in the Metashape software and to remove deformations in the sparse point cloud [44,45]. Coregistration was used to align the 2025 DTMs with the 2024 DTMs. Coregistration is a technique used for multitemporal topographic datasets, where the second DTM is aligned with the first DTM by identifying areas where no elevation change has taken place (e.g., rocks, roads) and using these areas as GCPs in the second DTM [44,46]. A Trimble Catalyst GPS was used to record the x, y, and z coordinates of the GCPs and several transects within the erosion plots with an average accuracy of 0.01 meters. In total, 12 transects were surveyed, one for each block of plots. As a measure of DTM model error, the Root Mean Square Error (RMSE) between the elevation of GPS points and the elevation of corresponding points in the DTMs was calculated, where the GPS data was used as ground truth data. Table 3 presents an overview of the drone and GPS surveys per study site.

2.4. Erosion Estimation from DTMs

DEMs of Difference (DoD) were calculated for each plot to estimate the erosion between the 2024 and 2025 drone flights. DoDs were calculated by subtracting the 2025 DTMs from the 2024 DTMs. The error propagation to the DoD was calculated as shown in formula 1 [47,48,49]. A “level of detail” (LoD) is used to include the DoD error in the erosion estimation. The LoD can be seen as a confidence interval for the change in elevation between the first and second DTM, where only elevation changes outside the confidence interval will be seen as actual topographic changes. Previous studies show that a 95% LoD is generally appropriate (formula 2) [47,50].

$${\sigma }_{DoD}=\pm \sqrt{{\sigma }_{DTM1}^2+{\sigma }_{DTM2}^2}$$

(1)

$$Lo{D}_{95\%}=1.96\times \sqrt{{\sigma }_{DTM1}^2+{\sigma }_{DTM2}^2}$$

(2)

where ${\sigma }_{DoD}$ is the propagate error, ${\sigma }_{DTM1}^2$ the error from the respective DTM, and $Lo{D}_{95\%}$ the 95% Level of Detail. The net eroded volume was calculated as the sum of all significant elevation differences in the DoD times the raster cell size. This volume was converted into eroded weight by multiplying by the plot average bulk density.

$${Weight}_{E_{net}}={\rho }_{plotaverage}*\sum \left({∆z}_{cell}*{Area}_{cell}\right)$$

(3)

2.5. Deriviation of Topographic Variables

Several variables were derived from the DTMs and orthomosaics, including plot area, length, and slope, as well as surface roughness and vegetation cover.

For the surface roughness calculations, the DTMs were first detrended to remove global effects such as slope that could affect the roughness calculations [24]. DTMs were detrended by fitting a linear, first-order polynomial to the DTM, and subtracting this polynomial from the original DTM. The Root Mean Squared Height (RMSH), also called the Random Roughness Index (RR) was calculated to quantify surface roughness [51].

$$RMSH=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^n(z-\bar{z}̄{)}^2}{n-1}}}$$

(4)

where z is the elevation in one cell, $\bar{z}$ the average elevation, and n the number of cells. Vegetation cover was derived from the orthomosaics using the VARI index, vegetation index that only uses the red, green, and blue wavelength bands, and does not require near-infrared data [52].

$$VARI={\displaystyle \frac{{\rho }_{Green}-{\rho }_{Red}}{{\rho }_{Green}+{\rho }_{Red}-{\rho }_{Blue}}}$$

(5)

where $\rho$equals the reflectance of the respective wavelength band. A threshold was used to determine vegetated areas from the resulting VARI index maps. All cells with a value above the threshold were classified as vegetation and below as bare soil. The threshold was selected by visual inspection of several threshold options in comparison with the original orthomosaics. The threshold differed slightly between the 2024 and 2025 drone flights and between areas but was generally between 0 and 0.1. The VARI maps were classified into vegetation and bare soil, and the vegetation cover was calculated as the percentage of vegetation cells divided by the total amount of cells per plot.

2.6. Statistical Analysis

All statistical analyses were performed in the R environment using the RStudio package. The accuracy of the net erosion estimates was assessed by calculating the RMSE and R2 metrics. The RMSE is a measure of model accuracy, with a lower RMSE meaning a more accurate estimation of erosion compared to the ground truth. The coefficient of determination (R2) is a measure of the goodness of fit for a model compared to the ground truth, with a higher value meaning a better fit.

$$RMSE=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^n(Y_{observed}-Y_{calculated}{)}^2}{n}}}$$

(6)

$$R^2={\displaystyle \frac{{cov}^2\left(Y_{observed},Y_{calculated}\right)}{var\left(Y_{observed}\right)\times var\left(Y_{calculated}\right)}}$$

(7)

Relative error in field measurements of sediment weights from erosion plots becomes larger as the actual weight of sediments in the sediment traps becomes smaller. In cases of very low erosion rates, this can significantly affect the accuracy assessment of erosion model estimations [53]. Nearing (2000) defined occurrence intervals between which the relative difference between actual and predicted erosion is deemed acceptable for different erosion rates [54]. The occurrence interval’s upper and lower boundaries of the occurrence interval are defined with the following equations:

$${Rdiff}_{occ}=\{\begin{array}{c}{Rdiff}_{low}=0.236{LOG}_{10}\left(M\right)-0.641\\ \\ {Rdiff}_{high}=-0.179{LOG}_{10}\left(M\right)+0.416\end{array}$$

(8)

As the erosion rates in some plots were very low during this study, Nearing’s effectiveness coefficient was calculated as an extra accuracy metric, as this metric can give a fairer assessment of model accuracy at very low erosion rates. The effectiveness coefficient is defined as the percentage of model predictions of which the relative difference with the field measurements falls within the 95% occurrence interval:

$$Rdiff={\displaystyle \frac{\left(Predicted-Measured\right)}{\left(Predicted+Measured\right)}}$$

(9)

$${ϵ}_{\left(\alpha =0.05\right)}={\displaystyle \frac{Numberofpredictionswithinoccurenceinterval}{Totalnumberofpredictions}}$$

(10)

The effects of the different treatments on measured erosion, soil roughness and vegetation cover parameters were statistically analyzed using one-way ANOVA, where the treatments and study area locations were the qualitative predictors. Logarithmic transformation of the data was performed in case the normality assumption was violated, which was checked with the Shapiro-Wilk test. If the ANOVA resulted in a significant outcome for the F-test, treatments were compared pairwise using Tukey’s procedure for pairwise comparisons.

A multiple linear regression was performed to investigate if the UAV-derived variables could be used to explain the erosion measured in the field. First, a correlation analysis was performed to identify variables that could be used in multiple linear regression. As several variables were not normally distributed, Spearman’s (non-parametric) rank correlation coefficient was calculated. Variables that had a significant correlation with the measured erosion were selected.

3. Results

3.1. DTM Accuracy Assessment

Table 4 shows summary statistics for the 2025 DTMs. The DTM resolution was lower than 0.02 meters for all areas. The number of coregistration points varied per area and was mainly dependent on the number of stable points that could be found between the 2024 and 2025 drone images. In general, RMSE and 95% LoD were higher in olive groves. The MOV1 and MOV2 areas were located next to each other, and both areas were surveyed in one continuous drone flight, resulting in one DTM containing both areas. However, separate GPS transects were surveyed for each area. This explains why the resolution and number of coregistration points are equal between MOV1 and MOV2, but the RMSE and 95% LoD differ.

Figure 2 shows graphs of the GPS transects per DTM, with the heights of the transect points on the y-axis and the distance along the transect on the x-axis. The blue line shows the corresponding height values of the DTMs along the transect. From Figure 2 it is visible that in the vineyards (FEA and MOV), and VFO, there was generally a good fit between GPS and DTM. The MOO DTMs had a much higher RMSE compared to the other areas, and Figure 2 shows that the GPS transects did not fit as well with the DTMs. For the 2025 MOO1 area, it was not possible to generate a correct DTM, explaining the misalignment with the GPS transect and the high RMSE.

Figure 3 compares both methods and their accuracy compared to the GPS transects. In general, the two methods performed similarly and the differences between the transects were small. In VFO3, the coregistration method performed much better, and Figure 3 shows that the transect line of the geo-referencing method follows a similar shape as the GPS and coregistration transects, however, there is a systematic difference in elevation between the GPS and geo-referencing lines.

3.3. Erosion Estimation from DTMs

Figure 4a shows the erosion predicted by the DoDs against the measured erosion. The DoDs could not accurately predict the erosion in most of the plots. In many of the plots, no erosion was detected at all. The DoDs showed a net deposition due to vegetation growth, which was not successfully filtered. After plots with net deposition were removed from the analysis, the final R2 was 0.1361, and the RMSE was 0.4377. Figure 4b shows the relative difference between the measured and predicted erosion for each plot. Relative differences that fall within the occurrence interval are deemed acceptable differences. The final effectiveness coefficient e was 0.55 for a 95% occurrence interval.

3.4. Patterns in DoD Elevation Changes

Elevation changes caused by agricultural management, external disturbances, and data processing were found to be much greater than changes caused by water-induced erosion. Figure 5 illustrates these changes, showing a selection of DoDs of the erosion plots. Elevation changes caused by agricultural management included compaction in ruts and raised sides of ruts caused by frequent tractor passages. Furthermore, in locations where vegetation was not successfully removed from the DTM, increases due to growth of vegetation strips were visible. Some areas of the plots were not covered by the drone surveys, usually due to high vegetation cover. These areas were interpolated during the creation of the DTMs, which introduced inaccurate elevation changes in the DoDs. Finally, some elevation changes were caused by soil samples taken between the 2024 and 2025 drone surveys.

3.5. Comparison Between Treatments and Areas

Figure 6 shows boxplots of measured erosion, random roughness, and vegetation cover. None of the investigated variables showed significant differences between treatments, but all variables showed significant differences between locations. Table 5 shows the outcomes of each ANOVA including p-values and post-hoc tests.

3.6. Multiple Linear Regression Analysis

Figure 7A shows a correlation matrix containing all correlations between the measured erosion and UAV-derived variables. Significant correlations were found between measured erosion, vegetation cover, and slope. No significant correlations were found between measured erosion, roughness, and plot length. As there was a high rate of collinearity between the 2024 and 2025 vegetation cover, only the initial (2024) vegetation cover was chosen for further analysis.

Figure 7B shows a forest plot displaying the results of the multiple linear regression model that explained the highest proportion of variance, while all predictors were significant. The model had an R2 of 0.72 and an RSE of 0.45. The slope and plot length predictors were not significant on their own, but the interaction between slope and plot length was significant and led to an increase in R2 and a decrease in RSE. Slope was not a significant predictor when added to the regression model. However, when an interaction term with plot length was added, the interaction between slope and plot length became a significant predictor. Model 2 shows that the slope had a regression coefficient of 0.95, and the plot length had a regression coefficient of 0.24. However, due to the interaction term, which had a regression coefficient of -0.02, the effect of each predictor is dependent on the other predictor. For example, a unit increase in slope would cause an increase in measured erosion of 0.95 - 0.02 * slope length.

4. Discussion

4.1. DTM Construction and Accuracy Assessment

All DTMs had a resolution higher than 0.02 meters. These results were similar to previous studies in agricultural orchards, although the resolution is highly dependent on several factors such as flight altitude, camera resolution, and overlap between pictures [55,56]. A comparison of the 2025 DTMs with the GPS transects showed large differences in RMSE between the different study areas. Furthermore, it showed that RMSE values were much lower in vineyards than in olive groves. This was an expected result, as areas where a larger percentage of the ground is obscured by vegetation canopy have been shown to lead to less accurate results for UAV photogrammetry [57,58]. There was a large difference in RMSE values between the two olive orchards, with MOO having much larger values than VFO. This outcome was counterintuitive, as the VFO area has a higher canopy cover and thus was expected to have higher RMSE values. However, this can be explained by the decrease in RMSE as the number of coregistration points decreased. Even so, it was not an obvious outcome that more coregistration points would lead to a lower RMSE, as the 2024 DTMs themselves were models and not true elevation data, and the GPS transects had been surveyed during the 2025 drone flights. However, what separated the 2024 DTMs from the 2025 DTMs was that the 2024 DTMs were directly georeferenced through ground control points, while the 2025 DTMS were indirectly georeferenced through alignment with the 2024 DTMs. This shows that georeferencing is crucial for accurate DTM generation, as was also concluded by previous studies [44,45,58,59]. Furthermore, it suggests that the coregistration process could be enhanced by georeferencing the 2025 DTMs with ground control points, concurrent with the alignment with the 2024 DTMs [60].

Comparing coregistration and georeferencing methods showed that, in general, coregistration did not produce better results than georeferencing. VFO 3 was the exception, as here there was a systematic error in the elevation of the georeferenced DTM. It is possible that this was caused by an error in the georeferencing process, but it could also have been caused by errors during the drone survey. Coregistration of multitemporal DTMs has mostly been used on larger scales [44,46] and is not often used in small-scale agricultural areas. This study did not show that coregistration was beneficial for erosion monitoring. The use of coregistration on a small scale in agricultural fields turned out to be limited, as it was difficult to identify sufficient areas that did not undergo any elevation change. Furthermore, whereas on larger scales, areas such as rocky outcrops can be used as coregistration points [46], coregistration in the study areas had to rely on singular rocks or other small features. It cannot be said with certainty that there was no elevation change in these features. If a coregistration point had an elevation change between the 2024 and 2025 drone flights, it could have introduced an extra error into the model. As such, the placement of permanent ground control points that remain in the study areas for the duration of monitoring could reduce the uncertainty of elevation changes in control points, reduce the error introduced by unknown elevation changes, and improve the accuracy of the DTMs [60,61].

4.2. DTM Erosion Estimation

It was not possible to accurately estimate the soil erosion in the plots, as the estimated erosion only had an R2 of 0.14 compared to the measured erosion. The accuracy in this study was much lower than other studies using a similar approach, although specific R2 values were not reported [62]. The accuracy of the DTMs and error propagation towards the DoDs could be largely held responsible for the low accuracy in erosion estimations. As in most plots, the total measured erosion did not exceed 0.2 tons per hectare; the actual elevation changes in the plots were far below the 95% Level of Detection threshold. This becomes clear when the measured erosion is calculated as depth of topsoil loss. For example, 207 grams of sediment was measured in plot FEA 1. This would equal a depth of 0.011 meters, which is significant as the 95% LoD for this plot was 0.117 meters. Similarly, plot MOV 9 saw 78 grams of erosion, equaling a depth of 0.005 meters which was smaller than the 95% LoD of 0.211 meters. This was also true for plots in the VFO area, where higher erosion rates were measured, but no erosion was found in the DoDs, as the 95% LoD was also larger in this area. This explains why no erosion was found in many of the plots. However, it does not explain why, in some plots, a net deposition was found, despite the plots being closed. Neither does it explain why the erosion of some plots was greatly overestimated.

The net deposition was caused by a combination of outliers in the DTMs and vegetation growth, which has been shown to significantly affect erosion estimates in previous studies [63]. While the filtering of vegetation with the VARI index was able to remove most vegetation, some ground vegetation was not detected and was included in the erosion estimations. If these areas with unfiltered vegetation grew during the study period, they would contribute to the deposition in the DoDs. As erosion rates were low in most plots, even small patches of unfiltered vegetation could cause a net deposition at the plot scale. The plots with greatly overestimated erosion were all found in the VFO study area. Here, a systematic difference between the two DTMs caused the 2025 DTM to be lower than the 2024 DTM, despite the coregistration. It is most likely that this systematic error was caused by the dense canopy and ground vegetation cover during the 2025 flights and ultimately shows the limitations of UAV photogrammetry in densely vegetated areas [57,58].

Another important factor explaining the low accuracy of the soil erosion estimations in this study is the natural seasonality of soil erosion in Mediterranean climates. Several studies have shown that under Mediterranean conditions, soil erosion varies strongly with season, with higher erosion occurring when vegetation cover is low and intense rainfall events occur after prolonged dry periods, and comparatively low erosion rates during wetter seasons with more vegetation cover. For instance, Ferreira and Panagopoulos (2014) used the Revised Universal Soil Loss Equation (RUSLE) approach and found that the highest predicted soil erosion rates occurred in autumn, contributing approximately 65% of the total annual erosion due to reduced vegetation cover after the dry summer, while predicted erosion in winter was much lower (≈ 7% of the total) despite relatively high rainfall erosivity during that season [64]. These findings suggest that the November-April period of the current study may be inherently less erosive in Mediterranean agricultural systems because the soil is wetter and less repellent and vegetation cover typically increases after early autumn rains, reducing the susceptibility of soil to detachment and transport [65].

While the erosion measured in the field was used as ground truth data in this study, it is important to realize that these measurements also contain a degree of error compared to the actual erosion. One source of error in measured erosion could be wind erosion, which can significantly influence soil loss in runoff plots [66]. Errors in sediment sampling, either due to equipment or human error, could also influence the error in the measured erosion [67]. Nearing et al. (1999) stated that there is a high variability in measured erosion for similar plots, and that this variability increases as erosion rates decrease. This variability will also affect the evaluation of soil erosion models [53]. As the acceptable relative difference between estimated and measured erosion becomes larger according to Nearing’s occurrence intervals, many of the relative differences between measured and predicted erosion in this study fell within the 95% occurrence intervals, resulting in an effectiveness coefficient of 0.55. Nearing’s occurrence intervals allow us to determine the performance of erosion estimations while taking this variability into account. As erosion rates become lower, the variability between erosion plots increases, which makes Nearing’s occurrence intervals especially relevant in studies with low erosion rates. This became evident from this study’s erosion estimations, as even though the R² of the DoD estimations was only 0.14, over half of the estimations could be considered acceptable according to the 95% occurrence intervals. This shows that the smaller the actual erosion is, the harder it becomes to estimate erosion through DoDs or erosion models.

4.3. Comparison Between Treatments and Areas

The analysis showed no significant differences between treatments and control. Application of mulch or mulch with biochar did not affect the measured variables at the time of this study. Previous studies have shown that the effectiveness of mulch and biochar changes over time [68,69,70], although the exact timeframe for which the application is effective is unknown. This could explain why there were no significant differences between treatments in this study. However, most studies on the long-term effects of mulch and biochar were focused on post-fire applications, and no literature for vineyards and olive groves was found, indicating a knowledge gap. As the treatments were applied one year before this study, it cannot be said if they had any effect closer to the time of application.

The ANOVA analysis showed some differences in the measured variables between some of the study areas. For the measured erosion, VFO had significant differences with all other areas, because VFO was the only area where higher erosion rates took place. For the 2024 vegetation cover, MOO had significant differences with all other areas. As MOO was the only organic orchard among the study areas, this explains the difference in vegetation cover, as organic farming has been shown to lead to higher winter vegetation cover in olive orchards compared to non-organic farming [71]. For the 2025 vegetation cover, MOO and VFO were significantly different from FEA and MOV, showing that olive groves generally have a higher vegetation cover [72].

4.4. Multiple Linear Regression Analysis

The multiple linear regression showed that a large amount of variability in measured erosion could be explained by the predictors. Compared to previous studies, the regression model performed similarly or better than comparable models in other studies [73,74,75], but did not perform as well as some studies using more advanced machine learning algorithms [74]. For example, Arnaez et al. (2007) achieved an R2 of 0.74 in vineyards, with other studies achieving similar scores through multiple regression [75]. However, Sahour et al. (2021) showed that it was possible to achieve better results using deep learning (R² = 0.95) or boosted regression trees (R² = 0.99) compared to multiple linear regression (R² = 0.77) [74]. It is important to note that the above-mentioned studies were conducted on a larger scale and time frame, and that no literature was found that used more complex machine learning techniques to investigate erosion for similar areas and time frames to this study.

The plot area was not important for the variation in measured erosion, even though it was a significant predictor in the regression. This is a logical outcome, as the measured erosion in the regression models was given as an erosion rate in tons/ha, so the plot area is already considered in the response variable. It is to be expected that if the response variable were given in total weight, the plot area would be a much more important predictor. In literature, some studies have reported a significant connection between plot area and erosion rates; larger plots could lead to a reduction in measured erosion rates due to more possibilities of sediment deposition within the plot [76,77], or erosion rates could increase due to more erosion processes taking place at larger scales [76]. However, the analysis shows that these factors did not affect erosion rates in this study. Scale effects for differences within the plot scale have been reported for smaller differences than those in this study [78]. The area may become a more relevant predictor for erosion, as the differences in area and subsequent scale effects become larger [79].

Higher initial ground cover led to lower measured erosion, which is in line with previous literature [80,81,82,83,84]. Model 4 confirmed that initial (2024) vegetation cover was an important predictor for erosion, as it shows that a model with just the initial vegetation cover can explain a good part of the variation in measured erosion. Interestingly, the initial groundcover was significant while the 2025 groundcover was not. The ground cover increased between 2024 and 2025 for the VFO study area, where the highest erosion rates were found. The MOO study area, where the least erosion took place, also saw a large increase in ground cover. High 2025 ground cover in both these areas explains why there was no correlation between the 2025 ground cover and measured erosion. Seasonal variations in ground vegetation cover could have played a role, as most of the study period took place during winter, meaning that the winter (2024) ground cover was present during most of the study period [85]. Furthermore, seasonal variations in vegetation cover have been shown to influence erosion rates [86].

As the slope and plot length interaction is comparable to the LS factor of the USLE model, and the ground cover can be seen as the cover management factor C, the regression analysis showed that erosion in the study areas was mostly impacted by the LS and C factors. These findings are in line with previous research, which states that the C and LS factors are the most important for the determination of soil loss potential [87,88]. This also becomes apparent from studies focusing on the Mediterranean [89]. López-Vicente & Navas (2009) state that “the C factor contributed most to variability; the LS and K factors contributed in a similar way” [89]. In this study, bulk density could be seen as a proxy for the erodibility (K) factor, as it has been shown that bulk density can partly explain variations in the K factor [90]. Interestingly, bulk density was not a significant predictor in the regression analysis, suggesting that the erodibility factor did not have a significant impact on erosion. It could be possible that bulk density was not a good proxy for the erodibility factor, as the high stoniness of Mediterranean soils has been shown to have a large impact on the erodibility [91]. The inclusion of the soil stone fraction in the regression analysis could have been a better representative of the K factor, especially considering that the high stoniness of Mediterranean soils is one of the main causes of low erosion rates in the Mediterranean [5,77].

4.5. Limitations

The erosion estimations were most likely influenced by disturbances that were specific to this study, including disturbance through tractor passages and disturbance by soil sampling. The tractor passages that occasionally took place compressed soil in the ruts, and pushing soil sideways caused raising of the edges of the ruts [92,93]. In the VFO area, rills were starting to form in some ruts during one field visit. However, during the next field visit, these rills disappeared due to compression caused by a tractor passage, and as such, the rills were no longer present during the 2025 drone flights and not visible on the 2025 DTMs. The tractor passages had another effect in the MOO area, where the soil was disturbed to such an extent that it was difficult to find any coregistration points. Furthermore, it was not known where and how often tractor passages occurred. In January, some students visited the study areas to collect soil samples for their research. This resulted in several holes per plot, which could have shown up in the 2025 DTMs as erosion and could have affected the erosion estimates.

This study used a Level of Detection threshold to differentiate between real elevation changes and noise. The downside of this methodology is that the smallest changes are discarded as noise, while they could have been actual elevation changes [49]. In this study, this would have limited the erosion estimations as there were mostly low elevation changes due to the low erosion rates. While many studies assume the same error for every cell of a DTM, in reality the error is spatially distributed and is dependent on topographic parameters such as slope and roughness, and data quality parameters such as point density [46]. Wheaton et al. (2010) proposed the use of a fuzzy inference system to make a spatial model of DTM error by using parameters that have been shown to influence the DTM error as inputs and showed that this led to an improvement in sediment budgets. In future studies, it could be beneficial to include spatially distributed errors in the erosion estimation, especially in areas with low erosion rates, to ensure that small, real elevation changes are not discarded as noise [49].

Vegetation cover was a limiting factor in the olive orchards, resulting in high RMSE values. The canopy cover of the olive trees caused the area below the trees to be obscured, which caused holes in the DTM under the trees, even though oblique photographs were added during the DTM generation process. In these areas, photogrammetry is the most suitable technique for DTM generation, and the use of a LIDAR scanner would lead to better results, as LIDAR scanners can penetrate vegetation and have multiple returns [58,94].

The high RMSE values of the DTMs were the factor that ultimately limited the DoD erosion estimations the most. Besides causing higher RMSE values, vegetation cover also limited the results by causing outliers in the DoDs, which further impacted the erosion estimates. The other limitations had less impact on the estimations, as even without disturbances, the erosion estimates would have been poor. If RMSE values had been much lower, disturbances caused by tractor passages and sampling would have had a larger impact on the results, especially since erosion rates were low. While the limitations discussed in this subchapter severely impacted the erosion estimation, they did not have a big impact on the results of the ANOVA analysis or the regression models.

5. Conclusions

This study investigated the use of UAVs for erosion monitoring in Mediterranean agricultural orchards by direct measurement of elevation differences through DTMs of Difference, and through the derivation of explanatory variables for erosion and multiple linear regression. Furthermore, this study investigated the use of coregistration to improve DoD estimates and reduce model error. Additionally, this study used UAV remote sensing to determine the effects of the application of mulch and mulch with biochar on erosion and other variables.

The main conclusions of this study are summarized below:

• Comparison between coregistration and georeferencing methods showed that coregistration did not lead to an improvement in DTM error compared to traditional georeferencing.
• The high RMSE values of the DTMs prevented accurate estimation of erosion with the DoDs, also because of low erosion rates during the study period.
• The ANOVA analysis showed no differences between the control, mulch, and mulch with biochar treatments. Only differences between areas were found, which could mostly be attributed to the type of agriculture and management in the study areas.
• The regression models were able to explain a large part of the variation in measured erosion and showed that the vegetation cover at the start of the rainy season contributed the most to the erosion, together with the slope length and steepness.