## Introduction

Soil erosion on cultivated lands has received much concern since it is considered to be one of the most critical forms of degradation.

Determining the removal of soil material and the deterioration of the soil system, soil erosion directly affects the quality of the soil, its agricultural productivity and its biological diversity. Moreover, accelerated erosion implies loss of water holding capacity of soils and turbidity issues through increased sediments in water. A European Commission analysis indicates that soil erosion rates continue to be more than soil formation rates across the European Union, but that the European agricultural policy is working to reduce this gap. Soil and sediments lost by water erosion in Europe are responsible for an estimated economic loss of about $20 billion per year, based on a restoration cost of $20 per tonne (Panagos *et al*., 2015).

Every erosion model must represent how climate, soil, topography and land use affect soil loss and related variables (Toy *et al*., 2002).

Many mathematical models have been developed to estimate plot soil erosion at different temporal scales. At present, attempts to improve empirical soil loss equations (Renard *et al*., 1997; Kinnell and Risse, 1998) are still carried out since these approaches continue to be attractive from a practical point of view (Cao *et al*., 2015) although this is not agreed by all scientists working on soil erosion. According to Boardman (2006), the continued use of the unmodified Universal Soil Loss Equation (USLE) (Wischmeier and Smith, 1978) could encourage a lack of interest in event driven erosion. On the other hand, Nearing (2013) recently stated that, at least for some purposes for which an erosion model is used, the USLE has been, and still can be, very successful. Moreover, the USLE, the revised USLE (RUSLE) (Renard *et al*., 1997) and the Water Erosion Prediction Project (WEPP) (Flanagan and Nearing, 1995) or other process-oriented models have to be considered as constituting a complementary suite of models to be chosen to meet the specific user need (Nearing, 2013). An implication of this reasoning is that improving empirical prediction of soil erosion continues to be sound from a scientific point of view. Obviously, empirical approaches must be physically plausible, which implies that empirical representation of the erosion phenomena should be in line with the physics of the erosion processes at the scales being investigated.

Complex erosion models, which can be solved by modern personal computers, typically require collecting spatially distributed and sometimes difficult to obtain input data. The performance of the model in terms of output quality does not always increase with its complexity, completeness and physical soundness.

The simplest process-based erosion model is divided into two components: one representing the interrill erosion processes and one to represent rill erosion. The main reason to use this scheme is that the soil erosion controlling variables affect interrill and rill erosion differently. As an example, slope length has no effect on interrill erosion whereas rill erosion increases with slope-length (Toy *et al*., 2002).

Any model should be validated using measured soil loss values and validation can be carried out fitting all measured values with the same error or to give priority to fitting the large values with minimal error. Erosion models used for soil conservation planning strategies should be most accurate in the range of 1-20 t ha^{–1} for obtaining a sufficiently accurate soil loss estimate useful to justify to a land user the need for reducing soil loss (Di Stefano and Ferro, 2016). Conversely, erosion models used for designing reservoirs should accurately estimate large values.

Fields plots are often used to obtain experimental data (soil loss values corresponding to different climate, soil, topography, crop and management conditions) for predicting and evaluating soil erosion and sediment yield. Plots are used to study physical phenomena affecting soil detachment and transport and their sizes are determined according to the experimental objectives and the type of data to be obtained. Studies on interrill erosion due to rainfall impact and overland flow need small plot widths (1-2 m) and lengths (<2 m), while studies on rill erosion require greater plot lengths (3-13 m).

An erosion model that fits the research database well does not assure that the model is adequate for all field conditions. Finally a model is characterised by a spatial (plot, hillslope, basin) and temporal (event, year) scale and its application out of its valid range determines a misuse. An example is a model deduced to estimate average annual soil loss, which is applied at the event scale.

In this paper, at first, the USLE and its revised versions are reviewed taking into account the methodologies developed to estimate the model’s factor and how some factors were modified both to become more physically-based and to improve the soil loss estimate. The WEPP, which represents a new process-oriented technology for soil erosion prediction at different spatial scales, is presented. The available criteria to discriminate between acceptable and unacceptable soil loss estimates are subsequently introduced. Finally, some research needs according to the authors’ points of view are delineated.

## The Universal Soil Loss Equation and its revised Versions

At present, the USLE (Wischmeier and Smith, 1978) is the most widely applied equation for estimating soil loss in the world. The USLE was developed at the National Runoff and Soil Loss Data Centre in cooperation between the United States Department of Agriculture (USDA) - Agricultural Research Service and Purdue University (Wischmeier and Smith, 1978) and it resulted from statistical analysis of more than 10,000 plot-years of basic runoff and soil loss data (Gilley and Flanagan, 2007).

For defining the mathematical structure of the USLE, a reference condition, named as the *unit plot*, was used. The unit plot was defined as a plot 22.1 m long, with a uniform 9% slope, maintained in a continuous regularly tilled fallow condition with up-and-down hill tillage. The unit plot was used to compare soil loss data collected on plots that had different slopes, lengths, cropping and management and conservation practices. The USLE was the result of the work of many individuals over a very long period of time (Laflen and Moldenhauer, 2003) and its deduction was an evolutionary process. Soil erosion data collected, assembled, summarised and statistically analysed by Wischmeier and co-workers (Gilley and Flanagan, 2007) allowed statistical separation of main factors affecting plot soil loss (rainfall, soil, morphology, crop cover and erosion control practices), detection of the representative variable of each factor, and definition of the model mathematical structure.

The USLE predicts long-term average annual erosion by water at an acceptable level of reliability (Larson *et al*., 1997) also in areas of western Europe (Bagarello *et al*., 2008, 2012). At present, the USLE is very widely applied in Europe and in many other Mediterranean countries and no scientific evidence exists that RUSLE (Renard *et al*., 1997) or RUSLE2 without adaption is better than the USLE in European and Mediterranean environments (Abu Hammad *et al*., 2005). In any case, the RUSLE equation conserves the same mathematical structure of the USLE, the revision being limited to improved estimation of some factors.

According to several authors, the USLE and RUSLE are still popular because they: i) combine acceptable accuracy with relative simplicity; ii) have the ability to use quite basic data; and iii) are the only models that can rely on a worldwide distributed dataset (Risse *et al*., 1993; Rosewell, 1993; Liu *et al*., 2001; Hann and Morgan, 2006; Salvador Sanchis *et al*., 2008).

The USLE/RUSLE model was originally designed to predict long-term average annual soil loss (Wischmeier and Smith, 1978), associated with interrill and rill erosion, and for this reason the model tends to over-predict small annual soil losses and under-predict large annual soil losses (Risse *et al*., 1993), although more process-oriented models like WEPP have the same performance (Tiwari *et al*., 2000; Kinnell, 2010). In particular for bare fallow plots in the USLE database, Kinnell (2010) showed that when soils have a low runoff coefficient, USLE over-estimates low event soil losses and under-estimates high event soil losses. The USLE/RUSLE model is based on the assumption that soil loss is not controlled by the runoff capacity to transport detached soil particles, *i.e*. USLE/RUSLE operates as a detachment limited model. In other words, the quantity of sediments leaving an area is determined by the amount of sediment made available by detachment processes. USLE does not take into account that sediment yield can be less than soil loss when the flow transport capacity is less than the detached soil to be transported. Therefore USLE does not estimate erosion for transport-limited cases, and thus when sediment yield has to be calculated then USLE has to be coupled with a mathematical operator simulating sediment delivery processes (sediment delivery ratio).

The USLE is given as:

where *A* (Mg ha^{–1}year^{–1}) is the mean annual soil loss per unit area, *R* (MJ mm ha^{–1}h^{–1}year^{–1}) is the rainfall/runoff erosivity factor, *K* (t ha h ha^{–1}MJ^{–1}mm^{–1}) is the soil erodibility factor, *L* is the slopelength factor, *S* is the slope-steepness factor, *C* is the cover and management factor, and *P* is the support practice factor.

Although the USLE/RUSLE was not originally designed to predict event soil loss, interest exists in determining soil loss at short temporal scales (year, season, event).

Ferro (2010) suggested that, according to the scheme of the plot erosion process and the reference condition used by Wischmeier and Smith (1978), the soil loss, *A _{e}* (kg m

^{–2}), for a rainfall event could be expressed by the following functional relationship:

