The experimental hillslope (Figure 1A) of the Leibniz-Centre of Agricultural Landscape Research (ZALF) in Müncheberg is located in the north-eastern part of Germany (52.6°N, 14.3°E; Deumlich et al., 2017). The site is characterized by an average annual precipitation of 547 mm (1992–2019) and an annual mean temperature of 9.3°C (https://open-research-data.zalf.de/default.aspx, DWD-ZALF Weather Station, March 2020). The soils along the hillslope are mostly Luvisols that developed from coarse-textured glacial sediments; the topsoil consists of loamy to silty sands with about 3% clay (<0.002 mm), 16% silt (0.002–0.063 mm), 81% sand (0.063–2 mm equivalent particle diameter) and about 6 g/kg of organic carbon (Deumlich et al., 2017). The arable field of the south-east exposed hillslope (length: 53.5 m, width; 6 m) ends at the footslope in tinplate funnel for runoff and sediment collection (Figures 1A,B). Cultivation was carried out together with the sowing of corn (Zea Maize, L.) with a grubber-drill combination machine on April 26, 2018 (row spacing was 0.75 m); fertilizer was mechanically applied 8 days later.

The automated runoff station at the footslope of the hillslope consists of a funnel-shaped runoff collector (Figures 1B(1)), a system of pipes and channels for distributing runoff water and sediments (Figures 1B(2)), a Coshocton-type sampler for splitting the runoff (Figures 1B(3)) with a subsurface installed automated sample collector with plastic bottles on a turntable, a runoff tank at ground level (Figures 1B(4)), for registration of the total amount of surface runoff, and a Hellmann rain gauge (Figures 1B(5)); the small tower meteorological tower was to measure wind speed, in 20, 50, 100, and 400 cm above the surface.

The plots for structure-from-motion (SfM) photogrammetry were installed on May 3, 2018, at the upper slope (US), middle slope (MS), and footslope (FS) in north-western direction on an identically-tilled area next to the large hillslope-plot (Figure 1A). Since soil erosion rates differ according to slope angle (e.g., Quan et al., 2020) these three plots were chosen for representing the different angles from 3° to 6° present at the hillslope. The potential flow lines at the soil surface of the hillslope runoff experiment determined from a digital elevation model (GlobalMapper 19.0, LiDAR, 2018; resolution: 1.2 cm × 1.2 cm, see Supplementary Appendix SA1) indicated that the SfM plots are not directly connected to the runoff collector at the footslope. The SfM plot size of 1 m length and 1.5 m width was selected such that all surface features (i.e., wheel track, non-compacted region, and 2 rows of corn) were included (Figure 1C) and the area was small enough to achieve mm-resolution due to SfM-processing. The distance between plots at footslope (FS) and MS was 16 m, and between plots at FS and US it was 38 m (Figure 1A). Replicates for the plots could not be identified at this field and were not required since the 3 plots at major hillslope positions could already sufficiently demonstrate the applicability of the SfM-technique and comparison of methods for soil settlement correction. The plots were marked by specially labelled sticks with black-and-white markers for ground control (GC) points (Figure 2). Sticks were driven into the ground down to at least 30 cm depth to ensure that their position was not affected by topsoil porosity changes. Each SfM plot received 15 sticks with GC markers, from which 3 or 4 were placed at the sides and 4 sticks were placed inside of the plot with the GC markers showing in different directions (Figure 1C). The local coordinates of the GC points were determined by using a ruler and a laser level (Einhell Bavaria BLW 400) relative to a reference point at the bottom left corner of each plot (Figures 1D(2)); UTM coordinates were obtained with GPS (Trimble Geo 7X, Handheld GNSS System, accuracy: 0.5–1.0 m) for reference points to determine the position of plots along the hillslope.

During the observation period from May 3 to May 16, 2018, two relevant rainfall events occurred on May 3 (5.8 mm) and May 15 (14.4 mm), the latter rain had the highest rain energy (311.5 J m⁻²) and erosivity (EI30) of 6 MJ mm ha⁻¹ h⁻¹ in terms of the maximal 30-min rain intensity (I30). The cumulative rainfall energy E of both events amounted to 406 J m⁻².

The bulk density, ρ _b (kg m⁻³) was determined gravimetrically using 100 cm³ intact soil cores (cylindrical steel cylinders of 5 cm height) by oven drying at 105°C for about 3 days. Samples were determined before the SfM measurements (May 3) and after the heavy rainfall (May 15) from the 1–6 cm soil depth related to the local surface elevation assuming that the value is valid initially after cultivation for most of the topsoil. The top 1 cm of soil could not be sampled without disturbance and was discarded. Since core sampling was destructive, we selected a region outside and downhill of the SfM plots for the sampling. Thus, the plot surface for SfM measurements remained intact and that the potential surface runoff from uphill was not affected by any disturbances of the soil surface. In each field campaign, 6 core samples were taken beneath each SfM plot, of which 3 samples were from the intact cultivated area and 3 soil cores from the area compacted by tractor wheels to capture the variability of soil bulk density related to visible soil structures of the plot (Figure 1C). The number of bulk density samples was limited because of limited soil area in the close vicinity of the SfM-plots that should remain intact for subsequent sampling and runoff observations. Note that each core sampling led to significant disturbance of the intact soil next to the plots. Also, the soil of the larger hillslope measurements should remain intact, thus only a relatively small area for bulk density sampling was available.

The bulk density after the rainfall event on May 15 was alternatively determined from the estimated porosity φ (t) (Linden and van Doren, 1987) as: φ = φ i − ( φ i − φ c ) ( 1 − e ( − a P − b E ) ) (1) where φ_i is the initial porosity, φ_c the final porosity of the re-consolidated soil, P [mm] is the amount of rainfall during the event and the cumulative rainfall energy E [J cm⁻²]. The bulk density ρ_b is obtained by ρ b = ( 1 − φ ) ∗ ρ s (2) where ρ_s is the density of the solid particles. The parameters for Eqs. 1, 2 are defined in Table 1. Parameters “a” and “b” in Eq. 1 were fitted manually. The optimization was based on the lowest root-mean-square-error (RMSE) in the mean between the measured and the calculated bulk densities inside and outside the tractor lane in the upper-, middle- and footslope.

Photo images of the plots for SfM-processing were taken on a daily basis with the compact digital camera SAMSUNG WB750 in approximately 1.5 m distance from the plot’s boundaries. The camera has a locked focal length, f, of 4 mm, a maximum aperture of f/3.2 (i.e., a maximum opening width of the objective lens) and a pixel size of 1.49 µm (Samsung, 2011). The DEMs of the soil surface were generated before and after two rainfall events. Each plot was photographed 30–50-times from different perspectives to ensure a spatial overlapping of the images of at least 60% as suggested previously (Westoby et al., 2012; Kaiser et al., 2015). There was no need to adjust the camera to similar perspectives or heights for subsequent photo-sessions at different days, since the camera positions were automatically determined during the image processing, and ground control points (GCP) ensured the georeferencing of the 3D-models. This was one of the major advantages of the SfM-photogrammetry in comparison to, e.g., laser scanning. A sensor size of 5 Mega Pixel (MP) (2592 × 1944 pixels) was used for the images during the initial period (May 2 to May 14) and a size of 12 MP (4096 × 3072 pixels) for the images taken until May 16. The 12 MP images were downscaled to a pixel sensor size of 5 MP before processing with Adobe Photoshop Elements (Adobe Systems, 2018 Adobe Photoshop Elements Version: 15.0 (20160905. m.97630) x64, operation system: Windows 8.1 64-Bit,Version: 8.1). The point cloud obtained with VisualSfM (Figure 2) appeared to have more evenly distributed points, as compared to a point cloud generated from original 5 MP images. This advantage of downscaling the pictures was found throughout the experiemental period, such that images were only taken in 12 MP resolution at the end of the observation period.

The image processing was carried out with the freely available software VisualSfM (Wu, 2011); the workflow (Figure 3) depicts the applied software for each of the subsequent steps, starting with the image alignment and the reconstruction of the sparse and the dense point clouds that were combined to 3D point clouds. The sparse point cloud contains all points that are found in three or more pictures (Westoby et al., 2012). The dense point cloud contains additional points that are reconstructed by the application of the CMVS- and PMVS2-algorithms. In VisualSfM, the Scale Invariant Feature Transform (SIFT) was used to identify common points and structures in the images independent of their size, illumination, and rotation. The Bundle Block Adjustment (BBA) carried out non-linear 3D spatial optimization of camera position related to GC points (Figure 1D) to find the common structures on images for generation of a condensed point cloud. The Clustering View for Multi-view Stereo (CMVS) routine divided data obtained with BBA-algorithm in smaller easier to handle point cloud clusters as a first step to aggregate the combined point cloud. Finally, the Patch-based Multi-view Stereo (PMVS2) routine was applied to independently reconstruct the 3D-spatial data clusters obtained with the CMVS as the second step in point cloud aggregation.

After the Visual SfM step (Figure 3), georeferencing of the 3D point clouds was carried out by assigning the measured local coordinates to 4 of the reconstructed GC points (Figure 1D) using MeshLab software (Cignoni et al., 2008; Cignoni, 2016). All points resulting from above-ground vegetation (i.e., maize plants) were manually removed in the May 16th surface models (i.e., end of the observation period). The Level of Detection (LoD) was then estimated using CloudCompare software (Cloud Compare, 2020 CloudCompare V2, EDF R&D/TELECOM ParisTech (ENST-TSI), Paris 2016) before exporting the DEMs derived from the point clouds to QGIS software (QGIS Development Team, 2018). The DEM generated from images after a rain event was subtracted from that derived from images before the event in QGIS to create a map of the pixel-based changes in soil surface micro-topography (for images of the workflow see Supplementary Appendix SA4). These DEMs of temporal Difference (DoD) were corrected for the uncertainty in the determination of the re-location of GC points by assigning a value of zero to all pixel values smaller than the LoD. Thus, we assumed that uncertainties caused by small differences in referencing the DEMs of two times were negligible. The plot-average changes in soil surface elevation, Δ z ¯ Δ t between two times, Δ t = t 2 − t 1 were obtained from the sum of DoD pixel values, Δ z i divided by the number of pixels, N _p = 1, … , i, as: Δ z ¯ Δ t = 1 N p ∑ i = 1 N p Δ z i . (3)

The LoD is defined as the smallest value of change in soil elevation that can actually be detected with SfM without being considered as noise. For the determination of the LoD only those GC points were deployed that were not used for georeferencing. For the LoD different definitions exist in the literature (e.g., Brasington et al., 2003). Here, the LoD between two measurements, Δ t , is calculated from the differences of the x-, y-, z-coordinates (isotropic) in all directions as L o D Δ t = 1 n ∗ k ∑ i = 1 n ∑ j = 1 k [ ( X P i , j t 2 − X P i , j t 1 ) 2 ] (4) where X denotes the k = 3 coordinates of the n = 4 GC points (P _i,j ) at times t ₁ and t ₂.

The component describing surface elevation changes due to consolidation and natural compaction (e.g., by rain impact) was considered by the mean surface elevation changes, Δ z ¯ c s Δ t (subscript cs denotes “compacted soil”), obtained from the original ( z ¯ t s   t 1 ) mean elevation (subscript ts means “tilled soil”) and the soil bulk density relations before ( ρ b t 1 ) and after the rain event ( ρ b t 2 ) (Hänsel et al., 2016) as: z ¯ c s t 2 = ( ρ b t 1 / ρ b t 2 ) ∗ z ¯ t s t 1 (5a) where the thickness of tilled soil before the rain at t ₁ was first assumed to correspond with the height of the soil core of 5 cm. The value of z ¯ c s t 2 (here in mm) depends on the thickness of the cultivated soil region, considered to be affected by consolidation; here we compare the effect of a thickness of 5 cm with that of when assuming a value of z ¯ t s t 1 of 10 cm. Other values were not considered because major elevation changes due to consolidation during single rain events are expected to occur within the uppermost 0–10 cm layer of the topsoil (Rousseva et al., 1988). The two different thickness values were assumed as possible range, since we did not know, whether the bulk density changes occurred in the upper 5 cm of the soil or reached down to 10 cm. Different consolidation was measured outside (2/3 of the plot) and inside the tractor lane (1/3 of the plot). This was accounted for by weighing the soil surface elevation change in the different parts of the plots as: z ¯ c s t 2 = ( ρ b , i t 1 / ρ b , i t 2 ) ∗ z ¯ t s t 1 ∗ 1 3 +   ( ρ b , o t 1 / ρ b , o t 2 ) ∗ z ¯ t s t 1 ∗ 2 3 (5b) The subscripts “i” and “o” denote “inside” and “outside” the tractor lane, respectively.

The plot-related mean settlement-induced component of the reduction in soil surface elevation, Δ z ¯ c s Δ t = ( z ¯ c s t 2 − z ¯ t s t 1 ) , was subtracted from the mean changes in surface elevation obtained from either the DoD maps between the two times, Δ z ¯ Δ t , Eq. 3, to yield a corrected mean value of surface elevation changes as: Δ z ¯ c o r Δ t = Δ z ¯ Δ t + Δ z ¯ c s Δ t (6) where subscript cor denotes “corrected”.

The fraction of the area of which the soil surface elevation increased, decreased or remained unchanged was calculated for the main soil surface structural areas of the plots. The surface structures inside the tractor lane (in TL), outside the tractor lane (out TL), and the seed row (SR) were defined and manually distinguished according to the visible structures in the DoDs. The DoDs of each plot were reduced to the individual soil surface structure and pixels were classified according to increase (+), decrease (−) and no change (0) in soil height and the number of pixels in each class was summed up. The area fraction of each class of soil surface structural feature was obtained (c.f., Eq. 3) from the sum of pixels divided by the total number of pixels.

Throughout the process of soil loss determination by SfM several errors accumulate: 1) GC points were manually levelled, 2) the generation of the dense point clouds in Visual SfM depends on the image quality and leads in case of low quality images to a lesser point density causing errors when creating surface models 3) georeferencing errors, and 4) errors in soil loss calculation from soil elevation changes due to natural consolidation.

