Comparison of the Applicability of Different Soil Erosion Models to Predict Soil Erodibility Factor and Event Soil Losses on Loess Slopes in Hungary

Climate change induces more extreme precipitation events, which increase the amount of soil loss. There are continuous requests from the decision-makers in the European Union to provide data on soil loss; the question is, which ones should we use? The paper presents the results of USLE (Universal Soil Loss Equation), RUSLE (Revised USLE), USLE-M (USLE-Modified) and EPIC (Erosion-Productivity Impact Calculator) modelling, based on rainfall simulations performed in the Koppány Valley, Hungary. Soil losses were measured during low-, moderate- and high-intensity rainfalls on cultivated soils formed on loess. The soil erodibility values were calculated by the equations of the applied soil erosion models and ranged from 0.0028 to 0.0087 t ha h ha−1 MJ−1 mm−1 for the USLE-related models. EPIC produced larger values. The coefficient of determination resulted in an acceptable correlation between the measured and calculated values only in the case of USLE-M. Based on other statistical indicators (e.g., NSEI, RMSE, PBIAS and relative error), RUSLE, USLE and USLE-M resulted in the best performance. Overall, regardless of being non-physically based models, USLE-type models seem to produce accurate soil erodibility values, thus modelling outputs.


Introduction
Water erosion is a worldwide problem that causes several economic and environmental impacts [1] due to deteriorating soil functions [2,3]. The negative impacts can occur on-site by reduced production and off-site by sediment transportation. Sediments can deliver pollutants into surface waters or onto other lands, threatening sustainable land use [4]. As Miller et al. [5] announced, nearly 90 percent of the total soil loss of nitrogen and phosphorus [6] can be delivered with the soil loss and later accumulated within the sedimented area. Recognizing the severity of accelerated soil erosion processes, a soil loss tolerance value was introduced [7,8] in the USA for soil conservation planning in the 1940s. Later, Mannering [9] and Skidmore [10] suggested the involvement of establishing tolerable rates into the regulation of farming. Morgan [11] defined the tolerable amount of soil loss as an amount that does not considerably decrease soil fertility; thus it is a key factor of sustainable soil management. The concept of tolerable soil loss was applied in soil conservation planning, as well as soil loss evaluation and mapping. The annual and long-term average soil loss can serve as the basis of soil erosion predictions and nonpoint source pollution control. With climate change, the occurrence of rare but more erosive rainfall events is going to increase [12]; therefore, studying the effect of these events on soil degradation is becoming more essential and crucial [13], because these events cause the larger part of annual soil erosion. Therefore, many scientists [14][15][16][17] emphasize that the conservation strategies should take into account large storms rather than average weather conditions.
With rising competition for the limited soil resources, soil erosion prediction has become an essential part of soil conservation planning and practices, as there is a need for simple and cost-effective soil erosion-and soil erodibility measurements [18][19][20][21][22]. There are several types of soil erosion models; the new ones require a lot of input parameters. This can be the reason why USLE (Universal Soil Loss Equation) [8] is still one of the most often used models. Compared to other models, USLE is still the most cited of all soil erosion models. USLE users do not need to calculate event rainfalls for the calculation of a yearly average [23][24][25]. Still, it has been recognised that there is a need to predict soil losses on a shorter time scale, which resulted in many studies with rainfall simulation measurements. Some studies [26][27][28][29] warn us that these simulations reflect single rainfall events which over-or underestimate average soil loss [30][31][32][33][34]. This phenomenon is common when the EI30 index is used in estimating erosion on bare fallow areas [35,36]. Boardman [37] stated that the application of USLE can assume unconcern in event-driven erosion, whereas it is unsatisfactory for estimating soil loss at event scale. Foster et al. [38] recommended that the erosivity factor, including rainfall amount, rainfall intensity and runoff amount, should be better than the single-storm erosion index (EI30) by Wischmeier and Smith [8]. On the other hand, Kinnel et al. [39] stated that USLE and similar models, such as RUSLE, can predict the event soil loss for runoff-producing events when we assume that the event soil loss is directly related to the product of event storm kinetic energy (E) and the maximum 30 min rainfall intensity (I30) [40]. This way, both USLE and RUSLE can model short-term soil losses. To provide better predictions, several similar but modified models were built, such as MUSLE [41], RUSLE [42], WEPP [43,44], USLE-M [45], etc.; while USLE and RUSLE are detachment limited models that hypothesise that soil loss is not affected by the runoff capacity of detached soil particles, USLE-M studies the runoff, operating, in this way, as a transport-limited model.
In practical soil conservation measures, the most important element of soil erosion modelling is soil erodibility or "K-factor", as it became known by the work of Wischmeier and Smith [8]. This factor is defined as the susceptibility of a soil to erode, representing the effect of soil properties on soil loss and the resistance against surface water runoff [42]. Unfortunately, the utilization of soil erodibility is not consistent in the literature [22,46,47], as it follows many inaccuracies and uncertainties within the soil erosion modelling soil erodibility determinations [48][49][50][51][52][53][54] and soil loss estimations [55]. Cassol et al. [56] also drew attention to the fact that the USLE K-factors can differ from annual estimation values in the case they originate from individual rainstorm events. Although there is a wellknown standard procedure for the USLE model [57] (22.13 m long with 9% uniform slope, natural and continuously tilled fallow field plots, long-term measurements), there are still several studies with different circumstances trying to simplify the observations [49]. Many pieces of literature exist, providing non-satisfying but detailed determination procedures; as a result, some uncertainty concerning the accuracy of the reported K-values exists. This means that different estimation methods behave differently and their results are relevant just under the given circumstances. The K-factor has high variability and mostly depends on the method of the experimental setup, the main soil properties and various and numerous input parameters [49]. Besides, in event-based, mostly rain-simulated erosional studies, erodibility is regarded as a dynamic variable due to the main exogenic forces [54], especially the rainfall and erosional processes driven by rain splash. To serve the purpose of simplifying the long-term measurements based on physical and chemical properties, several nomographs were created in the models (USLE, RUSLE and EPIC). At the same time, El-Swaify and Dangler [58] reported that erodibility estimation is quite risky, based on soil classification only. Based on the modelling research studies performed so far, we can conclude that there is a huge need for soil erodibility measurements [59].
Accordingly, the problem of how the soil erodibility values measured by event rainfall simulation can be comparable arises. However, to date, no comparison has been made relatively to the ability of the USLE, RUSLE and USLE-M models, nor of USLE and EPIC nomographs, in soil erodibility determination and soil loss estimations produced by individual high-intensity simulated rainfall events on different slope gradients and under different soil moisture conditions on soils formed on loess parent rock in the Koppány River Valley. Therefore, the specific objectives of this paper are the following: (a) to statistically test and compare the erodibility factors calculated by the USLE, RUSLE and USLE-M [35,36,45] models and the calculated values by the USLE and EPIC nomograph equations [60] based solely on soil physical and chemical properties; (b) to determine the difference between each predicted and measured soil loss based on high-intensity rainfall simulations; (c) to determine how accurate the event rainfall simulation results are, compared to other similar long term estimations; and (d) to establish how the different antecedent soil moisture content can affect the erodibility values and soil loss estimation results.