in which *F* is a functional symbol. The plot is subjected to a rainfall event having known values of both the specific kinetic energy, *E* (J m^{–2}), and the maximum rainfall intensity over a continuous 30-min period during the rainstorm, *I*_{30} (m s^{–1}). The rainfall erosivity is given by the product *EI*_{30} (J m^{–1} s^{–1}). The plot is also characterised by known values of slope length (), slope steepness (*s*), fraction of bare soil (*c*) and fraction of plot surface free from support practices (*p*). The soil covering the plot has a soil erodibility factor, *K* (s^{3}m^{–3}), which is defined as the soil loss per unit of erosivity. The units of the soil erodibility factor *K*, a derived physical variable, are obtained by the ratio between the soil loss (kg m^{–2}) occurring in a *reference condition*, to be defined, and the rainfall erosivity *E I*_{30} (kg m s^{–3}). The reference condition has pre-established values of slope length (_{o}), slope steepness (*s _{o}*), fraction of bare soil (

*c*) and fraction of plot surface free from support practices (

_{o}*p*). The mass depth per unit of time,

_{o}*m*(kg m

_{D}^{–2}s

^{–1}), is the product of soil particle density, ρ

_{s}(kg m

^{–3}), and linear soil depth (m) and it represents the layer of new soil generated in a given time period (s). In other words, the variables included in Eq. (2) are a set of independent quantities, which are necessary and sufficient for a complete definition of the studied physical phenomenon (plot soil loss according to the scheme of Wischmeier and Smith, 1978).

The functional relationship (Eq. 2) represents a physical phenomenon that does not depend on the choice of the measurement units and, according to the -Theorem of the dimensional analysis (Barenblatt, 1987), can be expressed using dimensionless groups. Ferro (2010) deduced the following explicit mathematical form of Eq. (2):

Since A vanishes when λ → 0, or s → 0, or c → 0, or p → 0, then the phenomenon is incompletely self-similar with respect to λ/λ_{o}, *s/s _{o}, c/c_{o}* and

*p/p*and the exact mathematical form of Eq. (3) is (Barenblatt, 1987):

_{o}in which *a*, *m*, *n*, *r* and *q* are numerical constants. In other words, the incompletely self-similar condition occurs because the soil loss has to be set equal to zero when the slope length or the slope steepness are equal to zero, no bare soil occurs in the plot, or support practices are effective enough to prevent any soil erosion of the plot area.

Using the reference condition suggested by Wischmeier and Smith (1978), *i.e*. *s _{o}=*9%,

_{o}=22.1 m,

*c*=1 (bare soil) and

_{o}*p*=1 (no support practice), Eq. (4) becomes:

_{o}Therefore the analysis by Ferro (2010) suggested that the multiplicative structure of the USLE can be obtained by dimensional analysis and self-similarity theory, using the soil erosion representative variables and the reference condition used by Wischmeier and Smith (1978).

One of the key elements of the USLE is the *K* factor, which is a measure of soil erodibility. The values of *K* were originally determined using, with reference to the unit plot condition, soil loss values measured from soil experiments under natural rainfall by the rainfall erosivity *E I*_{30} (Foster *et al*., 1981; Kinnell, 2010). A nomograph was proposed by Wischmeier *et al*. (1971) in which several measured soil properties are combined according to a pre-established scheme to determine *K*. The mathematical approximation of this nomograph, for the cases in which the percentage of silt + very fine sand particles (equivalent diameter ranging from 0.002 mm to 0.1 mm), *f*, does not exceed 70%, is:

in which *M*=*f* × (*f* + *g*), *g* is the percentage of coarse sand (0.1 to 2 mm), *OM* is the soil organic matter percentage, *SS* is the soil structure index and *PP* is the soil permeability index (Wischmeier and Smith, 1978; Foster *et al*., 1981).

The nomograph was developed with specific reference to U.S. soils and therefore it needs testing in other areas of the world, which requires the use of experimental data from those areas (Bagarello *et al*., 2012). A comparison between the experimental erodibility factor and the value calculated by the nomograph should be carried out for different soils, but even a check limited to a single soil may contribute to assess the applicability of the nomograph. This test is particularly onerous since a single measurement of *K* implies an experiment at plot scale prolonged for many years. This check is particularly important for clay soils, which were not included in the group of soils for which the nomograph is expected to be well suited (Romkens *et al*., 1997).

Bagarello *et al*. (2012) used a data set including surface soil samples collected at 1813 sampling points distributed throughout Sicily for developing a regional procedure to estimate the soil erodibility factor of the USLE, *K*, based only on soil textural data.

For each sampling point, the particle size distribution (PSD) was measured using conventional methods and particle size fraction data were classified according to the USDA standards. For each soil sample, *f* and *g* were also determined using the measured PSD. The geometric mean particle diameter, *D _{g}* (mm), was calculated according to Shirazi and Boersma (1984), using the clay, silt and sand particle size classes. The total organic carbon content,

*TOC*(%), was determined by the Walkley-Black method and the organic matter,

*OM*, content was estimated to be equal to 1.724 times the measured

*TOC*value. For each sample, the structure index,

*SS*, of the nomograph by Wischmeier

*et al*. (1971) was estimated using the available soil texture information and the classification reported in Figure 1. Table 1 was used for associating the permeability index,

*PP*, to each sampled soil. The soil erodibility factor,

*K*(t ha h ha

^{–1}MJ

^{–1}mm

^{–1}), and its first approximation,

*K’*(Wischmeier

*et al*., 1971), were calculated by the nomograph of Wischmeier

*et al*. (1971).

The nomograph *K* values were compared with the erodibility values, *K _{R}*

_{86}and

*K*

_{R}_{97}(both in SI units), estimated by the following relationships (Romkens

*et al*., 1997):

Eq. (7) was obtained by using data from 249 soils worldwide. Eq. (8), which was included in the RUSLE manual (Renard *et al*., 1997), was deduced on the basis of measured *K* values for 225 soils.

Table 2 provides the summary statistics of the data used by Bagarello *et al*. (2012) and also summarises both the first approximation, *K’*, and the final value, *K*, of the nomograph soil erodibility.

Eqs. (7) and (8) have the same mathematical form with different numerical coefficients. In particular, in the experimental range 0.002<*D _{g}*<0.99 mm for the Sicilian soils, Eq. (7) yielded different

*K*predictions as compared with Eq. (8) by a maximum of 45%. In addition, the maximum predicted soil erodibility factor (

*K*

_{R}_{86}=0.0423 for

*D*=0.030 mm and

_{g}*K*

_{R}_{97}=0.0439 for

*D*=0.022 mm) was appreciably lower than the maximum possible

_{g}*K*value, equal to 0.09 t ha h ha

^{–1}MJ

^{–1}mm

^{–1}(Wischmeier

*et al*., 1971; Foster

*et al*., 1981).

Figure 2 compares Eqs. (7) and (8) with the (*D _{g}*,

*K*of the nomograph) data pairs. A noticeable scattering of the data points was detected on the

*D*,

_{g}*K*plane and very low values of Nash- Sutcliffe efficiency index,

*NSEI*(0.028 for

*K*

_{R}_{86}and –0.004 for

*K*

_{R}_{97}) were obtained, showing that the proposed relationships cannot be used to obtain reliable estimates of the soil erodibility factor of the nomograph at a selected sampling point. This result may be explained by the origin of the tested equations, based on values of the erodibility factor for classes of soil. The generally high variation of

*K*for a given

*D*value (Figure 2) suggested that a point

_{g}*K*estimation cannot be obtained by empirically re-fitting a relationship of the form of Eqs. (7) and (8) to the Sicilian data. In the analysis of the medians for each log

*D*interval, both

_{g}*K*

_{R}_{86}and

*K*

_{R}_{97}were found to be significantly correlated with the nomograph

*K*(Figure 3). However, the two regression lines did not coincide with the 1:1 line. In other words, grouping soils by textural characteristics was not sufficient to obtain an estimate of the soil erodibility factor coinciding with the one obtained by the nomograph in the considered range of

*D*values.

_{g}For the Sicilian database, the relationship between *K’* and *M* (Wischmeier and Smith, 1978) can be expressed by the following equation:

which is characterised by a coefficient of determination *R*^{2}=0.67 and an exponent practically equal to one. The *K/K’ _{es}* ratios,

*K*being the soil nomograph erodibility and

*K’*the estimate of

_{es}*K’*obtained by Eq. (9), were grouped according to the value of the

*SS*×

*PP*index, and the median of

*K/K’*, denoted by the symbol

_{es}*MR*, was calculated for each

*SS*×

*PP*value. Plotting

*MR*against

*SS*×

*PP*showed that the two variables were related by the following linear relationship:

with *R*^{2}=0.94. Finally, the following relationship was obtained by combining Eqs. (9) and (10):

where *K _{es}* is the estimate of the soil erodibility factor, in SI units, based exclusively on soil textural data. The comparison between

