SMODERP2D—Sheet and Rill Runoff Routine Validation at Three Scale Levels

: Water erosion is the main cause of soil degradation in agricultural areas. Rill erosion can contribute vastly to the overall erosion rate. It is therefore crucial to identify areas prone to rill erosion in order to protect soil quality. Research on rainfall-runoff and subsequent sediment transport processes is often based on observing these processes at several scales, followed by a mathematical description of the observations. This paper presents the use of a combination of data obtained by different approaches at multiple scales to validate the SMODERP2D episodic hydrological-erosion model. This model describes infiltration, surface retention, surface runoff, and rill flow processes. In the model, the surface runoff generation is based on a water balance equation and is described by two separate processes: (a) for sheet flow, the model uses the kinematic wave approximation, which has been parameterized for individual soil textural classes using laboratory rainfall simulations, and (b) for rill flow, the Manning formula is used. Rill flow occurs if the critical water level of sheet flow is exceeded. The concept of model validation presented here uses datasets at different scales to study the surface runoff and erosion processes on the Býkovice agricultural catchment. The first dataset consisted of runoff generated by simulated rainfall on plots with dimensions of 2 × 8 m. The second dataset consisted of the runoff response to natural rainfall events obtained from long-term monitoring of 50 m 2 plots. These two datasets were used to validate and calibrate the sheet flow and infiltration parameters. The third dataset consisted of occurrence maps of rills formed during heavy rainfalls obtained using remote sensing methods on a field plot with an area of 36.6 ha. This last dataset was used to validate the threshold critical water level that is responsible in the model for rill flow initiation in the SMODERP2D model. The validation and the calibration of the surface runoff are performed well according to the Nash–Sutcliffe efficiency coefficient. The scale effect was evident in the 50 m 2 plots where parameters lower than the mean best fit the measured data. At the field plot scale, pixels with measured rills covered 5% of the total area. The best model solution achieved a similar rill cover for a vegetated soil surface. The model tended to overestimate the occurrence of rills in the case of simulations with bare soil. Although rills occurred both in the model and in the monitored data in many model runs, a spatial mismatch was often observed. This mismatch was caused by flow routing algorithm displacement of the runoff path. The suitability of the validation and calibration process at various spatial scales has been demonstrated. In a future study, data will be obtained from various localities with various land uses and meteorological conditions to confirm the transferability of the procedure.