Study Site and Soil Properties
The concerned area is situated in the north-eastern part of the South-Transdanubian Region, on the east side of Somogy county, in Hungary. Gerézdpuszta is situated in the Koppány valley, adjacent to the floodplain of the Koppány River ( Figure 1). This landscape is characterized by loess-covered asymmetric hilly areas, where the northern hillsides are short and steep, while the southern areas are longer and have gentle slopes [61]. Due to loess-like deposits, the area is prone to soil erosion; however, most of the arable lands are under conventional tillage where large-scale farming is typical. The less eroded hillslopes are characterized by Cambisols, while, in cultivated areas, Regosols and Colluviums can be found [62]. The climate is moderately warm and wet. The annual average temperature is between 10.0 • C and 10.2 • C and the average annual precipitation is between 605 mm and 700 mm. Most of the hillsides are characterized by agricultural cultivation and almost half of that is situated on slopes steeper than 12%; further, usually the farmers do not use any soil protection measures.
Based on the modelling research studies performed so far, we can conclude that there is a huge need for soil erodibility measurements [59].
Accordingly, the problem of how the soil erodibility values measured by event rainfall simulation can be comparable arises. However, to date, no comparison has been made relatively to the ability of the USLE, RUSLE and USLE-M models, nor of USLE and EPIC nomographs, in soil erodibility determination and soil loss estimations produced by individual high-intensity simulated rainfall events on different slope gradients and under different soil moisture conditions on soils formed on loess parent rock in the Koppány River Valley. Therefore, the specific objectives of this paper are the following: (a) to statistically test and compare the erodibility factors calculated by the USLE, RUSLE and USLE-M [35,36,45] models and the calculated values by the USLE and EPIC nomograph equations [60] based solely on soil physical and chemical properties; (b) to determine the difference between each predicted and measured soil loss based on high-intensity rainfall simulations; (c) to determine how accurate the event rainfall simulation results are, compared to other similar long term estimations; and (d) to establish how the different antecedent soil moisture content can affect the erodibility values and soil loss estimation results.