Here, the levelling the GC points with the laser level with an accuracy of ±1 mm. The standard deviation of georeferencing the point clouds in MeshLab and the calculation of the Level of Detection (LoD) in CloudCompare amounts to 0.7 mm. For the bulk density measurements used for the correction of soil settlement, a mean standard deviation for all plots of 70 kg m⁻³ was obtained. This value resulted in an error of 0.5 mm if a 5 cm topsoil layer z ¯ t s t 1 and of 1.1 mm if a 10 cm topsoil layer z ¯ t s t 1 is assumed (Eq. 5a). Thus, these errors add to a maximal value of either 2.2–2.8 mm if all sources of possible inaccuracies in data acquisition and processing are considered.

Soil surface roughness, as an important input for soil erosion models (e.g., Kaiser et al., 2015), was calculated in QGIS by employing the roughness algorithm derived from the GDAL DEM utility (QGIS Development Team 2014). This algorithm derives the roughness from the Terrain Ruggedness Index (TRI) according to Wilson et al. (2007) by averaging the absolute values of the differences in height between a pixel and its 8 neighbors. This algorithm was used because it provided a convenient and quick possibility to characterize morphologic soil surface changes within the plots. It was the only algorithm currently implemented in QGIS to determine soil roughness. In order to obtain an average roughness of the whole plot, the values of all pixels per plot were averaged.