Introduction
Soil loss due to water erosion has become an increasingly important topic with ongoing climate change, as a result of which more extreme precipitation events in central Europe are anticipated. The surface runoff from agricultural land is the main source of diffuse pollution (e.g., [1]), nutrient loss (e.g., [2]) and erosion (e.g., [3]). The amount of soil transported from an agricultural field can differ substantially depending on the intensity of rill erosion. Improved process-based surface runoff-erosion models and their usage for predicting extreme erosion events is crucial for the design of effective erosion control measures.
Empirical methods such as the Revised Universal Soil Loss Equation (RUSLE) are often used as a simple method for long-term modelling of the sediment yield [4]. Limitations of this method are that it is not able to estimate the soil loss from separate rainfall episodes and that it functions well only in areas without concentrated flow. These limitations make the application of RUSLE difficult, because extremely eroded areas must be excluded from the simulations. Physically based surface runoff and erosion models often calculate sheet and rill processes separately, and can therefore be used to identify rills.
The description of the surface runoff is often based on a kinematic or diffusion wave approximation of the Saint-Venant equation. For example, kinematic wave for solving surface runoff problems have been used, e.g., in [5,6], and in the HEC-1 model [7]. The diffusion wave approximation has also been widely used to solve surface runoff problems, e.g., in [8], and in software such as MIKE-SHE (DHI, Hørsholm, Denmark) [9], CASC2D (Colorado State University, Fort Collins, CO, USA) [10], WEPP (USDA, Washington, D.C., USA) [11], and SWAT (USDA-Agricultural Research Service, Temple, TX, USA) [12]. In Erosion3D [13], the model uses the Manning formula [14] for overland flow calculations. In addition, several studies have been published on quantifying the justification of a kinematic or diffusion wave approximation [15,16]. In this study, we utilize the SMODERP2D model, which is based on a kinematic approach to sheet and rill calculation. A simple comparison between the SMODERP2D model and the Erosion3D model was performed, focusing on the way in which these models and the RUSLE method identify extreme erosion [17].
The SMODERP2D model has been under development since the 1980s. The first versions of SMODERP1D worked as a 1D profile model, which calculated the tension along the slope caused by the overland flow to estimate the soil erosion and to design erosion control measures [18][19][20]. The key feature of the model was that it provided water flow parameters for various soil types, which were measured at the tilting flume. The user could choose water flow parameters from easy-to-measure soil characteristics, such as soil texture or clay content [18,21]. Further development of the model led to the inclusion of rill flow routines that could also estimate the rill development. In contrast with other erosion models, the rills are calculated explicitly and act as small water channels in DEM cells, in which the flow is governed by the Manning equation. Advances in software and hardware enabled the development of a 2D version of SMODERP1D. This model retained the same philosophy of rill flow calculation, and water routing was added in order to calculate the overland flow and the erosion in a complex landscape.
The application of physically based models requires calibration and validation which depends on measured data. Experimental field plots are most often used for erosion models, and RUSLE method [22] data can be used for runoff models. Plots of different length can be used [23], In [24], a plot length of 22.1 m was used. However, erosion events happen rather sparsely, and it is therefore difficult to obtain representative datasets for a study of processes that affected surface runoff.
Rainfall simulators (RS) can serve as a tool for obtaining data under controlled initial and boundary conditions [25,26]. RS allows the collection of precise replications with different rainfall [27]. The first rainfall simulators were designed in the 1930s [28]. The size of the irrigated plots can vary from less than one square meter to hundreds of square meters [29,30]. The generated raindrop sizes and their kinetic energy, usually measured with disdrometers, are crucial for rainfall-runoff and erosion experiments (e.g., [31,32]). The effect of soil freezing [33] or of the crop cover [34], both of which affect surface runoff generation and erosion, can also be studied using RS.
Photogrammetric methods are used to identify the impact of surface runoff, especially concentrated runoff in rills. Remote sensing data with the apparent consequences of a heavy rainstorm event are used to derive the pattern of erosion rills and ephemeral gullies. Aerial photographs [35] or UAV images [36] are used on field plots or small catchments. Supervised or unsupervised classification [37] or object-oriented analysis [38] can be used to quantify and assess erosion processes. A volumetric method [39] can be used for estimating the amount of eroded materials. The digital elevation model (DEM) can be obtained in high resolution from aerial photographs obtained from UAV data [36] or from LiDAR data, and can be used for rill erosion studies. For example, detailed changes in surface microtopography caused by runoff have been derived by Structure from Motion (SfM) methods [40], by terrestrial LIDAR [41], and by stereoscopic imagery [42]. High resolution DEMs have been used to study the impact of splash erosion [43].
The paper presented here describes methods for calibrating and for validating the SMODERP2D rainfall-runoff erosion model on an experimental field plot with three types of data at various scales. The derivation of the model runoff equation was based on laboratory rainfall simulator data. Laboratory conditions allow researchers to better control the conditions and to perform more repetitions, but they are carried out on reshaped soils. Although the topsoil is reshaped in the course of an agrotechnical year in the conditions of the Czech Republic, and the laboratory rainfall simulations are therefore to some extent representative, it is important to verify the model under natural conditions and with natural rains. Rainfall simulations on 16 m 2 plots and natural single events on plots 22.1 m × 2.27 m in size with a slope of 9% were used to calibrate and validate the flow routine of the SMODERP2D model. Remote sensing data, which provided the effects of a heavy rainstorm on rill development, was used to validate the rill flow routines in the model. The principles for validating the SMODERP2D model are shown for one field block.