Study Site and Soil Properties
The concerned area is situated in the north-eastern part of the South-Transdanubian Region, on the east side of Somogy county, in Hungary. Gerézdpuszta is situated in the Koppány valley, adjacent to the floodplain of the Koppány River ( Figure 1). This landscape is characterized by loess-covered asymmetric hilly areas, where the northern hillsides are short and steep, while the southern areas are longer and have gentle slopes [61]. Due to loess-like deposits, the area is prone to soil erosion; however, most of the arable lands are under conventional tillage where large-scale farming is typical. The less eroded hillslopes are characterized by Cambisols, while, in cultivated areas, Regosols and Colluviums can be found [62]. The climate is moderately warm and wet. The annual average temperature is between 10.0 °C and 10.2 °C and the average annual precipitation is between 605 mm and 700 mm. Most of the hillsides are characterized by agricultural cultivation and almost half of that is situated on slopes steeper than 12%; further, usually the farmers do not use any soil protection measures.  Table 1. The texture is sandy clay loam. The gentle slope section has low CaCO 3 content, while the steep slope has high CaCO 3 content; the soil organic carbon content is low in both sections. These characteristics indicate that the area has been severely prone to soil water erosion for a long time.
Plasticity indicates the amount of distilled water needed to reach the state of plasticity (or saturation). Soils are considered loam with values from 37 to 42 and clayey loam with values from 43 to 49 (accuracy is ±2). The pH (H 2 O) was measured with a WTW inolab pH7310 pH measuring device. SOC was measured with the wet combustion method. The particle sizes were measured with a Fritsch Laser analyser. Soil permeability was measured with a double-ring infiltrometer in the field. Soil structure was evaluated by experts on the field before the rainfall simulations.

Rainfall Simulator and Experimental Design
The Shower Power (SP) 02 is a portable field unit applicable on a plot scale with alternating nozzles which were designed and constructed by the Geographical Institute, Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences, in 2014. It is a portable field scale device with a measuring plot of 3 × 2 m; however-to exclude the border effects-the irrigated plot size was 12 m 2 . The drop-forming unit was an alternating axis equipped with two 80,100 Veejet nozzles placed next to each other with a distance of 2 m and 3 m in height. This way the overlapping between the nozzles was perfect, resulting in a uniform spatial intensity and drop-size distribution. The intensity can vary in the range 30-100 mm h −1 , depending on the axis alternation frequency. Using the V-jet nozzles, the rainfall energy was 77% of natural rainfall when the pressure and rainfall intensity of the simulator was constantly held at 0.41 KPa and 65 mm/h intensity rate [63]. The average kinetic energy of the simulator was 24 J m −2 mm −1 . The plot was fenced using metal sheets pushed into the soil to a depth of 15 cm. At the bottom of the plot, a pair of collector triangles gathered the runoff. During the rainfall simulation, the starting and ending time of simulation, ponding and runoff time were recorded. The volume of runoff was continuously measured by two scaled dishes as a function of time and the whole runoff was collected separately into buckets. Time and temporary runoff subsample values were measured and stored in controlling software. The subsample measurements and sediment concentration determination were made from the proportion of the total soil loss over the simulation. In the laboratory, the total amount of soil loss was oven-dried for 48 h at 60 • C and soil loss was measured afterward. Altogether, 18 simulations were conducted ( Table 2). Further details in Appendix A. Three rainfall simulation treatment rates were sequentially applied: (1) dry antecedent moisture content before the application of the low rainfall intensity rate (30 mm/h −1 ); (2) moist antecedent moisture after the low rainfall intensity rate (30 mm/h −1 ) and before the medium rainfall intensity (60 mm/h −1 ) (60 min after the previous treatment); and (3) wet antecedent moisture content after the medium rainfall intensity rate (60 mm/h −1 ) (after the ponds disappeared from the surface). In the first 8 experiments, the selected slope section was gentle, while the following 8 experiments were prepared on a steep slope. The real amount of rainfall was measured with the help of 0.5 L small buckets placed around the plot borders. The measured soil loss from the dry, moist and wet rainfall simulation treatments was used in this study. The total sum of the rainfall simulation experiments was denoted later as the pooled dataset.