Soil bulk density was initially higher in the tractor lanes (1,430–1,480 kg m⁻³) than in the soil regions between tractor lanes with 1,180–1,270 kg m⁻³ (Table 2, for statistical analyses see Supplementary Appendix SA2). Measured soil bulk density increased within the 20 days in most plots, except for the tractor lane regions on the middle and upper slopes; however, this effect of a decrease in soil bulk density was smaller than the standard deviation and negligible. The soil outside of the tractor lanes was generally more compacted at the second date indicated by a density increase of about 80–100 kg m⁻³ as compared to the first date and the soil within tractor lanes. On the plot at the footslope position, soil regions in and outside the tractor lanes were similarly more compacted as indicated by a density increase of 70–80 kg m⁻³. According to statistical analysis the bulk density differences between inside and outside the tractor lane were not significant, except for the middle and upper slope on May 2, 2018 (Supplementary Appendix SA2). The predicted bulk density for the May 22 was similar to the measured bulk density outside the tractor lane. However, modelling showed a stronger increase in soil bulk density inside the tractor lane than the measurements suggested (Table 2).

The final DoD-maps (Figure 4) show spatially-distributed patterns of increasing (green) and decreasing (red) soil surface elevations. The tractor lanes and the seed rows could be identified more clearly in the individual DEMs provided in the Appendix (Supplementary Appendix SA5). For the upper slope position, only relatively small changes in soil surface topography are noticed during the first 2 days after the installation of the plots (Figure 4A, top). Settlement of the soil can be observed more in the less compacted right part of the plot as compared to the more compacted seed rows and tractor lanes. After 12 days (Figure 4A, centre), the settling of the soil was more pronounced (i.e., more decreasing surface elevations) also in the plot region of the initially more compacted soil. Note that the first leaves of the maize plants could be identified as spots of larger elevation increase along the seed row in the middle of the plot. After the heavy rainfall event on May 15, the increased red spots in the lowest DoD map (Figure 4A, bottom) indicated a larger decrease in soil surface elevation esp. in the less and more compacted regions of the plot, where more than 80% of the area was subject to soil height loss (Figure 5).