*K*and

_{es}*K*of the nomograph suggested a minimal bias of the predictions and the discrepancies between

*K*and

_{es}*K*of the nomograph did not exceed a factor of two for 1711 data points, corresponding to 94.4% of the complete dataset, with a discrepancy of less than three in 99.2% of the cases. In terms of medians,

*K*was found to be significantly correlated with the nomograph

_{es}*K*(Figure 2), and the regression line was not significantly different from the 1:1 line because the 95% confidence intervals for the slope (0.86-1.41) and the intercept (–0.008-0.005) included one and zero, respectively.

Taking into account the multiplicative nature of the USLE, an approximation in the predicted erodibility by a factor of two or three implies an equal approximation in the predicted soil loss per unit area as compared with a prediction also using *OM* data. It seems reasonable to suggest that this level of approximation could be considered acceptable, at least when the order of magnitude of the soil erosion phenomenon has to be estimated.

The empirical origin of the USLE and its spatial (plot) and temporal (mean annual) scales of application have generated some criticism on the applicability of the equation (Bagarello *et al*., 2015a). For example, Morgan (2005) suggested that soil erosion cannot be adequately described by merely multiplying together six factors given that some interactions between the variables are ignored. Boardman (2006) underlined that, in the European context, the USLE is inappropriate because of differences in rainfall, dominant hydrological processes and landscape diversity as compared with the eastern USA. However, authors criticising the model also recognised its value to solve practical problems (Hann and Morgan, 2006).

A large number of simplified methods have been developed to estimate the model’s factors. Consequently, information on longterm rainfall amount (at daily, monthly or annual scales), soil texture, plot length and steepness, and land use could be enough to obtain some approximate estimate of soil loss. Due to the simplicity of the approach, particularly with reference to the required input data, it seems plausible to presume that the USLE will remain an attractive soil erosion prediction tool, especially for technicians and practitioners. This circumstance justifies efforts to continue improving the empirical equation, taking into account that our experimental and theoretical understanding of the soil erosion process continues to increase. The evolutionary concept in the development of the USLE should be rediscovered, also in terms of sharing data collected in different parts of the world and jointly cooperating to improve empirical soil erosion prediction. In other words, redoing now (with new data, knowledge and technologies) the work by Wischmeier and colleagues at a world scale represents an objective of important scientific and practical relevance, but may be too ambitious at present.

## Modified Universal Soil Loss Equation (MUSLE), USLE-M and USLE-MM

### Soil loss prediction

Foster *et al*. (1982) noticed that the USLE is somewhat unsatisfactory for estimating soil loss at the event scale and observed that erosivity factors that included rainfall amount, rainfall intensity and runoff amount were better predictors of erosivity than *EI*_{30}.

The USLE and its revised version (Renard *et al*., 1997) give no explicit consideration to runoff except when erosion results from snow melt although, for a given event, the soil loss per unit area, *A _{e}* (M L

^{–2}), is given by the product of the runoff amount,

*Q*(L

_{e}^{3}L

^{–2}), and the bulk sediment concentration,

*C*(M L

_{e}^{–3}).

Williams (1975) needed a function to determine sediment delivery from individual storms, so he replaced the rainfall energy factor in USLE with a runoff factor (volume of runoff × peak runoff rate for a storm):

where *SY* (tonnes) is sediment yield, *Q* (m^{3}) is volume of runoff, and m^{3}s^{–1}) and *q _{p}* is peak flow rate (m

^{3}s

^{–1}). Williams validated this modified USLE (MUSLE) function using data from 18 small watersheds in Texas and Nebraska, and found that it could explain 92% of the variation in individual storm sediment yields (Williams, 1975). MUSLE is one of the sediment generation equations currently available in the EPIC, APEX, and SWAT models.

Kinnell (2007), using runoff and soil loss measured in some plots of the USLE database, suggested that the sediment concentration for individual events is dependent on the event rainfall erosivity index, *EI*_{30} (MJ mm ha^{–1}h^{–1}), per unit quantity of rain, *P _{e}* (mm). Although the

*Q*

_{R}EI_{30}index,

*Q*(-) being the event runoff coefficient (=

_{R}*Q*/

_{e}*P*), has an empirical origin, it is based on the concept that the sediment discharged from an eroding area is given by the product of runoff and sediment concentration (Kinnell, 2007).

_{e}Using the *Q _{R}EI*

_{30}term as the erosivity index (Kinnell and Risse, 1998; Kinnell, 2007, 2010), the so-called USLE-M model was empirically derived:

where *A _{e}* (Mg ha

^{–1}) is the event soil loss per unit area,

*K*is the soil erodibility factor,

_{UM}*L*is the USLE slope-length factor,

*S*is the USLE slope-steepness factor,

*C*is the cover and management factor and

_{UM}*P*is the support practice factor (Kinnell and Risse, 1998; Kinnell, 2007). In other words, the factors of USLE-M are named

_{UM}*K*,

_{UM}*C*and

_{UM}*P*for underlining that changing the event erosivity factor from the original

_{UM}*EI*

_{30}to

*Q*

_{R}EI_{30}has impacts on the other factors of the model (Kinnell, 2010). Conceptually, the

*EI*

_{30}index accounts for the effect of runoff on erosion best when the soil surface is impervious,

*i.e*. when

*Q*approaches unity. Di Stefano

_{R}*et al*. (2017) demonstrated that the original structure of USLE-M can also be obtained by applying dimensional analysis and self-similarity theory and using the same soil erosion representative variables and the reference condition adopted by Wischmeier and Smith (1978). The USLE is known to describe a detachment limited erosion process (Wischmeier and Smith, 1978) while the USLE-M, including runoff as a term in the erosivity factor, should be more appropriate to describe a transport limited process. The analysis of soil loss measurements from plots having different lengths can allow one to establish the range of lengths for which the USLE can be applied,

*i.e*. assuring a consistency between the model and the physical process. This analysis can also allow one to quantify the improvement in soil loss prediction by the USLE-M as compared to the USLE. A comparison between the two models in different environments appears advisable to develop an analytical tool of a general validity that makes use of the most appropriate erosivity term for the particular situation.

At the Sparacia experimental station in Sicily (Italy), bare plot soil loss per unit area, *A _{e}* (Mg ha

^{–1}), was found to linearly increase with

*Q*

_{R}EI_{30}(MJ mm ha

^{–1}h

^{–1}) raised to an exponent

*b*

_{1}>1. The model was named USLE-MM (Bagarello

*et al*., 2008, 2010a, 2013a, 2013b, 2014, 2015a, 2015b) since it represents a modified version of the USLE-M, in which

*b*

_{1}=1 was assumed (Kinnell and Risse, 1998; Kinnell, 2007, 2010).

The first steps of the USLE-MM development began in 2008. Table 3 summarises the evolution of the investigations on this model, having the following general form (Bagarello *et al*., 2010a):

where *K _{MM}* (Mg ha

^{–1}per unit erosivity index) is the soil erodibility factor. The latest version of the model is (Bagarello

*et al*., 2015b):

where *s* (m m^{–1}) is the plot steepness, and *b*_{1}=1.37, *K _{MM}*=0.058 Mg ha

^{–1}per unit erosivity index, =0.86,

*a*=2.26 and

*b*=4.07 for the Sparacia site. Eq. (15) was developed by only imposing

*L*=1 for =22 m and

*S*=1 for

*s*=0.09. The model, that was calibrated by using a total of 492 (

*A*,

_{e}*Q*

_{R}EI_{30}) data pairs, had an NSEI of 0.80. The predicted and the measured soil loss values differed by a factor< 57. For the five highest measured values (

*A*>100 Mg ha

_{e}^{–1}), the factor of difference did not exceed 1.19.

It is clear that, currently, the USLE-MM cannot be proposed for routine use since there are several points needing clarification and further development. The first and perhaps most obvious uncertainty is that nothing can still be said about the general validity of the model, and the need to raise the *Q _{R}EI*

_{30}erosivity term to an exponent differing from (in particular, greater than) one could be viewed as a limitation of the approach. Values of

*b*

_{1}have been deduced for only a few datasets till now and we cannot state that these values can also be used in other environments. Therefore, there is the need to experimentally derive

*b*

_{1}at other sites, initially to establish if an exponent greater than one has to be used in general. Should this be the case,

*b*

_{1}estimation procedures need to be developed. Currently, we know that, in the USLE-MM,

*b*

_{1}is greater than 1 but it does not exceed approximately 1.6 in central and southern Italy locations (Bagarello

*et al*., 2008, 2010a, 2013a, 2013b, 2015a, 2015b) and it is equal to 1.55 at a Chinese experimental station (Gao

*et al*., 2012). Therefore, it appears plausible to suppose that a single