Soil Erosion Modelling
The evaluation of various soil erosion models was performed by comparing the measured and predicted soil losses for each rainfall simulation for the applied slopes and rainfall intensities. Since rainfall simulations were based on a specific design of rainfall events, the models had to be used based on the single rainfall event inputs. The event-based predictions were possible using the procedures for calculating individual design storm rainfall erosivity values provided by the rainfall simulation measurements.
K-factor values were determined by five direct methods using a rainfall simulator and by two indirect methods (based on soil physical properties).
The basis of the erodibility, K-factor calculation is provided by one of the several estimation methods, the widely spread and well-known USLE equation by Wischmeier and Smith [8]. The Revised Universal Soil Loss Equation (RUSLE) [64] computes with the same formula, as follows: where A (t ha −1 y −1 ) is the annual soil loss, R (MJ mm ha −1 h −1 y −1 ) is the rainfall erosivity factor, K (t ha h ha −1 MJ −1 mm −1 ) is the soil erodibility factor, L is the slope length (dimensionless), S is the slope gradient factor (dimensionless), C is the cropping cover management factor (dimensionless) and P is the agricultural (or support) practice factor (dimensionless). This way the soil erodibility can be calculated as A, the annual soil loss (t × ha −1 × y −1 ), is originated from the measured runoff in the field on a 6 m 2 plot area.
R, the rainfall erosivity index, is based on the equation by Foster et al. [38], where E is the kinetic energy (MJ ha −1 mm −1 ) and I 30 is the rainfall intensity (mm h −1 ). Kinetic energy can be obtained from the following expression: LS, the slope length and gradient factor in USLE, is based on the equation by Presbitero [65], where λ is the slope length; θ is the slope angle ( • ); m is equivalent to 0.5 for s > 5%, 0.4 for 3% < s ≤ 5%, 0.3 for 1% < s ≤ 3% and 0.2 for s ≤ 1% slope gradient. m, the exponent in RUSLE, based on Foster et al. [38], is defined as where β is a function of slope, moderately susceptibility for rill and inter-rill erosion [42], i.e., β = (sin θ/0.0896)/3(sin θ) 0.8 + 0.56 (6) where λ is the slope length, m is a variable length-slope exponent, β is a factor that varies with slope gradient and θ is the slope angle ( • ) for slopes 1-30 • . S, in RUSLE, based on McCool et al. [42], where the slopes are shorter than 4.6 m, is USLE-M can be used to evaluate parameters from single runoff events, using a combination of storm erosivity and runoff ratio, with Q R EI 30 term [35,36,45,66] as the erosivity index.
A e = Q R EI 30 K UM LS C UM P UM (8) where Q R (dimensionless) is the event runoff coefficient; EI 30 is the event erosivity factor. which has changed from the original EI 30 to Q R EI 30 factor; K UM is the soil erodibility factor (t ha h ha −1 MJ −1 mm −1 ); L is the USLE slope-length factor; S is the USLE slope-steepness factor; C UM is the crop factor; and P UM is crop management and support practice factor. The subscript U represents the factor associated with the USLE which has to be modified (subscript M).
The EPIC model k-value estimation was performed with the following formula by Williams et al. [41]: where SAN refers to the content of sand (0.1 mm size < 2 mm) in %; SIL refers to the silt particle content (0.002 mm size < 0.1 mm) in %; CLA refers to the clay particle size (<0.002 mm) in %; C refers to the organic carbon content in %; and SN1 = −SAN/100. The unit of the soil erodibility k-value is t ha h ha −1 MJ −1 mm −1 .
The event soil loss and runoff values, erosivity indices and other input parameters of different soil loss estimation methods were calculated by adding the measured variables at the event scale. These latter gave the basis of K-factor determinations which were inserted into the selected models. Soil loss predictions in the case of different models were obtained by the substitution of the mean, measured soil-and model-specific K-factors.