SMODERP2D
The SMODERP2D episodic rainfall-runoff/erosion model was used for the investigation presented here. The 2D model is based on a 1D profile version, in which the surface runoff and the erosion were typically calculated in several 1D profiles representing the main flow path in the hillslope [20]. The current generation of the SMODERP2D model is pixel distributed, and is implemented in python in order to be compatible with most GIS software. The development is presented on the github platform at github.com/storm-fsvcvut/smoderp2d (accessed on 1 January 2022).
The SMODERP2D model is primarily designed for surface runoff and erosion computation. The surface flow routing in the model is based on the digital elevation model (DEM). DEM also controls the spatial discretization of the model. The principle of the model is the cell-by-cell mass balance calculated in each time step. The change in the water level of the shear flow in any cell in time is controlled using the equation: where h is the surface water level (L), qin and qout are the sheet overland inflow and outflow into and out of a given raster cell (L.t −1 ), ep is the effective precipitation intensity (the potential precipitation reduced by the interception zone and the surface retention) (L.t −1 ), and inf is the infiltration rate (L.t −1 ). The kinematic wave approach is used in the calculation of the overland flow. The momentum is therefore expressed by the power-law: where a and b are power-law parameters. Equation (2) can be expressed in the form of the Manning-Strickler formula where X and Y are empirical parameters and s is the surface slope (L.L −1 ). The infiltration component of Equation (1) is calculated by the Philip infiltration equation [44] where S stands for sorptivity (L.t 1/2 ) and Ks stands for saturation hydraulic conductivity (L.t −1 ). The SMODERP2D model is subjected to uniform rainfall. The potential precipitation is reduced due to surface retention. The surface retention is the storage that needs to be filled before surface runoff can occur.
The flow routing of the surface runoff is based on the D8 one-direction flow algorithm [45]. The inflow to cell i is defined as the sum of the sheet outflows from the adjacent cells, as: where j is the index of the adjacent up-slope cells identified by the D8 flow algorithm. The time derivative in Equation (1) is calculated using an explicit method. The computation is therefore sensitive to the size of the time step. The size of the time step is controlled by the Courant criterion, which needs to be kept below the theoretical maximum value of 1, while the maximum value in practise is lower than 1 [46,47].
The sheet flow water level of the next time t + 1 step in Equation (1) which incorporates the sum (5) is calculated with the explicit time discretisation scheme for cell i as:

Rill Flow
The rill flow is also included in the calculation. In SMODERP2D, rill flow in a cell occurs if h > h , where h is the critical water level. The water flow corresponding to the water level above the critical water level has enough energy to start to carry the soil particles and to create a rill.
The critical water level h is calculated as: where h is the water corresponding to the critical velocity, and h is the water level corresponding to the critical shear stress. As shown in Formula (7), h uses several values obtained with a different approach. This approach is adopted in order to remain on the safe side of the emergence of a rill, since h is more sensitive to the sheet flow velocity and h is more sensitive to the slope of the soil surface. When the condition h > h is fulfilled, a rill starts to develop in a given cell and h = h . In SMODERP2D, the rill is a dynamic component and can increase as the rill flow increases. The rill volume is controlled by the volume of water corresponding to the rill water level h . This volume is calculated as: where: P stands for the size of the raster cell. The rill is simplified as a small channel at the soil surface with a rectangular cross section. The rectangle has width b and rill height y = 0.7b . The rill flow is as calculated with the Manning equation: where A is the cross section of the rill, n is the roughness in the rill, s is the surface slope, and is the hydraulic radius of the rill. As stated above, the size of the rill changes as ℎ increases. The scheme of this change is shown in Figure 1. The change in the rill flow varies with the change in in Equation (10), since ℎ increases, and therefore and also increase. The for an increasing rill is calculated as: During the recession limb of the hydrograph, the rill size "locks", and ℎ decreases until the rill is empty. The scheme of the emptying rill and the rill flow is shown in Figure  2. In this case, is calculated from fixed and decreasing ℎ . for decreasing rill flow is calculated as: The total water balance in cell i, where a rill is developed, is calculated as: where: and: The rill water level is recalculated to cover the whole cell and not just the bottom of the rill, as shown in Figures 1 and 2.