*b*

_{1}value could be used in general for practical purposes. This reasoning is not novel in development of empirical soil loss prediction methodologies. For example, the exponent of the plot length factor of the USLE (Wischmeier and Smith, 1978) varied widely from year to year and average values for locations varied between 0 and 0.9 (Laflen and Moldenhauer, 2003). Notwithstanding that, an exponent of 0.5 was suggested for most practical applications.

The soil erodibility factor has units that are specific to the value of *b*_{1}. This circumstance could be viewed as another problem limiting the interest for the model, given that two sites differing by the *b*_{1} value will also have non-comparable *K _{MM}* values. However, a common erodibility value was found to be usable for different plots established at Sparacia, which suggested that the soil erodibility factor of the USLE-MM is expressive of an intrinsic soil property (Bagarello

*et al*., 2010a). Hoping that a constant

*b*

_{1}value greater than 1 will be found to be a plausible assumption, experimentally determining

*K*at other sites according to Foster

_{MM}*et al*. (1981) is necessary to develop estimation criteria that should maintain the simplicity that characterises the USLE soil erodibility estimation procedure (Wischmeier

*et al*., 1971).

Bagarello *et al*. (2015b) showed that the developed equation for estimating *S* was in close agreement with the *S* predictions obtained according to Nearing (1997). Moreover, similar results were obtained between the new predictive relationship for *L* and the one used in the RUSLE (Renard *et al*., 1997). In other terms, the lack of specifically calibrated relationships for *L* and *S* does not seem to preclude a general use of the USLE-MM since alternative, and largely validated, relationships can be used to predict event soil loss with the USLE-MM. Evidently, confirming this finding in other environments is advisable to give generality to the conclusion by Bagarello *et al*. (2015b).

### Prediction of event runoff as a plot erosivity factor

Including a runoff term in an empirical plot soil loss model allows interpreting and predicting a soil loss per unit plot area that does not increase with . This is a common occurrence at the Sparacia site (Figure 4), and it depends on the fact that a decrease in runoff with (Figure 5) is offset by an increase in sediment concentration (Bagarello *et al*., 2015a). For the Sparacia area, rain intermittency was suggested to be a factor determining less runoff on longer plots (Bagarello *et al*., 2015a). An inverse relationship between *Q _{R}* and was detected in other investigations, including those by Joel

*et al*. (2002), showing decreasing runoff coefficients as the plot area increased from 0.25 to 50 m

^{2}, and Parsons

*et al*. (2006), working on plots varying in length from 2 to 28 m. Therefore, runoff-driven models such as the USLE-M or the USLE-MM, providing new support to the finding that runoff is an important predictor of soil erosion (Foster

*et al*., 1982), are particularly attractive since they take into account a phenomenon that is not considered in the USLE rainfall erosivity factor and the consequence is that a large amount of total variability in bare plot soil loss can be explained with these approaches. In other words, runoff represents an advance for erosion prediction (Nearing, 2000a) which justifies efforts to develop simple methods to predict plot runoff at the event temporal scale. A variety of approaches, briefly reviewed below, can be found in the literature.

Surface runoff prediction at the plot spatial scale and the erosive event temporal scale was the focus of the investigation by Yu *et al*. (1997), who considered 1-min rainfall and runoff rates from tropical and subtropical regions of Australia and south-east Asia. The model attempted to capture two important features in the runoff hydrographs and particularly that, even at the plot scale, the lag between runoff and rainfall can be significant and, due to the enormous spatial variation of soil hydrodynamic parameters, the apparent infiltration rate, *i.e*. the difference between rainfall and runoff, is strongly linked to rainfall intensity. Indeed, as rainfall intensity increases, the proportion of the surface for which rainfall intensity is greater than the maximum infiltration rate will increase, which implies an increase in rainfall excess and runoff rate. The great variability of infiltration rates at the plot scale was documented, as an example, by Bagarello *et al*. (2010b, 2013c), since the intensive sampling (176 measurement points) of 11×4 m^{2} plots established on the clay soil of Sparacia yielded point measurements of saturated soil hydraulic conductivity, *K _{s}*, varying by more than three orders of magnitude.

The model by Yu *et al*. (1997) has two components. The first component addresses infiltration, and therefore rainfall excess, and the second one deals with runoff routing down the slope length. The model assumes that a one-parameter exponential distribution can be used to characterise the spatial variation of the maximum infiltration rate and also that rainfall excess produced anywhere within the plot becomes measurable runoff at the plot outlet. Yu *et al*. (1997) minimised the sum of the squared errors between observed and predicted runoff rates to estimate the three parameters of the model, namely the initial infiltration amount before runoff occurs, a spatially averaged maximum infiltration rate which could be achieved if the entire plot produces runoff, and a dimensionless routing parameter. The optimised parameters were found to be usable for predicting total runoff for other runoff producing events not included in the optimisation procedure. Operating the model in a continuous mode was considered advisable by Yu *et al*. (1997) to isolate the effects of antecedent soil moisture conditions and soil crusting development on the event by event variation in the estimated parameter values and hence to develop estimation relationships for the model’s parameters. The approach by Yu *et al*. (1997) yields more information (runoff rates) than that strictly necessary to apply the USLE-M or the USLEMM (total runoff).

Especially in arid and semi-arid environments, the infiltration capacity of the soil and, consequently, runoff are influenced by sealing or crusting, *i.e*. by the very thin surface soil layer, having hydraulic parameters that differ, even substantially, from those of the undisturbed zone below (Chen *et al*., 2013). Vandervaere *et al*. (1998) suggested a conceptually simple approach to predict crustinduced surface runoff from a plot at the event temporal scale, making use exclusively of measured soil data. In particular, the starting point of the developed method was that the only measurable data, in addition to the initial (θ_{i}) and saturated (θ_{s}) soil water content, were the parameters of the Gardner (1958) hydraulic conductivity function for both the crust and the sub-soil. These last data can be obtained by tension infiltrometer experiments according to Vandervaere *et al*. (1997). Vertical infiltration is described by the Green-Ampt model since the water retention curve does not need to be known in this case. The soil is described as a two-layer medium with a homogeneous surface crust and semi-infinite homogeneous subsoil. The crust is saturated over the wetted depth and it has a hydraulic conductivity that coincides with the saturated hydraulic conductivity of the crust. Due to the impeding effect of the upper layer, the subsoil does not saturate during infiltration. Consequently, both the soil water content and the hydraulic conductivity of the subsoil above the wetting front do not coincide with those corresponding to fully saturated conditions. Determination of the pressure head at the wetting front in the subsoil takes into account non-saturation of this layer during infiltration. Runoff is simply given by the difference between rainfall and infiltration. An iterative scheme is used for calculations to allow infiltration from rainfall of variable intensity to be modelled. In particular, phases with and without runoff may occur successively. Vandervaere *et al*. (1998) tested their approach against total runoff from 26×5 m^{2} plots established in Niger. Rainfall was measured every minute and simulations were performed at 1 s time steps. The simulations allowed the authors to identify the substantial effect of the crust on soil infiltrability, *i.e*. earlier occurrence of runoff and higher runoff volumes as compared with the non-crusted case. In the approach by Vandervaere *et al*. (1998), surface detention and redistribution of water into the soil during the event are not taken into account and this represents a possible limitation needing further developments. Moreover, the field measured crust conductivity values used by Vandervaere *et al*. (1998) were characterised by a very large variability (Vandervaere *et al*., 1997). Therefore, there is a problem in establishing how the plot can be hydraulically parameterised. Another point to be considered is that plot length effects on runoff are not explicitly taken into account. Finally, a problem inherent in the Green-Ampt type models (Vandervaere *et al*., 1998) is the inconsistency between the crusted and non-crusted cases. Indeed, a crust with an infinite depth (onelayer condition) and a crust with a finite depth overlying a subsoil with characteristics identical to those of the crust (two-layer condition) yield different results although the configurations are, in fact, the same.

A modified Soil Conservation Service - Curve Number (SCSCN) model was applied by Gao *et al*. (2012) to obtain an estimate of the plot runoff coefficient specifically usable with a USLE-MM type soil erosion model. In short, the event runoff, *Q _{e}* (mm), is predicted by the following relationship:

where *P* (mm) is total precipitation depth, _{c} is the initial abstraction coefficient, representing losses due to interception, surface storage and infiltration, *S* (mm) is the potential maximum retention, and *M _{s}* (mm) is the moisture of the soil profile before the start of the storm. The parameters

*S*and