Statistical Evaluation of the Model Performance
The Pearson correlation coefficient (r), coefficient of determination (r 2 ), relative error (RE) and Nash-Sutcliffe model efficiency index (NSEI) were adopted to evaluate the models' performance.
The Pearson correlation coefficient (r) and the coefficient of determination (r 2 ) describe the degree of collinearity of predicted and measured data. The correlation coefficient ranges between −1 and 1 and indicates the linear relationship between measured and calculated data. R 2 ranges from 0 to 1, where the higher values indicate less error variance and above 0.5 can be considered acceptable.
The NSEI determines the relative magnitude of the residual variance compared to the measured data variance. The NSEI indicates how well the plot of measured versus predicted data fits the 1:1 line. It is calculated with the following relationship (Equation (11)): where A e,measured is the measured soil loss, A e,calculated is the soil loss value calculated by the models and A e,mean is the mean of the measured values. When the model accounts for all the variation in the measurements (A e,calculated = A e,measured ), the index has a value of 1.0; NSEI ranges between −∞ and 1.0. Values between 0.0 and 1.0 mean acceptable levels of performance. Zero indicates that the model predictions are as accurate as of the mean of the measured values. A negative NSEI indicates that the measured mean is better than the predictions. Based on Ahmad et al.'s [68] study on erosion prediction evaluation, the model performance criteria is NSEI > 0.4. Some of the commonly used error indices were applied. The relative error (RE) and the root-mean-square error (RMSE) were calculated to indicate the error between calculated and measured models. In the case of these error indices, the value of 0 indicates a perfect fit.
The percent bias (PBIAS) measures the average tendency of the predicted data to be larger or smaller than their measured counterparts [69]. The optimal value is 0.0, with low-magnitude values indicating an accurate model simulation, while positive values indicate under-estimation of bias and negative values indicate over-estimation of bias (Equation (14)). The PBIAS can be considered very good if PBIAS < ± 15%; good, if 15% ≤ PBIAS < ±30%; satisfactory, if 30% ≤ PBIAS < ±55%; and unsatisfactory, if PBIAS ≥ ±55%.
The Mann-Whitney U test was used to determine whether the estimates from the models were significantly different from the measured data (p < 0.05 means significant difference). The applied software was GraphPad Prism.

Evaluation of the Erodibility Estimated by the Applied Models
The first evaluation of erodibility values was performed by comparing the calculated K-results based on the measured parameters by USLE, RUSLE and USLE-M and by the EPIC and USLE nomographs ( Figure 2). In the separated, calculated results, in the case of the USLE and RUSLE models, the wet runs resulted in higher K-values than the dry runs, while USLE-M indicated higher erodibility with dry runs. These results owe to the high variance within the measured hydrological parameters. There was a higher difference among different slope categories with USLE-M calculations resulting in higher values for steep slopes, while the USLE and RUSLE models showed lower variation. The pooled data indicate higher variability with a higher difference between the K-values calculated by each model. The K-factors obtained by the USLE nomograph, in most cases, were under the other calculated values, while the K-factors resulting from the EPIC nomograph were above the USLE-type models. In the case of all models, the variability within the measurements increased with the increase in soil moisture content and slope degree.
The mean values ranged from 0.0028 to 0.0087 for the USLE-type models (Table 3). Comparing the calculated values based on the measured rainfall simulation parameters with the results of soil texture-based nomographs, the EPIC model returned

Evaluation of Event Soil Loss Estimation
Soil erosion measured results confirmed that the higher rainfall intensities with higher soil moisture content (moist and wet runs) generated more runoff and soil loss during each event containing higher suspended sediment concentration (Table 4).