Study Site
The Býkovice catchment is located in the central region of the Czech Republic in the Benešov municipality ( Figure 3). The elevation of the catchment is 370 to 510 m.a.s.l. The study area was established in 2011, and it consists of a single agricultural field. The field is divided into two subplots with different crops, and there are no erosion control measures. The topography of the catchment consists of one major thalweg. The total area of the field is 39.8 ha. The dominant soil type is Cambisol, which is typical for agricultural areas of the Czech Republic. The spatially distributed hydraulic conductivity is shown in Figure 3b. The depth of the topsoil was estimated to be 40 cm. According to particle size analyses of the upper horizon, the topsoil texture consists of sandy loam with stoniness up to 25%. The average annual precipitation is 600-700 mm, and the average annual temperature is 7-8 °C [49]. According to Kašpar [50], the sub-daily precipitation with a 20-year return period is 52.8 mm. The characteristic shapes of six-hour precipitation [51] are 37% "uninterrupted episodes" and 39% "most concentrated episodes". The precipitation and the flow in the stream have been measured since 2005. Long-term monitoring of surface runoff and erosion has been carried out since 2009. In 2020, a meteorological station (precipitation, air temperature, wind speed) and soil moisture sensors at three depths were established in the area. Three Wischmeier plots (22.13 × 2.27 m; slope 9%) were established to monitor the overland flow and erosion from 2010 to 2014 with tipping buckets to monitor the overland flow.

Plot-Scale Sheet Flow Calibration and Validation
The surface runoff modelled by SMODERP2D was calibrated and validated (i) on measurements monitoring the surface runoff process under controlled conditions using a rainfall simulator, and (ii) on measurements of natural precipitation events on the monitored Wischmeier plots.
A set of experiments with a field rainfall simulator were carried out in 2013 and 2014. A rain simulator with a rained area of 8 × 2 m was used. A total of eight measurements on arable land were performed, four measurements on bare soil and four on a vegetated land surface. In total, two bare soil plots were kept bare through the season. One plot was cultivated before each measurement. A surface sealing caused by the rainfall was kept at the second plot. In total, two plots were vegetated with the crop seeded at the field plot (wheat and clover).
The intensity of the rainfall was set to 60 mm/hour on all plots. The duration of the measurement was 30 min from the beginning of surface runoff. The rainfall simulation measurements were used to calibrate the parameters of the model. Overall, three parameters were used: the measured discharge, precipitation, and the surface runoff velocity. Linear regression and the least squares method were used for calibration. The calibration was performed in two steps. In the first step, the water balance Equation (1) parameters were calibrated (infiltration processes, surface retention and effective precipitation). In the second step, the parameters of the flow Equation (3) of the surface runoff were calibrated. The two-step calibration was based on the following measurement procedure. Rainfall and surface runoff were directly measured during the rainfall simulations. The remaining water was a combination of infiltration and surface retention. The calibration thus corresponds to values measured in the following way. First, the infiltration and retention parameters, i.e., the parameters in the balance part of the model were calibrated. The key point is the beginning of the surface processes, which is another important calibration parameter. In the second step, the shape of the runoff curve was adjusted using parameters b, X, Y. Measurements on bare soil and on soil with vegetation were processed separately. The mean values and the standard deviations of the modelled parameters were determined for two variants: with vegetation and without vegetation.
The parameter validation was carried out using long-term monitoring of surface runoff data. For the validation phase, a natural rainfall event was used. A total of six rainfallrunoff events observed in 2010 and 2011 were selected. Rainfall data were collected from the tipping bucket in combination with adjusted rainfall data [52] in 10-min intervals. The events were modelled for three states (based on RS simulation): (i) MEAN-the mean value of the parameters derived from the rainfall simulations, (ii) the MEAN + STD the mean value of the parameter derived from the rainfall simulation increased by the standard deviation of the derived parameters, and (iii) the MEAN-STD the mean value of the parameter derived from the rainfall simulation decreased by the standard deviation of the derived parameters.
Additionally, a Monte Carlo method was used to perform a simple sensitivity analysis based on the range obtained during calibration of the model. A total of 250 simulations were run for each bare soil plot and for each vegetated plot.