For the plots at the middle- (Figure 4B) and the footslope (Figure 4C) positions, the changes in soil surface topography are relatively similar during the first period between May 2, 4, and 14 (upper and central rows of the maps in Figure 4). An exception is the plot at the middle slope: the decrease in surface elevation was stronger in the more compacted tractor lane than in the looser region of the plot (Figure 4B, top). Twelve days later, the situation has changed completely: now the loose area shows higher settlement than the tractor lane (90% of the area outside the tractor lane was subject to soil height loss, Figure 5). This might be attributed to the different values in the level of detection, LoD, assigned to the DoDs. Changes in soil surface elevation could be masked by a higher limit used as LoD. On May 16, due to rainfall the soil surface elevation is more strongly decreasing (red spots) than increasing (green spots) especially for the plot at the middle slope (Figure 4B, bottom). Deposition of soil can be observed at the middle slope plot especially at the lower end of the seed row and in the imprints of the tyres. For the plot at the footslope position, the regions with an increase in surface elevation or deposition are largest and oriented along the seed row along the central part (Figure 4C, bottom and Figure 5); for this plot at the footslope, a gradient in surface elevation changes was observed ranging from decreasing elevations (erosion) in the upper to increasing elevations (deposition) in the footslope regions. Red and green spots next to each other, especially for plots at the middle and the bottom slope position (Figures 4B,C, maps at the bottom), indicate a flattening of the soil surface topography due to deposition, of soil particles that were detached at the local peaks for example, within the wheel tire marks.