Evaluation of Event Soil Loss Estimation
Soil erosion measured results confirmed that the higher rainfall intensities with higher soil moisture content (moist and wet runs) generated more runoff and soil loss during each event containing higher suspended sediment concentration (Table 4). The highest rates of erosion were recorded in storms of 63.16-103. 48   Generally, for dry runs, the highest variability in soil loss estimation occurred in the case of USLE, followed by the EPIC nomograph estimation. The closest results were derived from USLE-M estimations. USLE, RUSLE and the nomographs overestimated the measured amounts. Based on moist runs, soil losses estimated by USLE, RUSLE and USLE-M were close to the measured ones, while the EPIC nomograph resulted in nearly threefold higher values (also with higher variability). The nomograph underestimated the measured soil loss. The wet rainfall simulation results indicate higher variance within the data which resulted in a higher discrepancy among the estimations of USLE and USLE-M, as well as the EPIC nomograph. The USLE nomograph underestimated the average soil loss, while EPIC overestimated it. Correspondence could be found among the USLE, RUSLE and USLE-M models. Separating the data by slope steepness, then, in the case of the gentle slope, high variability occurred among the results in the attendance of many outlier values. The estimation of the EPIC nomograph indicated almost three- Generally, for dry runs, the highest variability in soil loss estimation occurred in the case of USLE, followed by the EPIC nomograph estimation. The closest results were derived from USLE-M estimations. USLE, RUSLE and the nomographs overestimated the measured amounts. Based on moist runs, soil losses estimated by USLE, RUSLE and USLE-M were close to the measured ones, while the EPIC nomograph resulted in nearly threefold higher values (also with higher variability). The nomograph underestimated the measured soil loss. The wet rainfall simulation results indicate higher variance within the data which resulted in a higher discrepancy among the estimations of USLE and USLE-M, as well as the EPIC nomograph. The USLE nomograph underestimated the average soil loss, while EPIC overestimated it. Correspondence could be found among the USLE, RUSLE and USLE-M models. Separating the data by slope steepness, then, in the case of the gentle slope, high variability occurred among the results in the attendance of many outlier values. The estimation of the EPIC nomograph indicated almost threefold greater loss than the measured one. The nearest soil loss values were those of USLE and RUSLE.
In the case of the steep slope, the closest estimation was obtained by USLE-M, then RUSLE and USLE, while other models indicated two-and threefold amounts of soil loss. Based on all pooled data, the best estimations were obtained by RUSLE, USLE, then USLE-M, while the EPIC nomograph overestimated the soil loss against the measured losses and the USLE nomograph underestimated it. At the same time, soil loss estimations by RUSLE, USLE and the USLE nomograph had the lowest variability.
The coefficient of determination shows an acceptable correlation between the measured soil loss and the predicted ones by USLE-M (Table 5). Table 5. Soil loss calculations by different soil loss estimation models (A U -USLE; A RU -RUSLE; A U-M -USLE-M; A nomo -USLE nomograph; A EPIC -EPIC nomograph) and their Pearson correlation coefficient (r), coefficient of determination (r 2 ), Nash-Sutcliffe coefficient (NSEI), root-mean-square error (RMSE), relative error (RE), percent bias (PBIAS), relative difference (R diff ) and p-value (Mann-Whitney U test), used to test models' performance. The model efficiency of almost all models-except RUSLE-was negative in the simulations. The negative NSEI values indicated that the mean measured soil loss from the field was a better representation of soil erosion for the concerned simulations than the estimations of specific soil erosion models. Based on the NSEI, the RUSLE model resulted in the best prediction, followed by the USLE and USLE-M, while the nomograph EPIC estimation seemed the worst predictor. RMSE values indicated that the best model performance was obtained in the case of the RUSLE, USLE and USLE nomographs. Based on the relative error (RE), a satisfactory relationship occurred between the measured and predicted datasets in the case of RUSLE, USLE and the USLE-M, which is the third with a negative value owing to higher measured soil loss. PBIAS (PB) with negative values indicated over-estimation of bias almost in all cases, except the USLE nomograph, which under-estimated the observed soil loss (Figure 4) with satisfactory estimation (49%); in the case of USLE (−13%), RUSLE (−3%) and USLE-M (−50%), we obtained good and satisfactory estimations. The estimates from the models departed significantly from the measured data (p < 0.05) with the EPIC model. The accuracy of model predictions and the tendencies of the models to over-or underpredict soil losses are graphically visualized in Figure 4a-e. The 1:1 line represents the perfect match between the measured and predicted values. The data points fitting or positioned closer to the 1:1 line indicate accurate predictions, while the values situated above the 1:1 line reveal overestimation and the ones situated below refer to the data that underestimated the amount of soil losses. USLE and RUSLE overestimated soil loss by an average of 13.4% and 3.6%, respectively (Figure 5a,b), while the EPIC model showed the worse performance, with an average overestimation of 212.95% (Figure 5d,e). RUSLE model predictions were sensitive for the measured values, when this increased above 0.6 t/ha; then, the model predictions indicated underestimations. USLE-M and the USLE nomograph (Figure 5c,e) underestimated the measured values by an average rate of 17.7% and 49.7%, respectively. However, except for some outlier values, the USLE-M showed the best performance and fit among the existing models which were not sensitive to the changes within the measured values. Water 2021, 13, x FOR PEER REVIEW 12 of 20  ured values, when this increased above 0.6 t/ha; then, the model predictions indicated underestimations. USLE-M and the USLE nomograph (Figure 5c,e) underestimated the measured values by an average rate of 17.7% and 49.7%, respectively. However, except for some outlier values, the USLE-M showed the best performance and fit among the existing models which were not sensitive to the changes within the measured values. Moreover, the residuals, the differences between measured and predicted values, of the data were investigated ( Figure 5). In the case of RUSLE and the USLE nomograph, high-density of points are close to the origin. In the case of USLE and USLE-M, the scattering from the origin increases but these are still independent and normally distributed, while, for the EPIC nomograph, the residuals and their deviation increase with the increase in measured soil loss; therefore, these are not independent.