Field-Scale Rill Flow Validation-A Case Study
To verify the capability of the SMODERP2D model to estimate rill flows, the calculated rills were compared with rills that occurred at the study site under natural rainfall. In this validation procedure, the calculations were performed over the whole area of the study site. Aerial photographs of the study area taken after heavy rainstorm events were used to visually derive the pattern of erosion rills and ephemeral gullies, and these photographs were used for the comparison with the modelled rills.
Aerial photogrammetric data with spatial resolution of 0.25 m were used to validate the rill processes. The erosion damage published in [53] and additional more recent data were collected, and the positions of rills were visually determined. The detected rills were divided into significant rills (SR) and noticeable rills (NR), according to their width. SRs were rills of width greater than the spatial resolution. NRs were less wide than the spatial resolution. The axes of the rills were vectorized and were converted to a raster in the same resolution of the DTM [48] as was used in the SMODERP2D model (resolution 5 × 5 m). Pixel values 5, 9, and 20 were assigned to SR, NR, and to cells without rills for subsequent comparison with the modelled rills.
Casual rainfalls were selected for each of the heavy erosion events described in the previous paragraph, using adjusted precipitation radar data [52]. Then, all of the events were simulated with the SMODERP2D model with spatial resolution of 5 × 5 m. In total, two variants of the land-use properties were used: (i) vegetation cover and (ii) bare soil. The spatial distribution of the soil hydrological properties from the soil map was used (Figure 3b). The events were modelled in three states (MEAN, MEAN + STD, and MEAN − STD) The formation of erosion rills is controlled by the critical water level excess calculated for each pixel, using Equation (7). For the first evaluation of this parameter, the model was run without the rill runoff, and the maximum water levels and the maximum tangential stress were stored for further analysis. The critical water level was then verified using the visually inspected rills and the output of the model. Subsequently, the model was run, including the rill flow processes, to calculate the maximum rill width in each pixel, which was used for the final comparison with the visually identified rills.
The rills identified by the model were also converted to SR and NR classes according to the modelled width of the rills. The thresholds of the modelled rills were the same as for the visually detected rills from the airborne survey. The reclassified values were 1, 2, and 0 for SR, for NR, and for cells without modelled rills.
A comparison of the modelled and observed rills was then performed using the ratio of the number of pixels (see Table 1), in which the occurrence of measured erosion rills and modelled erosion rills coincided. If modelled rills were observed in more pixels than the measured rills, the model overestimates the occurrence of rills. If the model modelled the occurrence of SRs where there are real NRs, then the model makes a slight overestimation. Conversely, where NRs are measured and SRs are modelled, the model slightly underestimates the occurrence of rills. The areas where rills are not measured and are also not modelled can be considered as a good prediction by the model. The pixel-by-pixel comparison method is shown in Table 1. Figure 4 presents the results of the calibration procedure. In total, four measurements on bare soil and four measurements on a vegetated surface were used. All measurements were carried out with rainfall intensity of ca 60 mm/h ( Table 2).  The bare soil models exhibited better agreement with the measured data. The beginning of the runoff, the increasing limb of the hydrographs, and the approach to a steady state were well simulated in the bare soil plot. The model performed poorly only in the case of a very gradual increase in surface runoff (Figure 4). The model performed less well for the vegetated plots. The beginning of the runoff, and also the increasing limb of the hydrographs, were well simulated. However, the later stages of surface runoff were underestimated. The Nash-Sutcliffe efficiency [54] provided satisfactory values, especially for the bare soil plots.

Model Calibration and Sensitivity
Parameter sets were derived for bare surfaces and for vegetated soil surfaces ( Table  2). The main differences between the bare soil surfaces and the vegetated soil surface were in the infiltration equation parameters, where the bare soil exhibited lower saturated hydraulic conductivity and sorptivity values. The parameters of the surface flow equation (b, X, Y) (Equation (3)) were derived from the soil texture characteristics, which was the same on both plots. Therefore, b and Y parameters were kept fixed to find the best solution. The X parameter was derived for soil texture and partly surface microtopography. The best solution had to be compensated during one of the bare soil models run with the X parameter. Changing of this parameter did not provide a better solution on the area with vegetation. STD in Table 2 marked with * shows an unchanged parameter. The parameters of the surface flow equation were unchanged, since these parameters are derived from the soil texture, which was the same on both plots. The vegetation parameters are missing at the bare soil plot. Although rainfall with the same intensity was applied, the duration was much longer for the vegetated plot due to the longer surface runoff lag time.
The results of the sensitivity analysis are shown in Figure 5. The MEAN − STD to MEAN + STD parameter ranges shown in Table 2 were used as the margins for the Monte Carlo simulations. Parameters without STD were set to a fixed level. As expected, the bare soil model parameters exhibited less uncertainty. Moreover, the runoff lag time varied less at the bare soil plots. The quantile margins were also narrower for the later stages of runoff for the bare soil model. The greatest variation of the bare soil plot was observed in the rising limb of the hydrograph. In the case of the vegetated soil model, the most distinct feature the variable runoff lag time.