*M*are calculated by the following relationships:

where *CN* is the curve number, *i.e*. a dimensionless variable that depends on land use, hydrological soil group, hydrologic conditions and antecedent moisture conditions, and *P*_{5} (mm) is the amount of the 5-day antecedent precipitation depth. The model was applied on 5 to 13 m long plots established with a slope gradient of 19-25° on a soil rich in silt (71-73%) and poor in clay (approximately 4%). Performance was good for both small and large runoff events when: i) was set equal to 0.05 instead of the value assumed in the original development of the model (=0.2); and ii) the *CN* value was determined taking into account the effect of the antecedent moisture condition, slope gradient and initial abstraction ratio. In particular, the *CN* value corresponding to normal moisture conditions was determined from the USDA-NRCS (2004) handbook for each runoff plot. These values, that were valid for a 5% slope steepness, were adjusted to the actual values using a relationship developed by Huang *et al*. (2006) in a similar environment to that considered by Gao *et al*. (2012). The adjusted values were converted to dry and wet conditions, depending on the magnitude of *P*_{5}, according to Hawkins *et al*. (1985). The relationships by Hawkins *et al*. (2002) were finally applied to convert the *CN* based on =0.20 to the *CN* based on =0.05. The model applied by Gao *et al*. (2012) undoubtedly has limitations, including determination of antecedent moisture condition on the basis of *P*_{5} alone, inability to take plot length effects on surface runoff into account, lack of consideration of the impact of rainfall intensity and duration on runoff amount, and use of locally developed relationships. On the other hand, Gao *et al*. (2012) were able to show that the approach can yield satisfactory results if properly applied. Moreover, the practical simplicity of the modified SCS-CN model agrees with the simplicity of the USLE family of models. Therefore, testing the model in other environments and situations appears advisable.

Another attempt to predict the runoff coefficient with the specific objective to simplify application of the USLE-MM was carried out by Todisco *et al*. (2015). In particular, the continuous rainfall- runoff model MISDc (*Modello Idrologico Semi-Distribuito in continuo*; Brocca *et al*., 2011) was used to estimate runoff volumes. This model incorporates a limited number of parameters and it is characterised by low computational requirements. However, it needs parameterisation, which was carried out by maximising the Nash-Sutcliffe efficiency index between the estimated and the observed runoff volume values for a set of calibration events (Todisco *et al*., 2015). The semi-distributed MISDc model was originally developed for flood simulation over large areas (>100 km^{2}; Brocca *et al*., 2011), *i.e*. for a very different condition from that corresponding to event runoff estimation at the plot scale. Notwithstanding this, Todisco *et al*. (2015) concluded that the MISDc model provided fairly accurate results even at the smaller plot scale.

Approximately fifteen years ago, Nearing (2000a) reminded us that Walt Wischmeier, at one point in his career, attempted to publish a Universal Runoff Equation but was unsuccessful in doing so. Probably, we are still far from a universal equation for predicting runoff and, perhaps, such an equation cannot be developed at all due to the complexity of the process to be modelled. However, our inability to simply and confidently predict event runoff at the plot or field scales should be considered as a solicitation to proceed with determination towards that goal, possibly with a broad involvement of the international scientific community, also because the lack of plot runoff data represents the most common situation.

There are still several problems to be solved with reference to the different alternative approaches that can be found in the literature to estimate plot runoff and hence to make use of the USLE-M or USLE-MM models a practical possibility. One is the need to go more in depth on the adopted theoretical or empirical approaches. For example, the lack of consideration of redistribution processes in the model by Vandervaere *et al*. (1998) could suggest a better performance of this approach under more continuous rain showers. In arid and semi-arid climates, however, a short duration of effective rain showers, that is, shorter than the concentration time necessary for a continuous flow along the hillslope, is rather common (Yair and Raz-Yassif, 2004). Another problem is related to the model parameterisation which means that efforts should be made to establish how to measure or determine a given parameter but also how to obtain a reliable estimate of a parameter for the area of specific interest. Many measurement methods exist, for example, for infiltration rate or soil hydraulic conductivity that have a clear physical meaning, but the passage from the point information to the areal characterisation is not free from uncertainties. Other parameters, such as those used in the SCS-CN model, have a more empirical character and they are not directly measurable. The need to adapt the estimation procedure of these parameters to the particular situation considered for the simulations (Gao *et al*., 2012) implies that generalising estimation procedures is currently not possible.

Surely the existing literature helps us to make choices on the approaches to be developed. For example, an important problem in central Italy is that high soil losses can occur when the soil is initially dry but simulated runoff increases with the antecedent soil moisture conditions. Therefore, explicitly considering the soil’s response to wetting in these cases (*e.g*., crusting phenomena) appears necessary to obtain reliable *Q _{R}* predictions (Todisco

*et al*., 2015). With respect to this point, the runoff prediction method by Vandervaere

*et al*. (1998) appears to be a good candidate method deserving further investigation. Other approaches have been developed very recently (Chen

*et al*., 2013, 2016) and they certainly merit consideration in future research.

### An alternative to runoff coefficient in USLE-MM type Models

According to Todisco *et al*. (2015), the pre-event soil water content, θ, could be used instead of the runoff coefficient to correct the rainfall-runoff erosivity factor, because obtaining soil moisture data is easier and less expensive than estimating surface runoff. Todisco *et al*. (2015) tested their approach at Masse (Italy), considering both satellite-derived and modelled soil water content but not direct measurements of θ. The performance of the (θ *EI*_{30})^{b1} erosivity index varied with the considered period of the year. In particular, high intensity rainfall events occurring in summer on initially dry soil yielded high sediment losses. Therefore, a small erosivity factor was associated with high soil losses. Neglecting the dry period events improved the performance of the alternative criterion (Todisco *et al*., 2015). In any case, the general conclusion by the authors was that including soil moisture in the erosivity factor of the USLE type model enhanced the capability of the model to account for variations in event soil loss as compared with the USLE alone.

A monitoring of θ should be carried out at the plot scale to establish the potential to use these data for event soil loss prediction. Moreover, we need to determine what has to be done when soil loss is high but the soil is initially dry.

## The Water Erosion Prediction Project

Since 1985, the USDA has been developing the WEPP model, as a new generation technology for soil erosion prediction on hillslope profiles, fields, and small watersheds (Flanagan *et al*., 2007). WEPP is a process-based, distributed parameter, continuous simulation erosion prediction computer simulation model. It simulates the important physical processes that affect soil loss and sediment delivery, including climate (*e.g*., rainfall occurrence, depth, duration, intensity), infiltration, percolation, plant growth, residue decomposition, soil tillage disturbance, runoff, detachment by raindrop impact and shallow overland flow (interrill erosion), detachment by excess flow shear stress in small channels (rill erosion), sediment transport, and sediment deposition (Flanagan and Nearing, 1995; Flanagan *et al*., 2007). Additionally when simulating small watersheds, the model determines soil detachment, sediment transport and deposition in larger channels (ephemeral gullies, grass waterways, *etc*.), and sedimentation in impoundments such as culverts, filter fences, or farm ponds (Lindley *et al*., 1995).

WEPP can function on a single storm basis using observed or generated climate inputs, or on a continuous simulation basis, again using either observed or generated climate inputs. Typically for model applications for soil conservation planning or natural resource assessments, it is necessary to use generated climate. The WEPP modelling system includes a stochastic weather generator called CLIGEN (CLImate GENerator) which utilises long-term observed monthly weather statistics from over 2600 stations in the United States. Rainfall occurrence is determined using a second order Markov chain, and a double exponential storm shape is assumed.

Every day in a simulation, WEPP updates the status of the soil (moisture, infiltrability, erodibility, critical shear stress, *etc*.), the status of the live plants (biomass, leaf area index, canopy cover, canopy height, root mass, root depth, *etc*.), and the status of the dead plant residues (flat residue mass in 3 pools (current, last, and old), buried residue mass in 3 pools, dead root mass in 3 pools, standing residue, *etc*.). If a tillage operation is scheduled to occur on a day, it will affect the soil bulk density, porosity, effective hydraulic conductivity, and adjusted erodibility parameters (increase adjusted interrill and rill erodibilities and decrease adjusted critical hydraulic shear stress), as well as reduce standing and flat residue mass and cover values (Alberts *et al*., 1995). Harvest operations for an annual crop convert live biomass to dead residue, and remove a fraction of the biomass as crop yield (Arnold *et al*., 1995). Residue management operations (Stott *et al*., 1995) may decrease surface residue cover (burning, baling, removal), may increase surface residue cover (residue addition), or may convert standing residue to flat (shredding/cutting).