Discussion
USLE and RUSLE resulted in higher soil erodibility values in the case of wet runs, while, surprisingly, USLE-M resulted in higher K-values for dry runs. This is not consistent with some studies where increased soil moisture content resulted in decreased sediment concentration, thus decreased soil erodibility values [59,[70][71][72]. At the same time, we can conclude that the variability of the calculated erodibility values increased in parallel with the increase in soil moisture content and slope degree.
The effect of slope gradient on soil loss and soil erodibility is well known. In this study, USLE-M indicated higher erodibility in steep slopes with higher variation than USLE-and RUSLE-based erodibilities. Kinnel et al. [3] reported, in their study, that values of both USLE-M and RUSLE varied with the slope gradient, which is parallel with our findings. In the case of gentle slopes, the event soil loss estimation's results showed high variance, but the best estimations were obtained by RUSLE and USLE-M for both gentle and steep slopes. These results indicate that, in the case of event soil loss measurements, many influencing factors should be taken into account in order to conclude with a correct interpretation.
The mean of the calculated soil erodibility values obtained by different soil erosion models ranged from 0.0028 to 0.0087 for USLE-type models. The previously reported K-factor, according to Regosol in Hungarian literature (derived from in situ rainfall simulation experiments), is 0. [73], which is significantly higher than our calculated K-factors. This difference could be derived from the difference in physical properties of the investigated Regosol, from the different initial soil moisture content and the applied replicates. The different rainfall duration and intensities were eliminated during the calculations. Generally, the EPIC nomograph estimated larger erodibility factors, while USLE estimated smaller ones. The observation of the EPIC nomograph's overestimation is in agreement with the findings by Zhang et al. [74] and Wang et al. [59], who reported gross overprediction of soil loss in the case of nomograph estimators.
The event soil loss evaluation confirmed that the initial soil moisture content and rainfall intensity has great importance in runoff and soil loss generation [75,76]; therefore, it is important to emphasize this parameter when we publish results and draw conclusions from these results, regarding the accuracy of models. For the dry run, when the amount of soil loss and the rainfall intensity was the lowest, USLE-M obtained the best performance. In the case of moist antecedent soil moisture content, where the soil loss and intensity increased, the USLE-type models predicted soil loss rates close to the measured ones. When the antecedent soil moisture content and intensity were higher, the estimated soil loss values were more variable, where the best performance was again obtained by the USLE-type models.
Based on the whole dataset and the statistical evaluation, the coefficient of determination resulted in an acceptable correlation between the measured and predicted values only in the case of USLE-M. Based on other statistical indicators, such as NSEI, RMSE, PBIAS and relative error, RUSLE and USLE, followed by USLE-M, resulted in the best performance. PBIAS indicated overestimation of bias in almost all cases, except the USLE nomograph. Similar results could be inferred by examining the graphical visualisation of predicted and measured values fitted to the 1:1 line. The lowest overestimation occurred with the application of RUSLE (3.6%) and USLE (13.4%), while a slight underestimation was obtained in the case of USLE-M (17.7%). At the same time, RUSLE indicated sensitivity in prediction, whereas, above 0.6 t/ha, event soil loss underestimation occurred in the estimations, an observation similar to the findings by Kinnel et al. [39] and Nearing [77]. While the USLE nomograph resulted in constant underestimation, EPIC overestimated the soil losses. Risse et al. [29], Tiwari et al. [78] and Kinnel et al. [36] reported that USLE overestimates the small annual and event soil losses while underestimating the high annual and event soil losses. In this study, these findings do not occur, except for a general overprediction of soil losses. Di Stefano et al. [79] found that both USLE and USLE-M tended to overestimate low event soil losses, while, in our study, the opposite was observed. These observations are partly due to the size of the plot, which was a unit in their study, and the values of the measured soil losses originated from natural rainfall and not from simulated precipitation. Based on the residuals, the ranking is similar to the above-mentioned results concerning the USLE-type models. Kinnel et al. [39] emphasized that USLE and RUSLE can predict event soil losses, while USLE-M can predict event soil losses better than RUSLE. These statements are just in partial agreement with our findings, whereas the capability of USLE and RUSLE to predict event soil losses is confirmed; however, USLE-M did not perform better than USLE or RUSLE. The nomographs are mostly based on solely soil physical parameters; therefore, in this study, we confirmed that they are unsatisfactory for this kind of prediction and, applied alone, cannot be used to estimate individual soil loss events [58]. It seems that the more rainfall simulation data are analysed, the more we know about soil erodibility; however, due to the complex nature of soil, the data gathered and analysed so far did not provide enough spatial coverage and in-depth knowledge to cover all needs. Similar research studies have been conducted in the past decades, dealing with both estimations and evaluations of the K-factor under natural and artificial rainfall conditions [80]. Furthermore, similar experiments have been conducted in different settings, e.g., tropical climate [81], that also supported the idea that more experiments and evaluation are needed.

Conclusions
Using the measured event soil losses from six runoff plots under different antecedent soil moisture content and rainfall intensities in loess covered Koppány River Valley, this paper presents a set of measured K-values and soil loss data by a set of soil erosion modelling methods, namely, USLE, RUSLE and USLE-M, as well as USLE and EPIC nomographs. The differences in the measured and predicted soil losses were evaluated in the case of different initial moisture content, rainfall intensities and slope gradients. The results can serve as values for event rainfall simulations for areas without natural runoff plots. The main conclusions are as follows: (1) The K-values in Hungary for the loess-covered Koppány River Valley ranged from 0.0028 to 0.087 t ha h ha −1 MJ −1 mm −1 . The calculated erodibility results vary with initial soil moisture content and slope degree. Regarding the pooled data, the CV% indicated less relative variability in the case of USLE and RUSLE, followed by USLE-M. (2) In respect to the event soil loss estimation, the results indicate high variance with the increase in initial soil moisture content and slope steepness. Based on the pooled data, the best estimations were obtained by RUSLE, USLE and USLE-M. (3) Both soil erodibility and soil loss estimation draw attention to the importance of the effect of soil moisture conditions. (4) Investigating the measured and estimated soil losses statistically, we found that RUSLE resulted in the best performance in event soil loss estimation. At the same time, RUSLE indicated sensitivity for the measured values. (5) Within these experiments, the USLE-type models' capability for estimating event soil losses was confirmed. At the same time, we cannot emphasize enough that the non-standardized circumstances involving many affecting factors could result in high variance among the results. (6) This study points out the high variability of soil erodibility factors and soil loss estimation values in the case of event rainfall simulations. Besides, it draws attention to the importance of further studies on this topic, because rainfall simulation is necessary for the determination of each soil type's erodibility and short or long time soil loss estimation. That is why we need to develop an accurate, standardized method for event soil loss measurements, where the observable parameters are not solely the type of soil and other soil properties because these methods cannot give satisfactory results. (7) This study focused on soils formed on loess parent material. These soils cover a wide range of geographical regions. Thus, this study can provide useful information for many areas that struggle with similar water erosion problems, suggesting using erosion modelling as a tool for finding a solution.
Funding: This research received no external funding.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Acknowledgments:
The authors are grateful to the staff of Geographical Institute, Research Centre for Astronomy and Earth Sciences, Hungarian Academy of Sciences, for providing the rainfall simulator equipment for our measurements and to the laboratory of ELTE University for enabling the laboratory measurements. In addition, we thank all the colleagues who took part in the field measurements: Pálma Anna Varga-Kertész, Balázs Bolf, Bernadett Bolf, Igor Dekemati, Zsófia Dobó, Marianna Ringer, Norbert Keller. Besides, we are grateful to the field's owner, who contributed to the in situ measurements and water supply.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
Appendix A