Model Validation-Plot Scale Sheet Flow
The model was validated using the long-term monitoring Wischmmeier plots. In total, six surface runoff events with various surface covers, as shown in Figure 6, were selected for validation. The validation showed that the model fitted the measured data well during less intensive events. In particular, however, two high rainfall intensity events on 16 June 2011 were not represented well by the model, most likely due to misrepresentation of the initial conditions or the infiltration parameters. On the other hand, both the vegetated model and the bare soil model represented the measured data well during the rest of the episodes. Interestingly, the best model was the MEAN-STD parameters model ( Table 3). Moreover, a complex event on 7 August 2010 was simulated relatively successfully considering the simplicity of the model.

Model Validation-Field Scale Rill Flow
To compare the modelled occurrence and the position of the eroded rills with the measured data, 96 model runs were computed for rainfalls with rill detection in aerial photographs (basic information on the modelled rainfall in Table 4). The SMODERP2D model was run in the same way as was used for surface runoff validation in the three scenario parameter set variants: MEAN, MEAN + STD, and MEAN − STD (parameter values in Table 2). This variation around the mean parameter values corresponds to a dryer/wetter, rougher/less rough state of the soil. In first stage, the model runs were computed without a rill process. The critical tangential stress was not exceeded in the modelled area, and only the velocity limit for rill formation was exceeded. It was visually verified that the achieved velocity exceeded the expected threshold in a considerable part of the modelled area.
The rill process was implemented in the second stage of the modelling procedure. Manually assigned values of the critical velocity and of the critical tangential stress were tested in this stage. The evaluation of the occurrence of rills in the model and visual identification was performed by the pixel-by-pixel method. An overview is presented in Figure 7 in the form of bagplots.
Only 5% of the pixels contained measured rills. The differences in the correct identification of erosion rills are therefore less apparent in Figure 7. The modelled runs overestimated the occurrence of erosion rills in the case bare soil plots, which corresponds to the fact that on the unprotected surfaces the possibilities of rills occurrence is higher.
The results are visualized for one simulation, in the form of maps showing both measured rills and modelled rills. Figure 8a shows the measured situation (identification of the rills, vectorization and raster classification), and Figure 8b,c shows the modelled data.

Discussion
This paper has shown the use of a combination of calibration and validation methods of a rainfall-runoff erosion model on a specific agricultural plot. A physically based approach implemented to the SMODERP2D model for the solution of event-based surface processes was used. The concept of an event-based model neglects other components influencing the long-term water balance, such as evaporation and transpiration. The model uses the Philip infiltration equation, which is routinely applied in models and has available input parameters. For example, [55] compares the influence of infiltration methods. A disadvantage of the Philip equation is that it is not capable of modelling an interrupted rainfall.

Plot-Scale Sheet Flow Calibration and Validation
Episode models are sensitive to the initial condition settings. Variability of the initial conditions is also evident from the rainfall simulations and long-term measurements on the field plots that were used for calibrating the model. The variability of the experiments (four for bare soil and four for vegetated surfaces) with same precipitation clearly shows the variable runoff lag time of the more variable measurements with vegetated surfaces. The maximum differences at the beginning of surface runoff between the measurements is 6 min in the case of bare soil, and 75 min in the case of a vegetated surface. The surface runoff after initialization did not differ significantly among experiments with the same type of surface.
The MEAN − STD scenario exhibited the best correlation between the measured surface runoff and the modelled surface runoff during the plot scale validation. The MEAN and MEAN + STD variants did not produce any surface runoff in some cases. Different sizes of the calibration plots (16 m 2 ) and the validation plots (50.2 m 2 ) may have caused the MEAN − STD best fit, since a longer plot is more prone to surface runoff formation, which may be underestimated by the parameters calibrated for a shorter plot. Only in the case of simulation 7 August 2010 (validation measuring) the NSE exhibited below-zero values for all modelled scenarios. Statistical methods of comparing hydrographs, especially of multiple peak runoff, are not always appropriate description and is better to evaluate by an engineering approach [56].
The model performs best on bare soil, where the infiltration routine is easier to describe. Infiltration into a vegetated soil is affected by the more complex soil structure due to biomacropores. The root zone depth, distribution and macropore connectivity is not included in the model structures, and this causes large differences between the modelled values and the measured values. Moreover, in the bare soil model, the infiltration decreases rapidly due to sealing of the soil surface caused by the impact of the raindrops. Since the soil surface is protected by the vegetation on the vegetated plot, the infiltration rate decreases more slowly and, together with the vegetation interception, causes the greater difference in the runoff lag time as it was shown during the sensitivity analysis. The infiltration routines were also shown to cause differences in a comparison between the SMODERP2D model and the Erosion3D model [17].