After all soil, plant, and residue management operations and parameter updating has occurred, WEPP then determines if any rainfall, snowmelt, and/or irrigation has occurred on the day, and if it has, whether any surface runoff results. The model uses a modified Green-Ampt Mein-Larson infiltration equation to calculate excess rainfall through time sequences in each rain storm event, with consideration of depressional storage from soil roughness (Stone *et al*., 1995). Peak runoff rate is determined using a solution of the kinematic wave equation, and soil detachment is modelled as a quasi-steady-state process at this rate, for an effective duration that preserves the total storm volume previously determined (Foster *et al*., 1995; Stone *et al*., 1995).

Soil detachment on hillslope profiles is determined as the sum of detachment (or deposition) rate of sediment in rill channels and the interrill sediment delivered to these rills (laterally), with distance downslope. For a slope profile region having uniform soil and cropping/management [overland flow element (OFE)], the interrill delivery of sediment will be constant, while rill detachment or deposition will vary with slope gradient, flow shear stress, hydraulic roughness, sediment transport capacity, and the sediment load already present in the rill channel flow (Foster *et al*., 1995).

The steady-state sediment continuity equation used in WEPP is:

where *G* (kg s^{–1}m^{–1}) is the sediment load in the rill flow *x* (m) is the distance downslope, *D _{f}* (kg s

^{–1}m

^{–2}) is the rill erosion rate, and

*D*(kg s

_{i}^{–1}m

^{–2}) is the interrill sediment delivery to the rills. Interrill sediment delivery is always greater than or equal to zero, and is calculated in WEPP using:

where *K _{iadj}* (kg s m

^{–4}) is adjusted interrill erodibility,

*I*(m s

_{e}^{–1}) is effective rainfall intensity, σ

_{ir}(m s

^{–1}) is interrill runoff rate,

*SDR*is an interrill sediment delivery ratio that is a function of the soil random roughness, row sideslope, and sediment particle size distribution (Flanagan and Nearing, 2000),

_{RR}*F*is an adjustment factor to account for sprinkler irrigation nozzle impact energy variation,

_{nozzle}*R*(m) is the rill spacing, and

_{s}*w*(m) is the rill width. For croplands, the baseline condition for erodibilities is for a freshly tilled soil with no crop or residues present. Interrill (and rill) erodibilities are adjusted to lower values as a function of soil consolidation, growing plants, and residue cover (Alberts

*et al*., 1995).

Rill erosion rate is either positive in the case of detachment, or negative in the case of deposition. For the detachment case:

where *D _{c}* (kg s

^{–1}m

^{–2}) is the rill flow detachment capacity, and

*T*(kg s

_{c}^{–1}m

^{–1}) is the rill flow sediment transport capacity (Foster

*et al*., 1995). Sediment transport capacity down a hillslope profile in WEPP is computed using a simplified equation (Finkner

*et al*., 1989) whose coefficients are determined by utilising the full Yalin (1963) sediment transport equation at the end of the OFE. The expression in parentheses in the equation is a

*sediment-feedback*term that decreases rill detachment when the water at a point is already laden with sediment. Clear water has the maximum detachment capacity, which is calculated with:

where *K _{radj}* (s m

^{–1}) is the adjusted rill erodibility, τ (Pa) is the flow shear stress acting on the soil, and τ

_{cadj}(Pa) is the adjusted critical flow shear stress, below which no rill detachment will occur.

For the case where sediment load exceeds sediment transport capacity, deposition will be predicted using:

where β is a raindrop-induced turbulence factor (usually 0.5 under rainfall), *v _{f}* (m s

^{–1}) is an effective sediment particle fall velocity, and

*q*(m

^{2}s

^{–1}) is flow discharge per unit width. Since sediment load

*G*is greater than the sediment transport capacity

*T*, the

_{c}*D*values computed with this equation will be negative.

_{f}Erosion calculations are made at 100 evenly-spaced points on each OFE. For detachment a Runge-Kutta numerical technique is used to solve for sediment load at each point, while for deposition a closed-form equation solution is employed (Foster *et al*., 1995; Flanagan and Nearing, 2000).

Many more details on the hydrologic, erosion, and auxiliary processes, equations and logic are available in a multitude of WEPP publications, especially Nearing *et al*. (1989), Stone *et al*. (1992), Flanagan and Nearing (1995) and Flanagan *et al*. (2007). The auxiliary WEPP hillslope model components include water balance, plant growth, residue management and decomposition, soil tillage disturbance and consolidation, and irrigation (Flanagan and Nearing, 1995).

While the WEPP computational science model is coded in FORTRAN, a variety of user interface programs now exist that allow much easier application of the program. These include a Windows interface coded in C++ (Flanagan *et al*., 1998), an ArcView/ArcGIS extension for GIS application of the model (Renschler, 2003; Flanagan *et al*., 2013), web-based interfaces for hillslope simulations (Elliot, 2004), and web-based interfaces for watershed simulations (Flanagan *et al*., 2013).

The WEPP model has undergone considerable validation over the years, in particular using data from the USLE database (Zhang *et al*., 1996; Liu *et al*., 1997; Tiwari *et al*., 2000). Tiwari *et al*. (2000) used 1600 plot-years of USLE erosion plot data for verification and validation of uncalibrated WEPP model soil loss predictions, and also compared the WEPP results to those from USLE and RUSLE. They found that, for the 20 sites tested, the Nash and Sutcliffe model efficiency index (*NSEI*) was 0.71 for average annual soil loss for WEPP, and this was comparable to that from RUSLE (0.72) and USLE (0.80).

Pandey *et al*. (2008) applied WEPP to the 2793 ha Karso agricultural watershed in India, and examined its ability to predict observed runoff and sediment yield. After development of input files and calibration with data from 1996, validation simulations were conducted using data from 1992, 1993, 1995, 1997, and 2000. Annual *NSEI* ranged from 0.78 to 0.92 for the sediment yield events in these five years. Singh *et al*. (2011) also applied the WEPP model in the eastern Himalaya region of India to the 239 ha hilly Umroi watershed under bun agriculture. Calibration was performed using observed runoff and sediment loss data for 86 storms in 2003, and validation conducted using data from 98 storms in 2004. For the validation period storms, *NSEI* was 0.87 for runoff and 0.90 for sediment loss.

Detailed information on WEPP model application, calibration, and validation can be found in Flanagan *et al*. (2012). The article also include two case studies, one for a hillslope profile application of WEPP, and the other for a watershed application. Laflen *et al*. (2004) also provides additional information on WEPP model validation studies.

## Evaluating soil erosion models using measured plot data

Measurement is a sequence of steps or operations that produce a value representing the magnitude of a selected variable (Toy *et al*., 2002). Different experimental methods can be applied to collect plot soil erosion data but none has proven to be fully satisfactory under any circumstance (Bonilla *et al*., 2006). Each method has advantages and disadvantages and a good knowledge of the feasibility and limitations of a given method is necessary to properly plan experiments and also make a correct use of the data for calibrating and/or testing soil erosion models (Boix-Fayos *et al*., 2007). The technique for measuring soil loss can be based on two antithetic approaches: i) direct survey of the eroded area at given time intervals to detect the amount of soil removed from the area; ii) measurement of the amount of soil that, leaving the eroded area as sediment in the flow, accumulates in a collection device (*e.g*., sedimentation tank) where the measurement is carried out. The first approach can be applied using erosion pins, photogrammetery, or at-site measurements with radionuclide techniques. Since soil loss measurements from runoff plots are widely applied and were used for empirically deducing the USLE, errors associated with this technique will be discussed in the following section.

### Errors in soil loss measurement on equipped runoff plots

In a plot equipped for soil loss measurement, runoff from a bounded area is collected and carried to a sampling unit by a conveyance system (Zobisch *et al*., 1996; Bonilla *et al*., 2006). A simple method for making these measurements is to store all runoff in a single tank or a sequence of tanks. For example, runoff and associated sediments from plots (22×8 m to 44×8 m) established at Sparacia are collected into a storage system consisting of three tanks, each having a capacity of approximately 1 m^{3}, that are arranged in series at the base of each plot (Figure 6) (Bagarello *et al*., 2010a). The stored water volume is easily determined from the known geometry of the tank and a water depth measurement. Sediment concentration can be measured by either catching the whole sediment amount or collecting a sample of the mixed suspension (Lang, 1992; Bagarello and Ferro, 1998). As reported by Ciesiolka *et al*. (2006), the practice of stirring the suspension, taking a sample by immersing a collector to a substantial depth beneath the water surface in the tank, and oven-drying the sample provided the experimental data collected in the US and summarized in the USLE (Wischmeier and Smith, 1978) and it was also adopted worldwide (Rosewell, 1993).