Between May 02 and 16, the plot at the upper slope position experienced the smallest changes in mean surface elevation. This means that here a larger area of the plot remained at the same surface elevation than at the middle slope and footslope (Figure 5). As observed qualitatively the middle slope had the highest surface area fraction with soil height loss in all soil surface structural regions (>80%). For the footslope, the seed row experienced the most pronounced increase in soil height, since 40% of the seed row area fraction increased in soil height (Figure 5).

The plot-scale averages of the changes in soil surface elevation obtained from the DoD-maps between May 14 and May 16 were all larger as compared to the hillslope-related calculated soil loss of the sediment yield of the collector station (Table 3). The plot-scale mean elevation changes were all negative except for the plot at footslope position, if considering no soil consolidation and consolidation related to 5 cm topsoil (Figure 6). When consolidation was related to 10 cm topsoil, all plots (upper-, middle- and footslope) showed increase in soil surface elevation, except when no consolidation was accounted for (Figure 6).

The plot-scale averages of the changes in soil surface elevation depict values between −6 and −2 mm that are dependent on the slope (Figure 7A). One exception is the positive value found for the plot at the bottom slope position indicating predominately deposition during the rain event. The plot at the upper slope (Figure 7A) with the smallest inclination (6%) showed smallest elevation changes (−1.8 mm) but relatively large changes due to consolidation of approx. -5 mm (Table 3). The plot at the middle slope with 9% inclination experiences the most pronounced surface changes (−2.2 mm) and consolidation (up to −5.4 mm). Although the inclination was largest at the plot at the footslope (11%), only relatively small changes (−1.7 to −3.1 mm) in the mean surface elevation were observed before and deposition (+3 mm) was observed after the rainfall event on May 15 (Table 3). The plot at the middle slope shows the highest losses in surface elevation in contrast to the plots at the upper- and footslope, which corresponded with the loss in soil volume and soil mass (Table 3).