Field-Scale Rill Flow Validation-A Case Study
A comparison between monitored and modelled erosion rills on the study area was used for verifying the rill model routines. The comparison was made on a pixel-by-pixel basis. However, this method sets very strict criteria for the position of the rill on the soil surface. Three factors affect the results: (i) The accuracy of the DTM and the method for calculating the direction of the outflow. The monitored rill and the modelled rill may be in close proximity to each other, but may not fall into same pixel category. This effect is shown in Figure 9. It is clear that the DTM quality and the outflow routing is important for a description of a real situation or for mutual rill formation from DTMs [57]. (ii) The angle of the rill flow direction and the slope direction may differ due to the microtopography during the initial stages of rill development. This may also cause a difference between a modelled rill and a measured rill, although this is less pronounced for the later stages of rill network development [58,59]. It has also been shown that DTM in fine resolution can predict rill development well [60].
(iii) The calibration of the computational parameters of the model. The surface condition and the infiltration capacity have the main influence in the model. Sheet flow processes are calibrated in this study with data from RS. Bare soil significantly overestimates erosion rills in comparison with monitored rills. Model parameters such as infiltration and roughness affect the formation of erosion rills more than the exact calibration of the critical water level for the formation of erosion rills (see Figure 10b). The MEAN modelled situation with vegetation cover has the best agreement with measured data (Figure 10). The difference between the runoff and rill generation may also be affected by the soil data. More detailed soil properties could produce different results. The soil properties in the model were based on rainfall simulation and the use of the latest spatial soil data, that are based on digital mapping of the soil properties [61]. The area of the rills predicted by the model was 1838 m 2 for noticeable rills and 2351 m 2 for significant rills for the solution that best corresponded to the monitored rills. The rills from the corresponding aerial photographs have a total area of 6648 m 2 [53]. The boundaries of significant rills are more visible from the remote sensing data than noticeable rills, and are therefore more accurately vectorized. The width of a rill is in some places greater than the width of a pixel. In the SMODERP2D model, the width of the rill is fixed in proportion to its depth and volume (Equations (12) and (14)). This ratio was shown to be in agreement with other studies [62].

Conclusions
The validation and calibration procedure for the SMODERP2D rainfall-runoff erosion model was based on data obtained at three spatial scales. The naturally generated runoff validation showed that simulated rainfall measurements are a suitable method for surface runoff model calibration. Although it is complicated to perform, rainfall simulation is still an easier tool than long-term measurements on erosion plots for obtaining surface runoff data. This was also proven for plots of slightly different scales.
Field plot scale data were used for validating the rill calculation routine of the SMO-DERP2D model by a visual and pixel-by-pixel comparison between modelled rills and monitored rills created by natural heavy rainfalls. Some mismatch was observed between the modelled rills and the measured rills, but this mismatch was most likely caused by the displacement of the modelled rills due to the flow algorithm. Nevertheless, the concept was proven. The main goal of soil protection is to reduce soil loss, and rill detection is the first step on this path. A reliable modeling tool for predicting the development of rill networks can also be used to test new agricultural and technical soil loss reduction practices and measures at regional scale where there are multiple farmers operating the land.
More detailed data obtained, e.g., by the Lidar method or by the SfM method, can be used in the future to better identify the spatial distribution of the rill generation process. Data from different land-use and meteorological conditions can also be used to further investigate the transferability of the presented approach to other locations.