The suspended sediment concentration measured for a sample is representative of the whole collected suspension if the watersediment mixture is well mixed (complete-mixing condition) and the measured suspended concentration assumes the same value in each measurement point of the tank, equal to the actual one. Consequently, the sediment amount is calculated by multiplying the sample concentration by total runoff volume. However, Lang (1992) tested a bottle sampler for sampling clay soil-water mixtures and concluded that the actual suspended particle concentration was underestimated by a factor of two. Consequently, Lang (1992) threw doubt on the reliability of plot soil loss data collected by using a runoff sampling technique. According to Ciesiolka *et al*. (2006), the investigations on sampling techniques cast doubts on the correctness of the USLE soil loss database, that has been widely used, even to validate process-based erosion technologies (Zhang *et al*., 1996; Tiwari *et al*., 2000).

Bagarello and Ferro (1998) and Bagarello *et al*. (2004) evaluated factors affecting the measured sediment concentration and they derived calibration curves for sediment storage tanks. In a storage tank, the mean concentration *C _{m}* (g L

^{–1}) is measured by the concentration profile, unavoidably corresponding to an incomplete mixing condition, obtained by collecting samples of given volume at different depths. In particular, 10 sampling taps along a vertical line are used at Sparacia (Figure 6). For a mean concentration obtained by integrating the measured concentration profile along the axial vertical of a tank, Bagarello and Ferro (1998) theoretically deduced that the actual concentration

*C*(g L

^{–1}) is linked to

*C*by the following equation:

_{m}in which *b* is a coefficient that has to be experimentally determined. For the clay soil of the Sparacia area, *b* was found to be independent of both the field worker and the water depth, *h*, in the tank, and therefore it was assumed constant and equal to 4.12 (Bagarello *et al*., 2004).

Using the same procedure at the Masse experimental station (Umbria), in which the soil has a silty-loam texture, Todisco *et al*. (2012) confirmed the validity of Eq. (23) for a given water level in the tank but they also showed that *b* increased as *h* decreased. This last result was considered to be a consequence of the reduction in the total number of solid particles that are stored in the tank for lower *h* levels.

The validity of Eq. (23) for a given water level in the tank and a dependence of *b* on *h* were also detected for the clay soil of Sparacia when the sediment concentration, *C _{bt}*, obtained by dipping a bottle sampler into the tank to reach the bottom was considered instead of

*C*(Bagarello and Ferro, 1998). In particular, the slope of the

_{m}*C vs C*relationship increased with

_{bt}*h*, probably because higher

*h*levels made the mixing procedure more difficult (

*i.e*., more incomplete) due to the greater volume of stored water. According to Ciesiolka

*et al*. (2006), the soil loss data used to develop the USLE database were obtained under the assumption that

*C*=

_{bt}*C*,

*i.e*. that the sediment concentration in the suspension collected by a single bottle coincided with the actual sediment concentration in the tank.

Therefore, there are clear signs that the sediment concentration in the tank is underestimated if the collected field data are not corrected. The error is systematic, *i.e*. independent of *h*, if *b* in Eq. (23) does not vary with the water level in the tank, such as in the case of the Sparacia soil. The error may depend on *h* for the same soil (Sparacia) but another sampling procedure (*C _{bt}*) or the same sampling procedure (

*C*) for another soil (Masse).

_{m}Evidently, there is the need of additional experimental investigation on errors in soil loss measurement by storage tanks. These errors are not commonly taken into account although they can have a noticeable effect on the reliability of the soil loss data (Bagarello and Ferro, 2017).

### Physical model concept

Determining the quality of the predictions by a soil erosion model requires establishment of a criterion to discriminate between acceptable and unacceptable soil loss estimates. According to Nearing (1998), the best possible model to predict the erosion from an area of land is a physical model of the area that has similar soil type, land use, size, shape, slope and erosive inputs. On the basis of this reasoning, Nearing (2000b) suggested that a soil erosion model performs well if the difference between the model prediction and the measured plot data value lies within the population of differences between pairs of measured values.

Bagarello and Ferro (2012) tested the physical model concept with the soil loss data of Sparacia. Two data points were obtained from the data collected, for a given event, at two plots of given geometric characteristics. For the first data point, one value (A) of the pair was chosen to serve as the measured value of soil loss (kg m^{–2}) and the other (B) was considered to be the predicted value from the physical model. For the second data point, value (B) was used as the measured value and value (A) as the predicted one. The relative difference, *Rdiff*, was then calculated according to the following relationship (Nearing, 2000b):

where *P _{e}* is the predicted soil loss value from the physical model represented by the replicated plot and

*M*is the measured soil loss value. The Sparacia data set included many

_{e}*A*values lower that 0.01 kg m

_{e}^{–2}, that was the minimum soil loss in the US investigation by Nearing (2000b). The 95% occurrence interval for the data developed by this last author included approximately 88-89% of data collected at Sparacia, depending on the plots considered. Taking into account that this discrepancy was moderate,

*i.e*. a few percentage units, and considering that a large sample size and a wide variety of conditions were considered in the US study, the conclusion by Nearing (2000b) that the developed analysis should be usable for model validation studies in general appeared to be reasonable.

More recently, bare plot soil loss data collected at Sparacia and Masse were used by Bagarello *et al*. (2015c) to initially demonstrate that an intercept, *b _{i}*, closer to 0, a slope,

*b*, closer to 1 and a higher coefficient of determination,

_{s}*R*

^{2}, of the

*P*linear relationship have to be expected when the physical model is defined in terms of perfect planimetrical equivalence to the sampled plot. In other terms, the physical model of a plot is another plot having the same length and the same width of that plot.

_{e}vs M_{e}With reference to the developed Italian database (Masse + Sparacia, *N*=1468 data pairs), the percentage of *Rdiff* values falling within the 95% occurrence interval defined by Nearing (2000b) was equal to 89%. A similar result (88%) was obtained by only considering the *M _{e}* values greater than 0.01 kg m

^{–2}(

*N*=1194),

*i.e*. the minimum soil loss in the US investigation.

The premise of the analysis by Nearing (2000b) was that the measured data with greater erosion rates showed, on average, less relative differences between replicates. This tendency was not so clear with the Italian database (Figure 7). In particular, three regions were distinguishable on the *Rdiff vs M _{e}* plot, with approximate boundaries between adjacent regions at 0.01 and 1 kg m

^{–2}. Relative differences decreasing with an increase in

*M*were detected for

_{e}*M*>1 kg m

_{e}^{–2}. With reference to these data (

*N*=316), the percentage of

*Rdiff*values falling within the 95% occurrence interval was equal to 85%. The relative differences decreased with a decrease in

*M*for

_{e}*M*<0.01 kg m

_{e}^{–2}. This result was considered to be reasonable, since the differences between two replicated plots are expected to decrease when the event has a low erosive power. Two plots yield similar results because the rainfall-runoff event is able to detach and transport only a small amount of soil particles. For 0.01<

*M*<1 kg m

_{e}^{–2}, a relationship between

*Rdiff*and

*M*was not detectable and the data points practically occupied all the space of the graph, suggesting that plot heterogeneities have a more noticeable impact on soil loss for intermediate levels of the erosion phenomenon. Therefore, the premise by Nearing (2000b) was only confirmed with reference to the highest erosion rates (

_{e}*i.e*.,

*M*>1 kg m

_{e}^{–2}) but even in this particular case the correspondence of the predicted occurrence interval with the data was not fully satisfactory.

From a scientific point of view, the sign of the difference between *P _{e}* and

*M*has to be determined to establish how to improve a soil erosion predictive tool. From a practical point of view, however, the absolute

_{e}*P*difference is enough to establish the accuracy level of the predictions. The ⎜

_{e}– M_{e}*P*⎜values were found to increase with

_{e}– Me*M*according to the following relationship (Figure 8):

_{e}with an *R*^{2}=0.72 and a 95% confidence interval of the exponent of 0.88-0.94, denoting a non-linearity of the relationship. Eq. (25) can be viewed as an alternative approach for applying the physical model concept by Nearing (2000b) since it predicts, for a given soil loss value (*M _{e}*), what is the mean absolute difference associated with the sampling of another, identical plot. A soil loss prediction by a model is accurate enough if the absolute difference with the measured value does not exceed the ⎜

*P*⎜value calculated by Eq. (25). The least restrictive criterion, using a relationship enveloping all data points, could alternatively be proposed. An intermediate criterion between the suggested regression line and a data enveloping line could also be developed by carrying out a frequency analysis of the data divided into half log-cycle intervals, similar to the one carried out by Nearing (2000b) to estimate 95% occurrence intervals for the data.