The average soil roughness in terms of the Terrain Ruggedness Index (TRI) decreased at all plots from May 2 till May 16 (Figure 7B). The plot at the footslope had the highest initial soil roughness with 2.5 mm, whereas the upper slope had the smallest initial soil roughness (2.2 mm). Throughout the observation period, the footslope position showed the strongest decrease of TRI-values from 2.5 to 1.9 mm, which corresponds to the deposition and levelling of the soil surface that was observed for the plot surface at the footslope position. At middle and upper slope positions, TRI values decreased only by approximately 0.4 mm from 2.3 to 2.0 mm and 2.1 to 1.7 mm. The soil roughness of the soil surface at the footslope position most gradually decreased (Figure 7B), while at the middle and upper slopes, the surface roughness decreased initially more strongly from May 2 till May 4, remained at the same level till May 15, and decreased again after the rainfall on May 15.

The subplots for soil erosion assessment with SfM, installed at the upper, middle, and footslope position of the hillslope differed with respect to slope inclination. After correcting the mean soil surface elevation changes of these plots by using predicted bulk density changes, the elevation changes were generally smaller as compared to the changes obtained by using measured bulk density for the correction (Table 4). Both, measurements and predictions, resulted in an overall reduction in original mean surface elevation of the tilled soil ( z ¯ t s t 1 ) for plots at the upper and the middle slopes. However, an increase in soil surface elevation was found for the plot at footslope position when the settlement correction was related to a 5 cm soil layer. When bulk density changes were assumed to affect a soil layer of 10 cm thickness, a net deposition was obtained for plots at the upper, middle, and footslope positions. Here, the predictions of the soil settlement correction led to larger values of the surface elevation changes as compared to the measurements (Table 4).

Although the soil loss observed by SfM at the smaller plots cannot quantitatively be compared to soil loss by surface runoff and erosion at larger slopes (Parsons, 2019, see discussion 4.4), an extrapolation may be useful to check the plausibility range of the SfM data and evaluate ranges and relative soil masses quantified with the SfM method. Still, for comparing subplot information with those of the hillslope, the soil mass changes determined at the SfM-plots need to be weighted, here according to the upper-, middle-, and footslope area fractions (Table 4). Note that the typical relations between slope lengths and relative sizes of slope sub-divisional areas (Jha et al., 2015) was modified by increasing the length of the footslope sub-division from 3 to 10 m and proportionally reducing the length of the middle slope to account for the observation of deposition along this 10 m slope region (Supplementary Appendix SA3). The soil loss measured in the collector station was in the range of SfM-measurements corrected by using the predicted bulk density changes according to Linden and van Doren (1987) when assuming a topsoil height of 5 cm. With the correction based on measured bulk density values related to a 5 cm thick topsoil, the soil loss extrapolated by the SfM-technique for the total slope was 8- to 9-times that observed at the hillslope station (Table 4). The correction of the soil settlement applied to a topsoil height of 10 cm led to an increase in soil surface elevation after the erosion event for all plots (Figure 6 and Table 4).

Soil surface structural changes due to raindrop impact could be quantified with the SfM-technique at the three hillslope positions. The plot located near the footslope received more sediments than were eroded (Table 3) because in these regions, any surface runoff coming from upper slope regions is saturated with sediments and cannot take up more soil particles (Schmidt, 1996). For plots at middle and upper slope positions, surface runoff only locally affected the surface roughness (Figure 7B, Supplementary Appendix SA5) and probably generated not sufficient kinetic energy for initiating larger-scale erosion under the present conditions and slope angles. The observed internal distribution of the soil within one SfM-plot with similar slope angle (10%) was also found in the study by Quan et al. (2020) with soil elevation increase at the lower end of the plot. Kaiser et al. (2018) and Hänsel et al. (2016) reported similar local redistribution of soil from the exposed higher elevated plot regions into the depressions and lower elevated regions.

Between May 02 and 16, the plots show a decrease in soil surface elevation in all surface structural units (Figure 5). In the plot at upper slope position, larger area fractions remained unaffected from these changes than in the plots at middle and footslope positions. Since soil erosion is unlikely to occur so far upslope due to insufficient kinetic energy of the runoff (Schmidt, 1996) and according to visual inspection, soil consolidation must have been the main reason for soil surface elevation reduction. Note that the area fraction for which elevation changes were observed with SfM for the soil surface structural regions is not directly related to the amount of soil erosion in the plot. Relatively small areas of the plot with relatively large elevation changes could have affected the overall plot scale mass changes.

The decline in soil surface roughness (Figure 7B) could be attributed to the collapse of soil aggregates. Mechanisms for aggregate breakdown have been attributed to contractive forces of water in menisci between soil particles (Hartge et al., 2014), the decrease in structural stability of soil aggregates during wetting (Bergsma und Valenzuela, 1981), the destruction of soil aggregates by rain drop impact (Bolt and Koenigs, 1972) and the subsequent transport of smaller particles into larger pores (Schmidt, 1988). The soil surface maps at the footslope position with the strongest decrease in soil roughness supported the observation that soil particles were deposited in the local depressions leading to a levelling of the profile (Figures 4, 5, 7B). The small effect of slope inclination on soil surface dynamics was probably related to relatively high infiltration rates at the upper and middle slope positions, and the focused surface runoff on cultivation-induced features such as wheel tracks and plant rows (Supplementary Appendix SA3).

Soil loss estimated by SfM-photogrammetry was smaller when soil consolidation was accounted for either by bulk density measurements or predictions (Eq. 1) according to Linden and van Doren (1987) (Table 4). This seems more plausible because the overall soil loss of the total slope was significantly smaller than those derived from plot-scale data (Table 3). Soil consolidation needs to be considered as an important source for soil elevation decrease of a freshly cultivated soil in temperate climate zones (Schmidt, 1988).

When comparing the measurement with the prediction of soil settlement correction, similar trends for the overall change in soil height of the plots can be found: A decrease in soil height for the upper and middle slopes and increase (i.e., net deposition) for the footslope position (Table 4). The prediction method has the advantage that no elaborate bulk density measurements are needed after each precipitation event. Measurement errors of the soil bulk density can lead to larger errors, so a lot of soil samples need to bet taken and larger soil areas need to be disturbed for the sampling. In case of the bulk density prediction, however, site-specific characteristics cannot directly be included in the analysis. Additionally, data of bulk density changes as determined here by core sampling, cannot clearly distinguish between natural soil settlement and erosion or deposition (e.g., Knapen et al., 2008). Both processes might lead to a change in bulk density in addition to natural consolidation due to raindrop impact. An alternative approach in this case would be to determine the bulk density changes in a levelled area (i.e., control plot) that is not subject to soil erosion but only to natural soil consolidation via rainfall. However, this plot would have to be located in the vicinity to the sloped plot, to ensure comparable raindrop impact and soil conditions. To install such a control plot in the field, might be challenging; and could be only tested under simulated rainfall and in the laboratory (e.g., Kaiser et al., 2018). The settlement prediction considers bulk density changes only due to raindrop impact for levelled plots and thus predicts final bulk densities without the impact. The advantage of the prediction over the direct measurements is that the model could be calibrated by fitting parameters “a” and “b” to site-specific bulk density changes observed when erosion impact could be excluded.

Besides deviations caused by conditions in the field, uncertainties might occur also throughout the SfM-processing due to low precision in georeferencing of the 3D point clouds in MeshLab. Because of limited computational power, dense point clouds were not generated for every single part of the plots. In most cases, the centre of the GCPs was not exactly represented by a single point but was rather located in between two points. Consequently, throughout the georeferencing process only one of the points located a certain distance away from the actual GCP centre could be chosen for georeferencing leading to a deviation from the real coordinates. The described georeferencing error has been accounted for by considering a detection level, LoD (Figure 7, legends).