_{e}– M_{e}The approach by Bagarello *et al*. (2015c) appears promising because it was found to be valid for the entire range of the measured soil loss values. However, it cannot be suggested for general use because the analysis had an empirical character and data were only collected at two experimental stations. Developing a unique database by the contribution of authors working in different parts of the world appears advisable to develop the most robust criterion possible for plot soil erosion model validation studies.

## Research needs

Our knowledge of the soil erosion process is still incomplete, as indicated by the large number of papers on the dynamics of this process that continue to be published in international journals. Notwithstanding this, there is also the need to go ahead with the prediction of this process, to allow people to establish if and how the phenomenon has to be mitigated in an area of interest. In the following, a list of topics on soil erosion prediction deserving consideration, according to the point of view of the authors, is reported.

### USLE-MM: possibility to use a common exponent in different environments

Foster *et al*. (1982) noticed that the USLE can be unsatisfactory for estimating soil loss at the event scale and observed that erosivity factors that included rainfall amount, rainfall intensity and runoff amount were better predictors of erosivity than *EI*_{30}, in which *E* (MJ ha^{–1}) is the event kinetic energy and *I*_{30} (mm h^{–1}) is the maximum 30-min rainfall intensity observed during the event. The possibility to use a common exponent *b _{1}* of the

*Q*

_{R}EI_{30}index,

*i.e*. establishing if the validity of the model is not limited to a single station, has to be investigated. The circumstance that the USLE-M (

*b*

_{1}=1) was proposed with reference to the US database (Kinnell and Risse, 1998) could be viewed as a suggestion that the results obtained at Sparacia (

*b*>1) represented an exception and that

_{1}*b*

_{1}=1 should be assumed in general. However, the

*b*

_{1}=1 hypothesis was never tested with the US data to the best of our knowledge. On the other hand, there are signs in the literature that

*b*

_{1}>1 can also be detected in different environments since Gao

*et al*. (2012) obtained an exponent of 1.55 for a Chinese silt soil.

### Soil water content as an alternative to runoff coefficient for predictive purposes

The application of a rainfall-runoff driven equation for estimating soil loss needs the availability of measured or estimated runoff values. Soil moisture data are used as input for a soil water balance model or could be used for estimating a modified rainfallrunoff erosivity factor in which the soil water content pre-event is included. Soil moisture datasets, such as the ERA-Interim/Land reanalysis data of the European Centre for Medium-range Weather Forecasts (ECMWF) and satellite retrievals from the European Space Agency - Climate Change Initiative (ESA-CCI), are useful for this kind of analysis.

### Testing WEPP in other environments

The WEPP model has been fairly extensively tested and validated within the United States, but to a much lesser degree internationally. Additional model application and validation studies under a wider range of climatic, soil, and topographic conditions would help to better define input parameter refinements needed for greatly varying soil and vegetation conditions. Of particular concern are the baseline infiltration and erodibility inputs in the soil input file since the default parameterisation equations reported in the WEPP model documentation (Flanagan and Nearing, 1995; Alberts *et al*., 1995) were based upon results from rainfall simulation studies conducted on 33 important cropland soils in the United States. These may not provide satisfactory runoff and soil loss predictions when applying the WEPP model to greatly different soils (*e.g*., tropical soils, volcanic origin soils). Also, application to data from slope and climate conditions that are more extreme than those used in the model development may highlight shortcomings and potential areas in WEPP that could be improved in the future. The watershed components of WEPP have been much less tested than those for hillslope profile erosion computations. Additional studies using observed data from small watersheds that included ephemeral gullies would allow substantial improvement and confidence in applying WEPP for those conditions.

### Estimating event and annual soil loss of given return Period

Previous experimental investigations showed that a large proportion of total plot soil erosion over a long time period is generally due to a relatively few, large storms. Consequently, erosion models able to accurately predict the highest plot soil loss values have practical importance since they could allow improvement of the design of soil conservation practices in an area of interest. From an engineering point of view, a probabilistic analysis of the maximum annual soil loss data should be carried out using the historical series in order to estimate the soil loss of a known return period (Edwards and Owens, 1991; Hession *et al*., 1996; Larson *et al*., 1997; Baffaut *et al*., 1998; Mannaerts and Gabriels, 2000; Bagarello *et al*., 2010c, 2011a, 2011b). Long historical sequences of maximum annual soil loss data have to be used to develop a reliable frequency analysis and hence to estimate the soil erosion variable with a given return period (Baffaut *et al*., 1998; Mannaerts and Gabriels, 2000). Consequently, the design event soil loss is unpredictable if little or no erosion data exist at a site, which is a rather common occurrence. On the other hand, long sequences of rainfall data are generally available and they can be used to predict historical sequences of soil loss by a suitable algorithm.

If a USLE-type model can be calibrated using measured data (rainfall-runoff erosivity, erodibility, geometrical characteristics, vegetation cover) corresponding to the maximum annual event soil loss, this equation becomes a robust method to indirectly estimate the maximum annual soil loss at event scale in areas where soil loss is not measured. Then the values of the maximum annual event soil loss estimated by the calibrated USLE-type equation can be used to develop a frequency analysis useful for predicting the value corresponding to a given return period.

### Reliability of the soil erosion measurement methods

Sampling the collected suspension in a storage tank is a common procedure to obtain soil loss data and a calibration curve of the tank is required to obtain actual concentration values from those measured by sampling. However, literature suggests that using a tank calibration curve was not a common procedure in the past.

Many references indicate that a widespread practice for measuring plot soil loss is stirring the suspension stored inside a tank, taking a sample by immersing a collector to a substantial depth beneath the surface of the suspension inside the tank, and ovendrying the sample. The measured concentration is then assumed to coincide with the actual one. This practice provided the experimental data collected in the US and summarised in the USLE, and it was also adopted worldwide.

Some investigations cast doubt on the reliability of plot soil loss data collected using this simple runoff sampling technique due to the assumed coincidence between measured and actual concentration.

Specific tests should be carried out, with different soil types, since soil loss data collected by a sampling procedure can have a noticeable effect on the calibrated empirical models for soil loss prediction (Bagarello and Ferro, 2017).

### Rill and gully erosion measurement by ground monitoring techniques (drones, image-based techniques)

Several researchers have pointed out the importance of measuring and modelling ephemeral gully erosion (Casali *et al*., 1999; Nachtergaele and Poesen, 1999; Capra *et al*., 2009) and stressed the importance of this erosion type in the overall sediment yields of agricultural catchments. Accurate ground measurements of ephemeral gully erosion are frequently carried out using a microtopographic profiler, a tape or a ruler, to assess some channel cross-sections (Di Stefano and Ferro, 2011; Di Stefano *et al*., 2013). These ground measurements, which are simple and low-cost, allow measurement of the cross-sectional areas and the reach length and calculation of the ephemeral gully volume, even if, according to Casali *et al*. (2006), the limit of these simple methods is due to the lack of information on the errors associated with ground measurements of the cross-sections.

Traditional aerial photogrammetry was successfully applied for large-scale and long-term investigations (Ionita, 2006) and the need of increasing resolution for short-term applications has accelerated the use of aerial drones (Carollo *et al*., 2015).

The recent advances in automatic 3D-photo reconstruction techniques for oblique images from uncalibrated and non-metric cameras coupled with the availability of photogrammetric software packages encouraged the use of close-range photogrammetry to investigate soil erosion processes. Terrestrial digital photogrammetry is characterised by high spatial resolution (centimetres to millimetres) and minimal impact of the field activity on both the ground measurements and farming operations and it requires a time-consuming work for post-processing the digital photos. Image-based modelling creates a digital terrain model using a set of photographs taken from the same surface (Frankl *et al*., 2015). The procedure involves a solution with camera model parameters and scene geometry simultaneously using redundant information coming from oblique images and without using control points during the composition of the model (Gomez-Gutierrez *et al*., 2014).

Further investigation should be carried out for testing the applicability of an image-based technique for measuring rill erosion.

### Physical model

The performances of a soil erosion model have to be tested for establishing the expected reliability of the soil loss estimates. The quality of the predictions can be established only if a criterion to discriminate between *acceptable* and *unacceptable* soil loss estimates is available. The *physical model* represented by a replicated plot has to be considered the best possible, unbiased, real world model.

Different criteria (Eqs. 24 and 25) should be tested for comparing the performances of the selected models. In particular the approach of Bagarello *et al*. (2015c) should be tested using datasets from different parts of the world.