Between the points of the 3D point cloud, an interpolation was carried out in areas with a low point density during DEM generation in CloudCompare. This is the case especially in the regions close to the plot boundaries, where the coverage with images was lower than in the plot centre. For every time step, VisualSfM produced different point clouds depending on the photo images. This was also the case, when two 3D models of the same object were generated from a different set of pictures. Hence, for both 3D models, different point clouds existed as a template for the DEM generation so that interpolation between the points was different leading to differences in the DEMs of the same object. This interpolation error increased with the complexity of an object’s surface. Since soil surfaces were rather heterogeneous, this error was probably important. A possible solution could be to use pictures with a higher image resolution (Figure 2; i.e., from 5 to 12 MP). Unfortunately, VisualSfM software was unable to process such highly resolved data. By the use of downscaled pictures from 12 MP to 5 MP in Photoshop Elements, the point cloud density was not increased but the points were more evenly distributed throughout the point cloud. Other software such as PhotoScan (Agisoft, 2018) would be better able to handle a variable amount of data points (Jiang et al., 2020). However, this software was not available and required more computational power.

Challenges of small-scale erosion quantification by SfM and future needs.

The SfM-photogrammetry proved to be a useful tool to observe small scale soil surface micro-topography and structural changes at three plots or subplots along a hillslope. The advantage of our case study carried out in combination with the hillslope erosion experiment was that the same agricultural management was carried out uniformly over the whole field and that the basic conditions, soil, crop, tillage, and weather information could be directly used and compared with the complete hillslope. However, the soil loss found at the SfM-plots could not be related to that measured at the hillslope collector station for several reasons: For a start, the origin of the sediments collected at the footslope is uncertain, and according to the surface flow lines, sediments may have also passed the funneled collector (c.f. Supplementary Appendix SA1). It is not clear, where the sediment in the collector station might have come from. Travel distances of particles is finite and small (Parsons et al., 2010), thus the small plots can only estimate local redistribution. Also, the suggested approach to relate soil surface elevation changes of a smaller slope to the average soil surface elevation changes of a larger slope (Table 4) is dependent on the empirically adjusted slope length among other factors. Thus, this approach is site specific and cannot be transferred to other areas. Similar comparisons of smaller plots to larger slopes (Chaplot and Poesen, 2012) gave considerably higher sediment delivery rates from 1 m² plots as compared to hillslope-scale (899 versus 4.3 g m⁻² y⁻¹). These authors attributed this discrepancy to splash erosion being the dominant sediment detachment and transport mechanism at hillslopes. Martinez et al. (2017) also found lower sediment yields at larger plots (27 m²) as compared to smaller plots (0.7 m²). In contrast, Boix-Fayos et al. (2007) found higher sediment concentrations at larger plots (30 m²) than in 1 m²-sized plots.

Thus, the observed discrepancies between the different soil loss estimation techniques in this study can be attributed to the smaller size of the plots used for the SfM-measurements (1.5 m²) in contrast to the collector station that accumulates the eroded sediment from a 318 m² hillslope. The SfM-plots reveal the local deposition and erosion processes and do not allow estimating processes between plots (Parsons et al., 2010). Any comparison would improve, if DoD maps of the hillslope were generated provided the SfM-technique could be applied to the total area. Unfortunately, the resolution of the DoDs for a larger area would still be too coarse, the identification of effects of rain events on surface structure dynamics is limited (Kaiser et al., 2018). On the other hand, one could separate larger hillslopes into smaller areas (1–3 m²) that are each observed in detail with SfM and finally merged into a large DoD maps.

The SfM measurements basically provide quantitative and spatially-distributed information on the surface topography; it is not possible to distinguish between deposition of soil from uphill and erosion of soil that left the plot and the settlement. Furthermore, the change in the surface micro-topography includes the decline in surface roughness after rain. This may be considered as a kind of local erosion and deposition, which makes it difficult to separate between the deposition from uphill and local processes. The separation between input and output from changes in mean surface elevation requires additional assumptions that could be based on observations at the neighboring hillslope as follows:

Upper slope position: Based on observations it may be assumed that here the deposition from above was negligibly small such that the changes in surface elevation can be explained by runoff soil loss and by settlement.

The soil surface at the middle slope is in a through-flux position and has both deposition from above and soil loss towards downhill positions. Soil settlement could be the main unknown when assuming that lateral inputs equal outputs of soil mass. At the footslope position, there is clearly more deposition than erosion such it is assumed that the surface elevation changes account for net accumulation and some settlement.

Note that the observations do not allow to exactly quantify the rates of the different components of the soil mass changes but we can provide information on potentially relevant limits by making estimates when assuming minimal and maximal ranges limits from the comparison with the data obtained at the complete hillslope.